[Haskell-cafe] NumberTheory library
Daniel Fischer
daniel.is.fischer at web.de
Thu May 5 19:15:44 EDT 2005
Am Donnerstag, 5. Mai 2005 06:20 schrieb ajb at spamcop.net:
> G'day all.
>
> Quoting Bo Herlin <bo at gcab.net>:
> > 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 at . We check the analogous property for @n at .
-- Say, @m = 2^u*v@ with odd @v at . 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 at . 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 at . If that fails,
-- we pass the number to the rho-method using @\x -> x^2+3 at .
-- 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 at .
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
More information about the Haskell-Cafe
mailing list