{------------------------------------------------------------------------
            Factorizing Gaussian integers using Gofer               
------------------------------------------------------------------------}
factorise :: Gauss Int -> String
factorise = show . factors 

factors :: (Gauss Int) -> [Gauss Int]
factors z | z == zero  = []
              | norm z == 1  = [z]
              | otherwise    = w:factors (z/w) where
                 w | u `divides` z = u 
                    | otherwise     = u'
                           where
                            u:u':_ = irreducible p
                            p:_    = prime_divisors (norm z) 

prime_divisors :: Int -> [Int]
prime_divisors x = f (abs x) where
  f (2*n) = 2:f n                             -- deal with even case
  f (3*n) = 3:f n
  f (5*n) = 5:f n
  f (7*n) = 7:f n
  f n | n<2       = []
       | otherwise = d:f (n/d) where
            d                         = maybe 11
            maybe k | n `mod` k == 0  = k
                          | k*k > n         = n            -- n has to be prime
                          | otherwise       = maybe (k+2)  -- next odd candidate


sum_of_squares :: Int -> [Gauss Int]
sum_of_squares n = [ (a:+!b) | a<-[1..m], b<-[1..(n-a*a)], n == a*a+b*b]
            where m:_ = [ k-1 | k <- [1..], k*k > n ]

irreducible :: Int -> [Gauss Int]
irreducible p | p == 2         = [unit + i, unit - i]
                    | p `rem` 4 == 3 = p +! 0 
                    | otherwise      = sum_of_squares p

divides :: (Gauss Int) -> (Gauss Int) -> Bool
z `divides` z' = a `mod` d == 0 && b `mod` d == 0 where  
                a :+! b  = z' * (conjugate z)                   
                d       = norm z

----- End 
