(Spoilers!) Comments and discussion on Valentin's 5th "Then and Now" - Roots

[EdS2]

This thread is an offering to anyone and everyone who would like a place to discuss or digress on Valentin's latest challenge thread:
[VA] SRC #012e - Then and Now: Roots

Beware: there will probably be spoilers within. You are encouraged to read the challenge carefully, have a good think, and try to solve it before exposing yourself to the thoughts of others.

As Valentin does set some splendid challenges but strongly prefers a clean thread for later publication, here's a place where we can discuss our approaches, show any mathematical analysis, use any tools, calculators or computers which might help us to enjoy the challenge.

My aim is only to collect those posts which are not within Valentin's rules - please continue to post in the original thread, especially with your solutions as requested over there.

I do this in a spirit of friendly collaboration, hoping that one strict thread and one loose thread will be a good combination which can keep everyone happy. As such, please no argumentation about the way Valentin runs his challenges - let's stick to the calculations, the mathematics, the analysis.

Links to tangentially related papers, discussions, videos, are all very welcome.

[EdS2]

.
..
.
...
.
..
.
.... possible spoiler thoughts follow ...
.
..
.
...
.
..
.

First off, the article Valentin mentions

Quote: 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

is surely VA027 ("Small Fry - Primes A'counting") found on his articles page.

.
..
.
...
.
..
.

We can read more about Riemann's prime-counting function R(x) on Wikipedia here and in Wolfram's MathWorld here. Warning: at that second link you will find some information about the roots, which are the object of the challenge.

.
..
.
...
.
..
.

On this page you'll see more about this R(x) function, and in particular there's an animation showing how a smooth function progressively approximates the prime-counting stepwise function as more terms in the series are added.

.
..
.
...
.
..
.

I think this might be a clue: we have a wiggly function, and near the origin we might wonder if those wiggles take the function negative.

.
..
.
...
.
..
.

Perhaps by studying the function we can figure out how small the argument might need to be, and figure out how many terms we'll need to be convincing about the locations of the zeros, if we find them.

.
..
.
...
.
..
.

Perhaps the final code will look a little bit like the code found in VA027...

[PeterP]

Ok, short on time and once I saw the factorial I knew that the HP41 was going to be highly limited anyways. But here is how far I have gotten so far:
From the graph shared it is clear that any zero would be between ]0,1[. So to get rid of the logarithm, I made x->10^-y, y in the interval ]0, infinity[

For the denominator D i just noted that zeta(k+1) very quickly goes to 1 (k = 10 or so) within the accuracy of the hp41 or even the 71 I would guess. D is dominated by the factorial (which leaves the 10 digit accuracy of the HP41 past 13! (or k = 12)

Which makes R(y) = 1 + sum [(-y)^k/D] which is an alternating sum, suggesting a pairwise summation (given that we know it converges).

After a little bit of arithmetic (in this forum I'm ashamed to call it even this as most people would have been able to do this - and all of this, and probably have - in their head), I got the following:

R(y) = 1 + sum{n=1,infinity} [y^(2n-1) * (y*(2n-1)*Zeta(2n) - (2n)^2*Zeta(2n+1)) / [(2n-1)* (2n) * (2n)! * Zeta(2n) * Zeta(2n+1)

Where we shall call the new Denominator D'. It is dominated by (2n)!, with the multiplication of the two Zeta functions going even faster to 1. We also see that the HP41 could calculate terms only until about n = 4 for the denominator, unless y is larger and we do very fancy pairwise division. And even then...

so
for n = 1 the term is: y * [y*1*Zeta(2) - (2*1)^2*Zeta(3)] / D'
for n = 2 the term is: y^3 * [y*3*Zeta(4) - (2*2)^2*Zeta(5)] / D'
for n = 3 the term is: y^5 * [y*5*Zeta(6) - (2*3)^2*Zeta(3)] / D'
etc

only for large n relative to Y can the summation pairs get negative:
y(2n-1)Zeta(2n) - (2n)^2Zeta(2n+1)

If for simplicity and a rough aproxmiation we set the two Zeta functions to be the same within the accuracy of the calculator we get
y(2n-1) < (2n)^2

Which, for large n gets close to y < 2n

Given the accuracy of the HP41, I am limited to how far I can take n (because of the factorial in D'), something around 5 or maybe 6 before substantial trickery has to take place.
However, y < 12 does not give negative numbers meaning there is no zero between y = 1 (or x = 10^-1) and y = 12 (x = 10^-12).

I also suspect from looking at the graph and the rapid rise in size of D' that any values of R(y) are very small indeed, maybe even too small for the 10 digit accuracy of the calc. One would have to find a clever way to scale the whole calculation, like many of you came up with in the first problem which dramatically increased the accuracy.

Alas, I have run out of both time and ideas for the moment, so I thought I'd share it. So far, my code consisted of paper and pencil and looking up the values for Zeta (its actually a reMarkable2 with the pen set to pencil)

Would love to hear any corrections (likely many) and suggestions (if this makes any sense at all).

[EdS2]

It seems then, from findings posted, that the argument is going to be very small indeed... and it seems it's a good idea to change variables accordingly.

I found a very short paper which is interesting in its own right, but also sets about evaluating the largest of the roots: it's spoilerific, of course.
[link removed by request]

And a comment in the parent [delinked by request]:

Quote:Remark. The answer to this problem is truly shocking.

Edit: several days after I posted this comment, Valentin advised me in PM that it displeases him to see the above links in plain view, and so I have delinked them. I'll re-emphasise that they lead to material which is in the nature of spoilers and that people should not read the material if they wish to be fully challenged. The links do, of course, persist in archived versions of this comment.

[J-F Garnier]

Hi SRC's readers,

I found in a 1957 paper an asymptotic formula for the Exponential Integral that may help to create an approximation of the R(x) sum.
I have currently little time to spend on the problem but I will try to build a HP-71B solution in the next days, unless Valentin posts his solution in the meantime.
So I give this clue for others to try it.

J-F

[Valentin Albillo]

.
Hi, J-F,

(02-15-2023 04:42 PM)J-F Garnier Wrote:  I have currently little time to spend on the problem but I will try to build a HP-71B solution in the next days, unless Valentin posts his solution in the meantime.

Thanks a lot for your interest and perseverance, J-F, and don't worry, I won't post my original solution until some other has been posted. If after a week none are, I won't post mine either ("If I see interest ..." clause) but will simply make the solution+comments PDF available to those who contributed, which so far includes you, PeterP and EdS2.

Best regards.
V.

[Fernando del Rey]

Hi all,

As a fan of Valentín's challenges, I've also been silently working on this one, but have not been able to make much progress.

I have taken the code for the R(x) function from Valentín's Small Fry - Primes A'counting article.

You just need to modify line 170 to remove the rounding in the final result:

100 DEF FNZ(Z) @ IF Z=2 THEN FNZ=PI*PI/6 ELSE IF Z=3 THEN FNZ=1.20205690316
110 IF Z=4 THEN FNZ=1.08232323371 ELSE IF Z=5 THEN FNZ=1.03692775514
120 IF Z=6 THEN FNZ=1.01734306198 ELSE IF Z=7 THEN FNZ=1.00834927738
130 IF Z<8 THEN END ELSE S=1 @ T=0 @ N=2
140 S=S+N^(-Z) @ N=N+1 @ IF S<>T THEN T=S @ GOTO 140 ELSE FNZ=S
150 DEF FNR(N) @ J=LN(N) @ R=1 @ N=1 @ K=1
160 R=R+1/(K*FNZ(K+1))*J^K/FACT(K) @ IF R<>N THEN K=K+1 @ N=R @ GOTO 160
170 FNR=R @ END DEF

Playing with the FNR function, I see that there is a change from positive to negative values somewhere in the range x=4E-11 to x=5E-11. So there is at least a root inside this range.

I've tried using FNROOT to find the root in this range, and I get the result:

4.3618520797E-11

But FNR of this value is not zero, but something like 1.6E-4, so I don't trust that this is a real root.

Playing with other values of FNR, I see that for x=6E-11 the function is again negative. So there seems to be a another root between x=5E-11 and x=6E-11. Using FNROOT in that range, I get the root:

5.99999999973E-11

But then, curiously, I see these values for FNR:

FNR(5.99999999973E-11) = 3.26676946469E-4
FNR(6E-11) = -8.11848710502E-5

So, it appears that there is a root very close to 6E-11.

But maybe I'm only getting garbage results all the time, or the function is too tricky for FNROOT to handle.

I've tried a change of variables in the x axis. All you need to do is modify line 150:

150 DEF FNR(N) @ J=N @ R=1 @ N=1 @ K=1

And you see the R(x) function with a logarithmic scale in the x axis. But root finding with FNROOT does not seem to get any better by doing this.

This is how far I have gotten so far. I'll keep working on it, but I fear this challenge is utterly beyond my reach.

[Valentin Albillo]

.
Hi, Fernando,

(02-15-2023 09:31 PM)Fernando del Rey Wrote:  As a fan of Valentín's challenges, I've also been silently working on this one, but have not been able to make much progress. [...] This is how far I have gotten so far. I'll keep working on it, but I fear this challenge is utterly beyond my reach.

Thanks for your interest and for posting your sleuthing, Fernando, much appreciated.

Best of lucks tomorrow Thursday 16th ! ...
V.

[J-F Garnier]

(02-15-2023 09:31 PM)Fernando del Rey Wrote:  I have taken the code for the R(x) function from Valentín's Small Fry - Primes A'counting article.

You just need to modify line 170 to remove the rounding in the final result:
[...]
Playing with the FNR function, I see that there is a change from positive to negative values somewhere in the range x=4E-11 to x=5E-11. So there is at least a root inside this range.

Hi Fernando,
These 'roots' are just numeric garbage, as I also experimented with my slightly more accurate program using the Horner method to evaluate R(x), I found sign reversals starting with log(x)=32 i.e. x around 1E-14.

This comes from huge digit cancellation when summing the terms, for instance with log(x)=-30 we sum -30/1+900/4-27000/18+810000/96... until the k! takes precedence (°) and makes the terms tend to zero, and everything should sum to something close to 0.
So I believe the solution is with an approximation of R(x), but I'm afraid the idea I got yesterday is not practical.

Actually, we know, since the above post, that the roots are much much smaller, around log10(x)=-15000.
This doesn't spoil anything, because the challenge is to compute them on the 71B !

J-F

Edit:
(°) this happens at k~80: 30^80=1.5e118; 80! = 7e118

[J-F Garnier]

Well, it took me some efforts but I got an approximation of the R(X) function good enough locate the roots.
Warning: It's a final spoiling of the problem but I guess that it's no more important after more than one week with few posts.

The starting point was the link provided by EdS2. Following the references, I found the paper with the solution.
However porting it to the 71B was not easy (for me).
I'm not a regular user of WolframAlpha, I used it only a few times, for instance to check Free42 results. So I had to learn how to ask for a calculation in the right way to get the magic constants that the approximation is using.
So definitively I didn't follow the rules of the challenge, hence the post in this thread.

[Edit:]
I wrote:
> So definitively I didn't follow the rules of the challenge, hence the post in this thread.
So was my initial understanding, but a couple of feedbacks lead me to move my solution in the main thread where it seems to be more appropriate.

J-F