[Haskell-cafe] NumberTheory library

[Yitzchak Gale]

> As promised, here's the first attempt:
> 
>         darcs get http://andrew.bromage.org/darcs/numbertheory/

A few small comments about the function "factor" in
Prime.hs:

o I think you are testing w' * w' < n each time, even
   when you are repeating factors of the same prime p.
   You only need to do that when you move to the next p.

o You can use divMod (or quotRem) instead of using
   div and then multiplying. (That's an ancient trick
   from C. It may not make a difference on modern
   CPUs.)

o You are using the old functional programming
   trick of artificially creating state with
   extra function parameters. Nowadays we don't
   need that in Haskell - it is much simpler
   just to use the State monad.

o I personally would much prefer a function
   that only returns prime numbers, not -1, 0,
   or 1.

o I think that in most of the places where you
   fix the type as Int or Integer, you could leave
   it polymorphic and avoid a lot of coercing.

Below is a version of "factor" that
implements the above (except the last two,
so that I preserve your interface).

Regards,
Yitz

import Control.Monad.State

-- ...etc.

-- An infinite list of possible prime factors
wheel :: Integral a => [a]
wheel = iVerySmallPrimes ++
         [f + d | d <- [0,iWheelModulus..], f <- wheelSettings]
   where
     iWheelModulus = fromIntegral wheelModulus
     iVerySmallPrimes = map fromIntegral verySmallPrimes
     wheelSettings = [f | f <- [z,z+2..iWheelModulus+1],
              null [p | p <- iVerySmallPrimes, f `mod` p == 0]]
     z = last iVerySmallPrimes + 2

-- |Factorize a number.
factor :: Integral a => a -> [a]
factor 0     = [0]
factor 1     = [1]
factor n
  | n < 0     = -1 : factor (-n)
  | otherwise = evalState (factorS n) wheel

-- Factorize n with wheel as state
factorS :: Integral a => a -> State [a] [a]
factorS 1 = return []
factorS n = do
      p <- gets head
      if p * p > n
        then return [n]
        else factorPS n p

-- Factorize n at a given prime p
factorPS :: Integral a => a -> a -> State [a] [a]
factorPS n p
  | rem == 0  = do
      ps <- factorPS quot p
      return (p:ps)
  | otherwise = do
      modify tail
      factorS n
  where (quot, rem) = n `quotRem` p