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   = Node k [v] (Node k1 vs EmptyPQ EmptyPQ) (popOff n)
  | True   = f n v k where                                     -- (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)
      | k<k3  = Node k1 vs (insertWithPriority lt v k) rt
     f n@(Node k1 vs lt rt@EmptyPQ) v k
              = Node k1 vs lt (Node k [v] EmptyPQ EmptyPQ)
     f n@(Node k1 vs lt rt) v k
              = Node k1 vs (insertWithPriority lt v k) rt

removeAndReinsertWith f n@(Node k vs lt rt) = 
  let q= foldr (\(v,k') q->insertWithPriority 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] 

(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]