Prime Puzzle #190
 

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 [URL="http://www.primepuzzles.net/puzzles/puzz_190.htm"]Puzzle 190[/URL] 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...

[QUOTE=rogue;85057]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 [URL="http://www.primepuzzles.net/puzzles/puzz_190.htm"]Puzzle 190[/URL] 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...[/QUOTE]

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

concatenation
 

[QUOTE=drew;85058]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[/QUOTE]
:smile:
Something like that.
If it will be of any help as I'm also not into this sort.

[url]http://en.wikipedia.org/wiki/Concatenation[/url]

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

[QUOTE=drew;85058]I'm not familiar with these problems. What is meant by concatenation? Is it the concatenation of the decimal representations?

1+2+3 = 123[/QUOTE]

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

[QUOTE=rogue;85057]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.[/QUOTE]
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.
[CODE]
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
[/CODE]

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.

[QUOTE=R. Gerbicz;85103]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.[/QUOTE]

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.

[QUOTE=rogue;85113]
Here is the new output:
[/QUOTE]
p=2 is also a solution.
[QUOTE=rogue;85113]
The only thing I need to verify is that the values of none of your variables exceed 2^64.[/QUOTE]
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 )
[CODE]
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.
[/CODE]

[QUOTE=R. Gerbicz;85122]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 )
[/QUOTE]

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.

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.

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?

[QUOTE=rogue;85180]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?[/QUOTE]

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.

[QUOTE=rogue;85113]... my own PPC ASM version[/QUOTE]Seems you have a PowerPC version of your program. Can you make it available so that I see how fast it runs on Linux/Power5+ 1.65 GHz ?
T.

[QUOTE=T.Rex;85391]Seems you have a PowerPC version of your program. Can you make it available so that I see how fast it runs on Linux/Power5+ 1.65 GHz ?
T.[/QUOTE]

See attached. You need GMP installed and you need to compile with -m64.

BTW, searched to 2e11 with no new solutions. I stopped as I don't expect another solution anytime soon. I could run it to 1e13 in a couple of weeks, but why waste my resources. Maybe I'll come back to it in the future, but I've got too many things on my plate at the moment that require horsepower.

Thanks.
But ... how do you compile it ?
gcc -m64 consec.c expmod.s mulmod.s
shows problems with // comments in .s files and also when I've removed them (I'm not an expert in compiling assembler ...).
T.

With // comments:
expmod.s:1: Error: junk at end of line, first unrecognized character is `/'

Without // comments:
expmod.s:6: Error: unsupported relocation against r26

[QUOTE=T.Rex;85480]Thanks.
But ... how do you compile it ?
gcc -m64 consec.c expmod.s mulmod.s
shows problems with // comments in .s files and also when I've removed them (I'm not an expert in compiling assembler ...).
T.

With // comments:
expmod.s:1: Error: junk at end of line, first unrecognized character is `/'

Without // comments:
expmod.s:6: Error: unsupported relocation against r26[/QUOTE]

Try the compiler option:
-mregnames

After that, use google.

Uh, why do your programs list 7 and 11? :unsure: They are not solutions.
[QUOTE=R. Gerbicz;85103]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.

[/QUOTE]

[QUOTE=Maybeso;87224]Uh, why do your programs list 7 and 11? :unsure: They are not solutions.[/QUOTE]

I'm checking (S(p)*X2)%p, so if S(p) is divisible by p then I get that solution but we can get some false solutions (like 7 and 11), because it is possible that X2 is divisible by p, but S(p) isn't divisible by p. Note that there wasn't other false solutions up to 2e11 or something like that, what Rogue searched and it is easy to write a fast independent program to see if it is a good solution or not.

Ok, that makes sense. I was thinking that the magnitude of p compared to that of the lcm factors would prevent false hits for small p, but it seems to prevent them for large p.

Some of the lcm factors: [CODE]7 11 13 37 41 53 73 79 101 137 239 271 4649 9091 9901 21649 333667 513239 909091 265371653 ...[/CODE]Strange that only 7 and 11 create hits.

I have access to a PICK OS machine, which is string based. Programs, numbers and variables are all dynamic length strings. So I wrote a brute force program for this puzzle in compiled PICK BASIC, with the inner loop like this:[CODE] M = 1
FOR N = 2 TO P
M = MOD(M:N, P)
NEXT N
IF M ELSE PRINT P[/CODE]But the machine is an emulation running on linux, and it's multi-user, so although the code looks pretty, it doesn't run much faster than a program using normal integers and 'manual' concatenation. I do get arbitrary length numbers for free, but time is still the main obstacle.

ps. I could port over your ingenius method, but since it doesn't require concatenation, there would be no advantage.

Bruce