module Main where

main = do
  let m = 37000000000 
      ps = map (flip (-) m) $ primesFrom m   
      tw = filter (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)])
  putStr   $ "Using " ++ show(k+1) ++ " primes to test "
  putStrLn $ show (q*q-p*p-1) ++ " odds up to " ++ show (q*q) ++ "."
  putStrLn $ "Twin primes above " ++ show m ++ ":"
  print 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 :: [Int]
primes   = 2: 3: sieve 0 primes' 5 
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)

(192353,17358,[(192347,36997368409),(192373,37007371129)])
Using 17359 primes to test 10002719 odds up to 37007371129.
Twin primes above 37000000000:
[(649,651)
Timeout