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

[Valentin Albillo]

      
Hi, all,

Caveat lector:  this post is VERY long, so you might want to read it only when you can allocate enough free time for it.

Well, almost three weeks have elapsed since I posted Problem 5 and, as expected, its difficulty and advanced math subject matter resulted in very few posts and just a single solution (*) other than my original one. Back in ye goode olde forum, problems like this would've been made short work of ... Ah, those were the days ! ...

Again, no RPL (or RPN ) solutions at all, and this time the tireless contributors were Fernando del Rey and J-F Garnier (*), plus additional posts in a parallel thread by EdS2, J-F Garnier, Fernando del Rey and PeterP. Thank you very much to all of you for your interest and valuable contributions.

Now, this is my detailed sleuthing process and resulting original solution, plus additional comments:


My sleuthing process

First of all, we must decide which equivalent formula for R(x) is the most suitable for rootfinding. The canonical formula, namely:
      
is absolutely useless for the purpose, as μ(k), the Möbius function, is hard to compute for large integer k, while li(x) is time-consuming to evaluate accurately and last but not least, the convergence of the summation is extremely slow, thus this canonical formula is out of the question.

As for the second formula, aka the Gram series, which I used in my article, namely      
it does converge very nicely for large arguments x but not so much for arguments x << 1 because then ln(x) is large and negative and this causes the summation's terms to oscillate widely between positive and negative while initially increasing in magnitude enormously, as can be seen by executing this code from the command prompt (after running my solution program below at least once, to initialize all code and data):

>INPUT X @ SCI 4 @ FOR M=1 TO 10 @ LN(X)^M/(M*FACT(M)*FNZ(M+1)); @ NEXT M @@ STD

    X         Term 1   Term 2    Term 3   Term 4     Term 5    ...  Term 9    Term 10
--------------------------------------------------------------------------------------
? 1E-12  -> -1.6798E1 1.5878E2 -1.0828E3 5.8556E3  -2.6386E4 ...  -2.8717E6  7.1448E6
? 1E-100 -> -1.3998E2 1.1027E4 -6.2664E5 2.8239E7  -1.0604E9  ... -5.5655E14 1.1539E16
? 1E-499 -> -6.9850E2 2.7457E5 -7.7861E7 1.7508E10 -3.2807E12 ... -1.0676E21 1.1045E23

The terms ultimately tend to zero as the factorial in the denominator dominates over the powers of the (large) negative logarithm in the numerator, but adding and subtracting the initial very large terms to eventually get a very small result incurs in enormous rounding errors and/or would require heavy multiprecision computations, so ultimately Gram series is no good either.

What to do ?  We must search for another, better suited equivalent formula, that's what, and we can do it in just three logical steps, as follows:

1) In my OP, I mention an old article of mine which includes a prime-counting function, so a quick, easy search at my website section HP Calculator Articles, reveals it to be HP Article VA027 - Small Fry - Primes A'counting.

2) In the article, the FNR prime-counting function uses FNZ, which is identified there as the Riemann Zeta function, so doing a Google search on "Prime counting function Riemann zeta function" brings as one of the very first hits the link "Prime-counting function" - Wikipedia.

3) Clicking on that link, we find under the section Formulas for prime-counting functions the paragraph:

"Folkmar Bornemann proved [...] that
      
where ρ runs over the non-trivial zeros of the Riemann zeta function and t>0."

which includes a formula that's also equivalent to the previous ones but whose argument is e^{-2 ∏ t} and consisting of two infinite summations in terms of t, valid for t > 0. This means that t essentially acts as a logarithm and so the formula can be used to evaluate R(x) for incredibly small values of x, as needed.

Note: I noticed that the second summation, where ρ runs over all non-trivial complex zeros of ζ (which so far appear in conjugate pairs,) can be shortened to use just one zero from each pair, as the imaginary parts of both terms would cancel out, while the real parts are identical and would double in value, thus cancelling the 1/2 factor. This cuts the number of terms in half (thus doubling the speed) and the second summation gets simplified to this (the big R stands for real part and the big I stands for imaginary part):
            

Once we have this formula, solving Problem 5 is now just a matter of implementing it in any of the languages available in HP vintage calcs, which might seem difficult at first sight but it merely takes some work to either evaluate the special functions appearing in the formula or getting their values from references, without computation.

The special functions and values appearing in the formula are Riemann's ζ function for integer arguments, its non-trivial complex zeros and the complex values of its first derivative ζ' at those zeros. After a little experimentation, it turns out that to achieve 12-digit accuracy only the real values ζ(2) to ζ(8), the non-trivial complex zeros ρ₁ - ρ₆ of ζ and the corresponding complex values of ζ' at those zeros are needed and can be taken from references. Assuming your vintage HP calc has decent complex-math capabilities, the rest is easily computed.


My original solution

My original solution is this 15-line, 897-byte HP-71B program:   

 10 DESTROY ALL @ OPTION BASE 1 @ DIM Z(8),Z0(6),K,L,P,R,S,T,U,V,W,Y @ COMPLEX Z1(6),C

 20 DATA 0,PI^2/6,1.20205690316,PI^4/90,1.03692775514,PI^6/945,1.00834927738,PI^8/9450
 30 DATA 14.1347251417,21.0220396388,25.0108575801,30.4248761259,32.9350615877,37.5861781​588
 40 DATA (.783296511867,.124699829748),(1.10929556346,-.248729788516)
 50 DATA (1.29579560501,.450036709438),(1.12013084524,-.667509469349)
 60 DATA (1.16057006749,.750554150342),(1.85346624998,-.561004420496)

 70 READ Z,Z0,Z1 @ A=2 @ M=SGN(FNR(A)) @ FOR B=500 TO 20000 STEP 500 @ N=SGN(FNR(B))
 80 IF N#M THEN FIX 8 @ DISP "t:";FNROOT(A,B,FNR(FVAR));"-> X=";FNX$(RES,8);", R=";FVALUE @ M=N
 90 A=B @ NEXT B

100 DEF FNR(T) @ S=0 @ K=0 @ REPEAT @ V=S @ L=2*K+3 @ S=S+(-1)^K*T^(-L)/(L*FNZ(L))
110 K=K+1 @ UNTIL S=V @ R=S/PI @ S=0 @ K=0 @ REPEAT @ K=K+1 @ C=(.5,Z0(K)) @ V=S
120 S=S+REPT(T^(-C)/(C*COS(PI/2*C)*Z1(K))) @ UNTIL S=V @ FNR=R+S

130 DEF FNZ(N) @ IF N<9 THEN FNZ=Z(N) @ END ELSE IF N>37 THEN FNZ=1 @ END
140 Y=0 @ P=1 @ REPEAT @ P=P+1 @ W=Y @ Y=Y+P^(-N) @ UNTIL Y=W @ FNZ=Y+1

150 DEF FNX$(T,N) @ STD @ T=-2*PI*T/LN(10) @ FNX$=STR$(10^(FP(T)+1))[1,N]&"E"&STR$(INT(T))

Note: keywords FNROOT, FVAR, FVALUE, COMPLEX and REPT are from the Math ROM; REPEAT and UNTIL are from the JPC ROM.


Line 10  performs some initialization and defines all arrays and certain variables. In particular, DIM ...,K,L,...,W,Y and COMPLEX ...,C are necessary to create those variables here and not inside the DEF function definitions because of a known system bug. See this post for further details.

Lines 20-60  define all necessary data, which includes an initial dummy 0 and the values of ζ(2) .. ζ(8) in line 20, the imaginary parts of the first six zeros of ζ(z) in line 30, and the six complex values of ζ'(z) at those zeros in lines 40-60.

Lines 70-90  read all data and compute and output the seven roots x₁ .. x₇, as well as the corresponding t parameters and resulting values of R(x) at the computed roots, which are suitably near 0.

The roots are located by sweeping the interval t = [2, 500 .. 20000] using steps of 500 and when a sign change is detected, FNROOT is used to accurately compute the root using the subinterval's extremes as initial approximations, for fast convergence.

Lines 100-120  define R(x), computing each summation until the sums stop changing, then returns the final sum of both. The first summation calls FNZ to compute ζ(n), while the second summation uses the complex values of ζ'(z) from the array where they were stored and various complex functions and complex arithmetic operations. FNR can be called from the command prompt after the program is run.

Lines 130-140  define ζ(n) for integer n ≥ 2, returning the real values from the array where they were stored for 2 ≤ n ≤ 8, the constant 1 for n > 37, and otherwise it computes the summation until it stops changing and returns the final sum. FNZ can be called from the command prompt after the program is run and will quickly return values accurate to 12 digits for all integer arguments n ≥ 2.

Line 150  defines an utility string-valued function which converts the value of the parameter t to the corresponding value of x and returns it as an N-character mantissa plus its corresponding exponent, which can be >> 499.

Notes:

○  There's no need to specify RADIANS in the initialization because all complex trig functions always use radians regardless on the current angular mode.

○   It is possible to save 30 bytes of RAM and slightly reduce array Z's size by including just the values of ζ(n) which are actually used, so line 20 would look like this:

            20 DATA 0,0,1.20205690316,0,1.03692775514,0,1.00834927738

but this would break the ability to call FNZ from the prompt for certain values (2, 4, 6, 8) and would require index remapping so it isn't worth the trouble.

○  The INT function at line 150 can't be changed to IP lest you'd get E-14827 instead of the correct E-14828.

Let's run it !

>RUN

t: 5433.88846830 -> X=1.828642E-14828, R= 0
t: 5607.21036604 -> X=2.039534E-15301, R=-1.E-24
t: 7835.62497997 -> X=3.289421E-21382, R=-5.E-24
t: 9330.89271531 -> X=2.000957E-25462, R= 3.E-24
t: 11987.8440717 -> X=1.374136E-32712, R= 2.92E-23
t: 14739.1970154 -> X=2.378127E-40220, R= 2.48E-23
t: 18576.1969026 -> X=1.420375E-50690, R= 1.69E-23

      

For the first, largest root, the 46-digit value is:

1.828643269752522610409732527318069320008652918 * 10^-14828

so we've got 7 correct mantissa digits (save 1 ulp) plus the correct 5-digit exponent.

As the HP-71B is a 12-digit machine we can do no better (7d mantissa + 5d exponent = 12d in all). The other six roots should have about 6-7 correct mantissa digits as well (give or take a few ulp,) plus the correct 5-digit exponent, of course.

For timing, execute instead:


>SETTIME 0 @ CALL @ TIME

which tells us that all 7 roots are obtained in 4.14" in go71b at 128x, 0.55" in Emu71/Win at 972x, or just 8' 50" on a physical HP-71B.

Additional comments

●  This problem is easily solved directly from the command line using recent versions of Mathematica, e.g. to find the first root correct to 17 digits, simply execute:

      10^-t /. FindRoot[RiemannR[10^-t], {t, 14000}, PrecisionGoal → 15, WorkingPrecision → 21]

      1.8286432697525226 x 10^-14828

Then again, Mathematica is a multi-Gb software intended to run on powerful modern hardware, so it's quite remarkable that vintage HP calcs can stand their ground, evaluating the same formula which Mathematica also uses.

●  All seven roots' exponents are so small that they can't be represented in Free42 Decimal (which uses the Intel Decimal Floating-Point Math Library, i.e. it uses IEEE 754-2008 quadruple precision decimal floating-point, which has min. exponent -6,143,) but notice that the 7th root's exponent, -50,690, exceeds even the lower limit of the HP-71B's internal 15-form representation's exponent, which is -50,000, see HP Journal July 1984, p.34.

Further roots have even smaller exponents. e.g. beyond the 7 roots just computed there are 10 even smaller additional roots going down to 10^-500,000, as seen in the graph below.

●  Problem 5 asked for the first seven roots but it's very easy to compute additional roots (albeit with diminished accuracy), e.g. to compute all roots having 5-digit-long exponents (i.e. up to E-99999), we just need to raise the upper limit of the FOR loop in line 70 like this:

       70 ... @ FOR B=500 TO 37000 STEP 500 @ ...

>RUN
      (... the first seven roots, and then ...)
      t: 23080.0653626 -> X=1.618632E-62980, R= 1.156E-23
      t: 28910.7052546 -> X=6.835264E-78891, R=-1.176E-23
      t: 36044.9662268 -> X=1.587852E-98358, R= 8.28E-24

which runs in 5.99" in go71b at 128x, 0.79" in Emu71/Win at 972x, or just 12' 47" on a physical HP-71B.

When computing roots beyond the 3rd one, the STEP 500 at line 70 can be raised to 1000 for faster execution and even further as the spacing between then roots decreases more and more. For instance, using STEP 1000 the timings for the 10 roots above are 5.32", 0.70" and just 11' 21" on a physical HP-71B.

●  Once the roots have been computed, we can easily obtain additional interesting results, such as for instance the unique real positive value of x where R(x) tentatively attains its global minimum value, as seen in this graph:

      
      
To compute the value of such x and the global minimum, simply insert the following program line and run it (now the program will be 967 bytes long):

 95  D=FNROOT(11987,14739,FNR(FVAR+.5)-FNR(FVAR-.5)) @ DISP "Min: ";FNR(D);"at ";FNX$(D,6)

>RUN 95
            Min: -3.14748471571E-13 at 1.1088E-36168

Notes:

○  The expression computing x must be run as a program line because FNROOT can't be executed from the command prompt if it calls a user-defined function, FNR in this case. Attempting to do so results in an error issued by the Math ROM, namely Error 6,  Kybd FN in FNROOT/INTEGRAL.

○  The global minimum is located by finding a root of R'(x) between the limits discussed below, which is computed with good accuracy by using a very simple numerical approximation to the derivative. The resulting global minimum is accurate to 10-12 digits while the corresponding x is accurate to about 5 digits plus the fully correct 5-digit exponent.

○  From the above graph (or a tabulation of values) we notice that the global minimum is located between the 5th and 6th roots, so we use the integer parts of their t parameters (11987 and 14739, respectively) as the initial approximations for faster convergence.

Well, I hope you enjoyed my solution and comments to Problem 5, even if you didn't enjoy the problem itself, which I very much did.

Actually, I think it's indeed a problem with a most amazing, awesome result, where a simple-looking function in its canonical (and Gram series, too) form, with no ad-hoc constants or contrivances and whose graph is a dull logarithmic-like curve for macroscopic arguments, suddenly turns out wiggly and crosses the X-axis by the tiniest amount when evaluated at "nanoscopic" arguments (0.00000...{14,000+ zeros}...0000018286...), completely invisible and utterly unexpected.

In my humble opinion, this does add a measure of "magic" and mystery to Mathematics, and one wonders what other as-yet-undiscovered marvels might be lurking out there...

Next will be Problem 6, which will conclude my months-long six-pronged SCR #12 - Then and Now thread once its main point has been thoroughly made.

V.