[VA] SRC #012e - Then and Now: Roots

[Valentin Albillo]

  
Hi, all,

Well, today is Feb 7 and my 888^th post here, so Welcome to the 5^th part (next-to-last) of my ongoing SRC #012 - Then and Now, where I'm showing that advanced vintage HP calcs which were great problem-solvers back THEN in the 80's are NOW still perfectly capable of solving recent, highly non-trivial problems intended to be tackled using modern PCs, never mind ancient calcs.

In the next weeks I'm proposing six increasingly harder such problems for you to try and solve using your vintage HP calcs while abiding by the mandatory rules summarized here:

You must use VINTAGE HP CALCS (physical/virtual,) coding in either RPN/Mcode, RPL/SysRPL or HP-71B BASIC/FORTH/Assembler, and NO CODE PANELS.

Caveat: the reason for this rule is obvious: this Problem 5 has been tailored for and is challenging if and only if you attempt to solve it using a vintage calc. Tackling it using modern software (e.g. Mathematica, etc.) is meaningless and just succeeds in spoiling the fun, so please don't do it. Now back to business ..


Once P1 (Probability), P2 (Root), P3 (Sum) and P4 (Area) are over, now's the turn for P5 (Roots) which revisits rootfinding, only this time the function whose roots are to be found has its own peculiarities and the final result is both wholly unexpected and utterly amazing ! ... But first a relevant, hopefully interesting preamble:

Preamble

Some of you might remember an article where I computed an approximation to the ∏(x) function, which gives the number of primes up to some given limit (say, the number of primes up to 1,000, which is 168) by using an equivalent form of the R(x) function, which is canonically defined as
      
and μ(k) is the Möbius function. However, this form of R(x) is not convenient for computations because both li(x) and μ(k) are time-consuming to evaluate and the convergence is poor, so I used instead a more amenable, equivalent form valid for x > 0, namely:
      
Using this form we easily get R(10³) ~ 168 primes up to 1,000 (err = 0%), R(10⁶) ~ 78,527 primes up to one million (err ~ 0.037%) or R(10⁹) ~ 50,847,455 primes up to one billion (err ~ 0.00016%). Indeed, R(x) is an extremely good approximation, as seen in the figure below (left image) which compares the graphs of R(x) and ∏(x):



Now you may be wondering where the "rootfinding" fits in all this. If we look at the graphics above (notice the X-axis' logarithmic scale in the zoomed image), R(x) seems to be always positive for x > 0, as it certainly looks like it will never cross the X axis and thus will have no positive real roots at all ... but quoting Pink Floyd: "things are not what they seem", and so we have:


Problem 5:  Roots

Write a program to compute and output the 7 largest positive real roots of R(x), in decreasing order.

You can use any equivalent form of R(x) that suits you best, be it the canonical definition or the one I used in my article and in the examples here, or any other you deem appropriate.

Your program should have no inputs and must compute and output the 7 roots and end. You should strive to get at least 7 correct digits (give or take a few ulp) for each of the 7 roots and the faster the running time the better.

Note: The main goal here is to compute the 7 roots, but as for intermediate values (e.g. values of special functions,) you can choose whether to compute them on the fly or else to get them from references and include them as DATA statements or whatever. The former approach will result in a smaller but slower program while the latter will be faster but bigger, or you can mix both approaches as you see fit. Your choice.

If I see interest I'll post my own original solution for the HP-71B, a 15-line program which automatically does the job, plus comments. Meanwhile, let's see your very own clever solutions AND please remember the above rules.

V.