Re: [Haskell-cafe] NumberTheory library

Am Donnerstag, 5. Mai 2005 06:20 schrieb ajb@spamcop.net:

G'day all.

Quoting Bo Herlin :

...

But there are functions that I cant find and that I assume others before me must have missed and then perhaps also implemented them. Is there any standard library with functions like:

binomial isCatalan nthCatalan nextCatalan isPrime nthPrime nextPrime factorize isFibonacci nthFibonacci nextFibonacci

You might want to check out the archives of the mailing list, too. These sorts of problems occasionally get solved. For the record, here are a few of my attempts:

http://andrew.bromage.org/Fib.hs (Fairly fast Fibonacci numbers) http://andrew.bromage.org/WheelPrime.hs (Fast factorisation) http://andrew.bromage.org/Prime.hs (AKS primality testing)

Cheers, Andrew Bromage

The fibonacci numbers are really cool, for fib 100000 and fib 200000 I got the same performance as with William Lee Irwins 'fastest fibonacci numbers in the west' from a couple of months ago. The AKS primality testing is definitely not for Haskell, I'm afraid - or somebody would have to come up with a better way to model polynomials. Trying relatively small numbers: Prelude Prime> isPrime (2^19-1) True (2.76 secs, 154748560 bytes) Prelude Prime> isPrime (2^29-1) <interactive>: out of memory (requested 277872640 bytes) after a lot of swapping. I attempted to implement that algorithm, too, a while ago, your version is much better, but still not really useful, it seems. Thanks for WheelPrime, it's a good idea which I didn't know before. I incorporated it into my Primality-module, elementary factorization is now about twice as fast. However, if you want fast factorization, I recommend Pollard's Rho-Method - it's not guaranteed to succeed, but usually it does and it's pretty fast: Prelude Primality> factor (2^73-13) [23,337,238815641,5102337469] (42.75 secs, -2382486596 bytes) Prelude Primality> pollFact (2^73-13) [23,337,238815641,5102337469] (1.64 secs, 186986044 bytes) and easily implemented. But don't expect too much - factoring 2^137-1 took over 51 hours on my (old and slow) machine. For all who want it, the module is below. One odd thing, though. It was about 20% faster when compiled with ghc-6.2.2 rather than with ghc-6.4. Any Guru know why? If anybody has better algorithms, I'd be happy to know about them. Cheers, Daniel -- | This module provides tests for primality of Integers, safe and unsafe. -- Also some factorization methods and a process calculating all positive -- primes are given. module Primality ( -- * Primality tests -- ** Safe but slow isPrime, -- :: Integer -> Bool -- ** Unsafe but much faster isProbablePrime, -- :: Int -> Integer -> Bool isRatherProbablePrime, -- :: Integer -> Bool isVeryProbablePrime, -- :: Integer -> Bool -- ** Tests for (strong) pseudoprimality isPseudo, -- :: Integer -> Integer -> Bool isStrongPseudo, -- :: Integer -> Integer -> Bool -- * Factorization methods -- ** Safe but slow factor, -- :: Integer -> [Integer] -- ** Pollards rho-method, faster but it may fail pollFact, -- :: Integer -> [Integer] -- * List of primes primes -- :: [Integer] ) where import Data.List (insert) ---------------------------------------------------------------------- -- Exported functions -- ---------------------------------------------------------------------- -- A) Unsafe but relatively fast -- | This function checks if the second argument is a pseudoprime -- for the first argument (or indeed a prime). isPseudo :: Integer -> Integer -> Bool isPseudo a n = powerWithModulus n a (abs n - 1) == 1 -- | This function checks if the second argument is a prime or a -- strong pseudoprime for the first argument. isStrongPseudo :: Integer -> Integer -> Bool isStrongPseudo a n | n < 0 = isStrongPseudo a (-n) | n < 2 = False | n == 2 = True | even n = False | a < 0 = isStrongPseudo (-a) n | a < 2 = error "isStrongPseudo: illegitimate base" | otherwise = fst (checkStrong a n ((n-1) `quot` 2)) -- | The list of all positive primes. primes :: [Integer] primes = 2:3:5:7:11:13:17:19:23:[n | n <- [29,31 .. ], and [n `rem` p /= 0 | p <- takeWhile (squareLess n) primes]] -- | @isProbablePrime k n@ tests whether @n@ is (if not prime) a strong -- pseudoprime for the first @k@ primes. isProbablePrime :: Int -> Integer -> Bool isProbablePrime k n = and [isStrongPseudo p n | p <- take k (takeWhile (squareLess n) primes)] -- | Testing strong pseudoprimality for the first 200 primes. isRatherProbablePrime :: Integer -> Bool isRatherProbablePrime = isProbablePrime 200 -- | Testing for the first 2000 primes. isVeryProbablePrime :: Integer -> Bool isVeryProbablePrime = isProbablePrime 2000 -- B) Safe and slow -- | Safe primality test, based on trial division. isPrime :: Integer -> Bool isPrime n | n < 0 = isPrime (-n) | n < 2 = False | n < 4 = True | even n = False | otherwise = head (factor n) == n -- | Elementary factorization by trial division, tuned by wheel-technique. factor :: Integer -> [Integer] factor 0 = [0] factor 1 = [1] factor n | n < 0 = -1: factor' (-n) 0 smallPrimes | otherwise = factor' n 0 smallPrimes where factor' :: Integer -> Integer -> [Integer] -> [Integer] factor' 1 _ _ = [] factor' n w [] = let w' = wheelModulus + w in if w'^2 > n then [n] else factor' n w' $ map (w' +) wheelSettings factor' n w ps'@(p:ps) = case n `quotRem` p of (q,0) -> p:factor' q w ps' _ -> factor' n w ps -- C) Pollard's rho-method -- | Factorization by Pollard's rho-method after eliminating small -- prime factors. A \'veryProbablePrime\' is treated as prime, -- failure to factor a number known as composite is signalled by -- including a @1@ in the list of factors. pollFact :: Integer -> [Integer] pollFact n | n < 0 = (-1):pollFact (-n) | n == 0 = [0] | n == 1 = [] | even n = 2:pollFact (n `quot` 2) | otherwise = smallFacts n 0 smallPrimes ---------------------------------------------------------------------- -- Helper functions -- ---------------------------------------------------------------------- -- calculate a power with respect to a modulus, first tests and -- forwarding to another helper powerWithModulus :: Integer -> Integer -> Integer -> Integer powerWithModulus mo n k | mo < 0 = powerWithModulus (-mo) n k | mo == 0 = n^k | k < 0 = error "powerWithModulus: negative exponent" | k == 0 = 1 | k == 1 = n `mod` mo | mo == 1 = 0 | odd k = powerMod mo n n (k-1) | otherwise = powerMod mo 1 n k -- calculate the power with auxiliary value to account for -- odd exponents on the way, we have @mo >= 2@ and @k >= 1@ powerMod :: Integer -> Integer -> Integer -> Integer -> Integer powerMod mo aux val k | k == 1 = (aux * val) `mod` mo | odd k = ((powerMod mo $! ((aux*val) `mod` mo)) $! (val^2 `mod` mo)) $! (k `quot` 2) | even k = (powerMod mo aux $! (val^2 `mod` mo)) $! (k `quot` 2) -- If a prime @p@ divides @a^(2*e)-1@, @p@ divides exactly one of the -- factors @a^e+1@ and @a^e-1@. We check the analogous property for @n@. -- Say, @m = 2^u*v@ with odd @v@. We recursively check whether @n@ -- divides one of the factors @a^(2^j*v)+1@ for @0 <= j <= m@ or -- the factor @a^v-1@. Success is propagated without further calculation, -- failure also returns the value @a^(2^j*v) `rem` n@ to the caller. checkStrong :: Integer -> Integer -> Integer -> (Bool,Integer) checkStrong a n m | even m = case checkStrong a n (m `quot` 2) of (True,_) -> (True,0) (False,k) -> (b,r) where r = (k^2) `rem` n b = (r+1) `rem` n == 0 | otherwise = ((k-1) `rem` n == 0 || (k+1) `rem` n == 0, k) where k = powerWithModulus n a m -- condition used for various tests squareLess :: Integer -> Integer -> Bool squareLess n k = k^2 <= n ---------------------------------------------------------------------- -- Wheel Factoring -- ---------------------------------------------------------------------- -- wheel size is set so that factoring is relatively fast, -- but finding smallPrimes and wheelSettings does not take too long wheelSize :: Int wheelSize = 7 verySmallPrimes :: [Integer] verySmallPrimes = take wheelSize primes wheelModulus :: Integer wheelModulus = product verySmallPrimes smallPrimes :: [Integer] smallPrimes = takeWhile (<= wheelModulus) primes wheelSettings :: [Integer] wheelSettings = [f | f <- [1 .. wheelModulus-1], null [p | p <- verySmallPrimes, f `rem` p == 0]] ---------------------------------------------------------------------- -- Helpers for Pollard's rho-method -- ---------------------------------------------------------------------- -- small factors by wheel-factoring, then pass to rho-method smallFacts :: Integer -> Integer -> [Integer] -> [Integer] smallFacts 1 _ _ = [] smallFacts n w [] | w'^2 > n = [n] | w' > 5*10^6 = rhoFact n | otherwise = smallFacts n w' $ map (w' +) wheelSettings where w' = wheelModulus + w smallFacts n w ps'@(p:ps) = case n `quotRem` p of (q,0) -> p:smallFacts q w ps' _ -> smallFacts n w ps -- Factorization by the rho-method, first we try the function -- @\x -> x^2+1@ with starting value @x1 = 2@. If that fails, -- we pass the number to the rho-method using @\x -> x^2+3@. -- If two nontrivial factors are found, they are factored and -- the lists of factors are merged. rhoFact :: Integer -> [Integer] rhoFact n | isVeryProbablePrime n = [n] | otherwise = case rho1 n of [k] -> rhoKFact 1 k [a,b] -> merge (rhoFact a) (rhoFact b) -- Factorization using @\x -> x^2+2*k+1@, if that fails, we try the -- next k until we give up when @k > 100@. rhoKFact :: Integer -> Integer -> [Integer] rhoKFact k n | isVeryProbablePrime n = [n] | k > 100 = [1,n] | otherwise = case rhoK k n of [m] -> rhoKFact (k+1) m [a,b] -> merge (rhoFact a) (rhoFact b) -- start the search for factors with @x1 = 2@ and @x2 = 5@ rho1 :: Integer -> [Integer] rho1 m = find1 m 2 5 -- With @x = xk@ and @y = x2k@, if @m@ divides @x2k-xk@, we can't find -- a nontrivial factor. If @m@ and @x2k-xk@ are coprime, we check the -- next k, otherwise we have found two nontrivial factors. find1 :: Integer -> Integer -> Integer -> [Integer] find1 m x y | g == m = [m] | g > 1 = [g, m `quot` g] | otherwise = find1 m (s1 x) (s1 (s1 y)) where g = gcd m (y-x) s1 :: Integer -> Integer s1 x = (x^2+1) `rem` m -- merge two sorted lists to one sorted list merge :: [Integer] -> [Integer] -> [Integer] merge = foldr insert -- start the search for factors using @\x -> x^2+2*k+1@ with @x1 = 2@ rhoK :: Integer -> Integer -> [Integer] rhoK k n = kFind s n 2 (4+s) where s = 2*k+1 -- analogous to find1 kFind :: Integer -> Integer -> Integer -> Integer -> [Integer] kFind k m x y | g == m = [m] | g > 1 = [g, m `quot` g] | otherwise = kFind k m (s x) (s (s y)) where g = gcd m (y-x) s :: Integer -> Integer s x = (x^2+k) `rem` m