[VA] SRC #011 - April 1st, 2022 Bizarro Special

[Valentin Albillo]

 
Hi, all,

Thanks for your interest in my SRC #11 and your valuable posts, much appreciated. Now let's conclude Rd. Albizarro's story, featuring my Original Solution, results and comments:


Last time on Bizarro World:

Rd. Albizarro came back from a quick dinner and found the computed value of the gravitational force F already waiting in the HP-10C's display, which was indeed correct to the required three decimal digits. How was the feat accomplished ?

And now, the conclusion:

Rd. Albizarro was aware that the 79 bytes of RAM available in the HP-10C wouldn't allow for any deterministic cubature methods, as even just storing each of the six integration variables in their own register would require 42 bytes, leaving only 37 program steps for the method's implementation and no memory at all for anything else. Also, most deterministic cubature methods suffer from the so-called curse of dimensionality, requiring a number of sample points (and thus running time) exponentially growing with the dimension, to achieve even modest accuracy.

What cubature method would do, then ? Only a non-deterministic one would do, such as an OEM method (OLRAC ET NOM), which are free from this "curse" (though slowly converging) but still quite challenging to implement on the HP-10C because evaluating the integrand would require generating six uniform (pseudo-)random numbers per evaluation and the HP-10C's instruction set doesn't include a random number generator (RNG).

Thus, some suitable RNG would have to be implemented in RPN user code as a subroutine to be called 6 times per integrand's evaluation, a brilliant plan except for the fact that the HP-10C doesn't have subroutine capability either, which meant that the RNG code would have to appear in-lined six times.

To know how many program steps would be available for this, Rd. Albizarro first made a quick estimation of the number of storage registers needed: one for storing the seed; another to count down the number of times the program would loop through to compute the integral; another to store the number of samples, needed for obtaining the final value; and a fourth to keep the running summation of the integrand's evaluations, so 4 storage registers in all which would require 28 bytes.

That left 51 steps for the program, with the integrand's evaluation requiring at least 22 of them (assuming six 1-step (nonexistent!) RAN# instructions) and leaving only 29 available for everything else. However, implementing the nonexistent 1-step RAN# instruction in RPN meant that the ersatz RNG code would had to be (29+6)/6 ~ 5.8, i.e. 5 steps long at most.

Could a decent RNG be implemented in only 5 program steps or less ? Yes !.   Rd. Albizarro had recently read about some RNG for HP RPN models which had been advocated by one ugly-Earth native Mr. Albillo because of its extreme simplicity and fairly reasonable behavior, who then went to publish it in some magazine called "PPC Calculator Journal". This RNG could be implemented in as little as 4 program steps, namely:

    RCL (seed), R-D, FRAC, STO (seed)

where R-D was a radians-to-degrees conversion, which most fortunately (thanks, Yhprum !) did exist in the HP-10C under the name ->DEG, so it could be used to generate the six random numbers, requiring 24 program steps in all. However, that would still leave only 11 steps for everything else, which included some initialization, updating the ongoing summation, checking for loop termination, computing the final average, ...

Fortunately, after a modicum of further reflection, Rd. Albizarro discovered that each pair of random numbers could be generated using just 7 steps instead of 8, resulting in 14 steps being available for the remaining operations, which was more than enough so the HP-10C program was indeed feasible and could be implemented in just 49 program steps, like this:

   01  STO 1    11  STO 3    21  ENTER    31  ENTER    41  /
   02  STO 2    12  -        22  ->DEG    32  ->DEG    42  STO+ 0
   03  CLX      13  1        23  FRAC     33  FRAC     43  RCL 1
   04  STO 0    14  STO- 1   24  STO 3    34  STO 3    44  X=0?
   05  RCL 3    15  +        25  -        35  -        45  GTO 47
   06  ->DEG    16  ENTER    26  X^2      36  X^2      46  GTO 05
   07  FRAC     17  X^2      27  +        37  +        47  RCL 0
   08  ENTER    18  RCL 3    28  RCL 3    38  ENTER    48  RCL 2
   09  ->DEG    19  ->DEG    29  ->DEG    39  SQRT     49  /
   10  FRAC     20  FRAC     30  FRAC     40  x          (GTO 00, default end of program)

Before running it, the initial seed had to be stored in R3 and the number of samples (pairs of points to generate and use in the integrand's evaluation) had to be specified, which Rd. Albizarro heuristically estimated as follows:

A good way to test a program is to run it against problems whose solutions are known. In this case, the gravitational force F between two contacting spheres (instead of cubes) is:


G = 1, m₁ = 1, m₂ = 1, d = 1, F = G*m₁*m₂/d² = 1.000,

so a version of this program particularized for spheres instead of cubes was executed, watching for the number of samples needed to achieve 3-digit accuracy, which was 687 samples, returning F ~ 0.999.

However, the volume of a sphere ("ball" would be more correct but whatever) of diameter 1 is 4/3*Pi*(1/2)³ = Pi/6 = 0.524+, significantly smaller than the volume of a cube of side 1, which is exactly 1, thus to maintain the same sample density (in order to achieve a comparable accuracy) the number of samples must be multiplied times the volumes' ratio, 6/Pi, so the estimated number of samples for the cubes case would be 687*6/Pi ~ 1,312 samples.

So, Rd. Albizarro stored some suitable seed (say, 1) in register R3, entered the number of samples to use (1312) in the display, set [FIX 3] and executed the program, like this:

    1 [STO 3] 1312 [R/S] -> 0.925  (0.925 4711044)

which returned the sought-after gravitational force F = 0.925 after ~ 72', nicely spent having a tasty quick dinner:

   

The theoretically correct 3-digit result is 0.926 (0.925 9812605, to 10 correct digits) so indeed 3 correct digits (save 1 ulp) were obtained, as required.

Having produced the desired result and met the strict deadline, the mission was fully accomplished but after some months had elapsed a much more relaxed Rd. Albizarro leisurely pondered whether the sextuple integral could perhaps be tackled symbolically, and after a few days finally succeeded in reducing it to a triple integral first, then to a double integral, and finally to a single-dimensional definite integral, which once evaluated resulted in this nice, exact symbolic value:

   

and after some rearranging, this longish expression would also exactly fit as a 79-step HP-10C program with no inputs required and no registers used (though alas, none were available as the program uses up all 79 bytes of RAM), like this:

   01   PI    17  x     33   √      49  0     65   LN
       3         -         1         x         6
       /         6         +         +         x
       2         √         LN        5         +
       LN        -         2         √         6
       +         LASTX     x         1         √
       2         4         2         +         2
       6         +         √         LN        x
       x         5         1         3         TAN^-1
       2         x         +         5         2
       √         LN        LN        x         2 
       7         -         +         -         x
       -         2         5         6         -
       3         x         √         √         3
       √         +         +         1     79   /
       2         3         1         +

Assuming RADians mode: [R/S] -> 0.9259812667 (exact: 0.9259812605...) in a few seconds, correct to 8 decimal digits, the last two being lost to rounding errors throughout.

Nothing else to do here, so as Superman's imperfect duplicate Bizarro #1 would say ...

                     

---

Well, this concludes Rd. Albizarro's story, and I still have a number of hopefully interesting Comments ready to post as an Epilogue of sorts, including some real-life applications, but first let's hear from you.

Regrettably, nobody posted a working HP-10C program and/or the sextuple integral's value correct to three digits, as required (although J-F Garnier came pretty close on both counts, many thanks for your efforts and great posts, J-F), but if any of you can produce working code for any other HP model (say, the HP-71B + Math ROM or some RPL model), with or without using symbolic manipulations (say, dimensionality reduction) and/or have some comments of your own, now that a correct solution has been provided and thus there's no spoiling the challenge for anyone, this is the time to post them and yes, I'm looking at you !

Regards.
V.