{-# LANGUAGE DataKinds, TypeSynonymInstances, FlexibleContexts, FlexibleInstances, BangPatterns, ScopedTypeVariables #-}
module Math.Algebra.Polynomial.Univariate.Cyclotomic
( cyclotomic
, cyclotomicMoebius
, cyclotomicNaive
)
where
import Data.Ratio
import Data.Semigroup
import Data.Monoid
import qualified Math.Algebra.Polynomial.FreeModule as ZMod
import Math.Algebra.Polynomial.FreeModule ( FreeMod , FreeModule(..) , ZMod , QMod )
import Math.Algebra.Polynomial.Univariate
import Math.Algebra.Polynomial.Class
import Math.Algebra.Polynomial.Pretty
import Math.Algebra.Polynomial.Misc
test_product :: Int -> Bool
test_product n = lhs == rhs where
lhs = product [ cyclotomic d | d <- divisors n ]
rhs = Uni $ ZMod.fromList [ (U n , 1) , (U 0 , -1) ]
cyclotomic :: Int -> Univariate Integer "x"
cyclotomic = cyclotomicMoebius
cyclotomicMoebius :: Int -> Univariate Integer "x"
cyclotomicMoebius = Uni . ZMod.unsafeZModFromQMod . unUni . cyclotomicMoebiusQ
cyclotomicMoebiusQ :: Int -> Univariate Rational "x"
cyclotomicMoebiusQ n = unsafeDivide (Uni numer) (Uni denom) where
sqfd = squareFreeDivisors n
numer = ZMod.product [ term d | (d,Plus ) <- sqfd ]
denom = ZMod.product [ term d | (d,Minus) <- sqfd ]
term d = ZMod.fromList [ (U (div n d) , 1) , (U 0 , -1) ]
polynomialLongDivision :: forall p. (Polynomial p, Fractional (CoeffP p)) => p -> p -> (p,p)
polynomialLongDivision p0 q = go zeroP p0 where
(bq,cq) = case findMaxTerm q of
Just bc -> bc
Nothing -> error "polynomialLongDivision: division by zero"
go !acc !p = case findMaxTerm p of
Nothing -> (acc,zeroP)
Just (bp,cp) -> case divM bp bq of
Nothing -> (acc,p)
Just br -> let cr = (cp/cq)
u = scaleP cr (mulByMonomP br q)
p' = p - u
acc' = (acc + monomP' br cr)
in go acc' p'
divideMaybe :: (Polynomial p, Fractional (CoeffP p)) => p -> p -> Maybe p
divideMaybe p q = case polynomialLongDivision p q of
(s,r) -> if isZeroP r then Just s else Nothing
unsafeDivide :: (Polynomial p, Fractional (CoeffP p)) => p -> p -> p
unsafeDivide p q = case divideMaybe p q of
Just s -> s
Nothing -> error "Polynomial/unsafeDivide: not divisible"
findMaxTerm :: FreeModule f => f -> Maybe (BaseF f, CoeffF f)
findMaxTerm = ZMod.findMaxTerm . toFreeModule
newtype E2PiI = E2PiI Rational deriving (Eq,Ord,Show)
instance Pretty E2PiI where
pretty (E2PiI y) = "e^{2*pi*i*" ++ show p ++ "/" ++ show q ++ "}" where
p = numerator y
q = denominator y
mod1 :: Rational -> Rational
mod1 x = x - fromInteger (floor x)
mulE2PiI :: E2PiI -> E2PiI -> E2PiI
mulE2PiI (E2PiI x) (E2PiI y) = E2PiI (mod1 $ x+y)
instance Semigroup E2PiI where
(<>) = mulE2PiI
instance Monoid E2PiI where
mempty = E2PiI 0
mappend = mulE2PiI
half :: Rational
half = 1/2
reduce :: ZMod E2PiI -> Either (ZMod E2PiI) Integer
reduce input =
case ZMod.findMaxTerm input of
Nothing -> Right 0
Just (E2PiI y, c)
| y == 0 -> Right c
| otherwise -> case minimalZeroSumCircle y of
Just circle -> reduce $ ZMod.sub input (ZMod.scale c circle)
Nothing -> Left input
minimalZeroSumCircle :: Rational -> Maybe (ZMod E2PiI)
minimalZeroSumCircle y
| y < half = Nothing
| r == 1 = Just $ ZMod.fromList [ (E2PiI (i % q) , 1) | i<-[0..q-1] ]
| otherwise = Nothing
where
p = numerator y
q = denominator y
r = q - p
type X = U "x"
type Pre = ZMod E2PiI
instance IsSigned Pre where
signOf _ = Just Plus
preCyclotomic :: Int -> FreeMod Pre X
preCyclotomic n = ZMod.product terms where
terms = [ ZMod.sub x (ZMod.konst (e k)) | k <- [1..n] , gcd k n == 1 ] where
x = ZMod.generator (U 1)
qn = fromIntegral n :: Rational
e k = ZMod.generator $ E2PiI (mod1 $ fromIntegral k / qn) :: ZMod E2PiI
cyclotomicNaive :: Int -> Univariate Integer "x"
cyclotomicNaive = Uni . ZMod.mapCoeff fun . preCyclotomic where
fun term = case reduce term of
Right c -> c
Left {} -> error "cyclotomicNaive: shouldn't happen"