Twin Primes
July 26, 2019
According to number theory:
m is the base of a twin-prime pair (m, m+2) if and only if 4((m−1)! + 1) == –m (mod m (m+2)).
Your task is to write a program that implements the criterion given above, then calculate the twin primes less than a thousand. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
@programmingpraxis: Not understanding the filter. According to your calculations 3 is a solution. So substituting it into the filter — 4((m−1)! + 1) == –m (mod m (m+2)) — I get 4((3−1)! + 1) == –3 (mod 3 (3+2)) —> 4(2! + 1) == -3 (mod 3 5) —> 12 == -3 (mod 3 5) —> 12 == -3 3 —> 12 == -9, which fails.
What am I missing?
@bookofstevegraham: The left-hand side of the congruence is 12, which you computed correctly. The right-hand side of the congruence is -3, which you also computed correctly. The modulus of the congruence is m * (m+2) = 3 * 5 = 15. Since 12 and -3 differ by 15, they are congruent modulo 15 and 3 is the base of the twin prime pair (3, 5).
Cannot seem to upload solution (Steve)
Klong version (Try #5)
3+!997 [3 4 5 6 7 8 9 ... 999] {&/x!:\2+!_x^1%2}'3+!997 :" Which of the numbers are prime?" [1 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 1 ... 0] a@&{&/x!:\2+!_x^1%2}'a::3+!997 :" Only keep the primes" [3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743 751 757 761 769 773 787 797 809 811 821 823 827 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941 947 953 967 971 977 983 991 997] ,:'a@&{&/x!:\2+!_x^1%2}'a::3+!997 :" Get each pair" [[3 5] [5 7] [7 11] [11 13] [13 17] [17 19] [19 23] [23 29] ... [991 997]] {2=(x@1)-x@0}'b::,:'a@&{&/x!:\2+!_x^1%2}'a::3+!997 :" Which pairs have a difference of 2" [1 1 0 1 0 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0 0 1 0 1 0 1 0 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] &{2=(x@1)-x@0}'b::,:'a@&{&/x!:\2+!_x^1%2}'a::3+!997 :" Only keep the ones with the proper difference - these are indices" [0 1 3 5 8 11 15 18 24 26 31 33 39 41 43 47 50 55 58 62 67 79 81 87 96 102 107 111 114 118 138 140 142 146 150] b::b@&{2=(x@1)-x@0}'b::,:'a@&{&/x!:\2+!_x^1%2}'a::3+!997 :" These are the pairs" [[3 5] [5 7] [11 13] [17 19] [29 31] [41 43] [59 61] [71 73] [101 103] [107 109] [137 139] [149 151] [179 181] [191 193] [197 199] [227 229] [239 241] [269 271] [281 283] [311 313] [347 349] [419 421] [431 433] [461 463] [521 523] [569 571] [599 601] [617 619] [641 643] [659 661] [809 811] [821 823] [827 829] [857 859] [881 883]] mod::{:[x>0; x!y; (x+y)!y]} :dyad mod(12; 15) 12 mod(-3; 15) 12 {[m1 m2]; m1::x@0; m2::x@1; .p($x," = ",$mod(4*1+*/1+!m1-1; m1*m2)=mod(-m1; m1*m2))}'b;"" :" All pairs are twin primes" [3 5 = 1] [5 7 = 1] [11 13 = 1] [17 19 = 1] [29 31 = 1] [41 43 = 1] [59 61 = 1] [71 73 = 1] [101 103 = 1] [107 109 = 1] [137 139 = 1] [149 151 = 1] [179 181 = 1] [191 193 = 1] [197 199 = 1] [227 229 = 1] [239 241 = 1] [269 271 = 1] [281 283 = 1] [311 313 = 1] [347 349 = 1] [419 421 = 1] [431 433 = 1] [461 463 = 1] [521 523 = 1] [569 571 = 1] [599 601 = 1] [617 619 = 1] [641 643 = 1] [659 661 = 1] [809 811 = 1] [821 823 = 1] [827 829 = 1] [857 859 = 1] [881 883 = 1] ""In Python, using a “twin-wheel”. A twin can only happen at 11, 17 and 29 (repeating every 30).
from itertools import accumulate, cycle, chain def fac_mod(n, m): """calculates n! (mod m)""" if n >= m: return 0 fac = 1 for i in range(2, n+1): fac = fac * i % m if fac == 0: break return fac def is_twin(m): M = m * (m + 2) return (4 * (fac_mod(m-1, M) + 1) + m) % M == 0 def gtwin(limit=None): for cand in accumulate(chain((3, 2, 6), cycle((6, 12, 12)))): if limit is not None and cand > limit: break if is_twin(cand): yield cand print(list(gtwin(1000)))And just for fun, here a “twin-wheel” for the primes 2,3,5 and 7. Using this wheel only one in seven odd numbers have to be tested.
def gtwin2(limit=None): for cand in accumulate(chain((3, 2, 6, 6), cycle((12, 12, 18, 12, 30, 6, 30, 12, 18, 12, 12, 6, 12, 12, 6)))): if limit is not None and cand > limit: break if is_twin(cand): yield cand