[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