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

[Valentin Albillo]

  
Hi, all,
                
Well, my HP calc site reached the 100,000 downloads mark a few days ago and on top of that today it's April, 1st, so let's celebrate the occasion ... (drum roll) ...

Welcome to my SRC #11 - April 1st, 2022 Bizarro Special

intended to once again put your brains and your HP calculators to work, this time featuring a bizarre challenge which nevertheless has important, useful real-life applications. But first, some exposition (Note: What follows are mostly my own ramblings © me, DC Comics are not to blame !):

The seasoned veterans among you might fondly remember the classic Superman comics of the '60s, and that within the Superman universe there was the Bizarro World (also known as Htrae), a fictional planet home to the eponymous people, whose society is ruled by the Bizarro Code, which states:

"Us do opposite of all Earthly things !  Us hate beauty!  Us love ugliness !  Is big crime to make anything perfect on Bizarro World !"

Originally a normal planet, the Bizarro World is now cube-shaped because the Bizarros couldn't stand the perfection of its spherical shape and Bizarroformed it like this:

In time, the Bizarro society thrived to the point where their cubic-shaped planet was quickly becoming overpopulated, so they set upon themselves the huge task of attaching face to face a second almost-identical planet to their original one, thus nearly duplicating the habitable surface and best of all, augmenting the imperfection of their homeworld's shape relative to perfect sphericity by making it prismatic, like this:

However, they first needed to make sure that the project was feasible, in particular that the gravitational force F between the planets when they were in contact would be manageable. Therefore, they commissioned Lebon Prize laureate scientist Rd. Nitnelav Albizarro #1 to carry out the computation, who immediately set to the task of finding out F when their respective centers were initially separated by an arbitrary distance d, which would then be shortened until the planets were in contact face to face.

To simplify matters and as the results could be easily rescaled afterwards, it was assumed that the gravitational constant G was 1 (in some units) and the planets were homogeneous cubes of side 1 (ditto) and mass 1 (ditto), initially placed like this:

Now, Rd. Albizarro considered a pair of sample points, (x₁, y₁, z₁) in Htrae 1 and (x₂, y₂, z₂) in Htrae 2, knowing that their contribution to the overall force would be:

and duly taking into account the symmetry (which nullifies the force exerted in the y and z directions) and integrating over all possible values for the respective point coordinates, quickly got this sextuple integral for the value of the force F between the planets:

to be numerically evaluated for the particular case when both planets are in contact, i.e. the distance d between their centers is 1.

As it happened, Rd. Albizarro was a Code-abiding citizen of the Bizarro society and wouldn't go for an unnecessarily accurate (ugh!) result but would instead be satisfied with obtaining about three correct digits. To that effect, and aware of the need to use some computing device in order to meet the assigned deadline, Rd. Albizarro promptly proceeded to check out the ones available in ugly human Earth at the time (ca. 1982) and saw that computers were so perfectly fit for the job that using one would be disgusting and could result in being severely reprimanded or worse, so the alternative was to get instead a programmable calculator (progcalc for short) as non-programmable ones certainly wouldn't do.

Restricting thus the search to progcalcs, it was soon apparent that there were essentially two main contenders, the ones branded TI (probably standing for "Totally Ideal", ugh !), and the ones branded HP (possibly standing for "Hardly Perfect"), the latter being just what was needed !

Now there was the question of selecting which particular HP model to use but that was easy-peasy, just a matter of looking at the specs and chosing the least capable one, which happened to be the HP-10C, the proverbial runt of the litter, a severely limited model having only 79 bytes of RAM available for storing programs and data, no subroutines, no flags, no indirection, no loop intructions, just two conditional tests, a meager function set, no I/O, and very slow to boot ... indeed, the least perfect model for the task at hand.

But Rd. Albizarro was unfazed, thinking that "If ugly Earth's Newton could do his gravitational calculations using this thing, it'll do for me as well", and against all odds actually succeeded by quickly writing a clumsy RPN (Really Perfect Not) program for the HP-10C, keying in the pertinent inputs, pressing the [R/S] (Run Slowly) key, going out to have a quick dinner and lo and behold, upon returning the computed value of F was already waiting in the display, which indeed was correct to three decimal places, as required:

---

Well, once the story's been told, it's your move:

Try to emulate Rd. Albizarro's achievement and write a program for the HP-10C which computes in a reasonable time a numerical value for the above sextuple integral correct to at least ~ three decimal places.

That failing (shame on you !), see if you can do it using other more perfect (ugh!) models, preferably HP and preferably vintage.

Note:  No cheating whatsoever allowed. Also, doing symbolic manipulations (transformations, dimensionality reduction, changes of variables), either by hand or using any CAS and/or giving math lectures is a big crime, you'll get arrested or worse, so stick to purely numerical computations. This is my last SRC for a looong while so don't spoil it for me, Ok ? Thanks.

I'll post my Original Solution for the HP-10C with results and extensive comments within a few days ... Or not, after all this might be an elaborate April Fools' Day practical joke !!   

V.

[Massimo Gnerucci]

Just to say that I love this one: Bizarro... cubic worlds... all the silliness of my youth's Superman!

[Ren]

I remember some of the Bizarro storyline, but I couldn't afford to buy comic books back then, so could only read what came my way by those who loaned comics to my older brothers.

I think a parallel can be seen in Star Trek: Lower Decks interactions with the Paklid.

[J-F Garnier]

A sextuple integral, that even a HP-71B Math ROM couldn't handle, and on a HP-10C, no less. That's a challenge!

As a very crude and ugly approximation, we can put all the mass at the center of each cube, and calculate the force F = Mm/D²
M=m=D=1 so F=1.

To improve this first approximation, we can then divide each cube in 8 smaller cubes of side 1/2, and calculate the force between each pair of small cubes (putting again the mass of each at their center), and sum them up.
Then we can divide the cubes further in 27, 64, etc parts and hope that this all converges to some value.

This would be with no problem, in principle, with enough computing resources, but here we have only a 10C (or another machine but using only 10C features and resources) so this seems impractical.

Another approach would be to use a statistical method by taking random points in each cube and calculating their contribution to the force.
With enough randomly distributed points, we can have a good approximation of the force, and the implementation may be easier.
We need to take 6 random values for the x1, x2, y1, y2, z1, z2 values, calculate the integrand and sum it.
... and here we realize that the 10C has no random function.
We could use the well-known user-code random generator, but we also realize that the 10C has no subroutine, duplicating the random generator code 6 times is just not possible.

However, we don't need a perfect random generator,
So it appears that the challenge is to find a poor-man random generator, compact enough to fit 6 times in the 10C program memory, yet not too bad.
Any tentative?

J-F

[ijabbott]

(04-02-2022 09:32 AM)J-F Garnier Wrote:  A sextuple integral, that even a HP-71B Math ROM couldn't handle, and on a HP-10C, no less. That's a challenge!

As a very crude and ugly approximation, we can put all the mass at the center of each cube, and calculate the force F = Mm/D²
M=m=D=1 so F=1.

TBH, I don't get why that is not the correct answer!

[Valentin Albillo]

(04-02-2022 08:31 PM)ijabbott Wrote:  

(04-02-2022 09:32 AM)J-F Garnier Wrote:  As a very crude and ugly approximation, we can put all the mass at the center of each cube, and calculate the force F = Mm/D², M=m=D=1 so F=1.

TBH, I don't get why that is not the correct answer!

Because a cube doesn't have spherical symmetry. cf. Wikipedia:

"[...] an object with a spherically symmetric distribution of mass exerts the same gravitational attraction on external bodies as if all the object's mass were concentrated at a point at its center. (This is not generally true for non-spherically-symmetrical bodies.)"

V.

[rawi]

Garnier wrote:

Quote:However, we don't need a perfect random generator,
So it appears that the challenge is to find a poor-man random generator, compact enough to fit 6 times in the 10C program memory, yet not too bad.
Any tentative?

I do not think that you need the random number generator 6 times. Why not generate one random number, store it in a register and use it six times?

[vaklaff]

(04-03-2022 06:21 AM)rawi Wrote:  I do not think that you need the random number generator 6 times. Why not generate one random number, store it in a register and use it six times?

Or just use 17 which is the best random anyway :-)

[ijabbott]

(04-03-2022 01:25 AM)Valentin Albillo Wrote:  

(04-02-2022 08:31 PM)ijabbott Wrote:  TBH, I don't get why that is not the correct answer!

Because a cube doesn't have spherical symmetry. cf. Wikipedia:

"[...] an object with a spherically symmetric distribution of mass exerts the same gravitational attraction on external bodies as if all the object's mass were concentrated at a point at its center. (This is not generally true for non-spherically-symmetrical bodies.)"

V.

Indeed, splitting each cube into 8 cubelets and summing the gravitation forces between each pair of cubelets (one from Htrae 1 and one from Htrae 2) using the point mass assumption produces an answer that differs from 1 (actually about 1.09217), which means that cubes of uniform density cannot be replaced with point masses.

EDIT: As pointed out to me by Albert Chan, I forgot to resolve the vectors, so my sum is wrong. It should be less than 1.

[rprosperi]

Though I don't have the time, and probably am not up to the analysis involved here, I wanted to add some comments about this entertaining and thought-provoking post:

• Thanks to Valentin for the huge effort and true creativity to provide this great post!
• This is a typical post from Valentin in the sense that, even if I am not able to contribute a possible solution, I am entertained but I also learn a lot, simply by reading the ensuing discussions, and I want to thank Valentin for that. I hope this is not the last SRC for some time...
• I think Bizarro was found in Heavy Metal, wasn't it? A true treasure for many years, it inspired so many other authors
• Has anyone else noted the important word in low-contrast, barely-readable white font near the bottom...
  

[J-F Garnier]

Below is my attempt to implement the sextuple integral computation on the 10C.

Unfortunately, I didn't succeed to fit into the 10C memory, so this is a 41C version but using only the instruction set of the 10C (except the LBLs needed by the 41C).
The programs expects the number of random samples in X as an input.

I used a simple random number generator of mine (unless I knew it from past readings, can't say). The code is documented and easy to understand.

Now, the results. I run it with different numbers of random samples, between 20 and 999:
#samples integral
20 0.628
50 0.823
100 0.898
200 0.910
500 0.981
999 0.967

So I hardy get 1 correct digit if I round all results in FIX 0.
I'm even not sure if the integral is smaller or larger than 1.
Is it really a surprise? Even with 999 samples, this means an average of about 3 samples per dimension of each cube.

But probably Rd. Albizarro preferred a wrong result with 3 figures rather than one reliable estimation that would be big crime !

At the end, did it make any difference for Bizarro World ?

J-F


01*LBL "SRC11"
STO 02 STO 01 ; #samples
PI STO 00 ; rnd seed
0 STO 03 ; init sum

08*LBL 00
RCL 00 PI * FRC STO 00 ; z2
RCL 00 PI * FRC STO 00 ; z1
- X^2 ; (z2-z1)²
RCL 00 PI * FRC STO 00 ; y2
RCL 00 PI * FRC STO 00 ; y1
- X^2 + ; (y2-y1)²+(z2-z1)²
RCL 00 PI * FRC STO 00
1 + ; x2
RCL 00 PI * FRC STO 00 ; x1
- X^2 ; (x2-x1)²
+ ; d²
LASTX SQRT X<>Y ; (x2-x1) d²
ENTER^ SQRT * ; d^3
/ ; (x2-x1)/d^3
ST+ 03 ; add to sum
1 ST- 02 ; decr counter
RCL 02 X<>Y X<=Y? GTO 00 ; loop

63*LBL 01
RCL 03 RCL 01 / ; result
END

[Valentin Albillo]

.
Hi, all,

Thanks for your interest in my SRC#11 and most definitely for your posts. I'll give here some comments. Let's begin ...

Massimo Gnerucci Wrote:Just to say that I love this one: Bizarro... cubic worlds... all the silliness of my youth's Superman!

I'm very glad for your appreciation of the Bizarro-themed story and thanks for telling me, I've always appreciated the concept myself so this was a nice opportunity to include it in this bookend SRC.

Ren Wrote:I think a parallel can be seen in Star Trek: Lower Decks interactions with the Paklid.

Thanks for the tip. I know nothing about the plethora of new ST series, I only watched the Original Series, Deep Space 9, Voyager and Enterprise (as well as all the movies), so no idea what Lower Decks is about or the aliens (I suppose they're aliens) you mention, sorry.

rawi Wrote:

J-F Garnier Wrote:We need to take 6 random values for the x1, x2, y1, y2, z1, z2 values, calculate the integrand and sum it.

I do not think that you need the random number generator 6 times. Why not generate one random number, store it in a register and use it six times?

Because if you do as J-F says, i.e. "calculate the integrand", you'll find that if all the variables use the same random number you'll get a nice but useless 0 divided by 0 value for the integrand and there's not much you or anyone else can do with that.

rprosperi Wrote:Thanks to Valentin for the huge effort and true creativity to provide this great post!

Thank you very much for your continued appreciation of my productions, Bob, you're always too kind but the thing I appreciate the most is that you do take the trouble and time to post it and let me know. Thanks again.

J-F Garnier Wrote:Now, the results. I run it with different numbers of random samples, between 20 and 999 [...]

I'm curious ... why 999 instead of 1,000 ? I've inspected your code and see no reason for stopping short of the more natural 1,000. You're not using ISG or something like that.

J-F Garnier Wrote:But probably Rd. Albizarro preferred a wrong result with 3 figures rather than one reliable estimation that would be big crime !

Nope, Rd. Albizarro sticks to Bizarro society Code but Mathematics is neither perfect nor imperfect, it just is.

J-F Garnier Wrote:At the end, did it make any difference for Bizarro World ?

Of course it did, they got 10 habitable faces where previously they only had six

I'll post the rest of Rd. Albizarro's story (which includes my Original Solution and extensive comments) next Saturday 9 so you've got 5 additional days to try and come up with your own Solution ... or at least die trying !

Best regards.
V.

[J-F Garnier]

(04-04-2022 10:34 PM)Valentin Albillo Wrote:  

rawi Wrote:I do not think that you need the random number generator 6 times. Why not generate one random number, store it in a register and use it six times?

Because if you do as J-F says, i.e. "calculate the integrand", you'll find that if all the variables use the same random number you'll get a nice but useless 0 divided by 0 value for the integrand and there's not much you or anyone else can do with that.

No, it's not as severe. If you use the same random value, and since x2=1+rnd you get exactly 1 for the integrand, and 1 as the result !

Quote:

J-F Garnier Wrote:Now, the results. I run it with different numbers of random samples, between 20 and 999 [...]

I'm curious ... why 999 instead of 1,000 ?

I was expecting someone would ask... and you did!
It's just that, as you may know, I'm not very good in post formatting, so I used 999 instead of 1000 to keep the figures well aligned :-)

J-F

[Valentin Albillo]

(04-05-2022 07:23 AM)J-F Garnier Wrote:  

(04-04-2022 10:34 PM)Valentin Albillo Wrote:  

Because if you do as J-F says, i.e. "calculate the integrand", you'll find that if all the variables use the same random number you'll get a nice but useless 0 divided by 0 value for the integrand and there's not much you or anyone else can do with that.

No, it's not as severe. If you use the same random value, and since x2=1+rnd you get exactly 1 for the integrand, and 1 as the result !

Sorry but it doesn't fly. That x2=1+rnd doesn't appear explicitly in the integrand expression as such, only x2-x1 in both numerator and denominator, and further it only appears when using your stochastic method, not if using other possible methods.

Thus, I don't think that rawi, when trying to compute the value of the integrand as it appears in my OP, would think or know about the need to add 1 to x2, which is only implicit in the [d, d+1] limits of integration and only when substituting d=1, thus he'd surely get 0 divided by 0 as I said he would.

J-F Garnier Wrote:

Valentin Albillo Wrote:

I'm curious ... why 999 instead of 1,000 ?

I was expecting someone would ask... and you did!
It's just that, as you may know, I'm not very good in post formatting, so I used 999 instead of 1000 to keep the figures well aligned :-)

It doesn't fly either, you could have used 1E3, which aligns properly with 100, 200 and 500, and furthermore you include 20 and 50 which don't align with them either and you didn't seem to mind.

Best regards.
V.

[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.

[Albert Chan]

It would be nice to have back of envelope calculation for F

We can keep planet intact, and split the other planet, to 2x2x2 = 8 cubes.
With symmetry, only 1 corner to consider

(3/4, 1/4, 1/4) → D1^2 = 11/16, cos(θ1) = 3/√11
(5/4, 1/4, 1/4) → D2^2 = 27/16, cos(θ2) = 5/√27

Force is proportional to 1/D^2
Force is an vector, we need cos(θ) correction, to sum only force between planets.

F = ((16/11)*(3/√11) + (16/27)*(5/√27)) / 2 ≈ 0.94295

Since 1 planet still intact, we would expect F to be even less, i.e. F < 0.94295

[Albert Chan]

Difference of two uniform distribution form triangular distribution.
This reduced sextuple to triple integral.

We just integrate ∫∫∫(f(x,y,z) dz dy dx), with all variables from 0 to 1

Let x = |x2-x1|, y = |y2-y1|, z = |z2-z1|, all with distribution ◣
Squaring is going to remove abs() anyway.

For x, we do need ◢◣ distribution, center at d, thus (d+x), (d-x) terms.
Because of the discontinuity at x=0, we fold the 2 integrals as 1.
This moved the discontinuity to the edge.

Here is HP71B equivlent code (smaller P actually make things worse)

10 P=.01 @ D=1 @ T=TIME
20 DEF FNZ(X,Y)=INTEGRAL(0,1,P,X/(X*X+Y*Y+(1-SQR(IVAR))^2)^1.5)
30 DEF FNY(X)=INTEGRAL(0,1,P,FNZ(X,1-SQR(IVAR)))
40 DEF FNX(X)=FNY(D+X)+FNY(D-X)
50 DISP INTEGRAL(0,1,P,FNX(1-SQR(IVAR)))/2, IBOUND, TIME-T

>RUN
 .92596918938      1.84480545818E-2      21.66

[Albert Chan]

This version improve accuracy by removing square roots. Let x=u^2, v=1-u

∫(f(1-√x), x=0..1) = ∫(f(1-u)*(2u), u=0..1) = 2*∫(f(v)*(1-v), v=0..1)

Transformed integral (with discontinuity at x=0, evaluate in 2 pieces)

\(\displaystyle I = 4 \int_{-1}^{1} \int_0^1 \int_0^1
\frac{(d+x)\,(1-|x|)\,(1-y)\,(1-z)}
{[(d+x)^2 + y^2 + z^2]^{3/2}}\;dz\;dy\;dx \)

10 P=.0001 @ D=1 @ T=TIME
20 DEF FNZ(X,Y)=INTEGRAL(0,1,P,X/(X*X+Y*Y+IVAR*IVAR)^1.5*(1-IVAR))
30 DEF FNY(X)=INTEGRAL(0,1,P,FNZ(X,IVAR)*(1-IVAR))
40 I1=INTEGRAL(0,1,P,FNY(D+IVAR)*(1-IVAR))*4 @ DISP I1,IBOUND*4,TIME-T
50 I2=INTEGRAL(0,1,P,FNY(D-IVAR)*(1-IVAR))*4 @ DISP I2,IBOUND*4,TIME-T
60 DISP I1+I2

>RUN
 .236198653356      2.35758565639E-5      128.35
 .689787720788      6.90121464808E-5      763.87
 .925986374144

1st integral is sum of force for mass separated by D or more.
2nd integral is sum of force for all mass separated by D or less.

2nd integral (for D=1) contributed about 3/4 of full force.
Its singularity (where planets are touching) make it harder to calculate.

[Albert Chan]

(04-09-2022 09:41 PM)Valentin Albillo Wrote:  

Trivia, from my thread, SOHCAHTOA, for arc-trig

atan(√(24/1))             // TOA, O=24, A=1
= acos(√(1/25))         // CAH, H = O+A = 25
= acos(1/5)

[Valentin Albillo]

(04-10-2022 12:37 AM)Albert Chan Wrote:  

(04-09-2022 09:41 PM)Valentin Albillo Wrote:  

Trivia, from my thread, SOHCAHTOA, for arc-trig

atan(√(24/1))             // TOA, O=24, A=1
= acos(√(1/25))         // CAH, H = O+A = 25
= acos(1/5)

That shaves off 2 steps from the program which evaluates the exact F's expression, leaving it at 77 steps (for now ...). Well done !

Regards.
V.