Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly
Sign up
Appearance settings

Commit d830e46

Browse files
committed
Chinese Remainder Theorem
1 parent b7d2a97 commit d830e46

File tree

2 files changed

+107
-1
lines changed

2 files changed

+107
-1
lines changed

‎lib/ModularArithmetic.hs‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -62,3 +62,9 @@ instance Fractional N where
6262
(/) = coerce divM
6363
recip = coerce recipM
6464
fromRational = undefined
65+
66+
{-
67+
import qualified Data.Vector.Unboxing as U
68+
instance U.Unboxable N where
69+
type Rep N = Int64
70+
-}

‎lib/ModularArithmetic_TypeNats.hs‎

Lines changed: 101 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,10 @@
66
{-# LANGUAGE NoStarIsType #-}
77
module ModularArithmetic_TypeNats where
88
import Data.Int
9-
import GHC.TypeNats (Nat, KnownNat, natVal, type (^), type (+))
9+
import GHC.TypeNats
10+
import qualified Test.QuickCheck as QC
11+
import Data.Proxy
12+
import Control.Exception (assert)
1013

1114
--
1215
-- Modular Arithmetic
@@ -68,3 +71,100 @@ instance KnownNat m => Fractional (IntMod m) where
6871
_ -> error "not invertible"
6972
where modulus = fromIntegral (natVal t)
7073
fromRational = undefined
74+
75+
{-
76+
import qualified Data.Vector.Unboxing as U
77+
instance U.Unboxable (IntMod m) where
78+
type Rep (IntMod m) = Int64
79+
-}
80+
81+
recipM :: (Eq a, Integral a, Show a) => a -> a -> a
82+
recipM !x modulo = case exEuclid x modulo of
83+
(1,a,_) -> a `mod` modulo
84+
(-1,a,_) -> (-a) `mod` modulo
85+
(g,a,b) -> error $ show x ++ "^(-1) mod " ++ show modulo ++ " failed: gcd=" ++ show g
86+
87+
-- |
88+
-- Assumption: @gcd m1 m2 == 1@
89+
--
90+
-- >>> crt 3 6 2 7
91+
-- 9
92+
-- >>> crt 2 5 3 9
93+
-- 12
94+
crt :: (Eq a, Integral a, Show a) => a -> a -> a -> a -> a
95+
crt !a1 !m1 !a2 !m2 = let !(s1,s2) = case exEuclid m2 m1 of
96+
(1,b,c) -> (b `mod` m1, c `mod` m2)
97+
(-1,b,c) -> ((-b) `mod` m1, (-c) `mod` m2)
98+
(g,a,b) -> error $ "CRT: " ++ show m1 ++ " and " ++ show m2 ++ " not coprime; gcd=" ++ show (abs g)
99+
!_ = assert (s1 == recipM m2 m1) ()
100+
!_ = assert (s2 == recipM m1 m2) ()
101+
c1 = ((a1 `mod` m1) * s1) `rem` m1
102+
c2 = ((a2 `mod` m2) * s2) `rem` m2
103+
m = m1 * m2
104+
result = c1 * m2 + c2 * m1
105+
in if result < m then
106+
result
107+
else
108+
result - m
109+
{-# SPECIALIZE crt :: Int64 -> Int64 -> Int64 -> Int64 -> Int64 #-}
110+
111+
-- |
112+
-- Assumption: @gcd m1 m2 == 1@
113+
crt' :: (KnownNat m1, KnownNat m2) => IntMod m1 -> IntMod m2 -> IntMod (m1 * m2)
114+
crt' x1@(IntMod !a1) x2@(IntMod !a2) = let !(s1,s2) = case exEuclid m2 m1 of
115+
(1,b,c) -> (b `mod` m1, c `mod` m2)
116+
(-1,b,c) -> ((-b) `mod` m1, (-c) `mod` m2)
117+
(g,a,b) -> error $ "CRT: " ++ show m1 ++ " and " ++ show m2 ++ " not coprime; gcd=" ++ show (abs g)
118+
!_ = assert (s1 == recipM m2 m1) ()
119+
!_ = assert (s2 == recipM m1 m2) ()
120+
c1 = (a1 * s1) `rem` m1
121+
c2 = (a2 * s2) `rem` m2
122+
result = c1 * m2 + c2 * m1
123+
in IntMod $ if result < m then
124+
result
125+
else
126+
result - m
127+
where
128+
m1 = fromIntegral (natVal x1)
129+
m2 = fromIntegral (natVal x2)
130+
m = m1 * m2
131+
132+
--
133+
-- Tests
134+
--
135+
136+
instance KnownNat m => QC.Arbitrary (IntMod m) where
137+
arbitrary = let result = IntMod <$> QC.choose (0, m - 1)
138+
m = fromIntegral (natVal (proxy result))
139+
in result
140+
where
141+
proxy :: f (IntMod m) -> Proxy m
142+
proxy _ = Proxy
143+
144+
runTests :: IO ()
145+
runTests = do
146+
QC.quickCheck $ QC.forAll (QC.choose (2, 10^9+9)) $ \m x -> prop_recipM x (QC.Positive m)
147+
QC.quickCheck $ QC.withMaxSuccess 10000 $ QC.forAll (QC.choose (2, 10^9+9)) $ \m1 -> QC.forAll (QC.choose (2, 10^9+9)) $ \m2 x y -> prop_crt x (QC.Positive m1) y (QC.Positive m2)
148+
QC.quickCheck $ QC.withMaxSuccess 10000 $ QC.forAll (QC.choose (2, 10^9+9)) $ \m1 -> QC.forAll (QC.choose (2, 10^9+9)) $ \m2 x y -> prop_crt' x (QC.Positive m1) y (QC.Positive m2)
149+
150+
prop_recipM :: Int64 -> QC.Positive Int64 -> QC.Property
151+
prop_recipM x (QC.Positive m) = gcd x m == 1 && m > 1 && m <= 10^9 + 9 QC.==>
152+
let y = recipM x m
153+
in 0 <= y QC..&&. y < m QC..&&. ((x `mod` m) * recipM x m) `mod` m QC.=== 1
154+
155+
prop_crt :: Int64 -> QC.Positive Int64 -> Int64 -> QC.Positive Int64 -> QC.Property
156+
prop_crt a1 (QC.Positive m1) a2 (QC.Positive m2)
157+
= gcd m1 m2 == 1 QC.==> let r = crt a1 m1 a2 m2
158+
in r `mod` m1 QC.=== a1 `mod` m1 QC..&&. r `mod` m2 QC.=== a2 `mod` m2
159+
160+
prop_crt' :: Int64 -> QC.Positive Int64 -> Int64 -> QC.Positive Int64 -> QC.Property
161+
prop_crt' a1 (QC.Positive m1) a2 (QC.Positive m2)
162+
= gcd m1 m2 == 1 QC.==> case (someNatVal (fromIntegral m1), someNatVal (fromIntegral m2)) of
163+
(SomeNat p1, SomeNat p2) ->
164+
let x = fromIntegral a1 `asIntModProxy` p1
165+
y = fromIntegral a2 `asIntModProxy` p2
166+
IntMod r = crt' x y
167+
in 0 <= r QC..&&. r < m1 * m2 QC..&&. r `mod` m1 QC.=== a1 `mod` m1 QC..&&. r `mod` m2 QC.=== a2 `mod` m2
168+
where
169+
asIntModProxy :: IntMod m -> Proxy m -> IntMod m
170+
asIntModProxy x _ = x

0 commit comments

Comments
(0)

AltStyle によって変換されたページ (->オリジナル) /