:: Extended Euclidean Algorithm and CRT Algorithm :: by Hiroyuki Okazaki , Yosiki Aoki and Yasunari Shidama :: :: Received February 8, 2012 :: Copyright (c) 2012-2018 Association of Mizar Users :: (Stowarzyszenie Uzytkownikow Mizara, Bialystok, Poland). :: This code can be distributed under the GNU General Public Licence :: version 3.0 or later, or the Creative Commons Attribution-ShareAlike :: License version 3.0 or later, subject to the binding interpretation :: detailed in file COPYING.interpretation. :: See COPYING.GPL and COPYING.CC-BY-SA for the full text of these :: licenses, or see http://www.gnu.org/licenses/gpl.html and :: http://creativecommons.org/licenses/by-sa/3.0/. environ vocabularies RECDEF_2, TARSKI, XBOOLE_0, FINSEQ_1, RELAT_1, FUNCT_1, COMPLEX1, NAT_1, MCART_1, ZFMISC_1, SUBSET_1, NUMBERS, INT_1, INT_2, CARD_1, ARYTM_1, XXREAL_0, ARYTM_3, ORDINAL4, NTALGO_1, CARD_3, VALUED_0; notations TARSKI, XBOOLE_0, ZFMISC_1, SUBSET_1, FUNCT_1, XTUPLE_0, MCART_1, FUNCT_2, ORDINAL1, NUMBERS, XCMPLX_0, XXREAL_0, INT_1, INT_2, NAT_1, NAT_D, FINSEQ_1, RVSUM_1, VALUED_0; constructors XXREAL_0, NAT_D, REAL_1, BINARITH, RVSUM_1, XTUPLE_0, FINSEQ_3, RELSET_1, VALUED_0; registrations XREAL_0, RELSET_1, ORDINAL1, INT_1, SUBSET_1, FINSEQ_1, NAT_1, XXREAL_0, XBOOLE_0, VALUED_0, CARD_1, NUMBERS, XTUPLE_0, VALUED_1, SEQM_3; requirements NUMERALS, BOOLE, SUBSET, REAL, ARITHM; begin :: Euclidean Algorithm theorem :: NTALGO_1:1 for x,p be Integer holds (x mod p) mod p = x mod p; definition let a,b be Integer; func ALGO_GCD(a,b) -> Element of NAT means :: NTALGO_1:def 1 ex A,B be sequence of NAT st A.0 = |.a.| & B.0 = |.b.| & (for i be Nat holds A.(i+1) = B.i & B.(i+1) = A.i mod B.i) & it = A. (min*{i where i is Nat: B.i = 0} ); end; theorem :: NTALGO_1:2 ::Th1: for a,b be Integer holds ALGO_GCD(a,b) = a gcd b; begin ::Extended Euclidean Algorithm ::========================================================================================== ::NZMATH source code :: :: def extgcd(x,y): :: """ :: Return a tuple (u,v,d); they are the greatest common divisor d :: of two integers x and y and u,v such that d = x * u + y * v. :: """ :: # Crandall & Pomerance "PRIME NUMBERS",Algorithm 2.1.4 :: a,b,g,u,v,w = 1,0,x,0,1,y :: while w: :: q,t = divmod(g,w) :: a,b,g,u,v,w = u,v,w,a-q*u,b-q*v,t :: if g >= 0: :: return (a,b,g) :: else: :: return (-a,-b,-g) ::========================================================================================== scheme :: NTALGO_1:sch 1 QuadChoiceRec { A,B,C,D() -> non empty set, A0() -> Element of A(),B0() -> Element of B(), C0() -> Element of C(),D0() -> Element of D(), P[object,object,object,object,object,object,object,object,object] }: ex f being sequence of A(),g being sequence of B(), h being sequence of C(),i being sequence of D() st f.0 = A0() & g.0 = B0() &h.0 = C0() & i.0 = D0() & for n being Nat holds P[n,f.n,g.n,h.n,i.n,f.(n+1),g.(n+1),h.(n+1),i.(n+1)] provided for n being Nat,x being Element of A(),y being Element of B(),z being Element of C(),w being Element of D() ex x1 being Element of A(),y1 being Element of B(), z1 being Element of C(),w1 being Element of D() st P[n,x,y,z,w,x1,y1,z1,w1]; definition let x,y be Element of INT; func ALGO_EXGCD(x,y) -> Element of [:INT,INT,INT:] means :: NTALGO_1:def 2 ex g,w,q,t be sequence of INT, a,b,v,u be sequence of INT, istop be Nat st a.0 = 1 & b.0 = 0 & g.0 = x & q.0 =0 & u.0 = 0 & v.0 = 1 & w.0 = y & t.0 =0 & (for i be Nat holds q.(i+1) = g.i div w.i & t.(i+1) = g.i mod w.i & a.(i+1) = u.i & b.(i+1) = v.i & g.(i+1) = w.i & u.(i+1) = a.i - q.(i+1)*u.i & v.(i+1) = b.i - q.(i+1)*v.i & w.(i+1) = t.(i+1)) & istop = min*{i where i is Nat: w.i = 0} & (0 <= g.istop implies it =[a.istop,b.istop,g.istop] ) & (g.istop < 0 implies it =[-(a.istop),-(b.istop),-(g.istop)] ); end; theorem :: NTALGO_1:3 for i2,i1 being Integer st i2 <= 0 holds i1 mod i2 <= 0; theorem :: NTALGO_1:4 for i2,i1 being Integer st i2 < 0 holds -(i1 mod i2) < -i2; theorem :: NTALGO_1:5 for x,y be Integer st |.y.| <> 0 holds |.x mod y.| < |.y.|; theorem :: NTALGO_1:6 for x,y be Element of INT holds ALGO_EXGCD(x,y)`3_3 = x gcd y & ALGO_EXGCD(x,y)`1_3 * x + ALGO_EXGCD(x,y)`2_3 * y = x gcd y; ::========================================================================================== ::NZMATH source code :: def inverse(x,p): :: """ :: This function returns inverse of x for modulo p. :: """ :: x = x % p :: y = gcd.extgcd(p,x) :: if y[2] == 1: :: if y[1] < 0: :: r = p + y[1] :: return r :: else: :: return y[1] :: raise ZeroDivisionError("There is no inverse for %d modulo %d." % (x,p)) ::======================================================================================= definition let x,p be Element of INT; func ALGO_INVERSE(x,p) -> Element of INT means :: NTALGO_1:def 3 for y be Element of INT st y = (x mod p) holds ( ALGO_EXGCD(p,y)`3_3 = 1 implies ( ( ALGO_EXGCD(p,y)`2_3 < 0) implies (ex z be Element of INT st z = ALGO_EXGCD(p,y)`2_3 & it = p + z )) & ( (0 <= ALGO_EXGCD(p,y)`2_3) implies it = ALGO_EXGCD(p,y)`2_3) ) & ( ALGO_EXGCD(p,y)`3_3 <> 1 implies it = {} ); end; theorem :: NTALGO_1:7 for x,p,y be Element of INT st y = x mod p & ALGO_EXGCD(p,y)`3_3 = 1 holds ( ALGO_INVERSE(x,p) * x ) mod p = 1 mod p; begin ::Algorithm for Chinese Remainder Theorem :: def ALGO_CRT(nlist): :: """ :: This function is Chinese Remainder Theorem using Algorithm 2.1.7 :: of C.Pomerance and R.Crandall's book. :: :: For example: :: >>> ALGO_CRT([(1,2),(2,3),(3,5)]) :: 23 :: """ :: r = len(nlist) :: if r == 1 : :: return nlist [ 0 ] [ 0 ] :: :: product = [] :: prodinv = [] :: m = 1 :: for i in range(1,r): :: m = m*nlist[i-1][1] :: c = inverse(m,nlist[i][1]) :: product.append(m) :: prodinv.append(c) :: :: M = product[r-2]*nlist[r-1][1] :: n = nlist[0][0] :: for i in range(1,r): :: u = ((nlist[i][0]-n)*prodinv[i-1]) % nlist[i][1] :: n += u*product[i-1] :: return n % M ::========================================================================================== definition let nlist be non empty FinSequence of [:INT,INT:]; func ALGO_CRT(nlist) -> Element of INT means :: NTALGO_1:def 4 ( len nlist = 1 implies it = (nlist.1)`1) & ( len nlist <> 1 implies ex m,n,prodc,prodi be FinSequence of INT, M0,M be Element of INT st len m = len nlist & len n = len nlist & len prodc = len nlist - 1 & len prodi = len nlist - 1 & m.1 =1 & (for i be Nat st 1<=i & i<=(len m) - 1 holds ex d,x,y be Element of INT st x = (nlist.i)`2 & m.(i+1) = m.i * x & y = m.(i+1) & d = (nlist.(i+1))`2 & prodi.i = ALGO_INVERSE(y,d) & prodc.i = y ) & M0 = (nlist.(len m))`2 & M = (prodc.((len m)-1))*M0 & n.1 = (nlist.1)`1 & (for i be Nat st 1<=i & i<=len m - 1 holds ex u,u0,u1 be Element of INT st u0 = ( nlist.(i+1))`1& u1 = ( nlist.(i+1))`2 & u = ( (u0-n.i)*(prodi.i)) mod u1 & n.(i+1) = n.i + u*(prodc.i)) & it = n.(len m) mod M ); end; theorem :: NTALGO_1:8 for a,b be Integer st b <> 0 holds (a mod b),a are_congruent_mod b; theorem :: NTALGO_1:9 for a,b be Integer st b <>0 holds (a mod b) gcd b = a gcd b; theorem :: NTALGO_1:10 for a,b,c be Integer st c <>0 & a = b mod c & b,c are_coprime holds a,c are_coprime; theorem :: NTALGO_1:11 for nlist be non empty FinSequence of [:INT,INT:], a,b be FinSequence of INT st len a = len b & len a = len nlist & (for i be Nat st i in Seg (len nlist) holds b.i <> 0 ) & (for i be Nat st i in Seg (len nlist) holds (nlist.i)`1 = a.i & (nlist.i)`2 = b.i ) & (for i,j be Nat st i in Seg (len nlist) & j in Seg (len nlist) & i <> j holds b.i,b.j are_coprime ) holds for i be Nat st i in Seg (len nlist) holds ALGO_CRT(nlist) mod b.i = a.i mod b.i; theorem :: NTALGO_1:12 for x,y be Integer, b,m be non empty FinSequence of INT st 2 <=len b & (for i,j be Nat st i in Seg (len b) & j in Seg len b & i <> j holds b.i,b.j are_coprime ) & (for i be Nat st i in Seg len b holds x mod b.i = y mod b.i ) & m.1 = 1 holds for k be Element of NAT st 1<= k & k <= len b & (for i be Nat st 1<=i & i <=k holds m.(i+1) = m.i * b.i ) holds x mod m.(k+1) = y mod m.(k+1); theorem :: NTALGO_1:13 for b be complex-valued FinSequence st len b = 1 holds Product b = b.1; theorem :: NTALGO_1:14 for b be FinSequence of INT ex m be non empty FinSequence of INT st len m = len b + 1 & m.1 = 1 & (for i be Nat st 1<=i & i <=len b holds m.(i+1) = m.i * b.i) & Product b = m.(len b + 1); theorem :: NTALGO_1:15 for nlist be non empty FinSequence of [:INT,INT:], a,b be non empty FinSequence of INT, x,y be Element of INT st len a = len b & len a = len nlist & (for i be Nat st i in Seg (len nlist) holds b.i <> 0 ) & (for i be Nat st i in Seg (len nlist) holds (nlist.i)`1 = a.i & (nlist.i)`2 =b.i ) & (for i,j be Nat st i in Seg (len nlist) & j in Seg (len nlist) & i <> j holds b.i,b.j are_coprime ) & (for i be Nat st i in Seg (len nlist) holds x mod b.i = a.i mod b.i ) & y = Product b holds ALGO_CRT(nlist) mod y = x mod y;