primzahlen = sieb [2..]
   where sieb (p:xs) = p : sieb [n | n <- xs, n `mod` p > 0]

fact n (p:ys)  | p*p > n = [n]
               | n `mod` p == 0  = p : fact (n `div` p) (p:ys)
               | otherwise       = fact n ys
primfaktor 1 = []
primfaktor n = fact n primzahlen