Prime Puzzle #190

[rogue]

I don't know how many of you follow Carlos Rivera's Prime Puzzles site, but I view it every other week or so looking for interesting programming challenges. I worked on Puzzle 190 about 4 years ago. I actually took the search up to 5e8 (without finding a new solution), which took about two months of CPU time when I ran it. Today I could search that same range in about one day. I consider the puzzle to be interesting not because the next solution has not been found, but because optimizing the algorithm is. What's funny is that I'm looking at the code again after four years and thanks to my lack of documenting the code, don't remember clearly how it works. All I know is that it does work (as I tested the results against code that uses a brute force method).

I have a question and a challenge. First, who is interested in pursuing this search further? Second, I challenge other coders out there to write the most optimal code for the search that they can.

For a base, my code can search a range of 1e6 starting from 1e9 in about 8 minutes. I haven't tried to optimize further (since I need to refresh my knowledge on the algorithm I used), but I assume that I and others could do much better given some time.

Have fun...

[drew]

Originally Posted by rogue

I don't know how many of you follow Carlos Rivera's Prime Puzzles site, but I view it every other week or so looking for interesting programming challenges. I worked on Puzzle 190 about 4 years ago. I actually took the search up to 5e8 (without finding a new solution), which took about two months of CPU time when I ran it. Today I could search that same range in about one day. I consider the puzzle to be interesting not because the next solution has not been found, but because optimizing the algorithm is. What's funny is that I'm looking at the code again after four years and thanks to my lack of documenting the code, don't remember clearly how it works. All I know is that it does work (as I tested the results against code that uses a brute force method).

I have a question and a challenge. First, who is interested in pursuing this search further? Second, I challenge other coders out there to write the most optimal code for the search that they can.

For a base, my code can search a range of 1e6 starting from 1e9 in about 8 minutes. I haven't tried to optimize further (since I need to refresh my knowledge on the algorithm I used), but I assume that I and others could do much better given some time.

Have fun...

I'm not familiar with these problems. What is meant by concatenation? Is it the concatenation of the decimal representations?

1+2+3 = 123

And for multi-digit numbers?

1+20 = 120

Drew

[mfgoode]

Originally Posted by drew

I'm not familiar with these problems. What is meant by concatenation? Is it the concatenation of the decimal representations?

1+2+3 = 123

And for multi-digit numbers?

1+20 = 120

Drew

Something like that.
If it will be of any help as I'm also not into this sort.

http://en.wikipedia.org/wiki/Concatenation

BTW: their down load I found is very slow but I have referenced the site in my Favourites.
Mally

[rogue]

Originally Posted by drew

I'm not familiar with these problems. What is meant by concatenation? Is it the concatenation of the decimal representations?

1+2+3 = 123

Yes. Using the website's notation S(11) = 1234567891011.

[R. Gerbicz]

Originally Posted by rogue

I actually took the search up to 5e8 (without finding a new solution), which took about two months of CPU time when I ran it. Today I could search that same range in about one day.

I've also written a program.
Rogue, the last time is the total computation time=687 seconds up to 2^29,
(this is about 5*10^8) on my machine.

C:\>p190
Searching for solutions up to 2^29
p=2 Time=0.000 sec.
p=3 Time=0.000 sec.
p=5 Time=0.000 sec.
p=7 Time=0.000 sec.
p=11 Time=0.000 sec.
p=8167 Time=0.015 sec.
p=371321 Time=0.468 sec.
Time=687.640 sec.

See the attachment for the gmp source. ( just modify limitbit to
search for all solutions up to 2^limitbit ), limitbit<=50 is good.
Or here it is an easy UBASIC implementation, but this is slower and limited to 2^32, because for this I've not written a prime generation routine.

   10   cls:clr time:point -2
   20   dim R(10),S(10),T(10),U(10),W(10),Y(10)
   30   P=1:D=1:B=1:C=10
   40   U(1)=10^9:T(1)=10^9
   50   for I=0 to 10:R(I)=10^I:next I
   60   for I=0 to 10:W(I)=1:next I
   70   while P<4294967291
   80   P=nxtprm(P)
   90   if P<C then goto 140
  100   D+=1:B*=10:C*=10:X=1
  110   for I=1 to D:X=lcm(X,R(I)-1):next I
  120   X2=X*X
  130   for I=1 to D:S(I)=X2\(R(I)-1):Y(I)=S(I)\(R(I)-1):next I
  140   L=P-B
  150   for H=2 to D-1
  160   U(H)=modpow(U(H-1),10,P)
  170   T(H)=(modpow(T(H-1),10,P)*U(H))@P
  180   next H
  190   for H=D-2 to 1 step -1
  200   W(H)=(W(H+1)*T(H+1))@P
  210   next H
  220   A=0
  230   for I=1 to D-1
  240   A+=((T(I)-1)*Y(I)+(T(I)*R(I-1)-1)*S(I)-X2)*W(I)
  250   next I
  260   Q=modpow(R(D),L+1,P)
  270   A*=Q
  280   A+=(Q-R(D))*Y(D)+S(D)*Q*R(D-1)
  290   if A@P=0 then print P;time
  300   wend
  310   print time
  320   end

It isn't very known but you can calculate the required mod very quickly.
For example the concenation of the 1 digit numbers is:
A= 123456789 then
10A=1234567890 so
9A=1111111101=10^2*(10^8-1)/(10^1-1)+1 from this
A=(10^2*(10^8-1)/(10^1-1)+1)/(10^1-1)

Using this method the concenation of all H digits numbers is ( using an easier form )
((10^(9*H*10^(H-1))-1)/(10^H-1)+10^(9*H*10^(H-1)+H-1)-1)/(10^H-1)-1

You can calculate not just all H digit numbers but also the concenation of the first
L+1 numbers ( which are D digits ) by the same method.

In the two codes:
R(H)=10^H
U(H)=10^(9*10^(H-1)) mod p
T(H)=10^(9*H*10^(H-1)) mod p
W(H)=10^((D-1)*9*10^(D-2)+(D-2)*9*10^(D-3)+...+(H+1)*9*10^H) mod p
etc.

We use the U table to calculate the T table. By recursion the computation is easy.
To avoid the modular divisions I multiply each fraction by the least common multiplies
of the denominators. This is lcm(10^1-1,10^2-1,...,10^D-1)^2 if p has got D digits in
base 10. This is the reason why we get some possible false p values. S and Y arrays
contain the required precalculated multiplies for each fractions.

[rogue]

Originally Posted by R. Gerbicz

I've also written a program.
Rogue, the last time is the total computation time=687 seconds up to 2^29,
(this is about 5*10^8) on my machine.

C:\>p190
Searching for solutions up to 2^29
p=2 Time=0.000 sec.
p=3 Time=0.000 sec.
p=5 Time=0.000 sec.
p=7 Time=0.000 sec.
p=11 Time=0.000 sec.
p=8167 Time=0.015 sec.
p=371321 Time=0.468 sec.
Time=687.640 sec.

It isn't very known but you can calculate the required mod very quickly.
For example the concenation of the 1 digit numbers is:
A= 123456789 then
10A=1234567890 so
9A=1111111101=10^2*(10^8-1)/(10^1-1)+1 from this
A=(10^2*(10^8-1)/(10^1-1)+1)/(10^1-1)

Using this method the concenation of all H digits numbers is ( using an easier form )
((10^(9*H*10^(H-1))-1)/(10^H-1)+10^(9*H*10^(H-1)+H-1)-1)/(10^H-1)-1

You can calculate not just all H digit numbers but also the concenation of the first
L+1 numbers ( which are D digits ) by the same method.

In the two codes:
R(H)=10^H
U(H)=10^(9*10^(H-1)) mod p
T(H)=10^(9*H*10^(H-1)) mod p
W(H)=10^((D-1)*9*10^(D-2)+(D-2)*9*10^(D-3)+...+(H+1)*9*10^H) mod p
etc.

We use the U table to calculate the T table. By recursion the computation is easy.
To avoid the modular divisions I multiply each fraction by the least common multiplies
of the denominators. This is lcm(10^1-1,10^2-1,...,10^D-1)^2 if p has got D digits in
base 10. This is the reason why we get some possible false p values. S and Y arrays
contain the required precalculated multiplies for each fractions.

Brilliant! I've taken your code and made a few enhancements (such as replacing GMPs powm with my own PPC ASM version) and am using my own sieve, which allows me to input the start and stop range and gives no false primes.

Here is the new output:

p=3 Time=0.000 sec.
p=5 Time=0.000 sec.
p=7 Time=0.000 sec.
p=11 Time=0.000 sec.
p=8167 Time=0.000 sec.
p=371321 Time=0.080 sec.
The search completed in 117.37 seconds

and from 5e8 to 1e9 it took an additional 113.65 seconds.

I am now pushing the limits up to 1e11 or even higher. The only thing I need to verify is that the values of none of your variables exceed 2^64. If so, I will probably replace the rest of the GMP code and improve the speed even more.

[R. Gerbicz]

Originally Posted by rogue

Here is the new output:

p=2 is also a solution.

Originally Posted by rogue

The only thing I need to verify is that the values of none of your variables exceed 2^64.

The S and Y arrays that store the required multipliers for each fractions has got about 180 bits numbers if p is about 10^9. Because lcm(10^1-1,10^2-1,...,10^9-1) has got 94 bits. But it is possible to avoid these large numbers by a binary tree method: if D=9, modifying the definition of S and Y:
S[i]=((10^1-1)*(10^2-1)*...*(10^9-1))^2/(10^i-1) mod p
Y[i]=((10^1-1)*(10^2-1)*...*(10^9-1))^2/(10^i-1)^2 mod p
( so we multiple the fractions by ((10^1-1)*(10^2-1)*...*(10^9-1))^2 )

Precompute RR[i]=10^i-1 in the beginning of the computation.

V[1]=RR[1]*RR[2]*RR[3]*RR[4] mod p
V[2]=V[1]*RR[5]*RR[6] mod p
V[3]=V[2]*RR[7] mod p
V[4]=V[3]*RR[8] mod p
V[5]=V[3]*RR[9] mod p
V[6]=V[2]*RR[8]*RR[9] mod p
V[7]=V[1]*RR[7]*RR[8]*RR[9] mod p
V[8]=V[7]*RR[5] mod p
V[9]=V[7]*RR[6] mod p
V[10]=RR[5]*RR[6]*RR[7]*RR[8]*RR[9] mod p
V[11]=V[10]*RR[1]*RR[2] mod p
V[12]=V[11]*RR[3] mod p
V[13]=V[11]*RR[4] mod p
V[14]=V[10]*RR[3]*RR[4] mod p
V[15]=V[14]*RR[1] mod p
V[16]=V[14]*RR[2] mod p

Y[1]=V[16]^2 mod p
S[1]=Y[1]*RR[1] mod p
Y[2]=V[15]^2 mod p
S[2]=Y[2]*RR[2] mod p
Y[3]=V[13]^2 mod p
S[3]=Y[3]*RR[3] mod p
Y[4]=V[12]^2 mod p
S[4]=Y[4]*RR[4] mod p
Y[5]=V[9]^2 mod p
S[5]=Y[5]*RR[5] mod p
Y[6]=V[8]^2 mod p
S[6]=Y[6]*RR[6] mod p
Y[7]=V[6]^2 mod p
S[7]=Y[7]*RR[7] mod p
Y[8]=V[5]^2 mod p
S[8]=Y[8]*RR[8] mod p
Y[9]=V[4]^2 mod p
S[9]=Y[9]*RR[9] mod p

Requiring only 36 mulmod and 9 squaremod operations per p prime
if D=9 to construct this modified S and Y arrays.

[rogue]

Originally Posted by R. Gerbicz

p=2 is also a solution.

The S and Y arrays that store the required multipliers for each fractions has got about 180 bits numbers if p is about 10^9. Because lcm(10^1-1,10^2-1,...,10^9-1) has got 94 bits. But it is possible to avoid these large numbers by a binary tree method: if D=9, modifying the definition of S and Y:
S[i]=((10^1-1)*(10^2-1)*...*(10^9-1))^2/(10^i-1) mod p
Y[i]=((10^1-1)*(10^2-1)*...*(10^9-1))^2/(10^i-1)^2 mod p
( so we multiple the fractions by ((10^1-1)*(10^2-1)*...*(10^9-1))^2 )

My sieve doesn't spit out 2. I normally handle it as a special case. 2 is already a known solution, so I ignored it.

Right now, I only use my mulmod/expmod on the U, T, and R arrays and the value of l. AFAICT, only the upper values of the R array do that, but I only get to about R[10] in my testing.

[R. Gerbicz]

I think that S(p) mod p is a random value, so the probability that S(p) is divisible by p is 1/p. By Mertens theorem sumprime(p=2,n,1/p) is about log(log(n)), here the base is e. So there is about log(log(n)) solutions up to n, that's the reason why there is only few known solutions so far.

[rogue]

I've gone up to 1.6e11 (so far) with no new solutions. I have a question. Clearly the value of a can exceed p, but is a % p = S(p) % p for all p or only when S(p) % p = 0?

[R. Gerbicz]

Originally Posted by rogue

I've gone up to 1.6e11 (so far) with no new solutions. I have a question. Clearly the value of a can exceed p, but is a % p = S(p) % p for all p or only when S(p) % p = 0?

Actually I multiple the fractions by lcm(10^1-1,...,10^D-1)^2, this is X2 in the code, so a%p=(S(p)*X2)%p, it means that if S(p)%p=0 then clearly a%p=0, but it's not true that a%p=S(p)%p for all p.