Add ChineseRemainderS functions

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