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

[Valentin Albillo]

      
Hi, all,

First of all, thanks to all of you who were interested in this April 1st thread and in particular to those who posted code, results and/or comments, namely J-F Garnier, Juan14, C.Ret, John Keith, Gerson W. Barbosa and JoJo1973. Thank you very much for your interest and continued appreciation.

However, I feel somewhat let down by the fact that only 6 people 6 took the trouble to post anything at all (out of the ~9,400 members currently registered,) but that's life and now it's time for my original solutions and comments, which for sheer message-length reasons I'll post two LOL at a time, so let's get the party started with the first two, i.e. LOL the First: Squares and LOL the Second: GCD ...


1. LOL the First:  Squares

This mini-challenge asks you to calculate the following 2² squares. You must program your own multiprecision squaring routine in your vintage HP calc and use it to get the results.

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

My original solution is this 10-line, 415-byte HP-71B user-defined function, which accepts the number to be squared as a string and returns its square as another string:

  1  DEF FNS$[200](P$) @ STD @ DIM M,L,P,I,J,U,A$[206]
  2  OPTION BASE 1 @ A$=FNL$(P$) @ U=LEN(A$) DIV 6 @ K=10^6
  3  DIM A(U),C(2*U),C$[12*U] @ MAT A=ZER @ MAT C=ZER
  4  FOR I=1 TO U @ A(U+1-I)=VAL(A$[I*6-5,I*6]) @ NEXT I @ FOR I=1 TO U
  5  M=A(I) @ L=I @ FOR J=1 TO U @ P=M*A(J) @ C(L)=C(L)+P @ P=RES
  6  IF P>=K THEN C(L)=RMD(P,K) @ C(L+1)=C(L+1)+P DIV K
  7  L=L+1 @ NEXT J @ NEXT I @ I=2*U+1 @ REPEAT @ I=I-1 @ UNTIL C(I)
  8  C$=STR$(C(I)) @ FOR I=I-1 TO 1 STEP -1
  9  C$=C$&FNL$(STR$(C(I))) @ NEXT I @ FNS$=C$ @ END DEF
 10  DEF FNL$[206](A$) @ P=RMD(LEN(A$),6) @ FNL$=RPT$("0",6*(P#0)-P)&A$

Note: Branching isn't needed thanks to the keywords provided by the JPC ROM (REPEAT..UNTIL). It also uses MAT..ZER from the Math ROM and RPT$ from the STRNGLEX LEX file for speed and convenience, all of them replaceable by slower BASIC code if necessary.

• Lines 1-3 do the initialization.
• Line 4 splits the input string into an array of 6-digit numeric values.
• Lines 5-7 perform the multiprecision squaring and locate the leftmost nonzero element.
• Lines 8-9 conjoin the result array elements into a result string and return it.
• Line 10 is an auxiliary user-defined function which adds needed leftmost 0s, if any.
  

My program chops the input into 6-digit values (multiplying two 6-digit numbers never exceeds the 12-digit range of the HP-71B) instead of doing the squaring digit by digit, which would be about 36x more processing and so ~30x slower.

Also, using strings for input/ouput is extremely convenient for multiprecision values, as the user simply enters the value between quotes and the result can be directly stored in a string variable for further processing e.g. combined with other similar multiprecision UDFs. For instance, raising a multiprecision value to the fourth power is as simple as FNS$(FNS$(value))

Finally, there's no limit to the size of the result squares other than maximum string length (65,500+ characters i.e. digits). As listed above, it allows for up to 200-digit squares but that limit can be changed by simply replacing the constants 200 and 206 by the desired maximum size.

Now let's use FNS$ from the command line to quickly compute the requested squares:

>DESTROY ALL
>FNS$("130834904430015239")

      17117772217211221211117217772227121

>FNS$("471287714788971663493899")

      222112110111011100020110111110102200012010222201

>FNS$("25781108305591628417975738")

      664665545464645645665646644665564654546645556644644

>FNS$("141422082876067219949805050005")

      20000205525005225202000505202222520222202205205000550500025

As you can see, these beautiful large squares contain solely the digits (1,2,7), (0,1,2), (4,5,6) and (0,2,5), respectively. Alas, they aren't the only such squares, see if you can find some more !  

Additional comments:  No one posted all four squares so three of them were left missing (though once you've created your multiprecision squaring program it's all too easy to use it to compute the squares and post them) and no one posted any more examples of this pattern so here you are, a few more:

      77470059130002034719700749² =

            6001610061606011616611060006010661000616100111161001

      942575429577943326987798² =

            888448440444044400080440444440408800048040888804


2. LOL the Second:  GCD

You may remember the GCD function. 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 ?


In the question above there's no mention of the GCD keyword from the JPC ROM, which is only available for the HP-71B and is limited to integer arguments, while this mini-challenge (which can be solved using almost any programmable scientific HP calc, from the HP-19C/29C upwards) deals with the GCD (Greatest Common Divisor) math function.

Mathematically, the GCD of two integers is the largest positive integer that divides each of them. But as it happens with the factorial, which can be extended to non-integer arguments via the \(\Gamma\) function, or with the Harmonic series, for which the digamma function does likewise, the GCD function can also be extended to non-integer arguments, and even to complex ones. This is done using this simple formula (which you can find in Wikipedia, no need to hunt obscure scholar papers for it):

      

where n has to be odd and ≥ 1 but m can be complex.

This 2-line, 89-byte HP-71B user-defined function accepts a complex argument K and a positive odd integer argument N and returns GCD(K,N) as per the above formula:

5  DEF FNG(K,N) @ COMPLEX P,R @ P=1 @ R=-2*PI*K/N @ FOR M=0 TO N-1
6  P=P*(1+RECT((1,R*M))) @ NEXT M @ FNG=LOG(P)/LOG(2) @ END DEF

We can check that it works fine by using this 4-line, 106-byte driver program to test it for 100 pairs of random integer arguments, excluding co-prime pairs from the listing. The first integer is passed as a complex value and the second integer is forced to be odd:

1  DESTROY ALL @ RANDOMIZE 1 @ RADIANS
2  FOR I=1 TO 100 @ A=INT(RND*50)+1 @ B=2*INT(RND*25)+1
3  G=IROUND(REPT(FNG((A,0),B))) @ IF G#1 THEN DISP USING "4X,4D,4D";A,B,G
4  NEXT I @ END

>RUN
       6   3      3
      38  19     19
      22  33     11
       9  21      3
      30  45     15
        ...

so it seems to be working Ok. Of course we can also call the UDF from the keyboard, e.g.:

>FNG((30,0),45)  ->  (15.000000002, 3.81114434334E-9), i.e. 15
>FNG((1/2,0),3)  ->  (1.79248125035, -2.26618007092)

Now, for real arguments we can use  cos x = ½·(e^ix + e^-ix)  to reduce the formula to the real product:

      

The user-defined function now becomes simpler (71 bytes, no comples variables or functions) and faster:

5  DEF FNG(N,M) @ P=1 @ R=PI*M/N @ FOR K=1 TO (N-1)/2
6  P=P*COS(R*K) @ NEXT K @ FNG=N+LOG2(P*P) @ END DEF

      Note: do not "optimize" LOG2(P*P) to 2*LOG2(P)

Let's apply it to this mini-challenge's question. First, let's evaluate GCD(15, x) for x in [4..5] at steps of 0.2:

>DESTROY ALL
>FOR I=4 TO 5 STEP .2 @ DISP USING "3X,D.D,2X,S2D.7D";I,FNG(15,I) @ NEXT I

  x    GCD(15, x) 
------------------
 4.0   +1.0000000  <- GCD(15,4) = 1
 4.2   +3.0660268
 4.4   +1.6288312
 4.6   +2.2305063
 4.8   +5.0211098
 5.0   +5.0000000  <- GCD(15,5) = 5

and a cursory inspection reveals that there are at least three solutions to the original question "For what value of x is  GCD(15, x) = 2 ?". Let's find one of them using this 1-line main program together with the UDF:

1  DESTROY ALL @ RADIANS @ STD @ DISP FNROOT(4,5,FNG(15,FVAR)-2)
5  DEF FNG(N,M) @ P=1 @ R=PI*M/N @ FOR K=1 TO (N-1)/2
6  P=P*COS(R*K) @ NEXT K @ FNG=N+LOG2(P*P) @ END DEF

>RUN
     x₃ =  4.59247568122

Now we can check that this value of x really makes  GCD(15, x) = 2 :

>FNG(15,RES)  ->  2

In the table above we can see that there are likely two other solutions for x, one in [4.0, 4.2] and another in [4.2, 4.4], namely:

FNROOT(4.0, 4.2 ...) ->  x₁ = 4.06191388928, FNG(15,RES) = 2.0000000002
FNROOT(4.2, 4.4 ...) ->  x₂ = 4.38207921034, FNG(15,RES) = 1.9999999996

and that completes all three solutions to this mini-challenge in the interval [4, 5], which you can see listed with 20 decimal digits in the comments below.

Additionally, they are duly located and labeled in my graph of  y = GCD(15, x)  right below, where you can also see the perhaps unexpected fact that  GCD(15, 4.5)  tends to -\( \infty \) because in this case the cosine product includes one cosine which is 0, thus so is the product itself and its log₂ tends to -\( \infty \). Notice theres's a fourth solution just outside the [4..5] interval:

      

Additional comments:  These are the three solutions computed using the HP-42S Solver via the 34-digit Free42 Decimal simulator:

01  LBL "GCD15"  13   x         01  LBL "LOL2"       
02  PI           14  DSE 00     02  MVAR "X"
03   x           15  GTO 00 ►   03  RCL "X"
04  15           16  X^2        04  XEQ "GCD15"
05   ÷           17  LOG        05   2
06   7           18   2         06   -
07  STO 00       19  LOG        07  END
08  SIGN         20   ÷
09 ►LBL 00       21  15
10  RCL 00       22   +
11  RCLx ST Z    23  END
12  COS

[SOLVER]  ->  Select Solve Program
Touch [LOL2], 4.0, [X], 5.0, [X], [X]  ->  4.59247568122 (...121 559541826...)
              4.0, [X], 4.2, [X], [X]  ->  4.06191388926 (...926 216056306...)
              4.2, [X], 4.4, [X], [X]  ->  4.38207921034 (...034 278186774...)

You can press [SHOW] after each solution to see up to 34 digits (at least 30 will be correct.) The fourth "solution" just outside the [4..5] interval is at x₄ = 5.22459905977 (222552325...).

---

That's all for now. I'll post the next two LOLs in a couple' days or so. Meanwhile, let's see your comments.

Regards.
V.