[VA] SRC #017 - April 1st, 2024 Spring Special

[Valentin Albillo]

      
Hi, all,

Welcome to my new SRC #017 - April 1st, 2024 Spring Special

Once again April 1st is here and I want to celebrate both it and the new season with this Spring Special where I'm proposing a number of mini-challenges for you to tackle with your favorite vintage HP calc, plus interesting facts not widely known (if at all).

Note: No hard rules, post here whatever you want as long as it's on topic and absolutely NO CODE PANELS (lest I'll consider you a troll only fit for my Ignore list,) but I'd appreciate it if you'd use vintage HP calcs (physical/virtual,) unless you're clueless as far as programming them is concerned.

As April 1st is also referred to as April Fool day, I'll take the last letter from April (L) and the two last ones from Fool (OL) to title each separate Section with a jocular prefix i.e. LOL !


1. LOL the First:  Squares

Let's start nice and easy. Square numbers have been a source of beauty and admiration since Pythagoras did his thing with them millennia ago. Just look at these small (related) beauties !

      375,501² = 141,001,001,001
      751,002² = 564,004,004,004

Now, if these are cute just consider the sheer amounts of beauty you'll discover when dealing with their bigger relatives, so this mini-challenge asks you to calculate the following 2² squares:

      130834904430015239² = ?
      471287714788971663493899² = ?
      25781108305591628417975738² = ?
      141422082876067219949805050005² = ?

I know you can simply paste them in Wolfram Alpha or use the multiprecision library in your RPL model or mutiprecision canned software in any device or web site and get done with it, but it would be so lame that you'll risk ridicule. What I'm asking you to do is to program your own multiprecision squaring routine in your vintage HP calc and use it to get the results. It's not that difficult at all, you know ... that is, if you've got what it takes ..

I'll post my own squaring routine that I wrote for the HP-71B from scratch, just a mere 483 bytes long and taking the form of a user-defined function so that it can be called right from the command line. It produces results like this, which you can profitably use for testing your code:

      3147926784726726563780030042374237187751² =

          9909443041999946686063013207932562462851613614976494122324802303779777224438001

Remember: if you decide to solve this mini-challenge, you must post both the (beautiful) results and the code which produces them, do not post results alone. Comments are most welcome, plus you might try getting more results like these.


2. LOL the Second:  GCD

OK, full throttle now. You may remember the GCD function appearing in my recent SRC #016 - Pi Day 2024 Special, where it was used in the Second appearance to count the number of co-prime pairs of random integers in a range and then this count was used to compute an approximation to \(\pi\).

     Sample GCD values:  GCD(35, 13) = 1,  GCD(35, 14) = 7 and GCD(35, 15) = 5

In this mini-challenge you must write code for any HP calc to find out the answer to this simple question:

      As we have that  GCD(15, 4) = 1  and  GCD(15, 5) = 5,  for what value of x is  GCD(15, x) = 2 ?

As always, you must post both results and code. Comments also most welcome. I'll post mine next Sat/Sun.


3. LOL the Third:  Random

You may be aware that advanced 12-digit HP models of old incorporated an excellent RNG (pseudo-Random Number Generator) which could generate a trillion (i.e. 10¹²) full 12-digit real (pseudo-)random numbers in the interval [0-1) before repeating. This RNG passes the Spectral test so it's extremely reliable for use in simulations and other advanced topics (e.g. multidimensional integration, Monte Carlo algorithms, etc.) without fear of any bias or short period degrading the results.

As far as I know, the implementations for the HP-71B, the RPL models, the HP-42S and Free42 (and perhaps many other HP models) are functionally identical so they produce exactly the same sequence of RNDs from the same seed. Now, this mini-challenge's question is:

Two consecutive RNDs cannot be exactly equal or the generation of subsequent RNDs would be stuck in a loop. But how close can two consecutive RNDs actually be ?

Use RANDOMIZE 1 or the equivalent instruction to specify the value 1 as the seed and write code to output the list of ever closer consecutive RNDs and their index N in the sequence, like this:

>RUN
      N       First RND         Next RND          |Difference|
 ----------------------------------------------------------------
      2      .731362440213     .77207218067      .040709740457
     13      5.64471991805E-2  6.30768172146E-2  6.6296180341E-3
    125      .805774019056     .803607575861     .002166443195
    ...

If you intend to go very deep into the sequence you'd better use a fast virtual model such as Emu71/Win or Free42 for faster results and/or greater depths. You might even boldly examine the entire trillion-strong sequence to find out the absolute closest pair !

I'll post my 3-line, ~90 byte solution which produces the above, plus a 4-line, 149-byte variant which stops right before a given pair of consecutive RNDs is generated so that you can generate them manually right from the command line (i.e. >RND;RND) and see for yourself how close they indeed are.

And remember, post code and results. What's your record ? Me, I've found a pair of consecutive RNDs which are only 0.00000000001 apart ! Can you do better ?


4. LOL the Fourth:  Logs

This is not a mini-challenge but a somewhat unexpected fact I've found, so you don't need to do anything but read on and eventually comment.

I've always thought that the two logarithmic functions available in the 10-digit HP-15C (LOG and LN, base 10 and base e respectively) would have essentially the same accuracy overall, but lo and behold, I've found that LOG seems to be substantially more accurate in some cases than LN.

For instance, let's consider computing log_b(5⁷) / log_b(5), which should return exactly 7 for any base b ≥ 2. But the actual results are:

     LN(5⁷) / LN(5)      ->  7.000000004      { 4 ulp from the exact value }
     LOG(5⁷) / LOG(5)  ->  7                       { exact }

and the same happens with other powers N of 5, e.g. for N = 2 to 9 we have:

     N      LOG          LN        
 ---------------------------------
     2    2.00000001   2.00000001
     3    3            3.00000001
     4    4            4.00000001
     5    5.00000001   5.00000001
     6    6            6.00000002
     7    7            7.00000004
     8    8.00000001   8.00000002
     9    9            9.00000001
-----------------------------------
 Total ulps:    3           13               { more than 400% larger error overall }

This much larger error puzzled me no end but I thought that perhaps the limited 10-digit accuracy (13-digit internally) might be playing a role and decided to test the same computations using the 12-digit (15-digit internally) HP-71B, which has LN, LOG10 and LOG2 (respectively base e, base 10 and base 2.) These were the results for power 8:

     >STD
     >LN(5^8) / LN(5)              ->  8.00000000004      { 4 ulp from the exact value }
     >LOG10(5^8) / L0G10(5)   ->  8                           { exact }
     >LOG2(5^8) / LOG2(5)      ->  7.99999999999      { 1 ulp from the exact value }

so we see that all three results differ among them, being in error by 4 ulp, 0 ulp and 1 ulp, respectively, and again the same happens with other powers N of 5, e.g. for N = 2 to 9 we have:

     N      LOG10           LN        
 ---------------------------------------
     2    2               2.00000000001
     3    3               3.00000000001
     4    3.99999999999   4.00000000001
     5    5               5.00000000001
     6    6.00000000001   6.00000000001
     7    7               6.99999999999
     8    8               8.00000000004
     9    8.99999999999   9.00000000002
----------------------------------------
  Total ulps:      3           12               { 400% larger error overall }

and again the LN error is much larger, so having two extra digits didn't help at all. For LOG2 (not shown in the table above) the Total ulps for the range is 6 ulps, half LN's total error but twice LOG10's total error. In the end, LOG₁₀ seems to be the most accurate logarithm available.

If you want to pursue the matter, you can try the above examples in your HP models, both 10-digit and 12-digit, from the first HP-35 up to the latest RPL models, to see if LOG₁₀ and LN differ that much in their respective errors. Testing other range of values (here from 5² = 25 to 5⁹ = 1,953,125) could be revealing as well. Are there arguments with even larger errors ?

At any rate, comments and discussion would be most welcome.


5. LOL the Fifth:  Gamma

Here's something peculiar I noticed while playing around with my HP-71B several decades ago, in Experimental Mathematics fashion.

Back then, I executed this loop directly from the command line to list the values of \(\Gamma\)(10^-1), \(\Gamma\)(10^-2), ..., \(\Gamma\)(10^-10), and got the following unexpectedly peculiar results:

   >DESTROY ALL @ FOR N=1 TO 10 @ N;GAMMA(10^(-N))@ NEXT N

      N     \(\Gamma\)(10^-N)
    -------------------
      1  9.51350769867
      2  99.4325851191
      3  999.423772485
      4  9999.42288323
      5  99999.4227942
      6  999999.422785
      7  9999999.42278
      8  99999999.4228
      9  999999999.423
     10  9999999999.42

where the evergrowing integer part of each value is clearly 10^N - 1, while the fractional part seems to be quickly converging to some limit around ~0.42 but the 12-digit HP-71B lacks the accuracy needed to refine it further, so this mini-challenge is:

You must write code to find this limit to much greater accuracy (say 10-12 digits or more,) perhaps by simply using Free42 Decimal to compute the values. Once done, answer these two questions:

• Can you estimate the most accurate value you got ? Because, as the integer part is exponentially growing, even the 34-digit Free42 Decimal will soon begin to lose digits in the fractional part, as it happened with the HP-71B results above. In that case, an RPL model with arbitrary multiprecision capability would surely help.
        
• Can you identify the symbolic, closed form of that numeric value ? My IDENTIFY program v2.0 (unreleased) sure can.
  

As always, post here your results and code (comments also welcome,) and please try to steer clear of references such as Wolfram, MathWorld, Wikipedia, ... as you'd just ruin the pleasure of doing a little Experimental Mathematics all on your own.


6. LOL the Sixth:  Miscellanea

1. [HP-71B specific] Try and deduce what result will be output by this expression without actually executing it: some regular numeric value ? Perhaps Inf or NaN ? An error message ? A long-running or indefinite internal loop ?
   
         NEIGHBOR(INTEGRAL(FNROOT(EPS,EPS,EPS),EPS,EPS,EPS),GAMMA(EPS))
   
   Once you've thoughtfully made your deduction, execute the expression and see what comes out. Was your deduction correct ? Can you explain the result obtained ?
   
   Please post your deduction & comments on the result.
         
2. [HP-71B specific] I executed this command-line expression on my HP-71B but it just resulted in System Error. What did I do wrong ? Bad addresses ? Forbidden range for PEEK$ ? Can you find out what's wrong ?
   
         FOR I=64204 TO 64215 @ DISP CHR$(HTD(REV$(PEEK$(DTH$(2*I),2)))); @ NEXT I
         
3. Some nice results I got (you can check the √s with Free42 Decimal):
   
         √ 95888      = 309.657875727390470000000975517...
         √ 22008840 = 4691.3580123456789961013...      
   
         (HP-71B specific, RPL/12-digit models too, maybe other 12-digit models) In radians:
   
       >EXP(ACOS(430/433))  ->  1.12500000006     { ~  9/8 }
       >EXP(ACOS(538/541))  ->  1.11111111113     { ~ 10/9 }
   
   Notice that 433 = 430 + 3 and 541 = 538 + 3, i.e. both fractions are of the form N / (N+3) and produce results of the form (N+1) / N so perhaps there's some hidden pattern at large here. Or not !
         
4. [HP-71B specific] What does this HP-71B user-defined function compute ?
   
         10  DEF FNC(M,N)=M!/N!/(M-N)!
         
5. [HP-71B specific] The HP-71B does not allow for variable names beginning with a 2-letter or more prefix (e.g. AB, CD7, MYTAXES, etc.) as many other pocket computers do, but executing this assignment in a program or from the keyboard ...
   
         FORM=STOP
   
   ... doesn't result in a Syntax Error so what gives ?
         
6. Last but not least, you're probably acquainted with the number 42, not only because it's the model number of the renowned HP-42S, later simulated by the Free42 virtual version, but also because according to The Hitchhiker's Guide to the Galaxy 42 is "The Answer to the Ultimate Question of Life, the Universe, and Everything", which is saying a whole lot.
   
   However, this doesn't end 42's importance to this Universe and to prove the point you must write code to compute to full accuracy this simple function S(x) given the argument x:
   \[ S(x) = \frac{1}{\sum_{n=0}^{\infty} \binom {2 n}{n}^3 \frac{{nx \space + \space 5}}{2^{12 n \space + \space 4}}} \]Use your code to compute S(42). Astounding, isn't it ?  
   
   Do not use any built-in summation capability but write and use instead your own code and post here both code and result. I'll post my simple 3-line HP-71B user-defined function in a few days.

---

That's all. Waiting for your clever solutions and/or comments galore, I'll post mine next Saturday/Sunday so you've got plenty of time to concoct and polish your code. And no CODE pannels, please.

V.