diff --git a/Math/NumberTheory/Moduli/Chinese.hs b/Math/NumberTheory/Moduli/Chinese.hs index 685e818ef..5048c9613 100644 --- a/Math/NumberTheory/Moduli/Chinese.hs +++ b/Math/NumberTheory/Moduli/Chinese.hs @@ -6,19 +6,119 @@ -- -- Chinese remainder theorem -- - {-# LANGUAGE BangPatterns #-} module Math.NumberTheory.Moduli.Chinese - ( chineseRemainder + ( chineseRemainders + , chineseRemainders2 + , chineseRemainder , chineseRemainder2 - ) where + ) +where +import Data.Ratio (numerator, denominator) import Control.Monad (foldM) import Math.NumberTheory.Euclidean (extendedGCD) +import Math.NumberTheory.Moduli.Class +import Math.NumberTheory.Primes.Factorisation (factorise) import Math.NumberTheory.Utils (recipMod) +-- [Unsure where to put these. Does arithmoi keep QuickCheck tests in separate documents? +-- +-- The two tests are slightly different. The test for the binary function compares to the +-- result of an exhaustive search. The test for the list function merely tests that the results +-- are consistent with the requirements as exhaustive search over the LCM of the moduli of long +-- lists tends to take an unreasonable amount of time. +-- +-- import Test.QuickCheck +-- import Data.List (foldl') + +-- prop_chineseRemainders :: [(Integer, Positive Integer)] -> Bool +-- prop_chineseRemainders rms +-- | Just (SomeMod d) <- c = getMod d == l && all (\ (r, Positive m) -> (getVal d-r) `mod` m==0) rms +-- | otherwise = or [(r1-r2) `mod` g /= 0|((r1, Positive m1):rms2)<-tails rms, (r2, Positive m2)<-rms2, let g = gcd m1 m2] +-- where +-- l = foldl' lcm 1 $ map (getPositive.snd) rms +-- c = chineseRemainders $ [r `modulo` fromIntegral m|(r,Positive m)<-rms] + +-- prop_chineseRemainders2 :: Integer -> Positive Integer -> Integer -> Positive Integer -> Bool +-- prop_chineseRemainders2 xv (Positive xm) yv (Positive ym) +-- | Just (SomeMod d) <- c = getMod d == l && [getVal d] == sols +-- | otherwise = null sols +-- where +-- l = lcm xm ym +-- c = chineseRemainders2 (xv `modulo` fromIntegral xm) (yv `modulo` fromIntegral ym) +-- sols = [i|i<-[0..l-1],(i-xv) `mod` xm==0,(i-yv) `mod` ym==0] + + + + +{-# DEPRECATED chineseRemainder, chineseRemainder2 "Consider switching to the more general and correct *remainders* (note the terminal S) functions" #-} + +-- | Given a list @[r_1 `modulo` m_1, ..., r_n `modulo` m_n)]@ of @SomeMod@ +-- pairs, @chineseRemainders@ calculates the intersection between all the +-- congruence classes represented by the @SomeMod@. +-- +-- This result may be another congruence class @Just (r `modulo` n)@ if all +-- congruence classes have a non-empty intersection, @Nothing@ if they do not +-- (such as may, or may not, happen if the moduli are not pairwise co-prime), +-- or a specific rational number @Just (InfMod r)@, if it was one of the parameters +-- and is also a member of all of the other congruence classes. +-- +-- n.b. The result will always be @Nothing@ if there are two distinct @InfMod@ +-- parameters (they do not intersect) or if a parameter is a non-integral @InfMod@ +-- and at least one other is a @SomeMod@ (they too do not intersect). +-- +-- On an empty parameter list, @chineseRemainders@ returns @Just (0 `modulo` 1)@, the +-- congruence class of all integers. + +chineseRemainders :: [SomeMod] -> Maybe SomeMod +chineseRemainders (x:xs) = foldM chineseRemainders2 x xs +chineseRemainders _ = Just (0 `modulo` 1) + +-- | Given a pair of @SomeMod@, @chineseRemainders2@ determines their intersection. +-- if any. This function is the underlying worker for @chineseRemainders@. + +chineseRemainders2 :: SomeMod -> SomeMod -> Maybe SomeMod +chineseRemainders2 xsm@(SomeMod x) ysm@(SomeMod y) + | xm == 1 = Just ysm + | ym == 1 = Just xsm + | Just (SomeMod j) <- i = Just $ (xv+(yv-xv)*xm*getVal j) `modulo` (xn*yn) + | (xv-yv) `mod` gm == 0 = chineseRemainders2 (xv `modulo` xq) (yv `modulo` yq) + | otherwise = Nothing + where + xv = getVal x + xm = getMod x + xn = getNatMod x + + yv = getVal y + ym = getMod y + yn = getNatMod y + + i = invertSomeMod (xm `modulo` yn) + + gn = gcd xn yn + gm = fromIntegral gn + + (xq, yq) = foldr distribute (xn `div` gn, yn `div` gn) $ factorise gm + distribute (p',a) (xo, yo) + | xo `mod` p == 0 = (xo*p^a, yo) + | otherwise = (xo, yo*p^a) + where + p = fromIntegral p' + +chineseRemainders2 x@(SomeMod x') y@(InfMod y') + | denominator y' /= 1 = Nothing + | numerator y' `mod` getMod x' == getVal x' = Just y + | otherwise = Nothing + +chineseRemainders2 x@(InfMod _) y@(SomeMod _) = chineseRemainders2 y x + +chineseRemainders2 x y + | x == y = Just x + | otherwise = Nothing + -- | Given a list @[(r_1,m_1), ..., (r_n,m_n)]@ of @(residue,modulus)@ -- pairs, @chineseRemainder@ calculates the solution to the simultaneous -- congruences @@ -30,6 +130,9 @@ import Math.NumberTheory.Utils (recipMod) -- if all moduli are positive and pairwise coprime. Otherwise -- the result is @Nothing@ regardless of whether -- a solution exists. +-- +-- n.b. the @chineseRemainders@ and @chineseRemainders2@ (note the terminal S) functions +-- will find a solution in this case, if one exist. chineseRemainder :: [(Integer,Integer)] -> Maybe Integer chineseRemainder remainders = foldM addRem 0 remainders where