[Haskell-cafe] factorising to prime numbers
Mikael Johansson
mikael at johanssons.org
Fri Feb 9 10:00:58 EST 2007
On Fri, 9 Feb 2007, Dougal Stanton wrote:
> Hi folks,
>
> I recently read in my copy of Concrete Mathematics the relationship
> between prime factors powers and lcm/gcd functions. So I decided to
> reimplement gcd and lcm the long way, for no other reason than because
> I could.
>
> If you look at the definition of 'powers' you'll note it's infinite. So
> there's no easy way to take the product of this list, if I don't know
> how many items to take from it.
>
> Is there a better way to turn an integer N and a list of primes
> [p1,p2,p3,...] into powers [c1,c2,c3,...] such that
>
> N = product [p1^c1, p2^c2, p3^c3, ...]
>
> If I'm missing something really obvious I'll be very grateful. I can't
> really work out what kind of structure it should be. A map? fold?
>
So... Thinking out loud here. You'll want to do something that's sort of
like a fold is my feeling. A first naïve implementation, taking care of
all the things I think you should watch for would be
powers n = let
powersAcc acc 1 _ = acc
powersAcc acc N (p:ps) =
powersAcc (acc ++ [f N p]) (N`div`(p^(f N p)) ps
in powersAcc [] n primes
But this should, is my feeling, be possible to do using some sufficiently
funky function in the middle of a fold or other.
Basically, the core function will eat one prime off the list, do things
with it, and spit out one more component of a new list. So we'll want
powers to have the signature
Integer -> [Integer]
or even, taking the implied primes into account, we'll want
Integer -> [Integer] -> [Integer]
such that it's inverse to
(foldr (*) 1) . (zipWith (^))
If we, thus, want to do this with a foldr of some sort, we need to match
the foldr signature
foldr :: (a -> b -> b) -> b -> [a] -> b
so that we process the [Integer] containing our primes, and spit out a
[Integer] containing the powers.
So a = Integer, and b = [Integer], and we get our foldr signature
foldr :: (Integer -> [Integer] -> [Integer]) -> [Integer] -> [Integer] ->
[Integer]
and the trick would seem to lie in writing this
[Integer] -> [Integer] -> [Integer].
Now, if we use the first element of a partially constructed list of
exponents to signify what's left to work out, we could write something
along the lines of
nextStep :: Integer -> [Integer] -> [Integer] -- I know, bad name...
nextStep p (n:exps) =
let c = f n p
n' = n `div` (p^c)
in n':exps ++ [c]
and we use this by
powers n = foldr nextStep [n] primes
Of course, this also suffers from halting problems, just as the map in the
original code. To make this really work out, we'd want some fold that can
be told that after this point, nothing will ever change the accumulated
value.
Maybe this is a point where foldM needs to be used, maybe with the Maybe
monad, in order to terminate the fold at some point.
> D.
>
>
> -- Concrete Mathematics
> -- Graham, Knuth & Patashnuk
>
> module Concrete where
>
> import Data.List
>
> -- the sieve of eratosthenes is a fairly simple way
> -- to create a list of prime numbers
> primes =
> let primes' (n:ns) = n : primes' (filter (\v -> v `mod` n /= 0) ns)
> in primes' [2..]
>
> -- how many of the prime p are in the unique factorisation
> -- of the integer n?
> f 0 _ = 0
> f n p | n `mod` p == 0 = 1 + f (n `div` p) p
> | otherwise = 0
>
> powers n = map (f n) primes
>
> --gcd :: Integer -> Integer -> Integer
> --gcd = f . map (uncurry min)
>
>
--
Mikael Johansson | To see the world in a grain of sand
mikael at johanssons.org | And heaven in a wild flower
http://www.mikael.johanssons.org | To hold infinity in the palm of your hand
| And eternity for an hour
More information about the Haskell-Cafe
mailing list