63 lines
2.2 KiB
Haskell
63 lines
2.2 KiB
Haskell
import Control.Concurrent (setNumCapabilities)
|
|
import Control.Monad.Par (runPar, get, spawnP)
|
|
import Control.Monad (join, (>=>))
|
|
import Data.List.Split (chunksOf)
|
|
import Data.List (intercalate, mapAccumL, genericTake, genericDrop)
|
|
import Data.Bifunctor (bimap)
|
|
import GHC.Conc (getNumProcessors)
|
|
import Math.NumberTheory.Primes (factorise, unPrime)
|
|
import Text.Printf (printf)
|
|
|
|
chowla :: Word -> Word
|
|
chowla 1 = 0
|
|
chowla n = f n
|
|
where
|
|
f = (-) =<< pred . product . fmap sumFactor . factorise
|
|
sumFactor (n, e) = foldr (\p s -> s + unPrime n^p) 1 [1..e]
|
|
|
|
chowlas :: [Word] -> [(Word, Word)]
|
|
chowlas [] = []
|
|
chowlas xs = runPar $ join <$>
|
|
(mapM (spawnP . fmap ((,) <*> chowla)) >=> mapM get) (chunksOf (10^6) xs)
|
|
|
|
chowlaPrimes :: [(Word, Word)] -> (Word, Word) -> (Word, Word)
|
|
chowlaPrimes chowlas range = (count chowlas, snd range)
|
|
where
|
|
isPrime (1, n) = False
|
|
isPrime (_, n) = n == 0
|
|
count = fromIntegral . length . filter isPrime . between range
|
|
between (min, max) = genericTake (max - pred min) . genericDrop (pred min)
|
|
|
|
chowlaPerfects :: [(Word, Word)] -> [Word]
|
|
chowlaPerfects = fmap fst . filter isPerfect
|
|
where
|
|
isPerfect (1, _) = False
|
|
isPerfect (n, c) = c == pred n
|
|
|
|
commas :: (Show a, Integral a) => a -> String
|
|
commas = reverse . intercalate "," . chunksOf 3 . reverse . show
|
|
|
|
main :: IO ()
|
|
main = do
|
|
cores <- getNumProcessors
|
|
setNumCapabilities cores
|
|
printf "Using %d cores\n" cores
|
|
|
|
mapM_ (uncurry (printf "chowla(%2d) = %d\n")) $ take 37 allChowlas
|
|
mapM_ (uncurry (printf "There are %8s primes < %10s\n"))
|
|
(chowlaP
|
|
[ (1, 10^2)
|
|
, (succ $ 10^2, 10^3)
|
|
, (succ $ 10^3, 10^4)
|
|
, (succ $ 10^4, 10^5)
|
|
, (succ $ 10^5, 10^6)
|
|
, (succ $ 10^6, 10^7) ])
|
|
|
|
mapM_ (printf "%10s is a perfect number.\n" . commas) perfects
|
|
printf "There are %2d perfect numbers < 35,000,000\n" $ length perfects
|
|
where
|
|
chowlaP = fmap (bimap commas commas) . snd
|
|
. mapAccumL (\total (count, max) -> (total + count, (total + count, max))) 0
|
|
. fmap (chowlaPrimes $ take (10^7) allChowlas)
|
|
perfects = chowlaPerfects allChowlas
|
|
allChowlas = chowlas [1..35*10^6]
|