Legendre prime counting function: Difference between revisions

Content added Content deleted
(→‎{{header|Haskell}}: added solution)
(→‎{{header|Haskell}}: simplified memoization)
Line 360: Line 360:


Memoization utilities:
Memoization utilities:
<lang haskell>type Memo2 a = Memo (Memo a)
<lang haskell>data Memo a = Node a (Memo a) (Memo a)

data Memo a = Node a (Memo a) (Memo a)
deriving Functor
deriving Functor


Line 370: Line 368:
| odd n = memo l (n `div` 2)
| odd n = memo l (n `div` 2)
| otherwise = memo r (n `div` 2 - 1)
| otherwise = memo r (n `div` 2 - 1)

memo2 :: Integral a => Memo2 b -> a -> a -> b
memo2 f = memo . memo f


nats :: Integral a => Memo a
nats :: Integral a => Memo a
nats = Node 0 ((+1).(*2) <$> nats) ((*2).(+1) <$> nats)
nats = Node 0 ((+1).(*2) <$> nats) ((*2).(+1) <$> nats)


memoize :: Integral a => (a -> b) -> Memo b
memoize :: Integral a => (a -> b) -> a -> b
memoize f = f <$> nats
memoize f = memo (f <$> nats)


memoize2 :: (Integral a, Integral b) => (a -> b -> c) -> Memo2 c
memoize2 :: (Integral a, Integral b) => (a -> b -> c) -> a -> b -> c
memoize2 f = memoize (memoize . f)
memoize2 f = memoize (memoize . f)


memoList :: [b] -> Memo b
memoList :: [b] -> Integer -> b
memoList (x:xs) = Node x (memoList l) (memoList r)
memoList = memo . mkList
where (l,r) = split xs
where
split [] = ([],[])
mkList (x:xs) = Node x (mkList l) (mkList r)
split [x] = ([x],[])
where (l,r) = split xs
split (x:y:xs) = let (l,r) = split xs in (x:l, y:r)</lang>
split [] = ([],[])
split [x] = ([x],[])
split (x:y:xs) = let (l,r) = split xs in (x:l, y:r)</lang>


Computation of Legendre function:
Computation of Legendre function:
Line 401: Line 398:
else go z (r `shiftR` 1) (q `shiftR` 2)
else go z (r `shiftR` 1) (q `shiftR` 2)


phi = memoize2 phiM
p :: Memo Integer
where
p = memoList (undefined : primes)
phiM x 0 = x

phiM x a = phi x (a-1) - phi (x `div` p a) (a - 1)
phi :: Integer -> Integer -> Integer
phi x 0 = x
p = memoList (undefined : primes)
phi x a = memo2 phiM x (a-1) - memo2 phiM (x `div` memo p a) (a - 1)

phiM :: Memo2 Integer
phiM = memoize2 phi


legendrePi :: Integer -> Integer
legendrePi :: Integer -> Integer
legendrePi n
legendrePi n
| n < 2 = 0
| n < 2 = 0
| otherwise = memo2 phiM n a + a - 1 where a = legendrePi (isqrt n)
| otherwise = phi n a + a - 1
where a = legendrePi (floor (sqrt (fromInteger n)))


main = mapM_ (\n -> putStrLn $ show n ++ "\t" ++ show (legendrePi (10^n))) [1..9]</lang>
main = mapM_ (\n -> putStrLn $ show n ++ "\t" ++ show (legendrePi (10^n))) [1..7]</lang>


<pre>λ> main
<pre>λ> main