Chinese remainder theorem: Difference between revisions
added R |
Added local variables |
||
Line 698:
' returns x where (a * x) % b == 1
_Mul_Inv ' ( a b -- n)
a = Pop()▼
y = 1▼
If b@ = 1 Then
Push 1
Return
EndIf
Do While a@ > 1
Loop
If
Push
Return
_Chinese_Remainder ' ( len -- n)
▲ s = 0
▲ a@ = Pop()
For i = 0 Step 1 While i < l▼
b@ = b@ * @(d@)
Next
For
Push
Next
Push (
Return
</lang>
{{out}}
<pre>
Line 744 ⟶ 749:
1000
0 OK, 0:
</pre>
|
Revision as of 16:46, 11 June 2014
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
You are encouraged to solve this task according to the task description, using any language you may know.
Suppose , , , are positive integers that are pairwise coprime. Then, for any given sequence of integers , , , , there exists an integer solving the following system of simultaneous congruences.
Furthermore, all solutions of this system are congruent modulo the product, .
Your task is to write a program to solve a system of linear congruences by applying the Chinese Remainder Theorem. If the system of equations cannot be solved, your program must somehow indicate this. (It may throw an exception or return a special false value.) Since there are infinitely many solutions, the program should return the unique solution where .
Show the functionality of this program by printing the result such that the 's are and the 's are .
Algorithm: The following algorithm only applies if the 's are pairwise coprime.
Suppose, as above, that a solution is required for the system of congruences:
Again, to begin, the product is defined. Then a solution can be found as follows.
For each , the integers and are coprime. Using the Extended Euclidean algorithm we can find integers and such that . Then, one solution to the system of simultaneous congruences is:
and the minimal solution,
- .
Ada
Using the package Mod_Inv from [[1]].
<lang Ada>with Ada.Text_IO, Mod_Inv;
procedure Chin_Rema is
N: array(Positive range <>) of Positive := (3, 5, 7); A: array(Positive range <>) of Positive := (2, 3, 2); Tmp: Positive; Prod: Positive := 1; Sum: Natural := 0;
begin
for I in N'Range loop Prod := Prod * N(I); end loop; for I in A'Range loop Tmp := Prod / N(I); Sum := Sum + A(I) * Mod_Inv.Inverse(Tmp, N(I)) * Tmp; end loop; Ada.Text_IO.Put_Line(Integer'Image(Sum mod Prod));
end Chin_Rema;</lang>
Bracmat
<lang bracmat>( ( mul-inv
= a b b0 q x0 x1 . !arg:(?a.?b:?b0) & ( !b:1 | 0:?x0 & 1:?x1 & whl ' ( !a:>1 & (!b.mod$(!a.!b):?q.!x1+-1*!q*!x0.!x0) : (?a.?b.?x0.?x1) ) & ( !x1:<0&!b0+!x1 | !x1 ) ) )
& ( chinese-remainder
= n a as p ns ni prod sum . !arg:(?n.?a) & 1:?prod & 0:?sum & !n:?ns & whl'(!ns:%?ni ?ns&!prod*!ni:?prod) & !n:?ns & !a:?as & whl ' ( !ns:%?ni ?ns & !as:%?ai ?as & div$(!prod.!ni):?p & !sum+!ai*mul-inv$(!p.!ni)*!p:?sum ) & mod$(!sum.!prod):?arg & !arg )
& 3 5 7:?n & 2 3 2:?a & put$(str$(chinese-remainder$(!n.!a) \n)) );</lang> Output:
23
C
When n are not pariwise coprime, the program crashes due to division by zero, which is one way to convey error. <lang c>#include <stdio.h>
// returns x where (a * x) % b == 1 int mul_inv(int a, int b) { int b0 = b, t, q; int x0 = 0, x1 = 1; if (b == 1) return 1; while (a > 1) { q = a / b; t = b, b = a % b, a = t; t = x0, x0 = x1 - q * x0, x1 = t; } if (x1 < 0) x1 += b0; return x1; }
int chinese_remainder(int *n, int *a, int len) { int p, i, prod = 1, sum = 0;
for (i = 0; i < len; i++) prod *= n[i];
for (i = 0; i < len; i++) { p = prod / n[i]; sum += a[i] * mul_inv(p, n[i]) * p; }
return sum % prod; }
int main(void) { int n[] = { 3, 5, 7 }; int a[] = { 2, 3, 2 };
printf("%d\n", chinese_remainder(n, a, sizeof(n)/sizeof(n[0]))); return 0; }</lang>
Erlang
<lang erlang>-module(crt). -import(lists, [zip/2, unzip/1, foldl/3, sum/1]). -export([egcd/2, mod/2, mod_inv/2, chinese_remainder/1]).
egcd(_, 0) -> {1, 0}; egcd(A, B) ->
{S, T} = egcd(B, A rem B), {T, S - (A div B)*T}.
mod_inv(A, B) ->
{X, Y} = egcd(A, B), if A*X + B*Y =:= 1 -> X; true -> undefined end.
mod(A, M) ->
X = A rem M, if X < 0 -> X + M; true -> X end.
calc_inverses([], []) -> []; calc_inverses([N | Ns], [M | Ms]) ->
case mod_inv(N, M) of undefined -> undefined; Inv -> [Inv | calc_inverses(Ns, Ms)] end.
chinese_remainder(Congruences) ->
{Residues, Modulii} = unzip(Congruences), ModPI = foldl(fun(A, B) -> A*B end, 1, Modulii), CRT_Modulii = [ModPI div M || M <- Modulii], case calc_inverses(CRT_Modulii, Modulii) of undefined -> undefined; Inverses -> Solution = sum([A*B || {A,B} <- zip(CRT_Modulii, [A*B || {A,B} <- zip(Residues, Inverses)])]), mod(Solution, ModPI) end.</lang>
- Output:
16> crt:chinese_remainder([{10,11}, {4,12}, {12,13}]). 1000 17> crt:chinese_remainder([{10,11}, {4,22}, {9,19}]). undefined 18> crt:chinese_remainder([{2,3}, {3,5}, {2,7}]). 23
Forth
Tested with GNU FORTH <lang forth>: egcd ( a b -- a b )
dup 0= IF 2drop 1 0 ELSE dup -rot /mod \ -- b r=a%b q=a/b -rot recurse \ -- q (s,t) = egcd(b, r) >r swap r@ * - r> swap \ -- t (s - q*t) THEN ;
- egcd>gcd ( a b x y -- n ) \ calculate gcd from egcd
rot * -rot * + ;
- mod-inv ( a m -- a' ) \ modular inverse with coprime check
2dup egcd over >r egcd>gcd r> swap 1 <> -24 and throw ;
- array-product ( adr count -- n )
1 -rot cells bounds ?DO i @ * cell +LOOP ;
- crt-from-array ( adr1 adr2 count -- n )
2dup array-product locals| M count m[] a[] | 0 \ result count 0 DO m[] i cells + @ dup M swap / dup rot mod-inv * a[] i cells + @ * + LOOP M mod ;
create crt-residues[] 10 cells allot create crt-moduli[] 10 cells allot
- crt ( .... n -- n ) \ takes pairs of "n (mod m)" from stack.
10 min locals| n | n 0 DO crt-moduli[] i cells + ! crt-residues[] i cells + ! LOOP crt-residues[] crt-moduli[] n crt-from-array ;
</lang>
- Output:
Gforth 0.7.2, Copyright (C) 1995-2008 Free Software Foundation, Inc. Gforth comes with ABSOLUTELY NO WARRANTY; for details type `license' Type `bye' to exit 10 11 4 12 12 13 3 crt . 1000 ok 10 11 4 22 9 19 3 crt . :2: Invalid numeric argument 10 11 4 22 9 19 3 >>>crt<<< .
Go
Go has the Extended Euclidean algorithm in the GCD function for big integers in the standard library. GCD will return 1 only if numbers are coprime, so a result != 1 indicates the error condition. <lang go>package main
import (
"fmt" "math/big"
)
var one = big.NewInt(1)
func crt(a, n []*big.Int) (*big.Int, error) {
p := new(big.Int).Set(n[0]) for _, n1 := range n[1:] { p.Mul(p, n1) } var x, q, s, z big.Int for i, n1 := range n { q.Div(p, n1) z.GCD(nil, &s, n1, &q) if z.Cmp(one) != 0 { return nil, fmt.Errorf("%d not coprime", n1) } x.Add(&x, s.Mul(a[i], s.Mul(&s, &q))) } return x.Mod(&x, p), nil
}
func main() {
n := []*big.Int{ big.NewInt(3), big.NewInt(5), big.NewInt(7), } a := []*big.Int{ big.NewInt(2), big.NewInt(3), big.NewInt(2), } fmt.Println(crt(a, n))
}</lang>
- Output:
Two values, the solution x and an error value.
23 <nil>
Haskell
<lang haskell>module CRT ( chineseRemainder ) where
egcd :: Integral a => a -> a -> (a, a) egcd _ 0 = (1, 0) egcd a b = (t, s - q * t)
where (s, t) = egcd b r (q, r) = a `quotRem` b
modInv :: Integral a => a -> a -> Maybe a modInv a b = case egcd a b of
(x, y) | a * x + b * y == 1 -> Just x | otherwise -> Nothing
chineseRemainder :: Integral a => [a] -> [a] -> Maybe a chineseRemainder residues modulii = do
inverses <- sequence $ zipWith modInv crtModulii modulii return . (`mod` modPI) . sum $ zipWith (*) crtModulii (zipWith (*) residues inverses) where modPI = product modulii crtModulii = map (modPI `div`) modulii
</lang>
- Output:
λ> chineseRemainder [10, 4, 12] [11, 12, 13] Just 1000 λ> chineseRemainder [10, 4, 9] [11, 22, 19] Nothing λ> chineseRemainder [2, 3, 2] [3, 5, 7] Just 23
Icon and Unicon
with error check added.
Works in both languages: <lang unicon>link numbers # for gcd()
procedure main()
write(cr([3,5,7],[2,3,2]) | "No solution!") write(cr([10,4,9],[11,22,19]) | "No solution!")
end
procedure cr(n,a)
if 1 ~= gcd(n[i := !*n],a[i]) then fail # Not pairwise coprime (prod := 1, sm := 0) every prod *:= !n every p := prod/(ni := n[i := !*n]) do sm +:= a[i] * mul_inv(p,ni) * p return sm%prod
end
procedure mul_inv(a,b)
if b = 1 then return 1 (b0 := b, x0 := 0, x1 := 1) while q := (1 < a)/b do { (t := a, a := b, b := t%b) (t := x0, x0 := x1-q*t, x1 := t) } return if x1 < 0 then x1+b0 else x1
end</lang>
Output:
->crt 23 No solution! ->
J
Solution (brute force):<lang j> crt =: (1 + ] - {:@:[ -: {.@:[ | ])^:_&0@:,:</lang> Example: <lang j> 3 5 7 crt 2 3 2 23
11 12 13 crt 10 4 12
1000</lang> Notes: This is a brute force approach and does not meet the requirement for explicit notification of an an unsolvable set of equations (it just spins forever). A much more thorough and educational approach can be found on the J wiki's Essay on the Chinese Remainder Thereom.
Mathematica
Very easy, because it is a built-in function: <lang Mathematica >ChineseRemainder[{2, 3, 2}, {3, 5, 7}] 23</lang>
OCaml
This is using the Jane Street Ocaml Core library. <lang ocaml>open Core.Std open Option.Monad_infix
let rec egcd a b =
if b = 0 then (1, 0) else let q = a/b and r = a mod b in let (s, t) = egcd b r in (t, s - q*t)
let mod_inv a b =
let (x, y) = egcd a b in if a*x + b*y = 1 then Some x else None
let calc_inverses ns ms =
let rec list_inverses ns ms l = match (ns, ms) with | ([], []) -> Some l | ([], _) | (_, []) -> assert false | (n::ns, m::ms) -> let inv = mod_inv n m in match inv with | None -> None | Some v -> list_inverses ns ms (v::l) in list_inverses ns ms [] >>= fun l -> Some (List.rev l)
let chinese_remainder congruences =
let (residues, modulii) = List.unzip congruences in let mod_pi = List.reduce_exn modulii ~f:( * ) in let crt_modulii = List.map modulii ~f:(fun m -> mod_pi / m) in calc_inverses crt_modulii modulii >>= fun inverses -> Some (List.map3_exn residues inverses crt_modulii ~f:(fun a b c -> a*b*c) |> List.reduce_exn ~f:(+) |> fun n -> let n' = n mod mod_pi in if n' < 0 then n' + mod_pi else n')
</lang>
- Output:
utop # chinese_remainder [(10, 11); (4, 12); (12, 13)];; - : int option = Some 1000 utop # chinese_remainder [(10, 11); (4, 22); (9, 19)];; - : int option = None
PARI/GP
<lang parigp>chivec(residues, moduli)={
my(m=Mod(0,1)); for(i=1,#residues, m=chinese(Mod(residues[i],moduli[i]),m) ); lift(m)
}; chivec([2,3,2], [3,5,7])</lang>
- Output:
23
Pari's chinese function takes a vector in the form [Mod(a1,n1), Mod(a2, n2), ...], so we can do this directly: <lang parigp>lift( chinese([Mod(2,3),Mod(3,5),Mod(2,7)]) )</lang> or to take the residue/moduli array as above: <lang parigp>chivec(residues,moduli)={
lift(chinese(vector(#residues,i,Mod(residues[i],moduli[i]))))
}</lang>
Perl 6
<lang perl6># returns x where (a * x) % b == 1 sub mul-inv($a is copy, $b is copy) {
return 1 if $b == 1; my ($b0, @x) = $b, 0, 1; ($a, $b, @x) = (
$b, $a % $b, @x[1] - ($a div $b)*@x[0], @x[0]
) while $a > 1; @x[1] += $b0 if @x[1] < 0; return @x[1];
}
sub chinese-remainder(*@n) {
my \N = [*] @n; -> *@a {
N R% [+] map { my \p = N div @n[$_]; @a[$_] * mul-inv(p, @n[$_]) * p }, ^@n
}
}
say chinese-remainder(3, 5, 7)(2, 3, 2);</lang>
- Output:
23
Python
<lang python>def mul_inv(a, b):
b0 = b x0, x1 = 0, 1 if b == 1: return 1 while a > 1: q = a / b a, b = b, a%b x0, x1 = x1 - q * x0, x0 if x1 < 0: x1 += b0 return x1
def chinese_remainder(n, a, lena):
p = i = prod = 1; sm = 0 for i in range(lena): prod *= n[i] for i in range(lena): p = prod / n[i] sm += a[i] * mul_inv(p, n[i]) * p return sm % prod
if __name__ == '__main__':
n = [3, 5, 7] a = [2, 3, 2] print(chinese_remainder(n, a, len(n)))</lang>
- Output:
23
R
<lang rsplus>mul_inv <- function(a, b) {
b0 <- b x0 <- 0L x1 <- 1L if (b == 1) return(1L) while(a > 1){ q <- as.integer(a/b) t <- b b <- a %% b a <- t t <- x0 x0 <- x1 - q*x0 x1 <- t } if (x1 < 0) x1 <- x1 + b0 return(x1)
}
chinese_remainder <- function(n, a) {
len <- length(n) prod <- 1L sum <- 0L for (i in 1:len) prod <- prod * n[i] for (i in 1:len){ p <- as.integer(prod / n[i]) sum <- sum + a[i] * mul_inv(p, n[i]) * p } return(sum %% prod)
}
n <- c(3L, 5L, 7L) a <- c(2L, 3L, 2L)
chinese_remainder(n, a)</lang>
- Output:
23
Racket
This is more of a demonstration of the built-in function "solve-chinese", than anything. A bit cheeky, I know... but if you've got a dog, why bark yourself?
Take a look in the "math/number-theory" package it's full of goodies! URL removed -- I can't be doing the Dutch recaptchas I'm getting. <lang racket>#lang racket (require (only-in math/number-theory solve-chinese)) (define as '(2 3 2)) (define ns '(3 5 7)) (solve-chinese as ns)</lang>
- Output:
23
REXX
algebraic
<lang rexx>/*REXX program uses the Chinese Remainder Theorem (Sun Tzu). */ parse arg Ns As . /*get optional arguments from CL.*/ if Ns== then Ns = '3,5,7' /*Ns not specified? Use default.*/ if As== then As = '2,3,2' /*As " " " " */ Ns=space(translate(Ns,,',')); #=words(Ns) /*elide superfluous blanks*/ As=space(translate(As,,',')); _=words(As) /* " " " */ if #\==_ then do; say "size of number sets don't match."; exit 131; end if #==0 then do; say "size of the N set isn't valid."; exit 132; end if #==0 then do; say "size of the A set isn't valid."; exit 133; end N=1 /*the product-to-be for prod(n.j)*/
do j=1 for # /*process each number for As, Ns.*/ n.j=word(Ns,j); N=N*n.j /*get an N.j and calculate prod*/ a.j=word(As,j) /* " " A.j from the As. */ end /*j*/
do x=1 for N /*use a simple algebraic method. */ do i=1 for # /*process each A.i number. */ if x//n.i\==a.i then iterate x /*is the modulus correct for F ? */ end /*i*/ /* [↑] limit solution to product*/ say 'found a solution with x=' x /*announce a possible solution. */ exit /*stick a fork in it, we're done.*/ end /*x*/
say 'no solution found.' /*oops, announce that ¬ found. */
/*stick a fork in it, we're done.*/</lang>
output
found a solution with x= 23
congruences sets
<lang rexx>/*REXX program uses the Chinese Remainder Theorem (Sun Tzu). */ parse arg Ns As . /*get optional arguments from CL.*/ if Ns== then Ns = '3,5,7' /*Ns not specified? Use default.*/ if As== then As = '2,3,2' /*As " " " " */ Ns=space(translate(Ns,,',')); #=words(Ns) /*elide superfluous blanks*/ As=space(translate(As,,',')); _=words(As) /* " " " */ if #\==_ then do; say "size of number sets don't match."; exit 131; end if #==0 then do; say "size of the N set isn't valid."; exit 132; end if #==0 then do; say "size of the A set isn't valid."; exit 133; end N=1 /*the product-to-be for prod(n.j)*/
do j=1 for # /*process each number for As, Ns.*/ n.j=word(Ns,j); N=N*n.j /*get an N.j and calculate prod*/ a.j=word(As,j) /* " " A.j from the As. */ end /*j*/
@.= /* [↓] converts congruences─►sets*/
do i=1 for #; _=a.i; @.i._=a.i; p=a.i do N; p=p+n.i; @.i.p=p; end /*build a list of modulo values. */ end /*i*/ /* [↓] find common number in sets*/ do x=1 for N; if @.1.x== then iterate /*find a number.*/ do v=2 to #; if @.v.x== then iterate x; end /*In all sets ? */ say 'found a solution with X=' x /*we found the lowest solution. */ exit /*stick a fork in it, we're done.*/ end /*x*/
say 'no solution found.' /*oops, there's not a solution. */
/*stick a fork in it, we're done.*/</lang>
output
found a solution with X= 23
Ruby
Brute-force. <lang ruby> def chinese_remainder(mods, remainders)
max = mods.inject( :* ) series = remainders.zip( mods ).map{|r,m| r.step( max, m ).to_a } series.inject( :& ).first #returns nil when empty
end
p chinese_remainder([3,5,7], [2,3,2]) #=> 23 p chinese_remainder([10,4,9], [11,22,19]) #=> nil </lang>
Tcl
<lang tcl>proc ::tcl::mathfunc::mulinv {a b} {
if {$b == 1} {return 1} set b0 $b; set x0 0; set x1 1 while {$a > 1} {
set x0 [expr {$x1 - ($a / $b) * [set x1 $x0]}] set b [expr {$a % [set a $b]}]
} incr x1 [expr {($x1 < 0) * $b0}]
} proc chineseRemainder {nList aList} {
set sum 0; set prod [::tcl::mathop::* {*}$nList] foreach n $nList a $aList {
set p [expr {$prod / $n}] incr sum [expr {$a * mulinv($p, $n) * $p}]
} expr {$sum % $prod}
} puts [chineseRemainder {3 5 7} {2 3 2}]</lang>
- Output:
23
uBasic/4tH
<lang>@(000) = 3 : @(001) = 5 : @(002) = 7 @(100) = 2 : @(101) = 3 : @(102) = 2
Push 3 : GoSub _Chinese_Remainder Print Pop()
' -------------------------------------
@(000) = 11 : @(001) = 12 : @(002) = 13 @(100) = 10 : @(101) = 04 : @(102) = 12
Push 3 : GoSub _Chinese_Remainder Print Pop()
' -------------------------------------
End
' returns x where (a * x) % b == 1
_Mul_Inv ' ( a b -- n)
Local (6)
b@ = Pop() a@ = Pop()
c@ = b@ d@ = 0 e@ = 1
If b@ = 1 Then Push 1 Return EndIf
Do While a@ > 1 f@ = a@ / b@ Push b@ : b@ = a@ % b@ : a@ = Pop() Push d@ : d@ = e@ - f@ * d@ : e@ = Pop() Loop
If e@ < 0 Then e@ = e@ + c@
Push e@
Return
_Chinese_Remainder ' ( len -- n)
Local (5)
a@ = Pop() b@ = 1 c@ = 0
For d@ = 0 Step 1 While d@ < a@ b@ = b@ * @(d@) Next
For d@ = 0 Step 1 While d@ < a@ e@ = b@ / @(d@) Push e@, @(d@) : GoSub _Mul_Inv c@ = c@ + (@(100 + d@) * Pop() * e@) Next
Push (c@ % b@)
Return </lang>
- Output:
23 1000 0 OK, 0:1054
zkl
Using the GMP library, gcdExt is the Extended Euclidean algorithm. <lang zkl>var BN=Import("zklBigNum"), one=BN(1);
fcn crt(xs,ys){
p:=xs.reduce('*,BN(1)); X:=BN(0); foreach x,y in (xs.zip(ys)){ q:=p/x; z,s,_:=q.gcdExt(x); if(z!=one) throw(Exception.ValueError("%d not coprime".fmt(x))); X+=y*s*q; } return(X % p);
}</lang> <lang zkl>println(crt(T(3,5,7), T(2,3,2))); //-->23 println(crt(T(11,12,13),T(10,4,12))); //-->1000 println(crt(T(11,22,19), T(10,4,9))); //-->ValueError: 11 not coprime</lang>