[VA] SRC #016 - Pi Day 2024 Special

[Valentin Albillo]

      
Hi, all,
     
Well, time to post my original solutions. As I said, my intent is to celebrate unusual \(\pi\) appearances by encouraging interested people to write simple code to put various incarnations of \(\pi\) right on the display. Thus, my solutions are pretty plain-vanilla code with no real regard for optimization or speed. That said, here we go:

1. First appearance

Solve this equation for x in [1, 6].

      
      
These two HP-71B command-line sequences will do the job right from the keyboard, namely:

   >DESTROY ALL @ STD @ I=INTEGRAL(0,1,1E-10,ATAN(ATANH(IVAR))/IVAR) -> 1.02576051092
   >FNROOT(1,6,FVAR*(LN(GAMMA(1/FVAR))-LN(GAMMA(1/2+1/FVAR))-LN(FVAR)/2)-I)

      3.14159265358

and, as \(\pi\) is 3.14159265359 when rounded to 12 digits, we've got the correct value save 1 ulp.

The value obtained is exact in the limit, not an approximation. The first command-line sequence computes the definite integral, then the second one solves the equation. Were it not for the 96-char command-line length restriction, both lines could be merged into a single one and then neither the DESTROY ALL statement nor the temporary I variable assignment would be needed.

I chose this particular expression because I liked the arctan(arctanh(...)) composite function on aesthetical grounds. Furthermore, this definite integral has a closed-form exact value, which is quite infrequent for integrals of transcendental composite functions, and it's not trivial to numerically compute to full accuracy so that the subsequent solve procedure can work with the best value obtainable.

2. Second appearance

You can get a nice approximation to the value of \(\pi\) (exact as the number of tries goes to infinity) by following these simple steps. First set Count to zero and select a number of tries N, then:

1. Choose two random integers A and B in [1, N].
2. Check if they are coprime, i.e. they have no common factors, in which case increment Count. Keep on performing steps 1 and 2 for N tries, then:
3. Output

This 3-line, 102-byte HP-71B program will output the value obtained by using 10, 100, 1000, ..., 1,000,000 tries:

   1  DESTROY ALL @ RANDOMIZE 1 @ FOR K=1 TO 6 @ N=10^K @ T=0
   2  FOR I=1 TO N @ T=T+(GCD(IP(RND*N+1),IP(RND*N+1))=1) @ NEXT I 
   3  DISP USING "8D,2X,D.4D";N,SQR(6*N/T) @ NEXT K

>RUN
        10   2.9277
       100   3.2733
      1000   3.1755
     10000   3.1448
    100000   3.1449
   1000000   3.1415  (1 ulp)  

The convergence to \(\pi\) is extremely slow but it's happening. Two integers A and B are coprime if their greatest common divisor (GCD) is 1. My solution above uses the JPC ROM's GCD function to check every pair A, B of randomly generated integers, while keeping a tally of how many coprimes were found.

Why does this procedure produce \(\pi\) ? Because the probability P_N that two randomly chosen integers in [1, ..., N] are coprime approaches 6/\(\pi\)² as N tends to infinity, and you can derive \(\pi\) from that value.

3. Third appearance

Solve this equation for x, where \(\phi\) is the Golden Ratio = (1+√5)/2.

     

This 2-line, 91-byte HP-71B program will output the root obtained by using 10, 100, 1000, ..., 100,000 terms of the infinite product. It does not use any solver but simply isolates x :

   1  DESTROY ALL @ FOR K=1 TO 5 @ P=1 @ FOR N=1 TO 10^K @ P=P*(1-1/(100*N*N)) @ NEXT N
   2  DISP USING "8D,2X,D.6D";N-1,5/(P*(1+SQR(5))/2) @ NEXT K

>RUN
        10   3.138604
       100   3.141280
      1000   3.141561
     10000   3.141590
    100000   3.141592

or simply from the command line, directly using 100,000 terms of the product:

   >DESTROY ALL @ FIX 6 @ K=(1+SQR(5))/2
   >P=1 @ FOR N=1 TO 10^5 @ P=P*(1-1/(100*N*N)) @ NEXT N @ 5/(P*K)

      3.141592

For 12-digit models such as the HP-71B, there's no point in using more than 100,000 terms of the product, as no additional accuracy will be obtained and in fact the rounding errors resulting from multiplying so many terms will begin to degrade the accuracy noticeably while also greatly increasing the run time.

I chose this particular appearance of \(\pi\) because I liked the fact that it produced the exact (in the limit) value of \(\pi\) from \(\phi\) and a short, unremarkable infinite product, as if \(\phi\) was the one which begat \(\pi\) ...

4. Fourth appearances

\(\pi\)'s prefixes of any length (first N digits of \(\pi\) for finite N) do actually appear inside irrational numbers, square roots for instance, like these:

● The 7-long prefix 3141592 appears in:

     √3485  =  59.033888...4733453141592349004..  at decimal 822
     √26401 = 162.483845...1905663141592268465..  at decimal  69
     √82777 = 287.709923...3523123141592695866..  at decimal  45   (8-long, actually !)
   
● The 10-long (rounded) prefix 3141592654 appears in:

     √2424609  = 1557.115602...6896933141592654600235..  at decimal 170
     √40850970 = 6391.476355...4838773141592654698929..  at decimal 112

See if you can find other appearances of assorted \(\pi\)'s prefixes in various irrationals. The smaller the argument and the earlier the appearance, the better.


I was hoping that interested readers would look into the matter and come up with some interesting appearances of \(\pi\)'s prefixes in other irrational values, perhaps within logs, trigs or exponentials but it wasn't to be, though I didn't require posting code so people could use more powerful harware/software, yet no one tried. A pity.

As for the appearances I posted, I believe them to be a fine compromise between smaller arguments and nearer locations for the first appearances inside the irrational value.

Of course you can always find prefixes of any length even for arguments as small as 2 (i.e. √2) but they'll first appear at tremendously deep locations, or else you can find the prefix beginning at the very first decimal place but for enormously large arguments. Thus the need to compromise.

5. Fifth appearances

Solve the following equations, either using a program or directly from the keyboard:

• Polynomial:  4 x³ - 22 x² + 29 x + 2 = 0  and  9 x⁴ - 19 x³ + 28 x² - 70 x - 344 = 0
• Fermat-like:  2063^x + 8093^x = 8128^x  and   1198^x + 4628^x = 4649^x
• Transcendental (where Γ is the Gamma function):   Γ ln(7 x⁵/19) = 16
  

See if you can find similar polynomial equations with roots even closer to \(\pi\). The smaller the degree and the size of the coefficients, the better.


These HP-71B command-line sequences will suffice, no programming needed:

● The first polynomial equation is solved by this sequence:

   >FIX 6 @ FNROOT(3,4,((4*FVAR-22)*FVAR+29)*FVAR+2)  ->  3.141593
      
● The second polynomial equation is solved by this sequence:

   >FIX 9 @ FNROOT(1,6,(((9*FVAR-19)*FVAR+28)*FVAR-70)*FVAR-344)  ->  3.141592654

● The first Fermat-like equation is solved by this sequence:

   >FIX 10 @ FNROOT(1,6,2063^FVAR+8093^FVAR-8128^FVAR)  ->  3.1415926535   (1 ulp)

● The second Fermat-like equation is solved by this sequence:

   >STD @ FNROOT(1,6,1198^FVAR+4628^FVAR-4649^FVAR) ->  3.14159265363   (4 ulp)

● The transcendental equation is solved by this sequence:

   >FIX 10 @ FNROOT(1,16,GAMMA(LN(7/19*FVAR^5))-16) ->  3.1415926535   (1 ulp)

All the roots are approximate values of \(\pi\), none are exact. The two polynomials are what I believe to be the minimal polynomials of degrees 3 and 4, i.e. the ones with the smallest coefficients for the achieved accuracy, and I found them with what I jokingly call my "Poor's Man PSLQ", which is actually a very simple search (not related to true PSLQ, mind you) with several tricks and optimizations applied to help run the search faster.

I used a similar ad-hoc approach to find the very short & simple transcendental equation, which I consider amazingly remarkable as it's got a root so close to \(\pi\) while featuring just 6 digits, 2 functions and 3 arithmetic operations in all. How's that for economy of resources ?

6. Sixth appearance

Given the following recurrence, with x₀ = 0 and { } denoting the fractional part function,

      

write a program to compute its succesive terms x₁, x₂, x₃, .. and for each term output in hexadecimal the value of d_n = IP(16 x_n ), where IP is the integer part function. Try and output as many correct hex digits as possible.


My original solution for the HP-71B is this 3-line, 112-byte plain-vanilla program:

1  DESTROY ALL @ X=0 @ FOR N=1 TO 9 @ M=4*N
2  X=FP(16*X+(30*M^2-89*M+64)/(8*M^4-64*M^3+178*M^2-206*M+84))
3  DISP BSTR$(IP(16*X),16);" "; @ NEXT N @ DISP

   >RUN  ->  2 4 3 F 6 A 8 8 8

and produces the first 9 (next is 2 instead of the correct 5) hexadecimal digits of \(\pi\)'s fractional part.

Using 12-digit arithmetic (HP-71B, HP-42S, etc.) that's the most correct digits we can get but using multiprecision arithmetic it's possible to obtain many more, e.g. using the 34-digit Free42 Decimal emulator with my 49-step, 89-byte program below:

      01  LBL "HEXPI"  13   x     25  X<>Y        37  RCLx 00     49  END
      02   0           14  64     26  X^2         38   +
      03  STO 00       15   -     27  30          39  FP
      04   1           16   x     28   x          40  STO 00
      05  STO 01       17  178    29  89          41  16
      06 ►LBL 00       18   +     30  RCLx ST T   42   x
      07   4           19   x     31   -          43  IP
      08  RCLx 01      20  206    32  64          44  HEXM
      09  ENTER        21   -     33   +          45  STOP
      10  ENTER        22   x     34  X<>Y        46  ISG 01
      11  ENTER        23  84     35   ÷          47  LBL 00
      12   8           24   +     36  16          48  GTO 00 ►

Executing XEQ "HEXPI" and pressing [R/S] after each digit shown we get:

      2 4 3 F 6 A 8 8 8 5 A 3 0 8 D 3 1 3 1 9 8 A 2 E 0 3

and the next hex digit is 6, which should be 7 so we've got 26 correct hex digits in all.

This recurrence has been checked up to 1,000,000 terms (i.e. hex digits, all of them correct) and it's known that there can only be a finite number of incorrect digits, if at all.

Furthermore, if the sequence of x_i happens to be uniformly distributed in [0-1] (as indeed it seems to be the case,) then this would imply the normality of \(\pi\) in base 2 (and thus in base 16 as well,) which would be a far-ranging, all-important achievement worthy of a Fields Medal or two, so perhaps you might want to try.

Last but not least, I was hoping from the get-go that some of you would use one of the powerful RPL models with multiprecision capabilities to get many more hex digits and luckily lo and behold, DavidM fulfilled my expectations. If any of you would like to try, these are the first 500 hex digits of \(\pi\) (ignore the 3.) for you to check your results:

   3.243F6A8885 A308D31319 8A2E037073 44A4093822 299F31D008 2EFA98EC4E 6C89452821 E638D01377 BE5466CF34 E90C6CC0AC
     29B7C97C50 DD3F84D5B5 B547091792 16D5D98979 FB1BD1310B A698DFB5AC 2FFD72DBD0 1ADFB7B8E1 AFED6A267E 96BA7C9045
     F12C7F9924 A19947B391 6CF70801F2 E2858EFC16 636920D871 574E69A458 FEA3F4933D 7E0D95748F 728EB65871 8BCD588215
     4AEE7B54A4 1DC25A59B5 9C30D5392A F26013C5D1 B023286085 F0CA417918 B8DB38EF8E 79DCB0603A 180E6C9E0E 8BB01E8A3E
     D71577C1BD 314B2778AF 2FDA55605C 60E65525F3 AA55AB9457 48986263E8 144055CA39 6A2AAB10B6 B4CC5C3411 41E8CEA154

7. Seventh appearance

Let x be the only real root of

      x³ - 6 x² + 4 x - 2 = 0

Compute ln(x²⁴ - 24) ÷ √H, where H is the number of hours in a week minus five.


This HP-71B command-line sequence will return the root, an approximation to \(\pi\) correct to a full 12 digits:

   >STD @ LN(FNROOT(1,10,FVAR^3-6*FVAR^2+4*FVAR-2)^24-24)/SQR(24*7-5) @ PI

      3.14159265359
      3.14159265359

but again, using the 34-digit Free42 Decimal arithmetic capabilities we can check whether there's additional accuracy to unfold. Using this trivial 15-step, 28-byte Solver program below, we can conduct the following session right from the keyboard:

   01  LBL "PIEQ"   09   x
   02  MVAR "X"     10   4
   03  RCL "X"      11   +
   04  ENTER        12   x
   05  ENTER        13   2
   06  ENTER        14   -
   07   6           15  END
   08   -

   - [SOLVER] -> Select Solve Program
   - Touch PIEQ, 0, touch X, touch X -> 5.31862821775
   - Press EXIT twice, 24, y^x, 24, -, LN, 24, ENTER, 7, x, 5, -, SQRT, ÷

           3.14159265359

   - Press and hold SHOW -> 3.141592653589793238462643383279504   (34 correct digits save 1 ulp)

thus this extremely simple procedure can produce ~ 34 correct digits of \(\pi\) with utmost, unexpected simplicity. Amazing indeed !

---

Well, that'll be all for now.

Thank you very much for your interest to all of you who posted solutions and comments, namely Gerson W. Barbosa, J-F Garnier, DavidM, C.Ret, ttw, Juan14, johnb and jonakeys, much appreciated and I'm glad you enjoyed it.

See you all next April 1^st.

V.