```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 ``` ```module Main where main = do let m = 31000000000 -- 31 bil on k', 11 bil on q' ps = map (flip (-) m) \$ primesFrom m tw = filter (uncurry ((==).(2+))) \$ zip ps (tail ps) (ph,(pt1:_)) = span (uncurry ((<).(2+))) \$ zip ps (tail ps) r = (floor.sqrt.fromIntegral) m (h,t) = span (<= r) primes' (k,p,q) = (length h,toInteger\$last h,toInteger\$head t) print (r,k,[(p,p*p),(q,q*q)]) -- putStrLn \$ "P[10,000]=" ++ show (primes!!10000) -- P[10,000]=104743 putStrLn \$ "Testing " ++ show (q*q-p*p-1) ++ " odds in the span up to " ++ show (q*q) ++ " with " ++ show(k+1) ++ " primes." putStrLn \$ "Primes above " ++ show m ++ " (up to the first twin pair):" print \$ map fst ph ++ [fst pt1,snd pt1] -- tw ------------------------------------------------------------- primesFrom m = sieveFrom (length h) ps \$ m`div`2*2+1 where (h,(_:ps)) = span (<= (floor.sqrt.fromIntegral) m) primes sieveFrom k (p:ps) x = let (h,t) = ([x,x+2..t-2],toInteger p^2) qs = map toInteger \$ take k primes' h' = filter (\x-> all ((/=0).(x`rem`)) qs) h in h' ++ sieveFrom (k+1) ps (t+2) ------------------------------------------------------------- (primes, primes') = (kprimes, kprimes') ------------------------------------------------------------- (kprimes, kprimes') = (primes, primes') where primes = 2: 3: sieve 0 primes' 5 :: [Int] primes' = tail primes {- sieve k (p:ps) x -- only 5 billion (?) = [x | x<-[x,x+2..p*p-2], and [(x`rem`p)/=0 | p<-take k primes']] ++ sieve (k+1) ps (p*p+2) -} sieve k (p:ps) x = let (h,t)= ([x,x+2..t-2],p^2) qs = take k primes' h' = filter (\x-> all ((/=0).(x`rem`)) qs) h in h' ++ sieve (k+1) ps (t+2) ------------------------------------------------------------- (qprimes, qprimes') = (primes, primes') where primes = 2: 3: sieve EmptyPQ primes' 5 :: [Int] primes' = tail primes sieve q (p:ps) x = h ++ sieve (addQ q' s (t+s)) ps (t+2) where (s,t) = (2*p,p*p) (h,q') = noComps q [x,x+2..t-2] noComps q xs = foldl f ([],q) xs where f (ps,q) x = maybe (ps++[x],q) ((,)ps) (consultQ q x) consultQ q x | isTopKey q x = Just \$ removeAndReinsertWith (\(v,k)->(v,k+v)) q | otherwise = Nothing -- prime addQ q v@step k@x = insertWithPriority q v k ------------------------------------------------------------- data PQ a = Node a [a] (PQ a) (PQ a) | EmptyPQ isTopKey (Node k _ _ _) x = k == x isTopKey _ _ = False instance (Show a,Integral a)=> Show (PQ a) where show n = showDep 0 n ++ "\n" where showDep d EmptyPQ = "" showDep d n@(Node k vs lt rt) = "" ++ showDep (d+2) lt ++ '\n' : replicate d ' ' ++ "{ " ++ show k ++ '-' : show (map (`div`2) vs) ++ showDep (d+2) rt popOff (Node k1 vs EmptyPQ rt ) = rt popOff (Node k1 vs lt EmptyPQ ) = lt popOff (Node k1 vs n2@(Node k2 vs2 _ _) n3@(Node k3 vs3 _ _)) | k2 < k3 = Node k2 vs2 (popOff n2) n3 | True = Node k3 vs3 n2 (popOff n3) insertWithPriority EmptyPQ v k = Node k [v] EmptyPQ EmptyPQ insertWithPriority n@(Node k1 vs lt rt) v k | k==k1 = Node k1 (v:vs) lt rt | k k1) f n@(Node k1 vs lt@(Node k2 vs2 lt2 rt2) rt@(Node k3 vs3 lt3 rt3)) v k | k==k2 = Node k1 vs (Node k2 (v:vs2) lt2 rt2) rt | k==k3 = Node k1 vs lt (Node k3 (v:vs3) lt3 rt3) | k>k3 = Node k1 vs lt (insertWithPriority rt v k) | kinsertWithPriority q v k') (popOff n) (map (f.flip(,)k) vs) in if isTopKey q k then removeAndReinsertWith f q else q -- 700,000,000 , both kprimes and qprimes the same output: -- [1,31,69,159,183,201,207,241,271,303,321,349,369,387,393,439,471,487,513,519,529,531] ```
 ```1 2 3 4 ``` ```(176068,15998,[(176063,30998179969),(176081,31004518561)]) Testing 6338591 odds in the span up to 31004518561 with 15999 primes. Primes above 31000000000 (up to the first twin pair): [27,67,133,153,159,189,219,237,249,307,309] ```