diff --git a/Math/NumberTheory/DirichletCharacters.hs b/Math/NumberTheory/DirichletCharacters.hs new file mode 100644 index 000000000..f4896ef72 --- /dev/null +++ b/Math/NumberTheory/DirichletCharacters.hs @@ -0,0 +1,614 @@ +-- | +-- Module: Math.NumberTheory.DirichletCharacters +-- Copyright: (c) 2018 Bhavik Mehta +-- Licence: MIT +-- Maintainer: Bhavik Mehta +-- +-- Implementation and enumeration of Dirichlet characters. +-- + +{-# LANGUAGE CPP #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE LambdaCase #-} +{-# LANGUAGE PatternSynonyms #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE ViewPatterns #-} + +module Math.NumberTheory.DirichletCharacters + ( + -- * Roots of unity + RootOfUnity + -- ** Conversions + , toRootOfUnity + , fromRootOfUnity + , toComplex + -- * An absorbing semigroup + , OrZero, pattern Zero, pattern NonZero + , orZeroToNum + -- * Dirichlet characters + , DirichletCharacter + -- ** Construction + , indexToChar + , indicesToChars + , characterNumber + , allChars + , fromTable + -- ** Evaluation + , eval + , evalGeneral + , evalAll + -- ** Special Dirichlet characters + , principalChar + , isPrincipal + , orderChar + -- ** Real Dirichlet characters + , RealCharacter + , isRealCharacter + , getRealChar + , toRealFunction + , jacobiCharacter + -- ** Primitive characters + , PrimitiveCharacter + , isPrimitive + , getPrimitiveChar + , induced + , makePrimitive + , WithNat(..) + -- * Debugging + , validChar + ) where + +#if !MIN_VERSION_base(4,12,0) +import Control.Applicative (liftA2) +#endif +import Data.Bits (Bits(..)) +import Data.Complex (Complex(..), cis) +import Data.Foldable (for_) +import Data.Functor.Identity (Identity(..)) +import Data.List (mapAccumL, foldl', sort, find, unfoldr) +import Data.Maybe (mapMaybe, fromJust, fromMaybe) +#if MIN_VERSION_base(4,12,0) +import Data.Monoid (Ap(..)) +#endif +import Data.Proxy (Proxy(..)) +import Data.Ratio ((%), numerator, denominator) +import Data.Semigroup (Semigroup(..), Product(..)) +import qualified Data.Vector as V +import qualified Data.Vector.Mutable as MV +import Data.Vector (Vector, (!)) +import GHC.TypeNats.Compat (Nat, SomeNat(..), natVal, someNatVal) +import Numeric.Natural (Natural) + +import Math.NumberTheory.ArithmeticFunctions (totient) +import Math.NumberTheory.Moduli.Chinese +import Math.NumberTheory.Moduli.Class (KnownNat, Mod, getVal) +import Math.NumberTheory.Moduli.Internal (isPrimitiveRoot', discreteLogarithmPP) +import Math.NumberTheory.Moduli.Multiplicative (MultMod(..), isMultElement) +import Math.NumberTheory.Moduli.Singleton (Some(..), cyclicGroupFromFactors) +import Math.NumberTheory.Powers.Modular (powMod) +import Math.NumberTheory.Primes (Prime(..), UniqueFactorisation, factorise, nextPrime) +import Math.NumberTheory.Utils.FromIntegral (wordToInt) +import Math.NumberTheory.Utils + +-- | A Dirichlet character mod \(n\) is a group homomorphism from \((\mathbb{Z}/n\mathbb{Z})^*\) +-- to \(\mathbb{C}^*\), represented abstractly by `DirichletCharacter`. In particular, they take +-- values at roots of unity and can be evaluated using `eval`. +-- A Dirichlet character can be extended to a completely multiplicative function on \(\mathbb{Z}\) +-- by assigning the value 0 for \(a\) sharing a common factor with \(n\), using `evalGeneral`. +-- +-- There are finitely many possible Dirichlet characters for a given modulus, in particular there +-- are \(\phi(n)\) characters modulo \(n\), where \(\phi\) refers to Euler's `totient` function. +-- This gives rise to `Enum` and `Bounded` instances. +newtype DirichletCharacter (n :: Nat) = Generated [DirichletFactor] + +-- | The group (Z/nZ)^* decomposes to a product (Z/2^k0 Z)^* x (Z/p1^k1 Z)^* x ... x (Z/pi^ki Z)^* +-- where n = 2^k0 p1^k1 ... pi^ki, and the pj are odd primes, k0 possibly 0. Thus, a group +-- homomorphism from (Z/nZ)^* is characterised by group homomorphisms from each of these factor +-- groups. Furthermore, for odd p, we have (Z/p^k Z)^* isomorphic to Z / p^(k-1)*(p-1) Z, an +-- additive group, where an isomorphism is specified by a choice of primitive root. +-- Similarly, for k >= 2, (Z/2^k Z)^* is isomorphic to Z/2Z * (Z / 2^(k-2) Z) (and for k < 2 +-- it is trivial). (See `lambda` for this isomorphism). +-- Thus, to specify a Dirichlet character, it suffices to specify the value of generators +-- of each of these cyclic groups, when primitive roots are given. This data is given by a +-- DirichletFactor. +-- We have the invariant that the factors must be given in strictly increasing order, and the +-- generator is as given by `generator`, and are each non-trivial. These conditions are verified +-- using `validChar`. +data DirichletFactor = OddPrime { _getPrime :: Prime Natural + , _getPower :: Word + , _getGenerator :: Natural + , _getValue :: RootOfUnity + } + | TwoPower { _getPower2 :: Int -- this ought to be Word, but many applications + -- needed to use wordToInt, so Int is cleaner + -- Required to be >= 2 + , _getFirstValue :: RootOfUnity + , _getSecondValue :: RootOfUnity + } + | Two + +instance Eq (DirichletCharacter n) where + Generated a == Generated b = a == b + +instance Eq DirichletFactor where + TwoPower _ x1 x2 == TwoPower _ y1 y2 = x1 == y1 && x2 == y2 + OddPrime _ _ _ x == OddPrime _ _ _ y = x == y + Two == Two = True + _ == _ = False + +-- | A representation of : complex +-- numbers \(z\) for which there is \(n\) such that \(z^n=1\). +newtype RootOfUnity = + RootOfUnity { -- | Every root of unity can be expressed as \(e^{2 \pi i q}\) for some + -- rational \(q\) satisfying \(0 \leq q < 1\), this function extracts it. + fromRootOfUnity :: Rational } + deriving (Eq) + +instance Show RootOfUnity where + show (RootOfUnity q) + | n == 0 = "1" + | d == 1 = "-1" + | n == 1 = "e^(πi/" ++ show d ++ ")" + | otherwise = "e^(" ++ show n ++ "πi/" ++ show d ++ ")" + where n = numerator (2*q) + d = denominator (2*q) + +-- | Given a rational \(q\), produce the root of unity \(e^{2 \pi i q}\). +toRootOfUnity :: Rational -> RootOfUnity +toRootOfUnity q = RootOfUnity ((n `rem` d) % d) + where n = numerator q + d = denominator q + -- effectively q `mod` 1 + -- This smart constructor ensures that the rational is always in the range 0 <= q < 1. + +-- | This Semigroup is in fact a group, so @'stimes'@ can be called with a negative first argument. +instance Semigroup RootOfUnity where + RootOfUnity q1 <> RootOfUnity q2 = toRootOfUnity (q1 + q2) + stimes k (RootOfUnity q) = toRootOfUnity (q * fromIntegral k) + +instance Monoid RootOfUnity where + mappend = (<>) + mempty = RootOfUnity 0 + +-- | Convert a root of unity into an inexact complex number. Due to floating point inaccuracies, +-- it is recommended to avoid use of this until the end of a calculation. Alternatively, with +-- the [cyclotomic](http://hackage.haskell.org/package/cyclotomic-0.5.1) package, one can use +-- @[polarRat](https://hackage.haskell.org/package/cyclotomic-0.5.1/docs/Data-Complex-Cyclotomic.html#v:polarRat) +-- 1 . @'fromRootOfUnity' to convert to a cyclotomic number. +toComplex :: Floating a => RootOfUnity -> Complex a +toComplex (RootOfUnity t) + | t == 1/2 = (-1) :+ 0 + | t == 1/4 = 0 :+ 1 + | t == 3/4 = 0 :+ (-1) + | otherwise = cis . (2*pi*) . fromRational $ t + +-- | For primes, define the canonical primitive root as the smallest such. For prime powers \(p^k\), +-- either the smallest primitive root \(g\) mod \(p\) works, or \(g+p\) works. +generator :: (Integral a, UniqueFactorisation a) => Prime a -> Word -> a +generator p k + | k == 1 = modP + | otherwise = if powMod modP (p'-1) (p'*p') == 1 then modP + p' else modP + where p' = unPrime p + modP = case cyclicGroupFromFactors [(p,k)] of + Just (Some cg) -> head $ filter (isPrimitiveRoot' cg) [2..p'-1] + _ -> error "illegal" + +-- | Implement the function \(\lambda\) from page 5 of +-- https://www2.eecs.berkeley.edu/Pubs/TechRpts/1984/CSD-84-186.pdf +lambda :: Integer -> Int -> Integer +lambda x e = ((powMod x (2*modulus) largeMod - 1) `shiftR` (e+1)) .&. (modulus - 1) + where modulus = bit (e-2) + largeMod = bit (2*e - 1) + +-- | For elements of the multiplicative group \((\mathbb{Z}/n\mathbb{Z})^*\), a Dirichlet +-- character evaluates to a root of unity. +eval :: DirichletCharacter n -> MultMod n -> RootOfUnity +eval (Generated ds) m = foldMap (evalFactor m') ds + where m' = getVal $ multElement m + +-- | Evaluate each factor of the Dirichlet character. +evalFactor :: Integer -> DirichletFactor -> RootOfUnity +evalFactor m = + \case + OddPrime (toInteger . unPrime -> p) k (toInteger -> a) b -> + discreteLogarithmPP p k a (m `rem` p^k) `stimes` b + TwoPower k s b -> (if testBit m 1 then s else mempty) + <> lambda (thingy k m) k `stimes` b + Two -> mempty + +thingy :: (Bits p, Num p) => Int -> p -> p +thingy k m = if testBit m 1 + then bit k - m' + else m' + where m' = m .&. (bit k - 1) + +-- | A character can evaluate to a root of unity or zero: represented by @Nothing@. +evalGeneral :: KnownNat n => DirichletCharacter n -> Mod n -> OrZero RootOfUnity +evalGeneral chi t = case isMultElement t of + Nothing -> Zero + Just x -> NonZero $ eval chi x + +-- | Give the principal character for this modulus: a principal character mod \(n\) is 1 for +-- \(a\) coprime to \(n\), and 0 otherwise. +principalChar :: KnownNat n => DirichletCharacter n +principalChar = minBound + +mulChars :: DirichletCharacter n -> DirichletCharacter n -> DirichletCharacter n +mulChars (Generated x) (Generated y) = Generated (zipWith combine x y) + where combine :: DirichletFactor -> DirichletFactor -> DirichletFactor + combine Two Two = Two + combine (OddPrime p k g n) (OddPrime _ _ _ m) = + OddPrime p k g (n <> m) + combine (TwoPower k a n) (TwoPower _ b m) = + TwoPower k (a <> b) (n <> m) + combine _ _ = error "internal error: malformed DirichletCharacter" + +-- | This Semigroup is in fact a group, so @stimes@ can be called with a negative first argument. +instance Semigroup (DirichletCharacter n) where + (<>) = mulChars + stimes = stimesChar + +instance KnownNat n => Monoid (DirichletCharacter n) where + mempty = principalChar + mappend = (<>) + +stimesChar :: Integral a => a -> DirichletCharacter n -> DirichletCharacter n +stimesChar s (Generated xs) = Generated (map mult xs) + where mult :: DirichletFactor -> DirichletFactor + mult (OddPrime p k g n) = OddPrime p k g (s `stimes` n) + mult (TwoPower k a b) = TwoPower k (s `stimes` a) (s `stimes` b) + mult Two = Two + +-- | We define `succ` and `pred` with more efficient implementations than +-- @`toEnum` . (+1) . `fromEnum`@. +instance KnownNat n => Enum (DirichletCharacter n) where + toEnum = indexToChar . fromIntegral + fromEnum = fromIntegral . characterNumber + succ x = makeChar x (characterNumber x + 1) + pred x = makeChar x (characterNumber x - 1) + + enumFromTo x y = bulkMakeChars x [fromEnum x..fromEnum y] + enumFrom x = bulkMakeChars x [fromEnum x..] + enumFromThenTo x y z = bulkMakeChars x [fromEnum x, fromEnum y..fromEnum z] + enumFromThen x y = bulkMakeChars x [fromEnum x, fromEnum y..] + +instance KnownNat n => Bounded (DirichletCharacter n) where + minBound = indexToChar 0 + maxBound = indexToChar (totient n - 1) + where n = natVal (Proxy :: Proxy n) + +-- | We have a (non-canonical) enumeration of dirichlet characters. +characterNumber :: DirichletCharacter n -> Integer +characterNumber (Generated y) = foldl' go 0 y + where go x (OddPrime p k _ a) = x * m + numerator (fromRootOfUnity a * fromIntegral m) + where p' = fromIntegral (unPrime p) + m = p'^(k-1)*(p'-1) + go x (TwoPower k a b) = x' * 2 + numerator (fromRootOfUnity a * 2) + where m = bit (k-2) :: Integer + x' = x `shiftL` (k-2) + numerator (fromRootOfUnity b * fromIntegral m) + go x Two = x + +-- | Give the dirichlet character from its number. +-- Inverse of `characterNumber`. +indexToChar :: forall n. KnownNat n => Natural -> DirichletCharacter n +indexToChar = runIdentity . indicesToChars . Identity + +-- | Give a collection of dirichlet characters from their numbers. This may be more efficient than +-- `indexToChar` for multiple characters, as it prevents some internal recalculations. +indicesToChars :: forall n f. (KnownNat n, Functor f) => f Natural -> f (DirichletCharacter n) +indicesToChars = fmap (Generated . unroll t . (`mod` m)) + where n = natVal (Proxy :: Proxy n) + (Product m, t) = mkTemplate n + +-- | List all characters for the modulus. This is preferred to using @[minBound..maxBound]@. +allChars :: forall n. KnownNat n => [DirichletCharacter n] +allChars = indicesToChars [0..m-1] + where m = totient $ natVal (Proxy :: Proxy n) + +-- | The same as `indexToChar`, but if we're given a character we can create others more efficiently. +makeChar :: Integral a => DirichletCharacter n -> a -> DirichletCharacter n +makeChar x = runIdentity . bulkMakeChars x . Identity + +-- | Use one character to make many more: better than indicesToChars since it avoids recalculating +-- some primitive roots +bulkMakeChars :: (Integral a, Functor f) => DirichletCharacter n -> f a -> f (DirichletCharacter n) +bulkMakeChars x = fmap (Generated . unroll t . (`mod` m) . fromIntegral) + where (Product m, t) = templateFromCharacter x + +-- We assign each natural a unique Template, which can be decorated (eg in `unroll`) to +-- form a DirichletCharacter. A Template effectively holds the information carried around +-- in a DirichletFactor which depends only on the modulus of the character. +data Template = OddTemplate { _getPrime' :: Prime Natural + , _getPower' :: Word + , _getGenerator' :: !Natural + , _getModulus' :: !Natural + } + | TwoPTemplate { _getPower2' :: Int + , _getModulus' :: !Natural + } -- the modulus is derivable from the other values, but calculation + -- may be expensive, so we pre-calculate it + -- morally getModulus should be a prefactored but seems to be + -- pointless here + | TwoTemplate + +templateFromCharacter :: DirichletCharacter n -> (Product Natural, [Template]) +templateFromCharacter (Generated t) = traverse go t + where go (OddPrime p k g _) = (Product m, OddTemplate p k g m) + where p' = unPrime p + m = p'^(k-1)*(p'-1) + go (TwoPower k _ _) = (Product (2*m), TwoPTemplate k m) + where m = bit (k-2) + go Two = (Product 1, TwoTemplate) + +mkTemplate :: Natural -> (Product Natural, [Template]) +mkTemplate = go . sort . factorise + where go :: [(Prime Natural, Word)] -> (Product Natural, [Template]) + go ((unPrime -> 2, 1): xs) = (Product 1, [TwoTemplate]) <> traverse odds xs + go ((unPrime -> 2, wordToInt -> k): xs) = (Product (2*m), [TwoPTemplate k m]) <> traverse odds xs + where m = bit (k-2) + go xs = traverse odds xs + odds :: (Prime Natural, Word) -> (Product Natural, Template) + odds (p, k) = (Product m, OddTemplate p k (generator p k) m) + where p' = unPrime p + m = p'^(k-1)*(p'-1) + +-- the validity of the producted dirichletfactor list here requires the template to be valid +unroll :: [Template] -> Natural -> [DirichletFactor] +unroll t m = snd (mapAccumL func m t) + where func :: Natural -> Template -> (Natural, DirichletFactor) + func a (OddTemplate p k g n) = (a1, OddPrime p k g (toRootOfUnity $ (toInteger a2) % (toInteger n))) + where (a1,a2) = quotRem a n + func a (TwoPTemplate k n) = (b1, TwoPower k (toRootOfUnity $ (toInteger a2) % 2) (toRootOfUnity $ (toInteger b2) % (toInteger n))) + where (a1,a2) = quotRem a 2 + (b1,b2) = quotRem a1 n + func a TwoTemplate = (a, Two) + +-- | Test if a given Dirichlet character is prinicpal for its modulus: a principal character mod +-- \(n\) is 1 for \(a\) coprime to \(n\), and 0 otherwise. +isPrincipal :: DirichletCharacter n -> Bool +isPrincipal chi = characterNumber chi == 0 + +-- | Induce a Dirichlet character to a higher modulus. If \(d \mid n\), then \(a \bmod{n}\) can be +-- reduced to \(a \bmod{d}\). Thus, the multiplicative function on \(\mathbb{Z}/d\mathbb{Z}\) +-- induces a multiplicative function on \(\mathbb{Z}/n\mathbb{Z}\). +-- +-- >>> :set -XTypeApplications +-- >>> chi = indexToChar 5 :: DirichletCharacter 45 +-- >>> chi2 = induced @135 chi +-- >>> :t chi2 +-- Maybe (DirichletCharacter 135) +induced :: forall n d. (KnownNat d, KnownNat n) => DirichletCharacter d -> Maybe (DirichletCharacter n) +induced (Generated start) = if n `rem` d == 0 + then Just (Generated (combine (snd $ mkTemplate n) start)) + else Nothing + where n = natVal (Proxy :: Proxy n) + d = natVal (Proxy :: Proxy d) + combine :: [Template] -> [DirichletFactor] -> [DirichletFactor] + combine [] _ = [] + combine ts [] = map newFactor ts + combine (t:xs) (y:ys) = case (t,y) of + (TwoTemplate, Two) -> Two: combine xs ys + (TwoTemplate, _) -> Two: combine xs (y:ys) + (TwoPTemplate k _, Two) -> TwoPower k mempty mempty: combine xs ys + (TwoPTemplate k _, TwoPower _ a b) -> TwoPower k a b: combine xs ys + (TwoPTemplate k _, _) -> TwoPower k mempty mempty: combine xs (y:ys) + (OddTemplate p k _ _, OddPrime q _ g a) | p == q -> OddPrime p k g a: combine xs ys + (OddTemplate p k g _, OddPrime q _ _ _) | p < q -> OddPrime p k g mempty: combine xs (y:ys) + _ -> error "internal error in induced: please report this as a bug" + newFactor :: Template -> DirichletFactor + newFactor TwoTemplate = Two + newFactor (TwoPTemplate k _) = TwoPower k mempty mempty + newFactor (OddTemplate p k g _) = OddPrime p k g mempty + -- rest (p,k) = OddPrime p k (generator p k) mempty + +-- | The gives a real Dirichlet +-- character for odd moduli. +jacobiCharacter :: forall n. KnownNat n => Maybe (RealCharacter n) +jacobiCharacter = if odd n + then Just $ RealChar $ Generated $ map go $ snd $ mkTemplate n + else Nothing + where n = natVal (Proxy :: Proxy n) + go :: Template -> DirichletFactor + go (OddTemplate p k g _) = OddPrime p k g $ toRootOfUnity ((toInteger k) % 2) + -- jacobi symbol of a primitive root mod p over p is always -1 + go _ = error "internal error in jacobiCharacter: please report this as a bug" + -- every factor of n should be odd + +-- | A Dirichlet character is real if it is real-valued. +newtype RealCharacter n = RealChar { -- | Extract the character itself from a `RealCharacter`. + getRealChar :: DirichletCharacter n + } + deriving Eq + +-- | Test if a given `DirichletCharacter` is real, and if so give a `RealCharacter`. +isRealCharacter :: DirichletCharacter n -> Maybe (RealCharacter n) +isRealCharacter t@(Generated xs) = if all real xs then Just (RealChar t) else Nothing + where real :: DirichletFactor -> Bool + real (OddPrime _ _ _ a) = a <> a == mempty + real (TwoPower _ _ b) = b <> b == mempty + real Two = True + +-- TODO: it should be possible to calculate this without eval/evalGeneral +-- and thus avoid using discrete log calculations: consider the order of m +-- inside each of the factor groups? +-- | Evaluate a real Dirichlet character, which can only take values \(-1,0,1\). +toRealFunction :: KnownNat n => RealCharacter n -> Mod n -> Int +toRealFunction (RealChar chi) m = case evalGeneral chi m of + Zero -> 0 + NonZero t | t == mempty -> 1 + NonZero t | t == RootOfUnity (1 % 2) -> -1 + _ -> error "internal error in toRealFunction: please report this as a bug" + -- A real character should not be able to evaluate to + -- anything other than {-1,0,1}, so should not reach this branch + +-- | Test if the internal DirichletCharacter structure is valid. +validChar :: forall n. KnownNat n => DirichletCharacter n -> Bool +validChar (Generated xs) = correctDecomposition && all correctPrimitiveRoot xs && all validValued xs + where correctDecomposition = sort (factorise n) == map getPP xs + getPP (TwoPower k _ _) = (two, fromIntegral k) + getPP (OddPrime p k _ _) = (p, k) + getPP Two = (two,1) + correctPrimitiveRoot (OddPrime p k g _) = g == generator p k + correctPrimitiveRoot _ = True + validValued (TwoPower k a b) = a <> a == mempty && (bit (k-2) :: Integer) `stimes` b == mempty + validValued (OddPrime (unPrime -> p) k _ a) = (p^(k-1)*(p-1)) `stimes` a == mempty + validValued Two = True + n = natVal (Proxy :: Proxy n) + two = nextPrime 2 + +-- | Get the order of the Dirichlet Character. +orderChar :: DirichletCharacter n -> Integer +orderChar (Generated xs) = foldl' lcm 1 $ map orderFactor xs + where orderFactor (TwoPower _ (RootOfUnity a) (RootOfUnity b)) = denominator a `lcm` denominator b + orderFactor (OddPrime _ _ _ (RootOfUnity a)) = denominator a + orderFactor Two = 1 + +-- | Test if a Dirichlet character is . +isPrimitive :: DirichletCharacter n -> Maybe (PrimitiveCharacter n) +isPrimitive t@(Generated xs) = if all primitive xs then Just (PrimitiveCharacter t) else Nothing + where primitive :: DirichletFactor -> Bool + primitive Two = False + -- for odd p, we're testing if phi(p^(k-1)) `stimes` a is 1, since this means the + -- character can come from some the smaller modulus p^(k-1) + primitive (OddPrime _ 1 _ a) = a /= mempty + primitive (OddPrime (unPrime -> p) k _ a) = (p^(k-2)*(p-1)) `stimes` a /= mempty + primitive (TwoPower 2 a _) = a /= mempty + primitive (TwoPower k _ b) = (bit (k-3) :: Integer) `stimes` b /= mempty + +-- | A Dirichlet character is primitive if cannot be 'induced' from any character with +-- strictly smaller modulus. +newtype PrimitiveCharacter n = PrimitiveCharacter { -- | Extract the character itself from a `PrimitiveCharacter`. + getPrimitiveChar :: DirichletCharacter n + } + deriving Eq + +data WithNat (a :: Nat -> *) where + WithNat :: KnownNat m => a m -> WithNat a + +-- | This function also provides access to the new modulus on type level, with a KnownNat instance +makePrimitive :: DirichletCharacter n -> WithNat PrimitiveCharacter +makePrimitive (Generated xs) = + case someNatVal (product mods) of + SomeNat (Proxy :: Proxy m) -> WithNat (PrimitiveCharacter (Generated ys) :: PrimitiveCharacter m) + where (mods,ys) = unzip (mapMaybe prim xs) + prim :: DirichletFactor -> Maybe (Natural, DirichletFactor) + prim Two = Nothing + prim (OddPrime p' k g a) = case find works options of + Nothing -> error "invalid character" + Just (0,_) -> Nothing + Just (i,_) -> Just (p^i, OddPrime p' i g a) + where options = (0,1): [(i,p^(i-1)*(p-1)) | i <- [1..k]] + works (_,phi) = phi `stimes` a == mempty + p = unPrime p' + prim (TwoPower k a b) = case find worksb options of + Nothing -> error "invalid character" + Just (2,_) | a == mempty -> Nothing + Just (i,_) -> Just (bit i :: Natural, TwoPower i a b) + where options = [(i, bit (i-2) :: Natural) | i <- [2..k]] + worksb (_,phi) = phi `stimes` b == mempty + +#if !MIN_VERSION_base(4,12,0) +newtype Ap f a = Ap { getAp :: f a } + deriving (Eq, Functor, Applicative, Monad) + +instance (Applicative f, Semigroup a) => Semigroup (Ap f a) where + (<>) = liftA2 (<>) + +instance (Applicative f, Semigroup a, Monoid a) => Monoid (Ap f a) where + mempty = pure mempty + mappend = (<>) +#endif + +-- | Similar to Maybe, but with different Semigroup and Monoid instances. +type OrZero a = Ap Maybe a +pattern Zero :: OrZero a +pattern Zero = Ap Nothing + +pattern NonZero :: a -> OrZero a +pattern NonZero x = Ap (Just x) +{-# COMPLETE Zero, NonZero #-} + +-- | Interpret an `OrZero` as a number, taking the `Zero` case to be 0. +orZeroToNum :: Num a => (b -> a) -> OrZero b -> a +orZeroToNum _ Zero = 0 +orZeroToNum f (NonZero x) = f x + +-- | In general, evaluating a DirichletCharacter at a point involves solving the discrete logarithm +-- problem, which can be hard: the implementations here are around O(sqrt n). +-- However, evaluating a dirichlet character at every point amounts to solving the discrete +-- logarithm problem at every point also, which can be done together in O(n) time, better than +-- using a complex algorithm at each point separately. Thus, if a large number of evaluations +-- of a dirichlet character are required, `evalAll` will be better than `evalGeneral`, since +-- computations can be shared. +evalAll :: forall n. KnownNat n => DirichletCharacter n -> Vector (OrZero RootOfUnity) +evalAll (Generated xs) = V.generate (fromIntegral n) func + where n = natVal (Proxy :: Proxy n) + vectors = map mkVector xs + func :: Int -> OrZero RootOfUnity + func m = foldMap go vectors + where go :: (Int, Vector (OrZero RootOfUnity)) -> OrZero RootOfUnity + go (modulus,v) = v ! (m `mod` modulus) + mkVector :: DirichletFactor -> (Int, Vector (OrZero RootOfUnity)) + mkVector Two = (2, V.fromList [Zero, mempty]) + mkVector (OddPrime p k (fromIntegral -> g) a) = (modulus, w) + where + p' = unPrime p + modulus = fromIntegral (p'^k) :: Int + w = V.create $ do + v <- MV.replicate modulus Zero + -- TODO: we're in the ST monad here anyway, could be better to use STRefs to manage + -- this loop, the current implementation probably doesn't fuse well + let powers = iterateMaybe go (1,mempty) + go (m,x) = if m' > 1 + then Just (m', x<>a) + else Nothing + where m' = m*g `mod` modulus + for_ powers $ \(m,x) -> MV.unsafeWrite v m (NonZero x) + -- don't bother with bounds check since m was reduced mod p^k + return v + -- for powers of two we use lambda directly instead, since the generators of the cyclic + -- groups aren't obvious; it's possible to get them though: + -- 5^(lambda(5)^{-1} mod 2^(p-2)) mod 2^p + mkVector (TwoPower k a b) = (modulus, w) + where + modulus = bit k + w = V.generate modulus f + f m + | even m = Zero + | otherwise = NonZero ((if testBit m 1 then a else mempty) <> lambda (toInteger m'') k `stimes` b) + where m'' = thingy k m + +-- somewhere between unfoldr and iterate +iterateMaybe :: (a -> Maybe a) -> a -> [a] +iterateMaybe f x = unfoldr (fmap (\t -> (t, f t))) (Just x) + +-- | Attempt to construct a character from its table of values. +-- An inverse to `evalAll`, defined only on its image. +fromTable :: forall n. KnownNat n => Vector (OrZero RootOfUnity) -> Maybe (DirichletCharacter n) +fromTable v = if length v == fromIntegral n + then Generated <$> traverse makeFactor tmpl >>= check + else Nothing + where n = natVal (Proxy :: Proxy n) + n' = fromIntegral n :: Integer + tmpl = snd (mkTemplate n) + check :: DirichletCharacter n -> Maybe (DirichletCharacter n) + check chi = if evalAll chi == v then Just chi else Nothing + makeFactor :: Template -> Maybe DirichletFactor + makeFactor TwoTemplate = Just Two + makeFactor (TwoPTemplate k _) = TwoPower k <$> getValue (-1,bit k) <*> getValue (exp4 k, bit k) + makeFactor (OddTemplate p k g _) = OddPrime p k g <$> getValue (toInteger g, toInteger (unPrime p)^k) + getValue :: (Integer,Integer) -> Maybe RootOfUnity + getValue (g,m) = getAp (v ! fromInteger (fromJust (chineseCoprime (g,m) (1,n' `quot` m)) `mod` n')) + +exp4terms :: [Rational] +exp4terms = [4^k % product [1..k] | k <- [0..]] + +-- For reasons that aren't clear to me, `exp4` gives the inverse of 1 under lambda, so it gives the generator +-- This is the same as https://oeis.org/A320814 +-- In particular, lambda (exp4 n) n == 1 (for n >= 3) +-- I've verified this for 3 <= n <= 2000, so the reasoning in fromTable should be accurate for moduli below 2^2000 +exp4 :: Int -> Integer +exp4 n = (`mod` bit n) $ sum $ map (`mod` bit n) $ map (\q -> numerator q * fromMaybe (error "error in exp4") (recipMod (denominator q) (bit n))) $ take n $ exp4terms diff --git a/Math/NumberTheory/Moduli/Internal.hs b/Math/NumberTheory/Moduli/Internal.hs new file mode 100644 index 000000000..c2148b987 --- /dev/null +++ b/Math/NumberTheory/Moduli/Internal.hs @@ -0,0 +1,126 @@ +-- | +-- Module: Math.NumberTheory.Moduli.Internal +-- Copyright: (c) 2020 Bhavik Mehta +-- Licence: MIT +-- Maintainer: Bhavik Mehta +-- +-- Multiplicative groups of integers modulo m. +-- + +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module Math.NumberTheory.Moduli.Internal + ( isPrimitiveRoot' + , discreteLogarithmPP + ) where + +import qualified Data.Map as M +import Data.Maybe +import Data.Mod +import Data.Proxy +import GHC.Integer.GMP.Internals +import GHC.TypeNats.Compat +import Numeric.Natural + +import Math.NumberTheory.ArithmeticFunctions +import Math.NumberTheory.Moduli.Chinese +import Math.NumberTheory.Moduli.Equations +import Math.NumberTheory.Moduli.Singleton +import Math.NumberTheory.Primes +import Math.NumberTheory.Powers.Modular +import Math.NumberTheory.Roots + +-- https://en.wikipedia.org/wiki/Primitive_root_modulo_n#Finding_primitive_roots +isPrimitiveRoot' + :: (Integral a, UniqueFactorisation a) + => CyclicGroup a m + -> a + -> Bool +isPrimitiveRoot' cg r = + case cg of + CG2 -> r == 1 + CG4 -> r == 3 + CGOddPrimePower p k -> oddPrimePowerTest (unPrime p) k r + CGDoubleOddPrimePower p k -> doubleOddPrimePowerTest (unPrime p) k r + where + oddPrimeTest p g = let phi = totient p + pows = map (\pk -> phi `quot` unPrime (fst pk)) (factorise phi) + exps = map (\x -> powMod g x p) pows + in g /= 0 && gcd g p == 1 && notElem 1 exps + oddPrimePowerTest p 1 g = oddPrimeTest p (g `mod` p) + oddPrimePowerTest p _ g = oddPrimeTest p (g `mod` p) && powMod g (p-1) (p*p) /= 1 + doubleOddPrimePowerTest p k g = odd g && oddPrimePowerTest p k g + +-- Implementation of Bach reduction (https://www2.eecs.berkeley.edu/Pubs/TechRpts/1984/CSD-84-186.pdf) +{-# INLINE discreteLogarithmPP #-} +discreteLogarithmPP :: Integer -> Word -> Integer -> Integer -> Natural +discreteLogarithmPP p 1 a b = discreteLogarithmPrime p a b +discreteLogarithmPP p k a b = fromInteger $ if result < 0 then result + pkMinusPk1 else result + where + baseSol = toInteger $ discreteLogarithmPrime p (a `rem` p) (b `rem` p) + thetaA = theta p pkMinusOne a + thetaB = theta p pkMinusOne b + pkMinusOne = p^(k-1) + pkMinusPk1 = pkMinusOne * (p - 1) + c = (recipModInteger thetaA pkMinusOne * thetaB) `rem` pkMinusOne + result = fromJust $ chineseCoprime (baseSol, p-1) (c, pkMinusOne) + +-- compute the homomorphism theta given in https://math.stackexchange.com/a/1864495/418148 +{-# INLINE theta #-} +theta :: Integer -> Integer -> Integer -> Integer +theta p pkMinusOne a = (numerator `quot` pk) `rem` pkMinusOne + where + pk = pkMinusOne * p + p2kMinusOne = pkMinusOne * pk + numerator = (powModInteger a (pk - pkMinusOne) p2kMinusOne - 1) `rem` p2kMinusOne + +-- TODO: Use Pollig-Hellman to reduce the problem further into groups of prime order. +-- While Bach reduction simplifies the problem into groups of the form (Z/pZ)*, these +-- have non-prime order, and the Pollig-Hellman algorithm can reduce the problem into +-- smaller groups of prime order. +-- In addition, the gcd check before solveLinear is applied in Pollard below will be +-- made redundant, since n would be prime. +discreteLogarithmPrime :: Integer -> Integer -> Integer -> Natural +discreteLogarithmPrime p a b + | p < 100000000 = fromIntegral $ discreteLogarithmPrimeBSGS (fromInteger p) (fromInteger a) (fromInteger b) + | otherwise = discreteLogarithmPrimePollard p a b + +discreteLogarithmPrimeBSGS :: Int -> Int -> Int -> Int +discreteLogarithmPrimeBSGS p a b = head [i*m + j | (v,i) <- zip giants [0..m-1], j <- maybeToList (M.lookup v table)] + where + m = integerSquareRoot (p - 2) + 1 -- simple way of ceiling (sqrt (p-1)) + babies = iterate (.* a) 1 + table = M.fromList (zip babies [0..m-1]) + aInv = recipModInteger (toInteger a) (toInteger p) + bigGiant = fromInteger $ powModInteger aInv (toInteger m) (toInteger p) + giants = iterate (.* bigGiant) b + x .* y = x * y `rem` p + +-- TODO: Use more advanced walks, in order to reduce divisions, cf +-- https://maths-people.anu.edu.au/~brent/pd/rpb231.pdf +-- This will slightly improve the expected time to collision, and can reduce the +-- number of divisions performed. +discreteLogarithmPrimePollard :: Integer -> Integer -> Integer -> Natural +discreteLogarithmPrimePollard p a b = + case concatMap runPollard [(x,y) | x <- [0..n], y <- [0..n]] of + (t:_) -> fromInteger t + [] -> error ("discreteLogarithm: pollard's rho failed, please report this as a bug. inputs " ++ show [p,a,b]) + where + n = p-1 -- order of the cyclic group + halfN = n `quot` 2 + mul2 m = if m < halfN then m * 2 else m * 2 - n + sqrtN = integerSquareRoot n + step (xi,!ai,!bi) = case xi `rem` 3 of + 0 -> (xi*xi `rem` p, mul2 ai, mul2 bi) + 1 -> ( a*xi `rem` p, ai+1, bi) + _ -> ( b*xi `rem` p, ai, bi+1) + initialise (x,y) = (powModInteger a x n * powModInteger b y n `rem` n, x, y) + begin t = go (step t) (step (step t)) + check t = powModInteger a t p == b + go tort@(xi,ai,bi) hare@(x2i,a2i,b2i) + | xi == x2i, gcd (bi - b2i) n < sqrtN = case someNatVal (fromInteger n) of + SomeNat (Proxy :: Proxy n) -> map (toInteger . unMod) $ solveLinear (fromInteger (bi - b2i) :: Mod n) (fromInteger (ai - a2i)) + | xi == x2i = [] + | otherwise = go (step tort) (step (step hare)) + runPollard = filter check . begin . initialise diff --git a/Math/NumberTheory/Moduli/JacobiSymbol.hs b/Math/NumberTheory/Moduli/JacobiSymbol.hs index 6d7256cb9..9be421d6a 100644 --- a/Math/NumberTheory/Moduli/JacobiSymbol.hs +++ b/Math/NumberTheory/Moduli/JacobiSymbol.hs @@ -17,6 +17,7 @@ module Math.NumberTheory.Moduli.JacobiSymbol ( JacobiSymbol(..) , jacobi + , symbolToNum ) where import Data.Bits @@ -48,6 +49,18 @@ negJS = \case Zero -> Zero One -> MinusOne +{-# SPECIALISE symbolToNum :: JacobiSymbol -> Integer, + JacobiSymbol -> Int, + JacobiSymbol -> Word, + JacobiSymbol -> Natural + #-} +-- | Convenience function to convert out of a Jacobi symbol +symbolToNum :: Num a => JacobiSymbol -> a +symbolToNum = \case + Zero -> 0 + One -> 1 + MinusOne -> -1 + -- | of two arguments. -- The lower argument (\"denominator\") must be odd and positive, -- this condition is checked. diff --git a/Math/NumberTheory/Moduli/Multiplicative.hs b/Math/NumberTheory/Moduli/Multiplicative.hs index e028329ed..2160e1664 100644 --- a/Math/NumberTheory/Moduli/Multiplicative.hs +++ b/Math/NumberTheory/Moduli/Multiplicative.hs @@ -7,9 +7,9 @@ -- Multiplicative groups of integers modulo m. -- -{-# LANGUAGE BangPatterns #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE ViewPatterns #-} +{-# LANGUAGE PatternSynonyms #-} module Math.NumberTheory.Moduli.Multiplicative ( -- * Multiplicative group @@ -26,22 +26,14 @@ module Math.NumberTheory.Moduli.Multiplicative import Control.Monad import Data.Constraint -import qualified Data.Map as M -import Data.Maybe import Data.Mod -import Data.Proxy import Data.Semigroup -import GHC.Integer.GMP.Internals import GHC.TypeNats.Compat import Numeric.Natural -import Math.NumberTheory.ArithmeticFunctions -import Math.NumberTheory.Moduli.Chinese -import Math.NumberTheory.Moduli.Equations +import Math.NumberTheory.Moduli.Internal import Math.NumberTheory.Moduli.Singleton import Math.NumberTheory.Primes -import Math.NumberTheory.Powers.Modular -import Math.NumberTheory.Roots -- | This type represents elements of the multiplicative group mod m, i.e. -- those elements which are coprime to m. Use @toMultElement@ to construct. @@ -84,27 +76,6 @@ newtype PrimitiveRoot m = PrimitiveRoot } deriving (Eq, Show) --- https://en.wikipedia.org/wiki/Primitive_root_modulo_n#Finding_primitive_roots -isPrimitiveRoot' - :: (Integral a, UniqueFactorisation a) - => CyclicGroup a m - -> a - -> Bool -isPrimitiveRoot' cg r = - case cg of - CG2 -> r == 1 - CG4 -> r == 3 - CGOddPrimePower p k -> oddPrimePowerTest (unPrime p) k r - CGDoubleOddPrimePower p k -> doubleOddPrimePowerTest (unPrime p) k r - where - oddPrimeTest p g = let phi = totient p - pows = map (\pk -> phi `quot` unPrime (fst pk)) (factorise phi) - exps = map (\x -> powMod g x p) pows - in g /= 0 && gcd g p == 1 && all (/= 1) exps - oddPrimePowerTest p 1 g = oddPrimeTest p (g `mod` p) - oddPrimePowerTest p _ g = oddPrimeTest p (g `mod` p) && powMod g (p-1) (p*p) /= 1 - doubleOddPrimePowerTest p k g = odd g && oddPrimePowerTest p k g - -- | Check whether a given modular residue is -- a . -- @@ -148,76 +119,3 @@ discreteLogarithm cg (multElement . unPrimitiveRoot -> a) (multElement -> b) = c CGDoubleOddPrimePower (unPrime -> p) k -> discreteLogarithmPP p k (toInteger (unMod a) `rem` p^k) (toInteger (unMod b) `rem` p^k) -- we have the isomorphism t -> t `rem` p^k from (Z/2p^kZ)* -> (Z/p^kZ)* - --- Implementation of Bach reduction (https://www2.eecs.berkeley.edu/Pubs/TechRpts/1984/CSD-84-186.pdf) -{-# INLINE discreteLogarithmPP #-} -discreteLogarithmPP :: Integer -> Word -> Integer -> Integer -> Natural -discreteLogarithmPP p 1 a b = discreteLogarithmPrime p a b -discreteLogarithmPP p k a b = fromInteger $ if result < 0 then result + pkMinusPk1 else result - where - baseSol = toInteger $ discreteLogarithmPrime p (a `rem` p) (b `rem` p) - thetaA = theta p pkMinusOne a - thetaB = theta p pkMinusOne b - pkMinusOne = p^(k-1) - pkMinusPk1 = pkMinusOne * (p - 1) - c = (recipModInteger thetaA pkMinusOne * thetaB) `rem` pkMinusOne - result = fromJust $ chineseCoprime (baseSol, p-1) (c, pkMinusOne) - --- compute the homomorphism theta given in https://math.stackexchange.com/a/1864495/418148 -{-# INLINE theta #-} -theta :: Integer -> Integer -> Integer -> Integer -theta p pkMinusOne a = (numerator `quot` pk) `rem` pkMinusOne - where - pk = pkMinusOne * p - p2kMinusOne = pkMinusOne * pk - numerator = (powModInteger a (pk - pkMinusOne) p2kMinusOne - 1) `rem` p2kMinusOne - --- TODO: Use Pollig-Hellman to reduce the problem further into groups of prime order. --- While Bach reduction simplifies the problem into groups of the form (Z/pZ)*, these --- have non-prime order, and the Pollig-Hellman algorithm can reduce the problem into --- smaller groups of prime order. --- In addition, the gcd check before solveLinear is applied in Pollard below will be --- made redundant, since n would be prime. -discreteLogarithmPrime :: Integer -> Integer -> Integer -> Natural -discreteLogarithmPrime p a b - | p < 100000000 = fromIntegral $ discreteLogarithmPrimeBSGS (fromInteger p) (fromInteger a) (fromInteger b) - | otherwise = discreteLogarithmPrimePollard p a b - -discreteLogarithmPrimeBSGS :: Int -> Int -> Int -> Int -discreteLogarithmPrimeBSGS p a b = head [i*m + j | (v,i) <- zip giants [0..m-1], j <- maybeToList (M.lookup v table)] - where - m = integerSquareRoot (p - 2) + 1 -- simple way of ceiling (sqrt (p-1)) - babies = iterate (.* a) 1 - table = M.fromList (zip babies [0..m-1]) - aInv = recipModInteger (toInteger a) (toInteger p) - bigGiant = fromInteger $ powModInteger aInv (toInteger m) (toInteger p) - giants = iterate (.* bigGiant) b - x .* y = x * y `rem` p - --- TODO: Use more advanced walks, in order to reduce divisions, cf --- https://maths-people.anu.edu.au/~brent/pd/rpb231.pdf --- This will slightly improve the expected time to collision, and can reduce the --- number of divisions performed. -discreteLogarithmPrimePollard :: Integer -> Integer -> Integer -> Natural -discreteLogarithmPrimePollard p a b = - case concatMap runPollard [(x,y) | x <- [0..n], y <- [0..n]] of - (t:_) -> fromInteger t - [] -> error ("discreteLogarithm: pollard's rho failed, please report this as a bug. inputs " ++ show [p,a,b]) - where - n = p-1 -- order of the cyclic group - halfN = n `quot` 2 - mul2 m = if m < halfN then m * 2 else m * 2 - n - sqrtN = integerSquareRoot n - step (xi,!ai,!bi) = case xi `rem` 3 of - 0 -> (xi*xi `rem` p, mul2 ai, mul2 bi) - 1 -> ( a*xi `rem` p, ai+1, bi) - _ -> ( b*xi `rem` p, ai, bi+1) - initialise (x,y) = (powModInteger a x n * powModInteger b y n `rem` n, x, y) - begin t = go (step t) (step (step t)) - check t = powModInteger a t p == b - go tort@(xi,ai,bi) hare@(x2i,a2i,b2i) - | xi == x2i, gcd (bi - b2i) n < sqrtN = case someNatVal (fromInteger n) of - SomeNat (Proxy :: Proxy n) -> map (toInteger . unMod) $ solveLinear (fromInteger (bi - b2i) :: Mod n) (fromInteger (ai - a2i)) - | xi == x2i = [] - | otherwise = go (step tort) (step (step hare)) - runPollard = filter check . begin . initialise diff --git a/Math/NumberTheory/Moduli/Sqrt.hs b/Math/NumberTheory/Moduli/Sqrt.hs index 45347cba6..074189aba 100644 --- a/Math/NumberTheory/Moduli/Sqrt.hs +++ b/Math/NumberTheory/Moduli/Sqrt.hs @@ -21,6 +21,7 @@ module Math.NumberTheory.Moduli.Sqrt -- * Jacobi symbol , JacobiSymbol(..) , jacobi + , symbolToNum ) where import Control.Monad (liftM2) diff --git a/arithmoi.cabal b/arithmoi.cabal index 4a55f54d9..32a8a772f 100644 --- a/arithmoi.cabal +++ b/arithmoi.cabal @@ -59,6 +59,7 @@ library Math.NumberTheory.ArithmeticFunctions.Moebius Math.NumberTheory.ArithmeticFunctions.SieveBlock Math.NumberTheory.Curves.Montgomery + Math.NumberTheory.DirichletCharacters Math.NumberTheory.Euclidean Math.NumberTheory.Euclidean.Coprimes Math.NumberTheory.Moduli @@ -94,6 +95,7 @@ library other-modules: Math.NumberTheory.ArithmeticFunctions.Class Math.NumberTheory.ArithmeticFunctions.Standard + Math.NumberTheory.Moduli.Internal Math.NumberTheory.Moduli.JacobiSymbol Math.NumberTheory.Moduli.SomeMod Math.NumberTheory.Primes.Counting.Approximate @@ -151,6 +153,7 @@ test-suite spec Math.NumberTheory.ArithmeticFunctions.MertensTests Math.NumberTheory.ArithmeticFunctions.SieveBlockTests Math.NumberTheory.CurvesTests + Math.NumberTheory.DirichletCharactersTests Math.NumberTheory.EisensteinIntegersTests Math.NumberTheory.GaussianIntegersTests Math.NumberTheory.EuclideanTests diff --git a/test-suite/Math/NumberTheory/DirichletCharactersTests.hs b/test-suite/Math/NumberTheory/DirichletCharactersTests.hs new file mode 100644 index 000000000..5b23c3c0a --- /dev/null +++ b/test-suite/Math/NumberTheory/DirichletCharactersTests.hs @@ -0,0 +1,253 @@ +-- | +-- Module: Math.NumberTheory.Moduli.DiscreteLogarithm +-- Copyright: (c) 2018 Bhavik Mehta +-- License: MIT +-- Maintainer: Andrew Lelechenko +-- +-- Tests for Math.NumberTheory.DirichletCharacters +-- + +{-# LANGUAGE GADTs #-} +{-# LANGUAGE PatternSynonyms #-} +{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE ViewPatterns #-} + +module Math.NumberTheory.DirichletCharactersTests where + +import Test.Tasty + +import Data.Complex +import Data.List (genericLength) +import Data.Maybe (isJust, mapMaybe) +import Data.Proxy +import Data.Ratio +import Data.Semigroup +import qualified Data.Vector as V +import Numeric.Natural + +import GHC.TypeNats.Compat (SomeNat(..), someNatVal, KnownNat, natVal, sameNat) +import Data.Type.Equality + +import Math.NumberTheory.ArithmeticFunctions (totient, divisorsList) +import Math.NumberTheory.DirichletCharacters +import qualified Math.NumberTheory.Moduli.Sqrt as J +import Math.NumberTheory.Moduli.Class (SomeMod(..), modulo) +import Math.NumberTheory.TestUtils (testSmallAndQuick, Positive(..)) + +rootOfUnityTest :: Integer -> Positive Integer -> Bool +rootOfUnityTest n (Positive d) = toComplex ((d `div` gcd n d) `stimes` toRootOfUnity (n % d)) == (1 :: Complex Double) + +-- | This tests property 6 from https://en.wikipedia.org/wiki/Dirichlet_character#Axiomatic_definition +dirCharOrder :: forall n. KnownNat n => DirichletCharacter n -> Bool +dirCharOrder chi = isPrincipal (totient n `stimes` chi) + where n = natVal @n Proxy + +-- | Tests wikipedia's property 3 (note 1,2,5 are essentially enforced by the type system). +testMultiplicative :: KnownNat n => DirichletCharacter n -> Natural -> Natural -> Bool +testMultiplicative chi (fromIntegral -> a) (fromIntegral -> b) = chiAB == chiAchiB + where chi' = evalGeneral chi + chiAB = chi' (a*b) + chiAchiB = (<>) <$> chi' a <*> chi' b + +-- | Test property 4 from wikipedia +testAtOne :: KnownNat n => DirichletCharacter n -> Bool +testAtOne chi = eval chi mempty == mempty + +dirCharProperty :: (forall n. KnownNat n => DirichletCharacter n -> a) -> Positive Natural -> Natural -> a +dirCharProperty test (Positive n) i = + case someNatVal n of + SomeNat (Proxy :: Proxy n) -> test chi + where chi = indexToChar @n (i `mod` totient n) + +realCharProperty :: (forall n. KnownNat n => RealCharacter n -> a) -> Positive Natural -> Int -> a +realCharProperty test (Positive n) i = + case someNatVal n of + SomeNat (Proxy :: Proxy n) -> test chi + where chi = chars !! (i `mod` length chars) + chars = mapMaybe isRealCharacter [principalChar @n .. maxBound] + +-- | There should be totient(n) characters +countCharacters :: Positive Natural -> Bool +countCharacters (Positive n) = + case someNatVal n of + SomeNat (Proxy :: Proxy n) -> + genericLength (allChars @n) == totient n + +-- | The principal character should be 1 if gcd k n is 1 and 0 otherwise +principalCase :: Positive Natural -> Positive Integer -> Bool +principalCase (Positive n) (Positive k) = + case k `modulo` n of + SomeMod a -> evalGeneral chi a == if gcd k (fromIntegral n) > 1 + then Zero + else mempty + where chi = principalChar + InfMod{} -> False + +-- | Test the orthogonality relations https://en.wikipedia.org/wiki/Dirichlet_character#Character_orthogonality +orthogonality1 :: forall n. KnownNat n => DirichletCharacter n -> Bool +orthogonality1 chi = magnitude (total - correct) < (1e-13 :: Double) + where n = natVal @n Proxy + total = sum [orZeroToNum toComplex (evalGeneral chi a) | a <- [0 .. maxBound]] + correct = if isPrincipal chi + then fromIntegral $ totient n + else 0 + +orthogonality2 :: Positive Natural -> Integer -> Bool +orthogonality2 (Positive n) a = + case a `modulo` n of + SomeMod a' -> magnitude (total - correct) < (1e-13 :: Double) + where total = sum [orZeroToNum toComplex (evalGeneral chi a') | chi <- allChars] + correct = if a' == 1 + then fromIntegral $ totient n + else 0 + InfMod {} -> False + +-- | Manually confirm isRealCharacter is correct (in both directions) +realityCheck :: KnownNat n => DirichletCharacter n -> Bool +realityCheck chi = isJust (isRealCharacter chi) == isReal' + where isReal' = and [real (evalGeneral chi t) | t <- [minBound..maxBound]] + real Zero = True + real (NonZero t) = t <> t == mempty + +-- | Check real character evaluation matches normal evaluation +realEvalCheck :: KnownNat n => RealCharacter n -> Int -> Bool +realEvalCheck chi i' = fromIntegral (toRealFunction chi i) == orZeroToNum toComplex (evalGeneral (getRealChar chi) i) + where i = fromIntegral i' + +-- | The jacobi character agrees with the jacobi symbol +jacobiCheck :: Positive Natural -> Bool +jacobiCheck (Positive n) = + case someNatVal (2*n+1) of + SomeNat (Proxy :: Proxy n) -> + case jacobiCharacter @n of + Just chi -> and [toRealFunction chi (fromIntegral j) == J.symbolToNum (J.jacobi j (2*n+1)) | j <- [0..2*n]] + _ -> False + +-- | Bulk evaluation agrees with pointwise evaluation +evalAllCheck :: forall n. KnownNat n => DirichletCharacter n -> Bool +evalAllCheck chi = V.generate (fromIntegral $ natVal @n Proxy) (evalGeneral chi . fromIntegral) == evalAll chi + +-- | Induced characters agree with the original character. +-- (Except for when d=1, where chi(0) = 1, which is true for no other d) +inducedCheck :: forall d. KnownNat d => DirichletCharacter d -> Positive Natural -> Bool +inducedCheck chi (Positive k) = + case someNatVal (d*k) of + SomeNat (Proxy :: Proxy n) -> + case induced @n chi of + Just chi2 -> and (V.izipWith matchedValue (V.concat (replicate (fromIntegral k) (evalAll chi))) (evalAll chi2)) + Nothing -> False + where d = natVal @d Proxy + matchedValue i x1 x2 = if gcd (fromIntegral i) (d*k) > 1 + then x2 == Zero + else x2 == x1 + +-- | Primitive checker is correct (in both directions) +primitiveCheck :: forall n. KnownNat n => DirichletCharacter n -> Bool +primitiveCheck chi = isJust (isPrimitive chi) == isPrimitive' + where isPrimitive' = all testModulus possibleModuli + n = fromIntegral (natVal @n Proxy) :: Int + possibleModuli = init (divisorsList n) + table = evalAll chi + testModulus d = not $ null [a | a <- [1..n-1], gcd a n == 1, a `mod` d == 1 `mod` d, table V.! a /= mempty] + +-- | Ensure that makePrimitive gives primitive characters +makePrimitiveCheck :: DirichletCharacter n -> Bool +makePrimitiveCheck chi = case makePrimitive chi of + WithNat chi' -> isJust (isPrimitive (getPrimitiveChar chi')) + +-- | sameNat also ensures the two new moduli are the same +makePrimitiveIdem :: DirichletCharacter n -> Bool +makePrimitiveIdem chi = case makePrimitive chi of + WithNat (chi' :: PrimitiveCharacter n') -> + case makePrimitive (getPrimitiveChar chi') of + WithNat (chi'' :: PrimitiveCharacter n'') -> + case sameNat (Proxy :: Proxy n') (Proxy :: Proxy n'') of + Just Refl -> chi' == chi'' + Nothing -> False + +orderCheck :: DirichletCharacter n -> Bool +orderCheck chi = isPrincipal (n `stimes` chi) && and [not (isPrincipal (i `stimes` chi)) | i <- [1..n-1]] + where n = orderChar chi + +fromTableCheck :: forall n. KnownNat n => DirichletCharacter n -> Bool +fromTableCheck chi = isJust (fromTable @n (evalAll chi)) + +-- A bunch of functions making sure that every function which can produce a character (in +-- particular by fiddling internal representation) produces a valid character +indexToCharValid :: KnownNat n => DirichletCharacter n -> Bool +indexToCharValid = validChar + +principalCharValid :: Positive Natural -> Bool +principalCharValid (Positive n) = + case someNatVal n of + SomeNat (Proxy :: Proxy n) -> validChar (principalChar @n) + +mulCharsValid :: KnownNat n => DirichletCharacter n -> DirichletCharacter n -> Bool +mulCharsValid chi1 chi2 = validChar (chi1 <> chi2) + +mulCharsValid' :: Positive Natural -> Natural -> Natural -> Bool +mulCharsValid' (Positive n) i j = + case someNatVal n of + SomeNat (Proxy :: Proxy n) -> + mulCharsValid (indexToChar @n (i `mod` totient n)) (indexToChar @n (j `mod` totient n)) + +stimesCharValid :: KnownNat n => DirichletCharacter n -> Int -> Bool +stimesCharValid chi n = validChar (n `stimes` chi) + +succValid :: KnownNat n => DirichletCharacter n -> Bool +succValid = validChar . succ + +inducedValid :: forall d. KnownNat d => DirichletCharacter d -> Positive Natural -> Bool +inducedValid chi (Positive k) = + case someNatVal (d*k) of + SomeNat (Proxy :: Proxy n) -> + case induced @n chi of + Just chi2 -> validChar chi2 + Nothing -> False + where d = natVal @d Proxy + +jacobiValid :: Positive Natural -> Bool +jacobiValid (Positive n) = + case someNatVal (2*n+1) of + SomeNat (Proxy :: Proxy n) -> + case jacobiCharacter @n of + Just chi -> validChar (getRealChar chi) + _ -> False + +makePrimitiveValid :: DirichletCharacter n -> Bool +makePrimitiveValid chi = case makePrimitive chi of + WithNat chi' -> validChar (getPrimitiveChar chi') + +testSuite :: TestTree +testSuite = testGroup "DirichletCharacters" + [ testSmallAndQuick "RootOfUnity contains roots of unity" rootOfUnityTest + , testSmallAndQuick "Dirichlet characters divide the right order" (dirCharProperty dirCharOrder) + , testSmallAndQuick "Dirichlet characters are multiplicative" (dirCharProperty testMultiplicative) + , testSmallAndQuick "Dirichlet characters are 1 at 1" (dirCharProperty testAtOne) + , testSmallAndQuick "Right number of Dirichlet characters" countCharacters + , testSmallAndQuick "Principal character behaves as expected" principalCase + , testSmallAndQuick "Orthogonality relation 1" (dirCharProperty orthogonality1) + , testSmallAndQuick "Orthogonality relation 2" orthogonality2 + , testSmallAndQuick "Real character checking is correct" (dirCharProperty realityCheck) + , testSmallAndQuick "Real character evaluation is accurate" (realCharProperty realEvalCheck) + , testSmallAndQuick "Jacobi character matches symbol" jacobiCheck + , testSmallAndQuick "Bulk evaluation matches pointwise" (dirCharProperty evalAllCheck) + , testSmallAndQuick "Induced character is correct" (dirCharProperty inducedCheck) + , testSmallAndQuick "Primitive character checking is correct" (dirCharProperty primitiveCheck) + , testSmallAndQuick "makePrimitive produces primitive character" (dirCharProperty makePrimitiveCheck) + , testSmallAndQuick "makePrimitive is idempotent" (dirCharProperty makePrimitiveIdem) + , testSmallAndQuick "Calculates correct order" (dirCharProperty orderCheck) + , testSmallAndQuick "Can construct from table" (dirCharProperty fromTableCheck) + , testGroup "Creates valid characters" + [ testSmallAndQuick "indexToChar" (dirCharProperty indexToCharValid) + , testSmallAndQuick "principalChar" principalCharValid + , testSmallAndQuick "mulChars" mulCharsValid' + , testSmallAndQuick "stimesChar" (dirCharProperty stimesCharValid) + , testSmallAndQuick "succ" (dirCharProperty succValid) + , testSmallAndQuick "induced" (dirCharProperty inducedValid) + , testSmallAndQuick "jacobi" jacobiValid + , testSmallAndQuick "makePrimitive" (dirCharProperty makePrimitiveValid) + ] + ] diff --git a/test-suite/Test.hs b/test-suite/Test.hs index 7a6e31b20..a1f033df3 100644 --- a/test-suite/Test.hs +++ b/test-suite/Test.hs @@ -48,6 +48,8 @@ import qualified Math.NumberTheory.SmoothNumbersTests as SmoothNumbers import qualified Math.NumberTheory.Zeta.RiemannTests as Riemann import qualified Math.NumberTheory.Zeta.DirichletTests as Dirichlet +import qualified Math.NumberTheory.DirichletCharactersTests as DirichletChar + main :: IO () main = defaultMainWithRerun tests @@ -102,4 +104,5 @@ tests = testGroup "All" [ Riemann.testSuite , Dirichlet.testSuite ] + , DirichletChar.testSuite ]