fork download 
  1. {-# OPTIONS_GHC -O2 #-}
  2. module Main where
  3. import Data.List hiding (union)
  4. import Data.Array.Unboxed
  5.  
  6. primesSA :: [Int]
  7. primesSA = 2 : prs ()
  8.   where 
  9.     prs () = 3 : sieve (prs ()) 3 []
  10.     sieve (p:ps) x fs = [i*2 + x | (i,True) <- assocs a] 
  11.                         ++ sieve ps (p*p) ((p,0) : 
  12.                              [(s, rem (y-q) s) | (s,y) <- fs])
  13.      where
  14.       q = (p*p-x)`div`2
  15.       a :: UArray Int Bool
  16.       a = accumArray (\ b c -> False) True (1,q-1)
  17.                      [(i,()) | (s,y) <- fs, i <- [y+s, y+s+s..q]]
  18.  
  19. main = print $ 
  20.          --  bprimes !!   400000 -- 100k:1.04-200k:2.69=n^1.37  -- ghc 7.6.3 !!
  21.                                  --           400k:7.35=n^1.45  
  22.          --  primes  !!  1000000 -- 200k:0.32-400k:0.83=n^1.38  -- CONSTANT
  23.                                  --           1mln:2.51=n^1.21
  24.                                  --           2mln:5.72=n^1.19  --   MEMORY !!!
  25.             primesSA !! 10000000 -- 400k:0.43-1mln:1.12=n^1.04    
  26.                                  --           2mln:2.26=n^1.01      9.3M
  27.                                  --           4mln:5.10=n^1.17(1.09)    9.4M
  28.                                  --          10mn:13.01=n^1.02(1.09,1.07)  9.4M
  29. bprimes :: [Int]  -- Richard Bird's sieve
  30. bprimes = _Y $ (2:) . minus [3..] . 
  31.               foldr (\p r-> p*p : union [p*p+p, p*p+2*p..] r) []
  32.  
  33. primes :: [Int]
  34. primes = [2,3,5,7] ++ _Y ((11:) . tail . minus (scanl (+) 11 wh11) 
  35.                . foldi (\(x:xs) r -> x : union xs r)
  36.                . map (\(w,p)-> scanl (\c d-> c + p*d) (p*p) w)
  37.                . equalsBy snd (tails wh11 `zip` scanl (+) 11 wh11))
  38.  
  39. wh3  = 2:wh3                 --  ([3],2)                {1*2,2*3}
  40. wh5  = 2:4:wh5               --  ([5,7],6)                  {2*4,6*5}
  41. wh7  = 4:2:4:2:4:6:2:6:wh7   --  ([7,11,13,17,19,23,29,31],30)  {8*6,30*7}
  42. wh11 = 2:4:2:4:6:2:6:4:2:4:6:6:2:6:4:2:6:4:6:8:4:2:4:2:  
  43.        4:8:6:4:6:2:4:6:2:6:6:4:2:4:6:2:6:4:2:4:2:10:2:10:wh11
  44.  
  45. _Y g = g (_Y g)        -- multistage with non-sharing fixpoint combinator
  46.      -- = g (fix g)    -- two stages with sharing fixpoint combinator
  47.  
  48. foldi f (h:t) = f h . foldi f . unfoldr (\(a:b:c)->Just(f a b,c)) $ t
  49.  
  50. union      a b = ordzipBy id  (:)  (:)  (:)  a b
  51. minus      a b = ordzipBy id  (:) skip skip  a b
  52. equalsBy k a b = ordzipBy k   skip (:) skip  a b
  53. skip       a b = b                       -- skip a = [] ; emit a = [a]
  54.  
  55. ordzipBy k f g h a b = loop a b where    -- concat$unfoldr pull(a,b)
  56.   loop a@(x:t) b@(y:r) = case compare (k x) y of
  57.     LT -> f x (loop t b)                 -- Just(f x,(t,b))
  58.     EQ -> g x (loop t r)                 -- Just(g x,(t,r))
  59.     GT -> h y (loop a r)                 -- Just(h y,(a,r))
Success #stdin #stdout 12.88s 9368KB
comments (?)
stdin 
Standard input is empty
stdout 
179424691
https://ideone.com/rHJ9ub
multistage telescoping production - general code - various wheels (with traditional union/minus: nuoLUE, with fused gaps/hits: vkXCXt) - and now with segmented arrays! -- STUArray: j24jxV (5.9x faster)
language:
Haskell (ghc 8.4.4)
created:
9 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!