/* Prime conspiracy (Rosetta code) in Picat. http://rosettacode.org/wiki/Prime_conspiracy """ A recent discovery, quoted from Quantamagazine (March 13, 2016): Two mathematicians have uncovered a simple, previously unnoticed property of prime numbers — those numbers that are divisible only by 1 and themselves. Prime numbers, it seems, have decided preferences about the final digits of the primes that immediately follow them. and This conspiracy among prime numbers seems, at first glance, to violate a longstanding assumption in number theory: that prime numbers behave much like random numbers. ─── (original authors from Stanford University): ─── Kannan Soundararajan and Robert Lemke Oliver The task is to check this assertion, modulo 10. * Lets call i -> j a transition if i is the last decimal digit of a prime, and j the last decimal digit of the following prime. Task Considering the first one million primes. Count, for any pair of successive primes, the number of transitions i -> j and print them along with their relative frequency, sorted by i . You can see that, for a given i, frequencies are not evenly distributed. Observation (Modulo 10), primes whose last digit is 9 "prefer" the digit 1 to the digit 9, as its following prime. Extra credit Do the same for one hundred million primes. Example for 10,000 primes 10000 first primes. Transitions prime % 10 → next-prime % 10. 1 → 1 count: 365 frequency: 3.65 % 1 → 3 count: 833 frequency: 8.33 % 1 → 7 count: 889 frequency: 8.89 % 1 → 9 count: 397 frequency: 3.97 % 2 → 3 count: 1 frequency: 0.01 % 3 → 1 count: 529 frequency: 5.29 % 3 → 3 count: 324 frequency: 3.24 % 3 → 5 count: 1 frequency: 0.01 % 3 → 7 count: 754 frequency: 7.54 % 3 → 9 count: 907 frequency: 9.07 % 5 → 7 count: 1 frequency: 0.01 % 7 → 1 count: 655 frequency: 6.55 % 7 → 3 count: 722 frequency: 7.22 % 7 → 7 count: 323 frequency: 3.23 % 7 → 9 count: 808 frequency: 8.08 % 9 → 1 count: 935 frequency: 9.35 % 9 → 3 count: 635 frequency: 6.35 % 9 → 7 count: 541 frequency: 5.41 % 9 → 9 count: 379 frequency: 3.79 % """ With arrays len = 10000 1 -> 1 count: 365 frequency: 3.65% 1 -> 3 count: 833 frequency: 8.33% 1 -> 7 count: 889 frequency: 8.89% 1 -> 9 count: 397 frequency: 3.97% 2 -> 3 count: 1 frequency: 0.01% 3 -> 1 count: 529 frequency: 5.29% 3 -> 3 count: 324 frequency: 3.24% 3 -> 5 count: 1 frequency: 0.01% 3 -> 7 count: 754 frequency: 7.54% 3 -> 9 count: 907 frequency: 9.07% 5 -> 7 count: 1 frequency: 0.01% 7 -> 1 count: 655 frequency: 6.55% 7 -> 3 count: 722 frequency: 7.22% 7 -> 7 count: 323 frequency: 3.23% 7 -> 9 count: 808 frequency: 8.08% 9 -> 1 count: 935 frequency: 9.35% 9 -> 3 count: 635 frequency: 6.35% 9 -> 7 count: 541 frequency: 5.41% 9 -> 9 count: 379 frequency: 3.79% CPU time 0.014 seconds. Backtracks: 0 len = 100000 1 -> 1 count: 4104 frequency: 4.10% 1 -> 3 count: 7961 frequency: 7.96% 1 -> 7 count: 8297 frequency: 8.30% 1 -> 9 count: 4605 frequency: 4.61% 2 -> 3 count: 1 frequency: 0.00% 3 -> 1 count: 5596 frequency: 5.60% 3 -> 3 count: 3604 frequency: 3.60% 3 -> 5 count: 1 frequency: 0.00% 3 -> 7 count: 7419 frequency: 7.42% 3 -> 9 count: 8387 frequency: 8.39% 5 -> 7 count: 1 frequency: 0.00% 7 -> 1 count: 6438 frequency: 6.44% 7 -> 3 count: 6928 frequency: 6.93% 7 -> 7 count: 3627 frequency: 3.63% 7 -> 9 count: 8022 frequency: 8.02% 9 -> 1 count: 8829 frequency: 8.83% 9 -> 3 count: 6513 frequency: 6.51% 9 -> 7 count: 5671 frequency: 5.67% 9 -> 9 count: 3995 frequency: 4.00% CPU time 0.203 seconds. Backtracks: 0 len = 1000000 1 -> 1 count: 42853 frequency: 4.29% 1 -> 3 count: 77475 frequency: 7.75% 1 -> 7 count: 79453 frequency: 7.95% 1 -> 9 count: 50153 frequency: 5.02% 2 -> 3 count: 1 frequency: 0.00% 3 -> 1 count: 58255 frequency: 5.83% 3 -> 3 count: 39668 frequency: 3.97% 3 -> 5 count: 1 frequency: 0.00% 3 -> 7 count: 72827 frequency: 7.28% 3 -> 9 count: 79358 frequency: 7.94% 5 -> 7 count: 1 frequency: 0.00% 7 -> 1 count: 64230 frequency: 6.42% 7 -> 3 count: 68595 frequency: 6.86% 7 -> 7 count: 39603 frequency: 3.96% 7 -> 9 count: 77586 frequency: 7.76% 9 -> 1 count: 84596 frequency: 8.46% 9 -> 3 count: 64371 frequency: 6.44% 9 -> 7 count: 58130 frequency: 5.81% 9 -> 9 count: 42843 frequency: 4.28% CPU time 3.439 seconds. Backtracks: 0 This Picat model was created by Hakan Kjellerstrand, hakank@gmail.com See also my Picat page: http://www.hakank.org/picat/ */ main => go. /* num_primes = 1000000 1 -> 1 count: 42853 frequency: 4.2853% 1 -> 3 count: 77475 frequency: 7.7475% 1 -> 7 count: 79453 frequency: 7.9453% 1 -> 9 count: 50153 frequency: 5.0153% 2 -> 3 count: 1 frequency: 0.0001% 3 -> 1 count: 58255 frequency: 5.8255% 3 -> 3 count: 39668 frequency: 3.9668% 3 -> 5 count: 1 frequency: 0.0001% 3 -> 7 count: 72827 frequency: 7.2827% 3 -> 9 count: 79358 frequency: 7.9358% 5 -> 7 count: 1 frequency: 0.0001% 7 -> 1 count: 64230 frequency: 6.4230% 7 -> 3 count: 68595 frequency: 6.8595% 7 -> 7 count: 39603 frequency: 3.9603% 7 -> 9 count: 77586 frequency: 7.7586% 9 -> 1 count: 84596 frequency: 8.4596% 9 -> 3 count: 64371 frequency: 6.4371% 9 -> 7 count: 58130 frequency: 5.8130% 9 -> 9 count: 42843 frequency: 4.2843% CPU time 3.191 seconds. Backtracks: 0 */ % Testing a matrix instead of a map, slightly faster go ?=> % N = 104_729, % 10_000 primes % N = 1_299_709, % 100_000 primes N = 15_485_863, % 1_000_000 primes % N = 2_038_074_743, % 100 000 000 prime % nope Primes = {P mod 10 : P in primes(N)}, Len = Primes.len, println(num_primes=Len), A = new_array(10,10), bind_vars(A,0), foreach(I in 2..Len) P1 = 1+Primes[I-1], % adjust for 1-based P2 = 1+Primes[I], A[P1,P2] := A[P1,P2] + 1 end, foreach(I in 0..9, J in 0..9, V = A[I+1,J+1], V > 0) printf("%d -> %d count: %5d frequency: %0.4f%%\n", I,J,V,100*V/Len) end, nl. go => true. go ?=> garbage_collect(300_000_000), % N = 104_729, % 10_000 primes % N = 1_299_709, % 100_000 primes N = 15_485_863, % 1_000_000 primes % N = 2_038_074_743, % 100 000 000 prime % nope Primes = {P mod 10 : P in primes(N)}, Len = Primes.len, println(num_primes=Len), A = new_array(10,10), bind_vars(A,0), foreach(I in 2..Len) P1 = 1+Primes[I-1], % adjust for 1-based P2 = 1+Primes[I], A[P1,P2] := A[P1,P2] + 1 end, foreach(I in 0..9, J in 0..9, V = A[I+1,J+1], V > 0) printf("%d -> %d count: %5d frequency: %0.4f%%\n", I,J,V,100*V/Len) end, nl. go => true. % Using sieve: Slightly slower than using primes/1 % We must preallocate memory otherwise it throws stack overflow... go1 ?=> garbage_collect(300_000_000), % N = 104_729, % 10_000 primes % N = 1_299_709, % 100_000 primes N = 15_485_863, % 1_000_000 primes % N = 2_038_074_743, % 100 000 000 prime % nope % Primes = {P mod 10 : P in primes(N)}, Primes = {P mod 10 : P in sieve(N)}, Len = Primes.len, println(num_primes=Len), A = new_array(10,10), bind_vars(A,0), foreach(I in 2..Len) P1 = 1+Primes[I-1], % adjust for 1-based P2 = 1+Primes[I], A[P1,P2] := A[P1,P2] + 1 end, foreach(I in 0..9, J in 0..9, V = A[I+1,J+1], V > 0) printf("%d -> %d count: %5d frequency: %0.4f%%\n", I,J,V,100*V/Len) end, nl. sieve(N) = Res => Sieve = {0 : I in 1..N}, foreach(I in 2..round(sqrt(N)), J in I*I..I..N) Sieve[J] := 1 end, Res := {I : I in 2..N, Sieve[I] == 0}. go1b => % N = 104_729, % 10_000 primes % N = 1_299_709, % 100_000 primes N = 15_485_863, % 1_000_000 primes % N = 2_038_074_743, % 100 000 000 prime % nope Primes = {P mod 10 : P in primes(N)}, Len = Primes.len, println(num_primes=Len), Map = new_map([[I,J]=0 : I in 0..9, J in 0..9]), foreach(I in 2..Len) P1 = Primes[I-1], P2 = Primes[I], Map.put([P1,P2],Map.get([P1,P2])+1) end, foreach([I,J] = V in Map.to_list().sort()) if V > 0 then printf("%d -> %d count: %5d frequency: %0.4f%%\n", I,J,V,100*V/Len) end end, nl. /* using next_prime/2 100_000 limit = 100000 [p1 = 1299689,p2 = 1299709,numPrimes = 100000] 1 -> 1 count: 4104 frequency: 4.1040% 1 -> 3 count: 7961 frequency: 7.9610% 1 -> 7 count: 8297 frequency: 8.2970% 1 -> 9 count: 4605 frequency: 4.6050% 2 -> 3 count: 1 frequency: 0.0010% 3 -> 1 count: 5596 frequency: 5.5960% 3 -> 3 count: 3604 frequency: 3.6040% 3 -> 5 count: 1 frequency: 0.0010% 3 -> 7 count: 7419 frequency: 7.4190% 3 -> 9 count: 8387 frequency: 8.3870% 5 -> 7 count: 1 frequency: 0.0010% 7 -> 1 count: 6438 frequency: 6.4380% 7 -> 3 count: 6928 frequency: 6.9280% 7 -> 7 count: 3627 frequency: 3.6270% 7 -> 9 count: 8022 frequency: 8.0220% 9 -> 1 count: 8829 frequency: 8.8290% 9 -> 3 count: 6513 frequency: 6.5130% 9 -> 7 count: 5671 frequency: 5.6710% 9 -> 9 count: 3994 frequency: 3.9940% CPU time 2.664 seconds. Backtracks: 0 % 1_000_000 is too slow: 88s limit = 1000000 [p1 = 15485857,p2 = 15485863,numPrimes = 1000000] 1 -> 1 count: 42853 frequency: 4.2853% 1 -> 3 count: 77475 frequency: 7.7475% 1 -> 7 count: 79453 frequency: 7.9453% 1 -> 9 count: 50153 frequency: 5.0153% 2 -> 3 count: 1 frequency: 0.0001% 3 -> 1 count: 58255 frequency: 5.8255% 3 -> 3 count: 39668 frequency: 3.9668% 3 -> 5 count: 1 frequency: 0.0001% 3 -> 7 count: 72827 frequency: 7.2827% 3 -> 9 count: 79358 frequency: 7.9358% 5 -> 7 count: 1 frequency: 0.0001% 7 -> 1 count: 64230 frequency: 6.4230% 7 -> 3 count: 68594 frequency: 6.8594% 7 -> 7 count: 39603 frequency: 3.9603% 7 -> 9 count: 77586 frequency: 7.7586% 9 -> 1 count: 84596 frequency: 8.4596% 9 -> 3 count: 64371 frequency: 6.4371% 9 -> 7 count: 58130 frequency: 5.8130% 9 -> 9 count: 42843 frequency: 4.2843% CPU time 88.393 seconds. Backtracks: 0 */ go4 => % Limit = 100_000, Limit = 1_000_000, println(limit=Limit), P1 = 2, P2 = 3, Map = new_map([[I,J]=0 : I in 0..9, J in 0..9]), NumPrimes = 2, while (NumPrimes < Limit) P1Mod = P1 mod 10, P2Mod = P2 mod 10, Map.put([P1Mod,P2Mod],Map.get([P1Mod,P2Mod])+1), P1 := P2, next_prime(P2,NextP), P2 := NextP, NumPrimes := NumPrimes + 1 end, println([p1=P1,p2=P2,numPrimes=NumPrimes]), foreach([I,J] = V in Map.to_list().sort()) if V > 0 then printf("%d -> %d count: %5d frequency: %0.4f%%\n", I,J,V,100*V/NumPrimes) end end, nl. % next prime next_prime(Num, P) :- Num2 = Num + 1, next_prime2(Num2, P). next_prime2(Num, P) :- prime(Num), P = Num. next_prime2(Num, P) :- Num2 = Num+1, next_prime2(Num2,P). % {{trans|Prolog}} % table of nth prime values (up to 100,000) go3 => conspiracy(6), nl. nthprime( 10, 29). nthprime( 100, 541). nthprime( 1000, 7919). nthprime( 10000, 104729). nthprime( 100000, 1299709). nthprime(1000000, 15485863). % added conspiracy(M) :- N is 10**M, nthprime(N, P), sieve(P, Ps), tally(Ps, Counts), Sorted = sort(Counts), show(Sorted). show(Results) :- foreach($tr(D1, D2, Count) in Results) bp.format("~d -> ~d: ~d~n", [D1, D2, Count]) end. % count results tally(L, R) :- tally(L, [], R). tally([_], T, T) :- !. tally([A|As], T0, R) :- [B|_] = As, Da is A mod 10, Db is B mod 10, count(Da, Db, T0, T1), tally(As, T1, R). count(D1, D2, [], [tr(D1, D2, 1)]) :- !. % count(D1, D2, [tr(D1, D2, N)|Ts], [tr(D1, D2, Sn)|Ts]) :- succ(N, Sn), !. count(D1, D2, [tr(D1, D2, N)|Ts], [tr(D1, D2, Sn)|Ts]) :- Sn = N+1, !. count(D1, D2, [T|Ts], [T|Us]) :- count(D1, D2, Ts, Us). % implement a prime sieve sieve(Limit, Ps) :- % numlist(2, Limit, Ns), Ns = 2..Limit, sieve(Limit, Ns, Ps). sieve(Limit, W, W) :- W = [P|_], P*P > Limit, !. sieve(Limit, [P|Xs], [P|Ys]) :- Q is P*P, remove_multiples(P, Q, Xs, R), sieve(Limit, R, Ys). remove_multiples(_, _, [], []) :- !. remove_multiples(N, M, [A|As], R) :- A =:= M, !, remove_multiples(N, M, As, R). remove_multiples(N, M, [A|As], [A|R]) :- A < M, !, remove_multiples(N, M, As, R). remove_multiples(N, M, L, R) :- % plus(M, N, M2), M2 = M+N, remove_multiples(N, M2, L, R).