|
| 1 | +module Juke ( |
| 2 | + fibonaccis, primes, primes', primes'', isPrime, isPrime', |
| 3 | + factors, factors', divides, primeFactors, multiFactors, multiProduct, |
| 4 | + coprime, totient, factorial, fibonacci, |
| 5 | + permutations, sublists, submultilists, partitions, partitionsx, partitionsx', |
| 6 | + fractionExpansion, |
| 7 | + ltoi, itol, |
| 8 | + isqrt |
| 9 | +) where |
| 10 | + |
| 11 | +import Data.Array (array, (!)) |
| 12 | +import Data.Ix (range) |
| 13 | +import Data.List (delete, foldl, foldl', foldl1, foldl1', group, nub, elemIndex) |
| 14 | +import Data.Maybe (fromJust) |
| 15 | + |
| 16 | +import qualified Rising |
| 17 | +import Sort |
| 18 | +import Base |
| 19 | + |
| 20 | +-- fibs |
| 21 | +fibonaccis :: [Integer] |
| 22 | +fibonaccis = 0 : 1 : zipWith (+) fibonaccis (tail fibonaccis) |
| 23 | + |
| 24 | +-- primes |
| 25 | + |
| 26 | +isPrime primes n = n >= head primes && (all notDivisor $ takeWhile small $ primes) where |
| 27 | + small p = p * p <= n |
| 28 | + notDivisor p = not $ p `divides` n |
| 29 | + |
| 30 | +isPrime' primes n = Rising.elem n primes |
| 31 | + |
| 32 | +primes :: (Integral a) => [a] |
| 33 | +primes = 2 : filter (isPrime primes) [3, 5 ..] |
| 34 | + |
| 35 | +primes' :: (Integral a) => [a] |
| 36 | +primes' = [2, 3, 5] ++ Rising.exclude [7, 9 ..] nonprimes where |
| 37 | + nonprimes :: (Integral a) => [a] |
| 38 | + nonprimes = foldr1 f $ map g $ tail $ primes' where |
| 39 | + f (x : xt) ys = x : Rising.merge xt ys |
| 40 | + g p = [p * p' | p' <- [p, p + 2 ..] ] |
| 41 | + |
| 42 | +data People a = VIP a (People a) | Crowd [a] |
| 43 | + |
| 44 | +primes'' :: (Integral a) => [a] |
| 45 | +primes'' = 2 : 3 : Rising.exclude [5, 7 ..] nonprimes where |
| 46 | + nonprimes :: (Integral a) => [a] |
| 47 | + nonprimes = serve $ foldTree mergeP $ map multiples $ tail $ primes'' where |
| 48 | + multiples p = vip [ p * k | k <- [p, p + 2 ..] ] |
| 49 | + vip (x : xs) = VIP x $ Crowd xs |
| 50 | + serve (VIP x xs) = x : serve xs |
| 51 | + serve (Crowd xs) = xs |
| 52 | + foldTree :: (a -> a -> a) -> [a] -> a |
| 53 | + foldTree f ~(x : xs) = f x $ foldTree f $ pairs $ xs where |
| 54 | + pairs ~(x : ~(y : ys)) = f x y : pairs ys |
| 55 | + mergeP :: (Ord a) => People a -> People a -> People a |
| 56 | + mergeP (VIP x xt) ys = VIP x $ mergeP xt ys |
| 57 | + mergeP (Crowd xs) (Crowd ys) = Crowd $ Rising.merge xs ys |
| 58 | + mergeP xs@(Crowd ~(x : xt)) ys@(VIP y yt) = |
| 59 | + case compare x y of |
| 60 | + LT -> VIP x $ mergeP (Crowd xt) ys |
| 61 | + EQ -> VIP x $ mergeP (Crowd xt) yt |
| 62 | + GT -> VIP y $ mergeP xs yt |
| 63 | + |
| 64 | +-- factors |
| 65 | +factors :: (Integral a) => a -> [a] |
| 66 | +factors n = [ x | x <- [1 .. n], 0 == rem n x ] |
| 67 | + |
| 68 | +factors' :: (Integral a) => a -> [a] |
| 69 | +factors' n | n < 1 = [] |
| 70 | +factors' n = quicksort $ map multiProduct $ submultilists $ multiFactors $ n |
| 71 | + |
| 72 | +k `divides` n = n `rem` k == 0 |
| 73 | + |
| 74 | +-- primeFactors |
| 75 | +primeFactors :: (Integral a) => a -> [a] |
| 76 | +primeFactors n = go n primes where |
| 77 | + go n ps@(p : pt) = |
| 78 | + if q < 1 then [] else |
| 79 | + if r == 0 then p : go q ps else |
| 80 | + go n pt |
| 81 | + where (q, r) = quotRem n p |
| 82 | + |
| 83 | +-- multiFactors |
| 84 | +multiFactors :: (Integral a) => a -> [(a, Int)] |
| 85 | +multiFactors n = [ (head xs, length xs) | xs <- group $ primeFactors $ n ] |
| 86 | + |
| 87 | +multiProduct :: (Integral a) => [(a, Int)] -> a |
| 88 | +multiProduct xs = product $ map (uncurry (^)) $ xs |
| 89 | + |
| 90 | +multiCombine f xs ys = filter k $ go f $ combine xs ys where |
| 91 | + k (n, m) = m /= 0 |
| 92 | + go f xs = map k xs where k (n, (m1, m2)) = (n, f m1 m2) |
| 93 | + combine [] [] = [] |
| 94 | + combine [] ys = map k ys where k (n, m) = (n, (0, m)) |
| 95 | + combine xs [] = map k xs where k (n, m) = (n, (m, 0)) |
| 96 | + combine xs@((xn, xm) : xt) ys@((yn, ym) : yt) = |
| 97 | + case compare xn yn of |
| 98 | + LT -> (xn, (xm, 0)) : combine xt ys |
| 99 | + GT -> (yn, (0, ym)) : combine xs yt |
| 100 | + EQ -> (xn, (xm, ym)) : combine xt yt |
| 101 | + |
| 102 | +wideLcm = multiProduct . foldl1 (multiCombine max) . map multiFactors |
| 103 | +wideGcd = multiProduct . foldl1 (multiCombine min) . map multiFactors |
| 104 | + |
| 105 | +-- Coprimes are numbers such that their gcd is 1. |
| 106 | +-- Coprimes have no factors in common, other than 1. |
| 107 | +coprime :: (Integral a) => Predicate a |
| 108 | +coprime a b = 1 == gcd a b |
| 109 | + |
| 110 | +-- Euler's "phi function" or "totient function" |
| 111 | +-- The number of units in Z_n. |
| 112 | +totient :: (Integral a) => a -> Int |
| 113 | +totient n = length [ x | x <- [0 .. n - 1], coprime n x ] |
| 114 | + |
| 115 | +factorial :: (Integral a) => a -> a |
| 116 | +factorial n = go 1 n where |
| 117 | + go k 0 = k |
| 118 | + go k n = go (k * n) (n - 1) |
| 119 | + |
| 120 | +fibonacci :: (Integral a) => a -> a |
| 121 | +fibonacci n = go 0 1 n where |
| 122 | + go a _ 0 = a |
| 123 | + go a b n = go b (a + b) (n - 1) |
| 124 | + |
| 125 | +-- examples of y-combinator and dp-combinator |
| 126 | +-- fib: exponential version |
| 127 | +fib 0 = 0 |
| 128 | +fib 1 = 1 |
| 129 | +fib n = fib (n-2) + fib (n-1) |
| 130 | +-- fib: version suitable for the y-combinator and dp-combinator: |
| 131 | +-- replace recursion with incomplete recursion |
| 132 | +fib' rec 0 = 0 |
| 133 | +fib' rec 1 = 1 |
| 134 | +fib' rec n = rec (n-2) + rec (n-1) |
| 135 | +-- y-combinator version |
| 136 | +fiby n = y fib' n |
| 137 | +-- dp-combinator version |
| 138 | +fiby' n = dp (0, n) fib' n |
| 139 | +-- the better algorithm |
| 140 | +fibg n = go 0 1 n where |
| 141 | + go a b 0 = a |
| 142 | + go a b n = go b (a + b) (n - 1) |
| 143 | +-- the better algorithm, list-based |
| 144 | +fiba n = fibonaccis !! n |
| 145 | +-- Dijkstra's recursion |
| 146 | +-- * F(2n) = (2 F(n-1) + F(n)) * F(n) |
| 147 | +-- * F(2n-1) = F(n-1)^2 + F(n)^2 |
| 148 | +fibd :: (Integral a) => a -> a |
| 149 | +fibd 0 = 0 |
| 150 | +fibd 1 = 1 |
| 151 | +fibd n = |
| 152 | + let (q, r) = quotRem n 2 in |
| 153 | + let [j, k, l] = map fibd $ map (+ q) $ [-1, 0, 1] in |
| 154 | + case r of |
| 155 | + 0 -> (2 * j + k) * k |
| 156 | + 1 -> k ^ 2 + l ^ 2 |
| 157 | +-- Exponentiation - Prelude already defines an efficient exponentiation for Num |
| 158 | +-- * |1 1|^n |F(n+1) F(n) | |
| 159 | +-- * |1 0| = | F(n) F(n-1)| |
| 160 | +data M2 a = M2 a a a a deriving (Eq, Show) |
| 161 | +instance (Num a, Ord a) => Num (M2 a) where |
| 162 | + (M2 a b c d) + (M2 a' b' c' d') = (M2 (a + a') (b + b') (c + c') (d + d')) |
| 163 | + (M2 a b c d) * (M2 a' b' c' d') = (M2 (a * a' + c * b') (b * a' + d * b') (a * c' + c * d') (b * c' + d * d')) |
| 164 | + negate (M2 a b c d) = M2 (-a) (-b) (-c) (-d) |
| 165 | + abs m = if signum m == -1 then negate m else m |
| 166 | + signum (M2 a b c d) = case compare (a * d) (b * c) of { LT -> -1 ; EQ -> 0 ; GT -> 1 } |
| 167 | + fromInteger n = M2 n' 0 0 n' where n' = fromInteger n |
| 168 | +fibx :: (Integral a) => a -> a |
| 169 | +fibx n = r where M2 _ r _ _ = (M2 1 1 1 0) ^ n |
| 170 | + |
| 171 | +-- get a list of all the permutations of a given list |
| 172 | +-- the size of the list of permutations, given a list of length n, is n! |
| 173 | +-- the permutations are given in the lexicographic order determined by the given list |
| 174 | + |
| 175 | +permutations :: (Eq a) => [a] -> [[a]] |
| 176 | +permutations [] = [[]] |
| 177 | +permutations xs = [ x : ys | x <- xs, ys <- permutations $ delete x $ xs ] |
| 178 | + |
| 179 | +sublists :: [a] -> [[a]] |
| 180 | +sublists [] = [[]] |
| 181 | +sublists (h : t) = ys ++ [ h : y | y <- ys ] where ys = sublists t |
| 182 | + |
| 183 | +submultilists :: (Integral n) => [(a, n)] -> [[(a, n)]] |
| 184 | +submultilists [] = [[]] |
| 185 | +submultilists ((h, hi) : t) = [ (h, hi') : ys | hi' <- [0, k hi .. hi], ys <- submultilists t ] where |
| 186 | + k x | x < 0 = -1 |
| 187 | + k x | otherwise = 1 |
| 188 | + |
| 189 | +multCombinations xs ys = filter k $ go xs ys where |
| 190 | + k (n, i) = i /= 0 |
| 191 | + go xs [] = xs |
| 192 | + go [] ys = ys |
| 193 | + go xs@(x@(xn, xi) : xt) ys@(y@(yn, yi) : yt) = |
| 194 | + case compare xn yn of |
| 195 | + LT -> x : go xt ys |
| 196 | + EQ -> (xn, xi + yi) : go xt yt |
| 197 | + GT -> y : go xs yt |
| 198 | + |
| 199 | +multCombinations' xs ys = filter k $ reverse $ go [] xs ys where |
| 200 | + k (n, i) = i /= 0 |
| 201 | + go a [] [] = a |
| 202 | + go a (x : xt) [] = go (x : a) xt [] |
| 203 | + go a [] (y : yt) = go (y : a) [] yt |
| 204 | + go a xs@(x@(xn, xi) : xt) ys@(y@(yn, yi) : yt) = |
| 205 | + case compare xn yn of |
| 206 | + LT -> go (x : a) xt ys |
| 207 | + GT -> go (y : a) xs yt |
| 208 | + EQ -> go ((xn, xi + yi) : a) xt yt |
| 209 | + |
| 210 | +quotCombinations xs ys = multCombinations xs (map k ys) where |
| 211 | + k (xn, xi) = (xn, -xi) |
| 212 | + |
| 213 | +quotCombinations' xs ys = multCombinations' xs (map k ys) where |
| 214 | + k (xn, xi) = (xn, -xi) |
| 215 | + |
| 216 | +prodCombinations xs = foldl multCombinations [] xs |
| 217 | + |
| 218 | +prodCombinations' xs = foldl' multCombinations' [] xs |
| 219 | + |
| 220 | +partitions :: (Integral a) => [a] -> a |
| 221 | +partitions xs = n `quot` d where |
| 222 | + n = factorial $ sum $ xs |
| 223 | + d = product $ map factorial $ xs |
| 224 | + |
| 225 | +partitionsx xs = multiProduct $ quotCombinations n d where |
| 226 | + n = factorial $ sum $ xs |
| 227 | + d = prodCombinations $ map factorial $ xs |
| 228 | + factorial n = prodCombinations $ map multiFactors $ [1 .. n] |
| 229 | + |
| 230 | +partitionsx' xs = multiProduct $ quotCombinations' n d where |
| 231 | + n = factorial $ sum $ xs |
| 232 | + d = prodCombinations $ map factorial $ xs |
| 233 | + factorial n = prodCombinations' $ map multiFactors $ [1 .. n] |
| 234 | + |
| 235 | +-- Fraction Expansion |
| 236 | +-- taken/adapted from http://xorlogic.blogspot.com/2007_09_01_archive.html |
| 237 | +-- (q, f, r) = fractionExpansion b n d |
| 238 | +-- q is the whole number part of n / d |
| 239 | +-- f is the non-repeating part of the decimal expansion (or the expansion in base b) of n / d |
| 240 | +-- r is the repeating part of the decimal expansion (or the expansion in base b) of n / d |
| 241 | +fractionExpansion :: (Integral a) => a -> a -> a -> (a, [a], [a]) |
| 242 | +fractionExpansion base num den = (q, fpart, rpart) where |
| 243 | + (q, r) = num `quotRem` den |
| 244 | + (fpart, rpart) = go num [] [r] |
| 245 | + go num fpart seen = |
| 246 | + let (q, r) = num `quotRem` den in |
| 247 | + if r == 0 then (reverse fpart, []) else |
| 248 | + let n = num `rem` den * base in |
| 249 | + let m = n `quot` den in |
| 250 | + let fpart' = m : fpart in |
| 251 | + let num' = n - den * m in |
| 252 | + case elemIndex num' seen of |
| 253 | + Just k -> splitAt k (reverse fpart') |
| 254 | + Nothing -> go num' fpart' (seen ++ [num']) |
| 255 | + |
| 256 | +ltoi b list = go 0 0 b (reverse list) where |
| 257 | +-- go (accum :: Integer) (place :: Integer) (base :: Integer) (revList :: Integer |
| 258 | + go a _ _ [] = a |
| 259 | + go a p b (x:xt) = go (a + b ^ p * x) (p + 1) b xt |
| 260 | + |
| 261 | +itol b n = go [] b n where |
| 262 | + go a b 0 = a |
| 263 | + go a b n = go (r : a) b q where |
| 264 | + (q, r) = quotRem n b |
| 265 | + |
| 266 | +-- Integers |
| 267 | +isqrt :: Integral a => a -> a |
| 268 | +isqrt = floor . sqrt . fromIntegral |
0 commit comments