HP49-50G Willan's formulae to find p(n), the n-th prime number

[Gil]

I tried to put, from Wikipedia, Willan's formulae to find p(n), the n-th prime number, on the HP50G.

Code:


\<< \-> n '1+\GS(ii=1,2^n,FLOOR((n/\GS(j=1,ii,FLOOR(COS(((j-1)!+1)*\pi/j)^2)))^(1/n)))' \->NUM
\>>

For n=1, I get 2, ok
For n=2, I get 3, ok
For n=3, I get 5, ok
For n=4, I get 7, ok
For n=5, I get 33!
For n=6, I get 65!

Could somebody tell me what is wrong in my code above?

NB
However, the following user RPN code works nicely with the loop instructions FOR... NEXT.

Code:


\<< \-> n
  \<< 1 1 2 n ^
    FOR ii 0 1 ii
      FOR j j 1 - ! 1 + j / \pi * COS SQ FLOOR +
      NEXT n SWAP / n INV ^ FLOOR +
    NEXT
  \>>
\>>

Thanks for your help.

Gil

[DavidM]

I believe there's a typo in your original code. The upper bound for the first Σ clause is listed as "2^n", but it should be "2*n".

Repairing that expression should fix the original symbolic version of your code.

[Michael Lopez]

Hi all,

Out of curiosity I changed Gil's RPL program as suggested by DavidM & found correct Prime numbers for n=5 & n=6. However, after that the Prime numbers no longer seem to be correct. Maybe some other error in algebraic formula.

Cheers,

Michael

[DavidM]

(05-19-2024 08:45 AM)Michael Lopez Wrote:  However, after that the Prime numbers no longer seem to be correct. Maybe some other error in algebraic formula.

You are correct, of course, Michael. I was the one with the typo, not Gil. I no longer have the scratchpad I made while checking this, but I'm sure I stopped comparing results as soon as I saw the correct result for p=5; I should have checked further before posting. Remind me not to post when tired.

Coding the Willans formula using the 50g's equation writer is messy, and manually constructing it in symbolic form is only marginally better. To add insult to injury, it's also very inefficient. It's far simpler (and much faster) to just set up an iterative loop using NEXTPRIME.

[Albert Chan]

Willans formula for primes.



Floored expression is simply (pi(i-1) < pi(i) < n), where pi = prime-counting function.
Or, if you prefer, (isprime(i) && pi(i)<n) = 0 or 1

2^n upper bound is thus overkill. All we need is upto p_n - 1

p_n = 1 + Σ((pi(i-1) < pi(i) < n), i = 1 .. p_n - 1)      // stop sum when pi(i) = n

This is still very inefficient, with many expensive pi(i)'s.
Instead, we can solve for pi(min(x)) = n, for x

Example, for x = 1000-th prime

pi(x) = n ≈ x/ln(x)   → Rough estimate, x ≈ n*ln(n) ≈ 6908

pi(2) = 1            // 1st prime = 2
pi(6908) = 888   // secant root, n=1000 --> x≈7780
pi(7780) = 985   // secant root, n=1000 --> x≈7915
pi(7915) = 999
pi(7917) = 999
pi(7919) = 1000 // x = p₁₀₀₀ = 7919

[Gil]

Thanks to all.

My first summation upper limit is indeed 2^n.

My second code in User RPN is OK.

Not very important for the first code as this method is incredibly slow.
I think however that my 1st code is correct, but that there is an internal bug in the equation writer of the HP50G, perhaps in relation with the double summation sign.

Regards

[Gil]

1) At the end, I put the square function just before the floor function as follows :

))^2))^(1/n)))'

Now, the algebraic version of the program looks like
\<< \-> n '1+\GS(ii=1,2^n,FLOOR((n/\GS(j=1,ii,FLOOR(COS(((j-1)!+1)*\pi/j))^2))^(1/n)))' \->NUM
\>>
and works perfectly well — but very slowly, even with EMU48 on the phone.

2) But I persist thinking that there is here a bug
as I wrote previously something of the kind
cos(A)², which is equal to (cos A)².

I presume that the equation writer might have understood
cos(A)² as perhaps, unduly, cos (A²).

3)The effort for correcting the algebraic expression in the equation writer was however worth : calculation of p(8) took 15.3 s, versus 86.1 s for the version with user RPN including the 2 FOR...NEXT instructions.

4) I tried then to speed up the user RPN version in approximate mode and adding —> NUM after the after the 2 FLOOR instructions. A surprisingly bad idea, at first glad, as now p(5) returns the number 33. Reason, now cos (n×pi), with pi approximation, will not give the expected integer.

Regards

[Albert Chan]

(05-19-2024 08:52 PM)Gil Wrote:  I presume that the equation writer might have understood
cos(A)² as perhaps, unduly, cos (A²).

I don't think so. cos(A)² should mean (cos(A))²

Why not make it unambiguous, as SQ(COS(A)) ?

[Albert Chan]

pi(x) = n ≈ x/ln(x) --> x has size around n*ln(n)

n^2 upper limit is more than enough. (there are more complicated, better one)

Code:

from gmpy2 import *
F = lambda j: (fac(j-1)+1) % j == 0
F = lambda x: 1 if x==1 else is_prime(x)    # equivalent, but faster

def prime(n):
    return 1 + sum(
      floor((n/sum(F(j) for j in range(1,i+1))) ** (1/n))
                        for i in range(1,n*n+1))

p2> for i in range(1,10): print i, prime(i)
...
1 2.0
2 3.0
3 5.0
4 7.0
5 11.0
6 13.0
7 17.0
8 19.0
9 23.0

[Gil]

I did, equivalent to your suggestion SQ (cos (A)), Albert,

(cos(A))².

[Albert Chan]

(05-18-2024 10:45 PM)Gil Wrote:  

Code:


\<< \-> n '1+\GS(ii=1,2^n,FLOOR((n/\GS(j=1,ii,FLOOR(COS(((j-1)!+1)*\pi/j)^2)))^(1/n)))' \->NUM
\>>

For n=1, I get 2, ok
For n=2, I get 3, ok
For n=3, I get 5, ok
For n=4, I get 7, ok
For n=5, I get 33!
For n=6, I get 65!

Could somebody tell me what is wrong in my code above?

p_n = 1 + Σ(G(ii), ii = 1 .. 2^n), for some complicated G(ii) = 0 or 1

For n=5 or 6, G(ii) never get to zero.
In other words, HP50g does not recognize cos(multiples of pi) = ±1

> 123456789 PI × 9 ÷ COS     → COS(123456789*PI/9)
> →NUM                                 → -.999999999999
> 123456789 9 ÷ PI × COS     → -1

Fix is probably multiply by PI last, like your RPN code.

[Gil]

Thanks, much more logical indeed to have the multiplication by pi at the end.

Unfortunately, the correct — according to me — algebraic expression with HP50G gives, for p(n=5), 33, a wrong result.

[Albert Chan]

(05-20-2024 09:56 PM)Gil Wrote:  Unfortunately, the correct — according to me — algebraic expression with HP50G gives, for p(n=5), 33, a wrong result.

I think floor expression was not simplified, but was calculated numerically.
This explained why RPN version (done symbolically) take so much longer.

Reminder: we wanted to simulate isprime(j), but with 1 considered a prime

lua> f = fn'j: floor(cos((tgamma(j)+1) / j * pi)^2)'
lua> g = fn'j: floor(cos((tgamma(j)+1) / j * pi))^2'
lua> for j=1,25 do print(j, f(j), g(j)) end
1      1      1
2      1      1
3      1      1
4      0      0
5      1      1
6      0      0
7      1      1
8      0      0
9      0      0
10    0      0
11    1      1
12    0      0
13    1      1
14    0      0
15    0      0
16    0      0
17    0      1
18    0      0
19    0      1
20    0      0
21    0      1
22    0      0
23    0      0
24    0      1
25    0      1

For 53 bits precision:
With f(j), p(n) would be correct for n ≤ pi(13) = 6
With g(j), p(n) would be correct for n ≤ pi(19) = 8