Skip to content

Instantly share code, notes, and snippets.

@Mercerenies
Created September 5, 2018 02:10
Show Gist options
  • Save Mercerenies/f9d462814b72bc8816226973e7d0140a to your computer and use it in GitHub Desktop.
Save Mercerenies/f9d462814b72bc8816226973e7d0140a to your computer and use it in GitHub Desktop.
Count for the number of primes and semiprimes up to a given quantity
-- Haskell program which generates the table
import Data.Ratio
import Data.List
divides :: Integral a => a -> a -> Bool
divides y x = x `mod` y == 0
primeFactors :: Integral a => a -> [a]
primeFactors n | n < 1 = error "invalid input"
primeFactors 1 = []
primeFactors n = let p = head $ filter (`divides` n) [2..n]
in p : primeFactors (n `div` p)
-- Short-circuits
lengthAtLeast :: Int -> [a] -> Bool
lengthAtLeast n _ | n <= 0 = True
lengthAtLeast n (_:xs) = lengthAtLeast (n - 1) xs
lengthAtLeast _ _ = False
-- Short-circuits
lengthIs :: Int -> [a] -> Bool
lengthIs n xs = lengthAtLeast n xs && not (lengthAtLeast (n + 1) xs)
isPrime :: Integral a => a -> Bool
isPrime = lengthIs 1 . primeFactors
isSemiPrime :: Integral a => a -> Bool
isSemiPrime = lengthIs 2 . primeFactors
matchesRule :: Integral a => a -> Bool
matchesRule x = isPrime x || isSemiPrime' x
where isSemiPrime' y = isSemiPrime x &&
length (nub $ primeFactors x) == 2
matchesRule' :: Integral a => a -> Bool
matchesRule' x = isPrime x || isSemiPrime x
count :: (a -> Bool) -> [a] -> Int
count f = length . filter f
generateStatistics :: Integral a => (a -> Bool) -> a -> Ratio Int
generateStatistics f n = count f [1..n^2] % (fromIntegral $ n^2)
showStats :: (Show a, Integral a) => (a -> Bool) -> a -> IO ()
showStats f n = print' n >> putStr "\t" >> print stats
where print' = putStr . show
stats :: Double
stats = realToFrac $ generateStatistics f n
firstException :: Integral a => a
firstException =
head . filter (\x -> generateStatistics matchesRule x < 0.5) $ [2..]
main :: IO ()
main = do
putStrLn "** Not counting squares of primes **"
mapM_ (showStats matchesRule) [1..50]
putStrLn ""
putStrLn "** Counting squares of primes **"
mapM_ (showStats matchesRule') [1..50]
** Not counting squares of primes **
1 0.0
2 0.5
3 0.5555555555555556
4 0.625
5 0.6
6 0.5833333333333334
7 0.5714285714285714
8 0.5625
9 0.5432098765432098
10 0.55
11 0.5371900826446281
12 0.5347222222222222
13 0.5266272189349113
14 0.5102040816326531
15 0.52
16 0.5078125
17 0.4982698961937716
18 0.5
19 0.4930747922437673
20 0.49
21 0.48072562358276644
22 0.48140495867768596
23 0.48015122873345933
24 0.4774305555555556
25 0.472
26 0.46597633136094674
27 0.4663923182441701
28 0.4642857142857143
29 0.4625445897740785
30 0.46
31 0.4578563995837669
32 0.4541015625
33 0.44995408631772266
34 0.44982698961937717
35 0.44816326530612244
36 0.44753086419753085
37 0.44631117604090575
38 0.4439058171745152
39 0.44115713346482577
40 0.438125
41 0.4366448542534206
42 0.43594104308390025
43 0.4337479718766901
44 0.43233471074380164
45 0.4311111111111111
46 0.42958412098298676
47 0.4291534631054776
48 0.4262152777777778
49 0.4248229904206581
50 0.4244
** Counting squares of primes **
1 0.0
2 0.75
3 0.7777777777777778
4 0.75
5 0.72
6 0.6666666666666666
7 0.6530612244897959
8 0.625
9 0.5925925925925926
10 0.59
11 0.5785123966942148
12 0.5694444444444444
13 0.5621301775147929
14 0.5408163265306123
15 0.5466666666666666
16 0.53125
17 0.5224913494809689
18 0.5216049382716049
19 0.5152354570637119
20 0.51
21 0.4988662131519274
22 0.49793388429752067
23 0.497164461247637
24 0.4930555555555556
25 0.4864
26 0.47928994082840237
27 0.4787379972565158
28 0.475765306122449
29 0.47443519619500596
30 0.4711111111111111
31 0.46930280957336107
32 0.46484375
33 0.46005509641873277
34 0.4593425605536332
35 0.45714285714285713
36 0.45601851851851855
37 0.4550766983199416
38 0.45221606648199447
39 0.44904667981591057
40 0.445625
41 0.4443783462224866
42 0.4433106575963719
43 0.4413196322336398
44 0.43956611570247933
45 0.4380246913580247
46 0.43620037807183365
47 0.4359438660027162
48 0.4327256944444444
49 0.4310703873386089
50 0.4304
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment