factorising to prime numbers

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? 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) -- Dougal Stanton

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.
So you need finite lists.
-- 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
Extremely high-level hint: Move from div to divMod (i.e., let f return a pair), and add List.unfoldr (or takeWhile) capacities to map. ;-) Wolfram

On Feb 9, 2007, at 9:20 AM, 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?
If I've understood correctly your list 'powers' will be all zeros after a certain point. Once that happens, you don't need to examine that part of the list anymore. This should at least occur as soon as the primes become larger than your number N (and probably sooner. sqrt(N) maybe? I forget). So, you should be able to only examine a prefix of the list 'primes'. The definition you have looks right, in that it correctly generates the correct list. If you want to test that its doing the right thing, you can just examine the prefix:
test n = product (zipWith (^) (takeWhile (
(untested, but I think it would work). or you can just create the portion of the powers list you need in the first place:
powersPrefix n = map (f n) (takeWhile (
(remember kids, a decidable problem is a semi-decidable problem where we can calculate a stopping condition).
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)
-- Dougal Stanton
Rob Dockins Speak softly and drive a Sherman tank. Laugh hard; it's a long way to the bank. -- TMBG

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@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

G'day all.
Quoting Dougal Stanton
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, ...]
It depends what you mean by "better". You can make it more concise by using zip or zipWith: product (zipWith (^) ps cs) If you're after performance, that's another issue entirely. For large integers, you're using entirely the wrong algorithm. Cheers, Andrew Bromage
participants (5)
-
ajb@spamcop.net
-
Dougal Stanton
-
kahl@cas.mcmaster.ca
-
Mikael Johansson
-
Robert Dockins