module MathObj.FactoredRational
(
T
, factorial
, eulerPhi
, divisors
) where
import qualified Algebra.Additive as Additive
import qualified Algebra.Ring as Ring
import qualified Algebra.Field as Field
import qualified Algebra.IntegralDomain as Integral
import qualified Algebra.ZeroTestable as ZeroTestable
import qualified Algebra.Real as Real
import qualified Algebra.ToRational as ToRational
import qualified Algebra.RealIntegral as RealIntegral
import qualified Algebra.ToInteger as ToInteger
import Data.Numbers.Primes
import PreludeBase
import NumericPrelude
data T = FQZero
| FQ Bool [Integer]
instance Show T where
show FQZero = "0"
show (FQ True pows) = "(-1)" ++ showPows pows
show (FQ False pows) = showPows pows
showPows :: [Integer] -> String
showPows pows = concat $ zipWith showPow primes pows
where showPow p 0 = ""
showPow p 1 = "(" ++ show p ++ ")"
showPow p n = "(" ++ show p ++ "^" ++ show n ++ ")"
instance Additive.C T where
zero = FQZero
FQZero + a = a
a + FQZero = a
x + y = fromRational' (toRational x + toRational y)
negate FQZero = FQZero
negate (FQ s e) = FQ (not s) e
instance Ring.C T where
FQZero * _ = FQZero
_ * FQZero = FQZero
(FQ s1 e1) * (FQ s2 e2) = FQ (s1 /= s2) (zipWithExt 0 0 (+) e1 e2)
fromInteger 0 = FQZero
fromInteger n | n < 0 = FQ True (factor (negate n))
| otherwise = FQ False (factor n)
_ ^ 0 = one
FQZero ^ _ = FQZero
(FQ s e) ^ n = FQ s (map (*n) e)
zipWithExt :: a -> b -> (a -> b -> c) -> [a] -> [b] -> [c]
zipWithExt da db f = zipWithExt'
where zipWithExt' [] bs = zipWith f (repeat da) bs
zipWithExt' as [] = zipWith f as (repeat db)
zipWithExt' (a:as) (b:bs) = f a b : zipWithExt' as bs
factor :: Integer -> [Integer]
factor n = factor' n primes
where
factor' 1 _ = []
factor' n (p:ps) = let (k,n') = logRem n p
in k : factor' n' ps
logRem :: Integer -> Integer -> (Integer, Integer)
logRem = logRem' 0
where logRem' k n p | n `mod` p == 0 = logRem' (k+1) (n `div` p) p
| otherwise = (k,n)
instance ZeroTestable.C T where
isZero FQZero = True
isZero _ = False
instance Eq T where
FQZero == FQZero = True
FQZero == (FQ _ _) = False
(FQ _ _) == FQZero = False
(FQ s1 e1) == (FQ s2 e2) = s1 == s2 && dropZeros e1 == dropZeros e2
where dropZeros = reverse . dropWhile (==0) . reverse
instance Ord T where
compare FQZero FQZero = EQ
compare FQZero (FQ False _) = LT
compare FQZero (FQ True _) = GT
compare (FQ False _) FQZero = GT
compare (FQ True _) FQZero = LT
compare (FQ False _) (FQ True _) = GT
compare (FQ True _) (FQ False _) = LT
compare fq1 fq2 = compare (toRational fq1) (toRational fq2)
instance Real.C T where
abs FQZero = FQZero
abs (FQ _ e) = FQ False e
instance ToRational.C T where
toRational FQZero = 0
toRational (FQ s e) = (if s then negate else id)
. product
. zipWith (^-) (map (%1) primes)
$ e
instance Integral.C T where
divMod a b =
if isZero b
then (undefined,a)
else (a/b,0)
instance RealIntegral.C T
instance ToInteger.C T where
toInteger FQZero = 0
toInteger (FQ s e) | any (<0) e = error "non-integer in FactoredRational.toInteger"
| otherwise = (if s then negate else id)
. product
. zipWith (^) primes
$ e
instance Field.C T where
recip FQZero = error "division by zero"
recip (FQ s e) = FQ s (map negate e)
factorial :: Integer -> T
factorial 0 = one
factorial 1 = one
factorial n = FQ False (takeWhile (>0) . map (factorialFactors n) $ primes)
factorialFactors :: Integer -> Integer -> Integer
factorialFactors n p = sum
. takeWhile (>0)
. map (n `div`)
$ iterate (*p) p
eulerPhi :: T -> Integer
eulerPhi FQZero = 1
eulerPhi (FQ _ pows) = product $ zipWith phiP primes pows
where phiP _ 0 = 1
phiP p a = p^(a1) * (p1)
divisors :: T -> [T]
divisors FQZero = [1]
divisors (FQ b pows) = map (FQ b) $ mapM (enumFromTo 0) pows