1.1 --- a/src/Tools/isac/Knowledge/Rational2.thy Thu Sep 20 10:07:02 2012 +0200
1.2 +++ b/src/Tools/isac/Knowledge/Rational2.thy Mon Sep 24 09:07:38 2012 +0200
1.3 @@ -1,6 +1,5 @@
1.4 (* rationals, i.e. fractions of multivariate polynomials over the real field
1.5 - author: Diana Meindl
1.6 - (c) Diana Meindl
1.7 + author: (c) Diana Meindl
1.8 Use is subject to license terms.
1.9 1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890
1.10 10 20 30 40 50 60 70 80 90 100
1.11 @@ -8,712 +7,660 @@
1.12
1.13 theory Rational2 imports Complex_Main begin
1.14
1.15 +chapter {* The Isabelle Theory *}
1.16 text {*
1.17 - Below we reference \cite{winkler96}
1.18 - F. Winkler, Polynomial Algorithms. ...
1.19 - by page. *}
1.20 + Below we reference \cite{winkler96} by page. From this book we adopt the
1.21 + algorithms as closely as possible. Isabelle coding standards \cite{TODO isarimpl}
1.22 + require to change the identifiers; further changes concern the transition from
1.23 + Winkler's imperative pseudocode to efficient functional code.
1.24
1.25 -subsection {* Foundamentals and auxiliary algorithms *}
1.26 -
1.27 -ML{*(* use modInvers(r,m) to get rinv= r^-1 mod m *)
1.28 -fun modi( r , rold , m , rinv )=
1.29 -if rinv < m
1.30 -then
1.31 - if r mod m = 1
1.32 - then ( rinv = rinv + 1 ;
1.33 - rinv )
1.34 - else modi ( rold * ( rinv + 1 ) , rold , m , rinv + 1 )
1.35 -else ( 0 )
1.36 -(*Errormeldung? es gibt kein Inverses sind nicht Teilerfremd!! schon zu Beginn darauf überprüfen!!*)
1.37 -
1.38 -fun modInvers ( r , m ) =
1.39 -modi ( r , r , m , 1 ) ;
1.40 -*}
1.41 -ML{*
1.42 -(* modInvers(r,m) --> rinv= r^-1 mod m *)
1.43 -*}
1.44 -ML{*
1.45 -modInvers(2,5);
1.46 -modInvers(3,5);
1.47 -modInvers(4,339);
1.48 -4*85
1.49 + However, in test/../reational.sml there is an almost original version of the
1.50 + algorithms, including the identifiers; this version was the first one
1.51 + developed in this thesis, the version is kept for simple reference to
1.52 + \cite{winkler96}.
1.53 *}
1.54
1.55 -ML{*(* r = r1 mod m1 and r= r2 mod m2 *)
1.56 -fun CRA_2 ( r1 : int , r2 : int , m1 : int , m2 : int ) =
1.57 -(*c=m1^-1 mod m2
1.58 - r1'=r1 mod m1
1.59 - s=(r2-r1')*c mod m2
1.60 - r= r1'+s*m1
1.61 - return r*)
1.62 - (r1 mod m1) + ((r2 - (r1 mod m1)) * (modInvers (m1 ,m2)) mod m2) * m1
1.63 +section {* Auxiliary Functions *}
1.64 +ML {*
1.65 + type unipoly = int list;
1.66 *}
1.67
1.68 +subsection {* modulo calculations etc. *}
1.69 +ML {*
1.70 + (*The "inverse modulo" functions: "mod_inv 5 7 = 3"
1.71 + because 5*3 mod 7 = 1, or more precisely 5^-1 mod 7 = 3.*)
1.72 + fun mod_inv r m =
1.73 + let
1.74 + fun modi (r, rold, m, rinv) =
1.75 + if rinv < m
1.76 + then
1.77 + if r mod m = 1
1.78 + then rinv
1.79 + else modi (rold * (rinv + 1), rold, m, rinv + 1)
1.80 + else 0
1.81 + in modi (r, r, m, 1) end
1.82
1.83 + fun leading_coeff (uvp: unipoly) =
1.84 + if nth uvp (length uvp - 1) <> 0
1.85 + then nth uvp (length uvp - 1)
1.86 + else leading_coeff (nth_drop (length uvp - 1) uvp);
1.87 +
1.88 + fun deg (uvp:unipoly) =
1.89 + if nth uvp (length uvp - 1) <> 0
1.90 + then length uvp - 1
1.91 + else deg (nth_drop (length uvp - 1) uvp)
1.92 +
1.93 + (* gcd is: greatest common divisor *)
1.94 + fun gcd a b = if b = 0 then a else gcd b (a mod b);
1.95 +
1.96 + (* Chinese remainder: given r1, r2, m1, m2
1.97 + find r such that r = r1 mod m1 and r = r2 mod m2 *)
1.98 + fun CRA_2 (r1: int, m1: int, r2: int, m2: int) =
1.99 + (r1 mod m1) + ((r2 - (r1 mod m1)) * (mod_inv m1 m2) mod m2) * m1
1.100 +*}
1.101 +
1.102 +subsection {* primes *}
1.103 +ML{*
1.104 + (* is n divisble by some number in primelist ? *)
1.105 + fun is_prime primelist n =
1.106 + if length primelist > 0
1.107 + then
1.108 + if (n mod (nth primelist 0)) = 0
1.109 + then false
1.110 + else is_prime (nth_drop 0 primelist) n
1.111 + else true
1.112 +
1.113 + (* sieve of erathostenes *)
1.114 + fun primeList p =
1.115 + let
1.116 + fun make_primelist list last p =
1.117 + if (nth list (length list - 1)) < p
1.118 + then
1.119 + if ( is_prime list (last + 2))
1.120 + then make_primelist (list @ [(last + 2)]) (last + 2) p
1.121 + else make_primelist list (last + 2) p
1.122 + else list
1.123 + in
1.124 + if p < 3
1.125 + then [2]
1.126 + else make_primelist [2,3] 3 p
1.127 + end
1.128 +*}
1.129 +text {*
1.130 +fun primeList' 2 =[2,3]|
1.131 + primeList' x= if x<2 then [2] else
1.132 + let
1.133 + fun make_primelist list last x =
1.134 + if (nth list (length list -1)) < x
1.135 + then
1.136 + if ( is_prime list (last + 2))
1.137 + then make_primelist (list @ [(last + 2)]) (last + 2) x
1.138 + else make_primelist list (last + 2) x
1.139 + else list
1.140 + in
1.141 + make_primelist [2,3] 3 x end
1.142 +*}
1.143 +
1.144 +text {*
1.145 + To find now a new prime number not dividing a number d, we find the next prime number higher than the old one with (primeList old+1) and test again if it divide our d. If not we found a new prime not dividing d, else we repeat it with the next prime.\\
1.146 +*}
1.147 +
1.148 +ML {*
1.149 +fun find_next_prime_not_divide a p =
1.150 + let
1.151 + val next = nth (primeList (p+1)) (length(primeList(p+1))-1) ;
1.152 + in
1.153 + if a mod next <> 0
1.154 + then next
1.155 + else find_next_prime_not_divide a next
1.156 +end
1.157 +*}
1.158 +
1.159 +text{*
1.160 +find_next_prime_not_divide 12 2;
1.161 +find_next_prime_not_divide 15 2;
1.162 +find_next_prime_not_divide 90 7;
1.163 +find_next_prime_not_divide 90 1;
1.164 +find_next_prime_not_divide 2 7;
1.165 +primeList 15;
1.166 +primeList 1;
1.167 +primeList 8;
1.168 +primeList 2;
1.169 +primeList 0;
1.170 +primeList' 15;
1.171 +primeList' 1;
1.172 +primeList' 8;
1.173 +primeList' 2;
1.174 +*}
1.175 +
1.176 +ML {*
1.177 +is_prime (primeList 7) 9;
1.178 +*}
1.179 +
1.180 +ML {*
1.181 +fun landau_mignotte_bound a b = 2603;
1.182 +(* 2^min(m,n)gcd(leading_coeff(a),leading_coeff(b)) min(1/(abs leading_coeff(a)) sqrt(sum(a_^2), 1/(abs leading_coeff(b)) sqrt(sum(b_i^2)) *)
1.183 +*}
1.184 +
1.185 +ML {*
1.186 +infix poly_mod
1.187 +fun uvp poly_mod m =
1.188 + let fun poly_mod' uvp m n =
1.189 + if n<length (uvp)
1.190 + then poly_mod' ((nth_drop 0 uvp)@[(nth uvp 0)mod m]) m ( n + 1 )
1.191 + else uvp (*end of poly_mod'*)
1.192 + in
1.193 + poly_mod' uvp m 0 end (*poly_mod*)
1.194 +*}
1.195 +
1.196 +ML {*
1.197 +infix %*
1.198 +fun (a:unipoly) %* (b:int) =
1.199 + let fun mul a b = a*b;
1.200 + in map2 mul a (replicate (length(a)) b)end
1.201 +
1.202 +infix %/
1.203 +fun (poly:unipoly) %/ (b:int) = (* =quotient*)
1.204 + let fun division poly b = poly div b;
1.205 + in map2 division poly (replicate (length(poly)) b) end
1.206 +*}
1.207 +
1.208 +
1.209 +
1.210 +ML {*
1.211 +
1.212 +infix %-%
1.213 +fun a %-% b =
1.214 + let fun minus a b = a-b;
1.215 + in map2 minus a b end
1.216 +
1.217 +infix %-
1.218 +fun p %- b =
1.219 + let fun minus p b = p-b;
1.220 +in map2 minus p (replicate (length(p)) b) end
1.221 +
1.222 +infix %+
1.223 +fun a %+ b =
1.224 + let fun plus a b = a+b;
1.225 + in map2 plus a b end
1.226 +*}
1.227
1.228
1.229 ML{*
1.230 -(*. univariate polynomials (uv) .*)
1.231 -(*. 5 * x^5 + 4 * x^3 + 2 * x^2 + x + 19 => [19,1,2,4,0,5] .*)
1.232 -type uv_poly = int list;
1.233 -(*. multivariate polynomials .*)
1.234 -(*. 5 * x^5 * y^3 + 4 * x^3 * z^2 + 2 * x^2 * y * z^3 - z - 19
1.235 - => [(5,[5,3,0]),(4,[3,0,2]),(2,[2,1,3]),(~1,[0,0,1]),(~19,[0,0,0])] .*)
1.236 -type mv_poly = ( int * int list ) list
1.237 +fun CRA_2_poly (ma, mb) (r1, r2) =
1.238 +let
1.239 + fun CRA_2' (ma, mb) (r1, r2) = CRA_2 (r1, ma,r2, mb)
1.240 +in map (CRA_2' (ma, mb)) (r1 ~~ r2)end
1.241 +*}
1.242 +
1.243 +ML {*(*a/b mod m*)
1.244 +fun mod_div (a:int) (b:int) (m:int) =
1.245 + a*(mod_inv b m) mod m
1.246 +
1.247 +fun poly_mod_div (a:unipoly) (b:int) m =
1.248 + let fun lm_div a b c m i =
1.249 + if length c = length a
1.250 + then c
1.251 + else lm_div a b (c @ [ mod_div (nth a i) b m ]) m (i+1)
1.252 + in lm_div a b [] m 0 end
1.253 +*}
1.254 +
1.255 +ML {*
1.256 +fun mod_mult a b m = (a * b) mod m;
1.257 +fun mod_minus a b m = (a - b) mod m;
1.258 +fun mod_plus a b m = (a + b) mod m;
1.259 *}
1.260
1.261 ML{*
1.262 -val a=[15,1,2,4];
1.263 -val b= [(5,[5,3,0]),(4,[3,0,2]),(2,[2,1,3]),(~1,[0,0,1]),(~19,[0,0,0])];
1.264 -val u:uv_poly=a;
1.265 +fun lc_of_list_not_0 [] = [](* and delete leading_coeff=0*)
1.266 +| lc_of_list_not_0 (uvp:unipoly) =
1.267 + if nth uvp (length uvp -1) =0
1.268 + then lc_of_list_not_0 (nth_drop (length uvp-1) uvp)
1.269 + else
1.270 + if nth uvp 0 = 0
1.271 + then lc_of_list_not_0 (nth_drop 0 uvp)
1.272 + else uvp;
1.273 *}
1.274
1.275 +ML {* (* if poly2|poly1 *)
1.276 +infix %/%
1.277 +fun [a:int] %/% [b:int] =
1.278 + if b<a andalso a mod b = 0
1.279 + then true
1.280 + else false
1.281 +| (poly1:unipoly) %/% (poly2:unipoly) =
1.282 + let
1.283 + val b = (replicate (length poly1 - length(lc_of_list_not_0 poly2)) 0) @ lc_of_list_not_0 poly2 ;
1.284 + val c = leading_coeff poly1 div leading_coeff b;
1.285 + val rest = lc_of_list_not_0 (poly1 %-% (b %* c));
1.286 + in
1.287 + if rest = []
1.288 + then true
1.289 + else
1.290 + if c<>0 andalso length rest >= length poly2
1.291 + then rest %/% poly2
1.292 + else false
1.293 + end
1.294 +*}
1.295 +
1.296 +
1.297 +ML {*
1.298 +fun mod_poly_gcd (upoly1:unipoly, upoly2:unipoly, m:int) =
1.299 + let
1.300 + val moda = upoly1 poly_mod m
1.301 + val modb = (replicate (length upoly1 - length(lc_of_list_not_0 upoly2)) 0)
1.302 + @ (lc_of_list_not_0 upoly2 poly_mod m) ;
1.303 + val c = mod_div (leading_coeff moda) (leading_coeff modb) m
1.304 + val rest = lc_of_list_not_0 ((moda %-% (modb %* c)) poly_mod m)
1.305 + in
1.306 + if rest = []
1.307 + then [upoly1, [0] ]
1.308 + else
1.309 + if length rest < length upoly2
1.310 + then mod_poly_gcd (upoly2, rest, m) (*[ rest,[c]@d]*)
1.311 + else mod_poly_gcd (rest, upoly2, m)
1.312 + end
1.313 +*}
1.314 +
1.315 +ML {*
1.316 +fun mod_polynom_gcd (poly1:unipoly) (poly2:unipoly) (m:int) =
1.317 + let
1.318 + val moda = poly1 poly_mod m;
1.319 + val modb = (replicate (length poly1 - length poly2) 0) @ (poly2 poly_mod m) ;
1.320 + val c = mod_div (leading_coeff moda) (leading_coeff modb) m;
1.321 + val rest = lc_of_list_not_0 (moda %-% (modb %* c) poly_mod m) poly_mod m;
1.322 + in (* nth_drop ( length ( upoly1 ) - 1 )*)
1.323 + if rest = []
1.324 + then [poly2, [0]]
1.325 + else
1.326 + if length rest < length poly2
1.327 + then [rest]
1.328 + else mod_poly_gcd (rest, poly2 , m)
1.329 + end
1.330 +*}
1.331
1.332 ML{*
1.333 -CRA_2(2,1,4,5);
1.334 +fun normirt_polynom (poly1:unipoly) (m:int) =
1.335 + let
1.336 + val poly1 = lc_of_list_not_0 poly1
1.337 + val lc_a=leading_coeff poly1;
1.338 + fun normirt poly1 b m lc_a i =
1.339 + if i=0
1.340 + then [mod_div (nth poly1 i) lc_a m]@b
1.341 + else normirt poly1 ( [mod_div(nth poly1 i) lc_a m]@b) m lc_a (i-1) ;
1.342 + in
1.343 + if length(poly1)=0
1.344 + then poly1
1.345 + else normirt poly1 [] m lc_a (length(poly1)-1)
1.346 +end
1.347 *}
1.348
1.349 +ML{*
1.350 +fun mod_gcd_p (poly1:unipoly) (poly2:unipoly) (p:int) =
1.351 + let
1.352 + val rest = mod_polynom_gcd poly1 poly2 p;
1.353 + in
1.354 + if length rest = 2
1.355 + then normirt_polynom (nth rest 0) p
1.356 + else mod_gcd_p poly2 (nth rest 0) p
1.357 +end
1.358 +*}
1.359 +
1.360
1.361 ML{*
1.362 -3 mod 2;
1.363 -((3+5)*5)mod 9;
1.364 -CRA_2(2,6,2,5);
1.365 -modInvers(2,5);
1.366 -modInvers(3,5);
1.367 -CRA_2(2,6,3,5);
1.368 -6 mod 5;
1.369 +fun gcd_n (list:int list) =
1.370 + if length list = 2
1.371 + then gcd (nth list 0)(nth list 1)
1.372 + else gcd_n ([gcd (nth list 0)(nth list 1)]@(nth_drop 0 (nth_drop 0 list)))
1.373 +*}
1.374 +
1.375 +ML {*
1.376 +fun pp (poly1:unipoly) =
1.377 + if (poly1 poly_mod (gcd_n poly1)) = (replicate (length (poly1)) 0)
1.378 + then poly1 %/ (gcd_n poly1)
1.379 + else poly1
1.380
1.381 *}
1.382
1.383 +ML {*
1.384 +pp [4,8,2,6];
1.385 +*}
1.386 +
1.387 +ML {*
1.388 +fun poly_centr (poly:unipoly) (m:int) =
1.389 + let
1.390 + fun midle m =
1.391 + let val mid = m div 2
1.392 + in
1.393 + if m mod mid = 0
1.394 + then mid
1.395 + else mid+1
1.396 + end
1.397 + fun centr a m mid =
1.398 + if a > mid
1.399 + then a - m
1.400 + else a
1.401 + fun polyCentr poly poly' m mid counter =
1.402 + if length(poly) > counter
1.403 + then polyCentr poly (poly' @ [centr (nth poly counter) m mid]) m mid (counter+1)
1.404 + else (poly':unipoly)
1.405 + in polyCentr poly [] m (midle m) 0
1.406 +end
1.407 +*}
1.408 +
1.409 +
1.410 +subsection {* GCD_MOD Algorgithmus *}
1.411 +
1.412 +ML {*
1.413 + fun GCD_MOD a b =
1.414 + let
1.415 +(*1*)val d = gcd (leading_coeff a) (leading_coeff b);
1.416 + val M = 2*d*landau_mignotte_bound a b;
1.417 + fun GOTO2 a b d M p = (*==============================*)
1.418 + let
1.419 +(*2*) val p = find_next_prime_not_divide d p
1.420 + val c_p = normirt_polynom ( mod_gcd_p a b p) p
1.421 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p
1.422 + fun GOTO3 a b d M p g_p = (*~~~~~~~~~~~~~~~~~~~~~~*)
1.423 +(*3*) if (deg g_p) = 0
1.424 + then [1]
1.425 + else
1.426 + let
1.427 + val P = p
1.428 + val g = g_p
1.429 + fun WHILE a b d M P g = (*------------------*)
1.430 +(*4*) if P > M
1.431 + then g
1.432 + else
1.433 + let
1.434 + val p = find_next_prime_not_divide d p
1.435 + val c_p = normirt_polynom ( mod_gcd_p a b p) p
1.436 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p
1.437 + in
1.438 + if deg g_p < deg g
1.439 + then GOTO3 a b d M p g_p
1.440 + else
1.441 + if deg g_p = deg g
1.442 + then
1.443 + let
1.444 + val g = CRA_2_poly (P,p)(g,g_p)
1.445 + val P = P*p
1.446 + val g = poly_centr(g poly_mod P)P
1.447 + in WHILE a b d M P g end
1.448 + else WHILE a b d M P g
1.449 + end (*----------------------------------*)
1.450 +
1.451 +
1.452 + val g = WHILE a b d M P g (* << 1.Mal -----*)
1.453 + in g end (*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
1.454 +
1.455 + val g = GOTO3 a b d M p g_p (* << 1.Mal ~~~~~~~~~~*)
1.456 +(*5*) val g = pp g
1.457 + in
1.458 + if (a %/% g) andalso (b %/% g)
1.459 + then g
1.460 + else GOTO2 a b d M p
1.461 + end (*==============================================*)
1.462 +
1.463 +
1.464 + val g = GOTO2 a b d (*p*)1 M (* << 1. Mal =============*)
1.465 + in g end;
1.466 +*}
1.467
1.468
1.469 ML{*
1.470 -fun lc uv =
1.471 -if nth uv (length uv -1)<>0
1.472 -then nth uv (length uv -1)
1.473 -else lc (nth_drop (length uv-1) uv);
1.474 -
1.475 -fun deg uv = length uv - 1;
1.476 -
1.477 + fun GCD_MOD' a b =
1.478 + let
1.479 +(*1*)val d = gcd (leading_coeff a) (leading_coeff b);
1.480 + val M = 2*d*landau_mignotte_bound a b;
1.481 + fun GOTO2 a b d M p = (*==============================*)
1.482 + let
1.483 +(*2*) val p = find_next_prime_not_divide d p
1.484 + val c_p = normirt_polynom ( mod_gcd_p a b p) p
1.485 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p
1.486 + fun WHILE a b d M P p g = (*------------------*)
1.487 +(*4*) if P > M
1.488 + then
1.489 + if a %/% (pp g) andalso b %/% (pp g)
1.490 + then pp g
1.491 + else GOTO2 a b d M p
1.492 + else
1.493 + let
1.494 + val p = find_next_prime_not_divide d p
1.495 + val c_p = normirt_polynom (mod_gcd_p a b p) p
1.496 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p
1.497 + in
1.498 + if deg g_p < deg g
1.499 + then
1.500 +(*3*) if deg g_p = 0
1.501 + then [1]
1.502 + else WHILE a b d M p p g_p
1.503 + else
1.504 + if deg g_p = deg g
1.505 + then
1.506 + WHILE a b d M (P*p) p ( poly_centr(CRA_2_poly (P,p)(g,g_p))(P*p))
1.507 + else WHILE a b d M P p g
1.508 + end (*--------------WHILE--------------------*)
1.509 + in
1.510 +(*3*) if deg g_p = 0 then [1]
1.511 + else WHILE a b d M p p g_p
1.512 + end (*========================GOTO2======================*)
1.513 +in GOTO2 a b d (*p*)1 M end
1.514 *}
1.515
1.516 -ML{*
1.517 -(* modInvers(r,m) --> rinv= r^-1 mod m *)
1.518 -(* CRA_2(r1,r2,m1,m2) --> r= r1 mod m1 = r2 mod m2 *)
1.519 -(* lv uv --> leading coefficient of univariate polynomials, deg uv *)
1.520 +subsection {* GCD_MODm Algorgithmus *}
1.521 +
1.522 +ML {*
1.523 +type monom = (int * int list);
1.524 +type multipoly = monom list;
1.525 +val mpoly1 = [(~3,[0,2]),(3,[1,0]),(1,[0,5]),(~6,[1,1]),(~1,[1,3]),(2,[1,4]),(1,[2,3]),(~1,[3,1]),(2,[3,2])];
1.526 +val mpoly2 = [(2,[3,1]),(~1,[3,0]),(~1,[2,2]),(1,[2,1]),(~1,[1,3]),(4,[1,1]),(~2,[1,0]),(2,[0,2])];
1.527 +val (test:monom) =(22,[5,6,7]);
1.528 +*}
1.529 +
1.530 +
1.531 +ML {*
1.532 +fun get_coef ((coef,var):monom) = coef;
1.533 +fun get_varlist ((coef,var):monom) = var;
1.534 +fun get_coeflist (poly:multipoly) =
1.535 + let
1.536 + fun get_coeflist' (poly:multipoly) list =
1.537 + if poly = []
1.538 + then
1.539 + list
1.540 + else
1.541 + get_coeflist' (nth_drop 0 poly) (list @ [get_coef(nth poly 0)])
1.542 + in
1.543 + get_coeflist' poly []
1.544 +end
1.545 +
1.546 +(* Eine Funktion für verschiedene Typen? *)
1.547 +fun lc_mv (poly:multipoly) =
1.548 + if get_coef (nth mpoly1 (length mpoly1 - 1)) <> 0
1.549 + then get_coef (nth poly (length poly-1))
1.550 + else lc_mv (nth_drop (length poly-1) poly);
1.551 +fun deg_mv (mvp:multipoly) =
1.552 + if get_coef (nth mvp (length mvp-1)) <> 0
1.553 + then nth (get_varlist (nth mvp (length mvp-1))) 0
1.554 + else deg_mv (nth_drop (length mvp-1) mvp)
1.555 +fun greater_deg (monom1:monom) (monom2:monom)=
1.556 + if deg_mv [monom1] > deg_mv [monom2]
1.557 + then 1
1.558 + else if deg_mv [monom1] < deg_mv [monom2]
1.559 + then 2
1.560 + else if length (get_varlist monom1) >0
1.561 + then
1.562 + greater_deg (1,(nth_drop 0 (get_varlist monom1))) (1,(nth_drop 0 (get_varlist monom2)))
1.563 + else 0
1.564 +
1.565 *}
1.566
1.567 ML {*
1.568 -lc a;
1.569 -lc [2,4,4,1,0,0];
1.570 +get_coef test;
1.571 +get_varlist test;
1.572 +val t = [[1,2,4,5],[4,5,6,7]];
1.573 +nth (nth t 1) 1;
1.574 +lc_mv mpoly1;
1.575 + deg_mv mpoly1;
1.576 + deg_mv [(5,[4,5,6])];
1.577 +get_coeflist mpoly1;
1.578 +val monom =(5,[4,5,6]) ;
1.579 + deg_mv [monom];
1.580 +get_varlist monom;
1.581 + [(1,(nth_drop 0 (get_varlist monom)))];
1.582 +greater_deg (5,[4,5,6]) (5,[4,5,5]);
1.583 *}
1.584
1.585 ML {*
1.586 -fun absolut a =
1.587 -if a <= 0
1.588 - then ~a
1.589 - else a;
1.590 +infix %%/
1.591 +fun (poly:multipoly) %%/ (b:int) = (* =quotient*)
1.592 + let fun division monom b = ((get_coef monom div b),get_varlist monom);
1.593 + in map2 division poly (replicate (length(poly)) b) end
1.594 *}
1.595
1.596 ML {*
1.597 -absolut(3);
1.598 -absolut ~3;
1.599 -5 div 3 ;
1.600 -4 div 10;
1.601 -
1.602 -*}
1.603 -ML{*
1.604 -(* modInvers(r,m) --> rinv= r^-1 mod m *)
1.605 -(* CRA_2(r1,r2,m1,m2) --> r= r1 mod m1 = r2 mod m2 *)
1.606 -(* lv uv --> leading coefficient of univariate polynomials, deg uv *)
1.607 -(* abs x --> |x| *)
1.608 +mpoly1;
1.609 +mpoly1 %%/ 2;
1.610 +~11 div 2;
1.611 +11 div 2;
1.612 +[3,4,~4,~3,~1] %/ 2;
1.613 *}
1.614
1.615 ML {*
1.616 -infix divide
1.617 -fun a divide b =
1.618 - if abs b > 0
1.619 - then abs b divide abs a - abs a div abs b * abs b
1.620 - else abs a
1.621 -fun INT_GCD a b = a divide b;
1.622 -(* funktoniert nicht wenn ich es mit INT_GCD rekursiev mache
1.623 - die Struktur passt ihm dann nicht *)
1.624 +fun cont (poly1:unipoly) = gcd_n poly1
1.625 +fun cont' (poly:multipoly) = gcd_n (get_coeflist poly)
1.626 +*}
1.627 +
1.628 +
1.629 +ML {*
1.630 +[3,4,5] = [3,4,5];
1.631 +[3,4,5]=[5,4,3];
1.632 +[3,4,5] @ [get_coef (nth mpoly1 (length mpoly1 - 1))
1.633 + - get_coef (nth mpoly2 (length mpoly2 - 1))];
1.634 +mpoly1 = [] orelse mpoly2 =[];
1.635 +get_varlist (nth mpoly1 (length mpoly1 - 1)) = get_varlist (nth mpoly2 (length mpoly2 - 1))
1.636 *}
1.637
1.638 ML {*
1.639 -12 divide 6;
1.640 -10 divide 3;
1.641 -15 divide 6;
1.642 -6 divide 15;
1.643 -INT_GCD 6 21;
1.644 +infix %%-
1.645 +fun (mpoly1:multipoly) %%- (mpoly2:multipoly) =
1.646 + let
1.647 + fun min' (mpoly1:multipoly) (mpoly2:multipoly) (result:multipoly) =
1.648 + if mpoly1 = [] andalso mpoly2 =[]
1.649 + then result
1.650 + else if mpoly1 = []
1.651 + then let
1.652 + val result = [(0 - get_coef (nth mpoly2 (length mpoly2 - 1)),
1.653 + get_varlist (nth mpoly2 (length mpoly2 - 1)))] @ result
1.654 + in
1.655 + min' mpoly1 (nth_drop (length mpoly2 - 1) mpoly2) result
1.656 + end
1.657 + else if mpoly2 = []
1.658 + then
1.659 + mpoly1 @ result
1.660 +
1.661 + else
1.662 + if get_varlist (nth mpoly1 (length mpoly1 - 1)) = get_varlist (nth mpoly2 (length mpoly2 - 1))
1.663 + then
1.664 + let
1.665 + val result = [(get_coef (nth mpoly1 (length mpoly1 - 1))
1.666 + - get_coef (nth mpoly2 (length mpoly2 - 1)),
1.667 + get_varlist (nth mpoly1 (length mpoly1 - 1)))] @ result
1.668 + in
1.669 + min' (nth_drop (length mpoly1 - 1) mpoly1) (nth_drop (length mpoly2 - 1) mpoly2) result
1.670 + end
1.671 + else
1.672 + if greater_deg (nth mpoly1 (length mpoly1 - 1)) (nth mpoly2 (length mpoly2 - 1)) = 1
1.673 + then
1.674 + min' (nth_drop (length mpoly1 - 1) mpoly1) mpoly2 ( [nth mpoly1 (length mpoly1 - 1)] @ result)
1.675 + else
1.676 + let
1.677 + val result = [(0 - get_coef (nth mpoly2 (length mpoly2 - 1)),
1.678 + get_varlist (nth mpoly2 (length mpoly2 - 1)))]@ result
1.679 + in
1.680 + min' mpoly1 (nth_drop (length mpoly2 - 1) mpoly2) result
1.681 + end
1.682 + in
1.683 + min' mpoly1 mpoly2 []
1.684 +end
1.685 +*}
1.686 +
1.687 +ML {*
1.688 +[(~6,[1,1]),(~1,[1,3])]%%-[(~4,[1,1]),(3,[1,3])];
1.689 +[(~6,[1,1]),(~1,[1,5])]%%-[(~4,[1,1]),(3,[2,3])];
1.690 +*}
1.691 +
1.692 +ML {*
1.693 +*}
1.694 +
1.695 +ML {*
1.696 +*}
1.697 +
1.698 +ML {*
1.699 +*}
1.700 +
1.701 +ML {*
1.702 +*}
1.703 +
1.704 +ML {*
1.705 +*}
1.706 +
1.707 +ML {*
1.708 +*}
1.709 +
1.710 +ML {*
1.711 +
1.712 +*}
1.713 +text {*
1.714 +
1.715 +fun GCD_MODm a b n s =
1.716 +(*0*) if s = 0
1.717 + then g = (gcd (cont a) (cont b)) %* (GCD_MOD pp(a) pp(b))
1.718 + else
1.719 + let
1.720 +(*1*) val M = 1 + Int.min (DEG_X a, DEG_X b);
1.721 +(*2*) fun GOTO2 a b n s x M =
1.722 + let
1.723 + val r = AN_Integer_DEG_H;
1.724 + val g_r = GCD_MODm (REDUCED a)(REDUCED b) n (s-1)
1.725 +(*3*) fun GOTO3 a b n s x M g_r r=
1.726 + let
1.727 + val m = 1;
1.728 + val g = g_r
1.729 + fun WHILE a b n s x M m r =
1.730 + if m > M
1.731 + then
1.732 + if g TEILT a andalso g TEILT a
1.733 + then g
1.734 + else GOTO2 a b n s x M
1.735 + else
1.736 + let
1.737 + val r = NEW_INTEGER;
1.738 + val g_r = GCD_MODm (REDUCED a)(REDUCED b) n (s-1);
1.739 + in
1.740 + if DEG_XN g_r < DEG_XN g
1.741 + then GOTO3 a b n s x M g_r r
1.742 + else if DEG_XN g_r < DEG_XN g
1.743 + then INCORPORATE g_r
1.744 + else WHILE a b n s x M (m+1) r
1.745 + end (* WHILE *)
1.746 + in (* GOTO3*)
1.747 + WHILE a b n s x M m r
1.748 + end (*GOTO3*)
1.749 + in (* GOTO2*)
1.750 + GOTO3 a b n s x M g_r r
1.751 + end (*GOTO2*)
1.752 + in
1.753 + GOTO2 a b n s x M
1.754 + end (*end 0*)
1.755
1.756 *}
1.757
1.758 ML {*
1.759 -val list = [ 2 , 3 ];
1.760 -length ( list );
1.761 -list @ [(4+2)];
1.762 + val a = [~18, ~15, ~20, 12, 20, ~13, 2];
1.763 + val b = [8, 28, 22, ~11, ~14, 1, 2];
1.764 + if GCD_MOD a b = [~2, ~1, 1] then ()
1.765 + else error "GCD_MOD a b = [~2, ~1, 1]";
1.766 *}
1.767
1.768 -ML{*
1.769 -(* modInvers(r,m) --> rinv= r^-1 mod m *)
1.770 -(* CRA_2(r1,r2,m1,m2) --> r= r1 mod m1 = r2 mod m2 *)
1.771 -(* lv uv --> leading coefficient of univariate polynomials, deg uv *)
1.772 -(* abs x --> |x| *)
1.773 -(* INT_GCD a b *)
1.774 -*}
1.775 -
1.776 -ML {*
1.777 -infix listMod;
1.778 -fun m listMod list =
1.779 -if length( list )>0 andalso m mod nth list 0 <> 0
1.780 -then m listMod nth_drop 0 list
1.781 -else length( list )
1.782 -
1.783 -
1.784 -fun makePrimeList list last x =
1.785 -if last < x
1.786 -then
1.787 -if (( last + 2) listMod list ) = 0
1.788 -then makePrimeList ( list @ [( last + 2)]) ( last + 2 ) x
1.789 -else makePrimeList list ( last + 2 ) x
1.790 -else list
1.791 -
1.792 -fun primeList x =
1.793 -if x<3
1.794 -then [2]
1.795 -else makePrimeList [2,3] 3 x
1.796 -*}
1.797 -
1.798 -
1.799 -
1.800 -
1.801 -ML {*
1.802 -makePrimeList [2,3] 3 7 ;
1.803 -makePrimeList [2,3] 3 5;
1.804 -makePrimeList [2,3] 3 3;
1.805 -makePrimeList [2,3] 3 29;
1.806 -primeList 29;
1.807 -primeList 2;
1.808 -primeList 28;
1.809 -*}
1.810 -
1.811 -ML{*
1.812 -(* modInvers(r,m) --> rinv= r^-1 mod m *)
1.813 -(* CRA_2(r1,r2,m1,m2) --> r= r1 mod m1 = r2 mod m2 *)
1.814 -(* lv uv --> leading coefficient of univariate polynomials, deg uv *)
1.815 -(* abs x --> |x| *)
1.816 -(* INT_GCD a b *)
1.817 -(* m listMod list --> new list where every entry mod m*)
1.818 -(*primeList x --> a list with primenumbers until >= x *)
1.819 -*}
1.820 -
1.821 -ML {*
1.822 -fun nextPrimeList p list list' =
1.823 -if length ( list ) = length ( list' )
1.824 - then
1.825 - if (( p + 2) listMod list ) = 0
1.826 - then nextPrimeList ( p + 2 ) (list @ [( p + 2)]) list'
1.827 - else nextPrimeList ( p +2) list list'
1.828 - else list
1.829 -
1.830 -
1.831 -fun findNextPrime p =
1.832 - let
1.833 - val list = primeList p;
1.834 - val list' = list;
1.835 - in
1.836 - if p = 2
1.837 - then nth (nextPrimeList ( p + 1 ) list list' ) (length( list ))
1.838 - else nth (nextPrimeList p list list') (length( list ))
1.839 - end
1.840 -
1.841 -fun findNextPrimeNotDivide x p =
1.842 -let
1.843 -val next = findNextPrime p ;
1.844 -in
1.845 -if x mod next <> 0
1.846 -then next
1.847 -else findNextPrimeNotDivide x next
1.848 end
1.849 -*}
1.850 -
1.851 -ML {*
1.852 -fun GCD_MOD a b =
1.853 -let
1.854 -val d = INT_GCD ( lc a) ( lc b);
1.855 -val m = 2 * d ;
1.856 -val list = primeList d;
1.857 -val p = findNextPrimeNotDivide d 1;
1.858 -
1.859 -in
1.860 -d * m * p
1.861 -end
1.862 -*}
1.863 -
1.864 -
1.865 -ML {*
1.866 -findNextPrime 33 ;
1.867 -findNextPrimeNotDivide 90 7;
1.868 -findNextPrimeNotDivide 90 1;
1.869 -
1.870 -*}
1.871 -
1.872 -ML{*
1.873 -(* modInvers(r,m) --> rinv= r^-1 mod m *)
1.874 -(* CRA_2(r1,r2,m1,m2) --> r= r1 mod m1 = r2 mod m2 *)
1.875 -(* lc uv --> leading coefficient of univariate polynomials, deg uv *)
1.876 -(* abs x --> |x| *)
1.877 -(* INT_GCD a b *)
1.878 -(* m listMod list --> new list where every entry mod m*)
1.879 -(*primeList x --> a list with primenumbers until >= x *)
1.880 -(*findNextPrime p *)
1.881 -(*findNextPrimeNotDievide x p --> find the primenumber after p wich not divide x*)
1.882 -*}
1.883 -
1.884 -ML {*
1.885 -fun min a b =
1.886 - if a < b
1.887 - then a
1.888 - else b
1.889 -*}
1.890 -
1.891 -ML {*
1.892 -min 3 5;
1.893 -min ~6 5 ;
1.894 -
1.895 -*}
1.896 -
1.897 -ML {*
1.898 -fun landau_mignotte_bound a b = 100000
1.899 -(* 2^min(m,n)gcd(lc(a),lc(b)) min(1/(abs lc(a)) sqrt(sum(a_^2), 1/(abs lc(b)) sqrt(sum(b_i^2)) *)
1.900 -
1.901 -*}
1.902 -
1.903 -ML {*
1.904 -val a=[2,~13,20,12,~20,~15,~18];
1.905 -val b=[2,1,~14,~11,22,28,8];
1.906 -
1.907 -*}
1.908 -
1.909 -
1.910 -ML {*
1.911 -val list = a;
1.912 -nth_drop 1 list;
1.913 -list;
1.914 -(nth_drop 0 list)@[(nth list 0)mod 5];
1.915 -*}
1.916 -ML{*
1.917 -(* modInvers(r,m) --> rinv= r^-1 mod m *)
1.918 -(* CRA_2(r1,r2,m1,m2) --> r= r1 mod m1 = r2 mod m2 *)
1.919 -(* lc uv --> leading coefficient of univariate polynomials, deg uv *)
1.920 -(* abs x --> |x| *)
1.921 -(* INT_GCD a b *)
1.922 -(* m listMod list --> new list where every entry mod m*)
1.923 -(*primeList x --> a list with primenumbers until >= x *)
1.924 -(*findNextPrime p *)
1.925 -(*findNextPrimeNotDievide x p --> find the primenumber after p wich not divide x*)
1.926 -*}
1.927 -ML {*
1.928 -fun list_mod' list m n =
1.929 -if n<length ( list )
1.930 -then list_mod' ((nth_drop 0 list)@[(nth list 0)mod m]) m ( n + 1 )
1.931 -else list
1.932 -
1.933 -infix list_mod
1.934 -fun list list_mod m =
1.935 -list_mod' list m 0
1.936 -
1.937 -fun mul a b = a*b;
1.938 -fun list_mult a b = map2 mul a (replicate (length(a)) b);
1.939 -
1.940 -fun minus a b = a-b;
1.941 -fun list_minus a b = map2 minus a b;
1.942 -
1.943 -fun plus a b = a+b;
1.944 -fun list_plus a b = map2 plus a b;
1.945 -
1.946 -*}
1.947 -ML {*
1.948 -4 div 2;
1.949 -a;
1.950 -b;
1.951 -a list_mod 5;
1.952 -([12,0,9]) list_mod 5;
1.953 -b list_mod 5;
1.954 -*}
1.955 -
1.956 -
1.957 -ML {*
1.958 -nth;
1.959 -val mv= [[1,2,3],[4,5,6]];
1.960 -map2 (nth) mv [0,2];
1.961 -fun modulo a b = a mod b;
1.962 -map2 (modulo) a [5,5,5,5,5,5,5];
1.963 -*}
1.964 -
1.965 -ML {*
1.966 -CRA_2 (4,3,5,5);
1.967 -*}
1.968 -ML{*
1.969 -(* modInvers(r,m) --> rinv= r^-1 mod m *)
1.970 -(* CRA_2(r1,r2,m1,m2) --> r= r1 mod m1 = r2 mod m2 *)
1.971 -(* lc uv --> leading coefficient of univariate polynomials, deg uv *)
1.972 -(* abs x --> |x| *)
1.973 -(* INT_GCD a b *)
1.974 -(* m listMod list --> new list where every entry mod m*)
1.975 -(*primeList x --> a list with primenumbers until >= x *)
1.976 -(*findNextPrime p *)
1.977 -(*findNextPrimeNotDievide x p --> find the primenumber after p wich not divide x*)
1.978 -*}
1.979 -
1.980 -ML{*
1.981 -fun CRA_2_poly a b ma mb =
1.982 -(* ( r1 mod m1 ) + ( (( r2-(r1 mod m1)) * ( modInvers ( m1 , m2 ))) mod m2 ) * m1*)
1.983 -list_plus (a list_mod ma) ((list_mult (list_mult (list_minus b (a list_mod ma))( (modInvers ( ma , mb ))mod mb)list_mod mb) ma) )
1.984 -*}
1.985 -
1.986 -ML {*
1.987 -fun m_div a b m=
1.988 -if a mod b =0
1.989 -then a div b
1.990 -else m_div ( a + m ) b m
1.991 -
1.992 -infix mod_div (* a/b mod m *)
1.993 -fun ( a , b ) mod_div m = m_div a b m ;
1.994 -
1.995 -fun lm_div a b c m i =
1.996 -if length c = length a
1.997 -then c
1.998 -else lm_div a b (c @ [((nth a i), b) mod_div m]) m (i+1)
1.999 -infix list_mod_div
1.1000 -fun ( a , b ) list_mod_div m = lm_div a b [] m 0;
1.1001 -(*fun minus a b = a-b;*)
1.1002 -fun mod_mult a b m = ( a * b ) mod m;
1.1003 -fun mod_minus a b m = ( a - b ) mod m;
1.1004 -fun mod_plus a b m = (a+b) mod m;
1.1005 -
1.1006 -*}
1.1007 -
1.1008 -ML {*
1.1009 -([1,2,3,4,5],3)list_mod_div 5;
1.1010 -m_div 3 2 5;
1.1011 -map2 minus [3,3,3] [2,2,2];
1.1012 -(1,4) mod_div 5;
1.1013 -(2,2) mod_div 5;
1.1014 -mod_mult 2 4 5;
1.1015 -a; b;
1.1016 -a list_mod 5;
1.1017 -b list_mod 5;
1.1018 -mod_minus 0 1 5;
1.1019 -replicate (length(a))
1.1020 -
1.1021 -*}
1.1022 -
1.1023 -ML{*
1.1024 -
1.1025 -(* modInvers(r,m) --> rinv= r^-1 mod m *)
1.1026 -(* CRA_2(r1,r2,m1,m2) --> r= r1 mod m1 = r2 mod m2 *)
1.1027 -(* lc uv --> leading coefficient of univariate polynomials, deg uv *)
1.1028 -(* abs x --> |x| *)
1.1029 -(* INT_GCD a b *)
1.1030 -(* m listMod list --> new list where every entry mod m*)
1.1031 -(*primeList x --> a list with primenumbers until >= x *)
1.1032 -(*findNextPrime p *)
1.1033 -(*findNextPrimeNotDievide x p --> find the primenumber after p wich not divide x*)
1.1034 -(*list_mod list m, list_mult list b*)
1.1035 -(*(a,b) mod_div m, mod_minus a b m, mod_mult a b m, mod_plus a b m*)
1.1036 -*}
1.1037 -ML{*
1.1038 -fun lc_of_list_not_0 [] = []
1.1039 -| lc_of_list_not_0 list =
1.1040 -if nth list (length list -1) =0
1.1041 -then lc_of_list_not_0 (nth_drop (length list-1) list)
1.1042 -else list;
1.1043 -
1.1044 -
1.1045 -*}
1.1046 -ML {*
1.1047 -
1.1048 -fun mod_poly_div (a,b,m,d)=
1.1049 -let
1.1050 -val moda= a list_mod m;
1.1051 -val modb= ( replicate ( length ( a ) - length ( b )) 0 ) @ (b list_mod m) ;
1.1052 -val c = ( lc ( moda ) , lc ( modb )) mod_div m;
1.1053 -val rest = lc_of_list_not_0 ( list_minus moda (( list_mult modb c ))list_mod m );
1.1054 -in
1.1055 -if rest= []
1.1056 -then [b, [0],[c]@d]
1.1057 -else
1.1058 - if length ( rest ) < length ( b )
1.1059 - then [ rest,[c]@d]
1.1060 - else ( mod_poly_div ( rest , b , m ,[c]@ d ))
1.1061 -end
1.1062 -
1.1063 -infix mod_polynom_div
1.1064 -fun ( a , b ) mod_polynom_div m =
1.1065 -let
1.1066 -val moda= a list_mod m;
1.1067 -val modb= ( replicate ( length ( a ) - length ( b )) 0 ) @ (b list_mod m) ;
1.1068 -val c = ( lc ( moda ) , lc ( modb )) mod_div m;
1.1069 -val rest =lc_of_list_not_0 ( list_minus moda ( list_mult modb c ) list_mod m )list_mod m;
1.1070 -in (* nth_drop ( length ( a ) - 1 )*)
1.1071 -if rest= []
1.1072 -then [b,[0],[c]]
1.1073 -else
1.1074 - if length ( rest ) < length ( b )
1.1075 - then [rest,[c]]
1.1076 - else mod_poly_div (rest, b , m , [c])
1.1077 -end
1.1078 -
1.1079 -*}
1.1080 -
1.1081 -ML {*
1.1082 -val m=5; val a= [2,1]; val b=[3,4];
1.1083 -
1.1084 -([2,1],[3,4]) mod_polynom_div 5;
1.1085 -
1.1086 -val er=([2,0,2],[3,4]) mod_polynom_div 5;
1.1087 -*}
1.1088 -
1.1089 -ML {*
1.1090 -(*mod_polynom_div 1000000000000000;*)
1.1091 -val b=[~18,~15,~20,12,20,~13,2];
1.1092 -val a=[8,28,22,~11,~14,1,2];
1.1093 -*}
1.1094 -
1.1095 -
1.1096 -ML{*
1.1097 -(* modInvers(r,m) --> rinv= r^-1 mod m *)
1.1098 -(* CRA_2(r1,r2,m1,m2) --> r= r1 mod m1 = r2 mod m2 *)
1.1099 -(* lc uv --> leading coefficient of univariate polynomials, deg uv *)
1.1100 -(* abs x --> |x| *)
1.1101 -(* INT_GCD a b *)
1.1102 -(* m listMod list --> new list where every entry mod m*)
1.1103 -(*primeList x --> a list with primenumbers until >= x *)
1.1104 -(*findNextPrime p *)
1.1105 -(*findNextPrimeNotDievide x p --> find the primenumber after p wich not divide x*)
1.1106 -(*list_mod list m, list_mult list b*)
1.1107 -(*mod_div a b m, mod_minus a b m, mod_mult a b m, mod_plus a b m*)
1.1108 -(*(a,b) mod_polynom_div m --> (a/b)mod m = [rest, quotient *)
1.1109 -*}
1.1110 -
1.1111 -ML{*
1.1112 -fun normirt a b m lc_a i =
1.1113 -if i=0
1.1114 -then [(((nth a i) ,lc_a)mod_div m)]@b
1.1115 -else normirt a ( [(((nth a i) ,lc_a)mod_div m)]@b) m lc_a (i-1) ;
1.1116 -
1.1117 -fun normirt_polynom a m =
1.1118 -let
1.1119 -val lc_a=lc a;
1.1120 -in
1.1121 -if length(a)=0
1.1122 -then a
1.1123 -else normirt a [] m lc_a (length(a)-1)
1.1124 -end
1.1125 -*}
1.1126 -
1.1127 -ML{*
1.1128 -lc [2,3,4,2];
1.1129 -val list_a=[3,3,3,3]; val lc_a =lc list_a;val m = 5;
1.1130 -normirt_polynom list_a 5;
1.1131 -*}
1.1132 -ML{*
1.1133 -(* modInvers(r,m) --> rinv= r^-1 mod m *)
1.1134 -(* CRA_2(r1,r2,m1,m2) --> r= r1 mod m1 = r2 mod m2 *)
1.1135 -(* lc uv --> leading coefficient of univariate polynomials, deg uv *)
1.1136 -(* abs x --> |x| *)
1.1137 -(* INT_GCD a b *)
1.1138 -(* m listMod list --> new list where every entry mod m*)
1.1139 -(*primeList x --> a list with primenumbers until >= x *)
1.1140 -(*findNextPrime p *)
1.1141 -(*findNextPrimeNotDievide x p --> find the primenumber after p wich not divide x*)
1.1142 -(*list_mod list m, list_mult list b*)
1.1143 -(*mod_div a b m, mod_minus a b m, mod_mult a b m, mod_plus a b m*)
1.1144 -(*(a,b) mod_polynom_div m --> (a/b)mod m = [rest, quotient *)
1.1145 -(*normirt_polynom a m*)
1.1146 -*}
1.1147 -
1.1148 -ML{*
1.1149 -infix modGCD_p;
1.1150 -fun (a , b) modGCD_p p =
1.1151 -let
1.1152 -val rest= ((a,b) mod_polynom_div p);
1.1153 -in
1.1154 -if length(rest)=3
1.1155 -then normirt_polynom (nth rest 0) m
1.1156 -else ( b ,nth rest 0 ) modGCD_p p
1.1157 -end
1.1158 -*}
1.1159 -
1.1160 -ML{*
1.1161 -(a,b) modGCD_p 5;
1.1162 -(a,b) modGCD_p 7;
1.1163 -(a,b) modGCD_p 11;
1.1164 -(a,b) modGCD_p 13;
1.1165 -(a,b) modGCD_p 17;
1.1166 -*}
1.1167 -
1.1168 -ML{*
1.1169 -(* modInvers(r,m) --> rinv= r^-1 mod m *)
1.1170 -(* CRA_2(r1,r2,m1,m2) --> r= r1 mod m1 = r2 mod m2 *)
1.1171 -(* lc uv --> leading coefficient of univariate polynomials, deg uv *)
1.1172 -(* abs x --> |x| *)
1.1173 -(* INT_GCD a b *)
1.1174 -(* m listMod list --> new list where every entry mod m*)
1.1175 -(*primeList x --> a list with primenumbers until >= x *)
1.1176 -(*findNextPrime p *)
1.1177 -(*findNextPrimeNotDievide x p --> find the primenumber after p wich not divide x*)
1.1178 -(*list_mod list m, list_mult list b*)
1.1179 -(*mod_div a b m, mod_minus a b m, mod_mult a b m, mod_plus a b m*)
1.1180 -(*(a,b) mod_polynom_div m --> (a/b)mod m = [rest, quotient *)
1.1181 -(*normirt_polynom a m*)
1.1182 -(* (a,b)modGCD_p p ggt(a_p,b_p)*)
1.1183 -*}
1.1184 -
1.1185 -
1.1186 -ML {*
1.1187 -5 mod 3;
1.1188 -*}
1.1189 -
1.1190 -
1.1191 -ML {*
1.1192 -fun gcd_mod' a b g p' d m =
1.1193 -let
1.1194 -val p=findNextPrimeNotDivide d p';
1.1195 -in
1.1196 -if p' <= m
1.1197 -then
1.1198 - ( if deg (list_mult ((a,b)modGCD_p p)(d mod p))< deg g
1.1199 - then gcd_mod' a b g p d m
1.1200 - else
1.1201 - ( if deg (list_mult ((a,b)modGCD_p p)(d mod p))= deg g
1.1202 - then gcd_mod' a b (CRA_2_poly g (list_mult ((a,b)modGCD_p p)(d mod p)) ((findNextPrime p')*p') (findNextPrime p') ) ((findNextPrime p')*p') d m
1.1203 - else list_mult ((a,b)modGCD_p p)(d mod p)
1.1204 - )
1.1205 - )
1.1206 -else g
1.1207 -end
1.1208 -*}
1.1209 -ML {* (*----------------------------------------------------------------*)
1.1210 -val g = [2,2,2,2];
1.1211 -val g_p = [3,2,3,2];
1.1212 -val P = 5;
1.1213 -val p = 7;
1.1214 -*}
1.1215 -ML {*
1.1216 -val (r1, r2, m1, m2) = (2, 3, 5, 7);
1.1217 -fun CRA_2' (m1, m2) (r1, r2) = CRA_2 (r1, r2, m1, m2)
1.1218 -*}
1.1219 -ML {*
1.1220 -map (CRA_2' (m1, m2)) (g ~~ g_p)
1.1221 -*}
1.1222 -ML {*
1.1223 -*}
1.1224 -ML {*
1.1225 -*}
1.1226 -ML {*
1.1227 -*}
1.1228 -ML {*
1.1229 - (*\cite{winkler96} p.94 ----------------------------------------------*)
1.1230 -val a = [~18, ~15, ~20, 12, 20, ~13, 2];
1.1231 -val b = [8, 28, 22, ~11, ~14, 1, 2];
1.1232 -*}
1.1233 -ML {*
1.1234 -fun aaaaa (a, b) d g P =
1.1235 - let
1.1236 - val p = findNextPrimeNotDivide d P;
1.1237 - val c_p = (a , b) modGCD_p p;
1.1238 - val g_p = list_mult c_p (d mod p);
1.1239 -(*---
1.1240 - val g =
1.1241 - if deg g_p = deg g then (*goto (3)*) [2]
1.1242 - else
1.1243 - let
1.1244 - val P = p
1.1245 - val g = g_p;
1.1246 - in g end;
1.1247 ----*)
1.1248 -
1.1249 -
1.1250 - in 123 end;
1.1251 -
1.1252 -*}
1.1253 -text {*
1.1254 -
1.1255 -(*start of fun GCD_MOD*)
1.1256 -lc a; lc b;
1.1257 -(*1*)
1.1258 -val d = INT_GCD (lc a) (lc b)
1.1259 -val M = 10412;
1.1260 -(*2*)
1.1261 -val p = findNextPrimeNotDivide d 3;
1.1262 -val c_p = (a , b) modGCD_p p;
1.1263 -val g_p = list_mult c_p (d mod p);
1.1264 -(*3*)
1.1265 -deg g_p;
1.1266 -fun bbbbb p 0 = [1]
1.1267 - | bbbbb p g_p =
1.1268 - let
1.1269 - val P = p
1.1270 - val g = g_p
1.1271 - in aaaaa (a, b) d g P
1.1272 -
1.1273 -
1.1274 -if deg g_p = 0 then [1]
1.1275 -else
1.1276 - let
1.1277 - val P = p
1.1278 - val g = g_p;
1.1279 - in g end;
1.1280 -*}
1.1281 -ML {*
1.1282 -(*4*)
1.1283 -*}
1.1284 -ML {*
1.1285 -
1.1286 -*}
1.1287 -ML {*
1.1288 -*}
1.1289 -ML {*
1.1290 -*}
1.1291 -ML {*
1.1292 -
1.1293 -fun gcd_mod a b p' m d = (*Fuer die WHILE-Schleife*)
1.1294 -let
1.1295 -val p=findNextPrimeNotDivide d p';
1.1296 -val g_p = list_mult ((a,b)modGCD_p p) (d mod p)
1.1297 -in
1.1298 -if deg g_p =0
1.1299 -then 1
1.1300 -else g=gcd_mod' a b g_p p m d;
1.1301 -if length((((a,b)modGCD_p p),a) mod_polynom_div m)=3 && length((((a,b)modGCD_p p),b) mod_polynom_div m)=3
1.1302 -then g
1.1303 -else gcd_mod a b p m d
1.1304 -end
1.1305 -
1.1306 -fun GCD_MOD a b =
1.1307 -let
1.1308 -val d = INT_GCD (lc a) (lc b);
1.1309 -val m = 2*d*landau_mignotte_bound a b;
1.1310 -val p'=1;
1.1311 -val p = findNextPrime 1;
1.1312 -(*val c_p = (a,b)modGCD_p p;*)
1.1313 -val g_p =list_mult ((a,b)modGCD_p p)(d mod p);
1.1314 -in
1.1315 -if deg g_p=0
1.1316 -then 1
1.1317 -else (*while-Schleife*)
1.1318 -gcd_mod a b g_p p m d
1.1319 -end
1.1320 -
1.1321 -*}
1.1322 -
1.1323 -ML{*
1.1324 -lc [1,2,3,4,5];
1.1325 -
1.1326 -*}
1.1327 -
1.1328 -
1.1329 -end
2.1 --- a/test/Tools/isac/Knowledge/rational2.sml Thu Sep 20 10:07:02 2012 +0200
2.2 +++ b/test/Tools/isac/Knowledge/rational2.sml Mon Sep 24 09:07:38 2012 +0200
2.3 @@ -4,17 +4,453 @@
2.4 12345678901234567890123456789012345678901234567890123456789012345678901234567890
2.5 10 20 30 40 50 60 70 80
2.6 *)
2.7 +"--------------------------------------------------------";
2.8 +"table of contents --------------------------------------";
2.9 +"--------------------------------------------------------";
2.10 +"----------- special case for partial fractions ---------";
2.11 +"----------- fun mod_invers -----------------------------";
2.12 +"----------- fun CRA_2 ----------------------------------";
2.13 +"----------- fun lc -------------------------------------";
2.14 +"----------- fun deg ------------------------------------";
2.15 +"----------- fun INT_GCD --------------------------------";
2.16 +"----------- fun primeList ------------------------------";
2.17 +"----------- fun find_next_prime_not_divide -------------";
2.18 +"----------- fun poly_mod -------------------------------";
2.19 +"----------- fun %* -------------------------------------";
2.20 +"----------- fun %/ -------------------------------------";
2.21 +"----------- fun %-% ------------------------------------";
2.22 +"----------- fun %+ -------------------------------------";
2.23 +"----------- fun CRA_2_poly -----------------------------";
2.24 +"----------- fun mod_div --------------------------------";
2.25 +"----------- fun poly_mod_div ---------------------------";
2.26 +"----------- fun mod_mult -------------------------------";
2.27 +"----------- fun mod_minus ------------------------------";
2.28 +"----------- fun mod_plus -------------------------------";
2.29 +"----------- fun lc_of_list_not_0 -----------------------";
2.30 +"----------- fun %/% ------------------------------------";
2.31 +"----------- fun mod_polynom_gcd ------------------------";
2.32 +"----------- fun normirt_polynom ------------------------";
2.33 +"----------- fun poly_centr -----------------------------";
2.34 +"----------- fun mod_gcd_p ------------------------------";
2.35 +"----------- fun gcd_n ----------------------------------";
2.36 +"----------- fun pp -------------------------------------";
2.37 +"----------- fun GCD_MOD --------------------------------";
2.38 +"--------------------------------------------------------";
2.39 +"--------------------------------------------------------";
2.40
2.41 -writeln "test/../rational2";
2.42 +"----------- values used for many tests -----------------";
2.43 + val a = [~18, ~15, ~20, 12, 20, ~13, 2]: unipoly;
2.44 + val b = [8, 28, 22, ~11, ~14, 1, 2]: unipoly;
2.45 +"--------------------------------------------------------";
2.46 +(*============ inhibit exn WN120920 ==============================================
2.47 +a %- 10;
2.48
2.49 +val list = mod_gcd_p a b 7;
2.50 +(list %* 2) poly_mod 7;
2.51 +val list = mod_gcd_p a b 11;
2.52 +(list %* 2) poly_mod 11;
2.53 +val list = mod_gcd_p a b 13;
2.54 +(list %* 2) poly_mod 13;
2.55 +val list = mod_gcd_p a b 17;
2.56 +(list %* 2) poly_mod 17;
2.57 +
2.58 +a;b;
2.59 +find_next_prime_not_divide 2 7;
2.60 +
2.61 +a;
2.62 +b;
2.63 +val gcd = GCD_MOD a b;
2.64 +a %/% gcd;
2.65 +b %/% gcd;
2.66 +
2.67 +a;b;
2.68 +val gcd = GCD_MOD' a b;
2.69 +a %/% gcd;
2.70 +b %/% gcd;
2.71 +
2.72 +(*poly :: 'a poly => 'a => 'a;*)
2.73 +(*[:4,0,1:];*)
2.74 +(*[:0, 1, 1:] =!= 0*)
2.75 +
2.76 +pp [4,8,12,~2,0,~4];
2.77 +cont a;
2.78 +cont [4,8,12,~2,0,~4];
2.79 +cont' mpoly1;
2.80 +gcd_n a;
2.81 +nth a 3;
2.82 +nth_drop 3 a;
2.83 +
2.84 +============ inhibit exn WN120920 ==============================================*)
2.85 +
2.86 +
2.87 +
2.88 +"----------- special case for partial fractions ---------";
2.89 +"----------- special case for partial fractions ---------";
2.90 +"----------- special case for partial fractions ---------";
2.91 (*specific simplification for partial fractions in Z-transform*)
2.92 -val ctxt = ProofContext.init_global @{theory Isac};
2.93 -val ctxt = declare_constraints' [@{term "z::real"}] ctxt;
2.94 +val ctxt = ProofContext.init_global @{theory Rational2};
2.95 +val ctxt = Variable.declare_constraints @{term "z::real"} ctxt;
2.96 +(* parseNEW is defined in Isac; wait for that
2.97 val SOME fract1 =
2.98 parseNEW ctxt "(z - 1 / 2) * (z - -1 / 4) * (A / (z - 1 / 2) + B / (z - -1 / 4))";
2.99 -val SOME (fract2,_) = rewrite_set_ @{theory Isac} false norm_Rational fract1;
2.100 +val SOME (fract2,_) = rewrite_set_ @{theory Rational2} false norm_Rational fract1;
2.101 +*)
2.102 (* TODO
2.103 if term2str fract2 = "A * (1 / 4 + z) + B * (-1 / 2 + z)" then ()
2.104 else error "specific simplification for partial fractions";
2.105 *)
2.106
2.107 +"----------- fun mod_invers -----------------------------";
2.108 +"----------- fun mod_invers -----------------------------";
2.109 +"----------- fun mod_invers -----------------------------";
2.110 +"~~~~~ fun mod_inv, args:"; val (r, m) = (5,7);
2.111 +"~~~~~ fun modi, args:"; val (r, rold, m, rinv) = (r, r, m, 1);
2.112 +rinv < m; (*=true*)
2.113 +r mod m = 1; (*=false*)
2.114 +"~~~~~ fun modi, args:"; val (r, rold, m, rinv) = (rold * (rinv + 1), rold, m, rinv + 1);
2.115 +rinv < m;(*=true*)
2.116 +r mod m = 1; (*=false*)
2.117 +"~~~~~ fun modi, args:"; val (r, rold, m, rinv) = (rold * (rinv + 1), rold, m, rinv + 1);
2.118 +rinv < m;(*=true*)
2.119 +r mod m = 1; (*=true*)
2.120 +"~~~~~ to mod_inv return val:"; val (rinv) = (rinv);
2.121 +
2.122 +if mod_inv 5 7 = 3 then () else error "mod_inv 5 7 = 3 changed";
2.123 +if mod_inv 3 7 = 5 then () else error "mod_inv 3 7 = 5 changed";
2.124 +if mod_inv 4 339 = 85 then () else error "mod_inv 4 339 = 85 changed";
2.125 +
2.126 +"----------- fun CRA_2 ----------------------------------";
2.127 +"----------- fun CRA_2 ----------------------------------";
2.128 +"----------- fun CRA_2 ----------------------------------";
2.129 +
2.130 +if CRA_2(17,9,3,4) = 35 then () else error "CRA_2(17,9,3,4)=35 changed";
2.131 +if CRA_2(7,2,6,11) = 17 then () else error "CRA_2(7,2,6,11) = 17 changed";
2.132 +
2.133 +"----------- fun lc -------------------------------------";
2.134 +"----------- fun lc -------------------------------------";
2.135 +"----------- fun lc -------------------------------------";
2.136 +if lc [3,4,5,6] = 6 then () else error "ls_mv (3,4,5,6) = 6 changed" ;
2.137 +if lc [3,4,5,6,0] = 6 then () else error "ls_mv (3,4,5,6,0) = 6 changed" ;
2.138 +
2.139 +"----------- fun deg ------------------------------------";
2.140 +"----------- fun deg ------------------------------------";
2.141 +"----------- fun deg ------------------------------------";
2.142 +if deg [3,4,5,6] = 3 then () else error "lc (3,4,5,6) = 3 changed" ;
2.143 +if deg [3,4,5,6,0] = 3 then () else error "lc (3,4,5,6) = 6 changed";
2.144 +
2.145 +"----------- fun INT_GCD --------------------------------";
2.146 +"----------- fun INT_GCD --------------------------------";
2.147 +"----------- fun INT_GCD --------------------------------";
2.148 +if INT_GCD 15 6 = 3 then () else error "INT_GCD 15 6 = 3 changed";
2.149 +if INT_GCD 10 3 =1 then () else error "INT_GCD 10 3 = 1 changed";
2.150 +if INT_GCD 6 24 = 6 then () else error "INT_GCD 6 24 = 6 changed";
2.151 +
2.152 +"----------- fun primeList ------------------------------";
2.153 +"----------- fun primeList ------------------------------";
2.154 +"----------- fun primeList ------------------------------";
2.155 +if primeList 1 = [2] then () else error " primeList 1 = [2] changed";
2.156 +if primeList 3 = [2,3] then () else error " primeList 3 = [2,3] changed";
2.157 +if primeList 6 = [2,3,5,7] then () else error " primeList 6 = [2,3,5,7] changed";
2.158 +if primeList 7 = [2,3,5,7] then () else error " primeList 7 = [2,3,5,7] changed";
2.159 +
2.160 +"----------- fun find_next_prime_not_divide -----------------";
2.161 +"----------- fun find_next_prime_not_divide -----------------";
2.162 +"----------- fun find_next_prime_not_divide -----------------";
2.163 +if find_next_prime_not_divide 15 2 = 7 then () else error "findNextPrimeNotdivide 15 2 = 7 changed";
2.164 +if find_next_prime_not_divide 90 1 = 7 then () else error "findNextPrimeNotdivide 90 1 = 7 changed";
2.165 +
2.166 +"----------- fun poly_mod -------------------------------";
2.167 +"----------- fun poly_mod -------------------------------";
2.168 +"----------- fun poly_mod -------------------------------";
2.169 +if ([5,4,7,8,1] poly_mod 5) = [0, 4, 2, 3, 1] then ()
2.170 +else error "[5,4,7,8,1] poly_mod 5 = [0, 4, 2, 3, 1] changed" ;
2.171 +if ([5,4,~7,8,~1] poly_mod 5) = [0, 4, 3, 3, 4] then ()
2.172 +else error "[5,4,~7,8,~1] poly_mod 5 = [0, 4, 3, 3, 4] changed" ;
2.173 +
2.174 +"----------- fun %* ------------------------------";
2.175 +"----------- fun %* ------------------------------";
2.176 +"----------- fun %* ------------------------------";
2.177 +if ([5,4,7,8,1] %* 5) = [25, 20, 35, 40, 5] then ()
2.178 +else error "[5,4,7,8,1] %* 5 = [25, 20, 35, 40, 5] changed";
2.179 +if ([5,4,~7,8,~1] %* 5) = [25, 20, ~35, 40, ~5] then ()
2.180 +else error "[5,4,~7,8,~1] %* 5 = [25, 20, ~35, 40, ~5] changed";
2.181 +
2.182 +"----------- fun %/ -------------------------------";
2.183 +"----------- fun %/ -------------------------------";
2.184 +"----------- fun %/ -------------------------------";
2.185 +if ([4,3,2,5,6] %/ 3) =[1, 1, 0, 1, 2] then ()
2.186 +else error "%/ [4,3,2,5,6] 3 =[1, 1, 0, 1, 2] changed";
2.187 +if ([4,3,2,0] %/ 3) =[1, 1, 0, 0] then ()
2.188 +else error "%/ [4,3,2,0] 3 =[1, 1, 0, 0] changed";
2.189 +
2.190 +"----------- fun %-% ------------------------";
2.191 +"----------- fun %-% ------------------------";
2.192 +"----------- fun %-% ------------------------";
2.193 +if ([8,~7,0,1] %-% [~2,2,3,0])=[10,~9,~3,1] then ()
2.194 +else error "[8,~7,0,1] %-% [~2,2,3,0]=[10,~9,~3,1] changed";
2.195 +if ([8,7,6,5,4] %-% [2,2,3,1,1])=[6,5,3,4,3] then ()
2.196 +else error "[8,7,6,5,4] %-% [2,2,3,1,1]=[6,5,3,4,3] changed";
2.197 +
2.198 +"----------- fun %- -----------------------------";
2.199 +"----------- fun %- -----------------------------";
2.200 +"----------- fun %- -----------------------------";
2.201 +if ([8,~7,0,1] %- 8 )=[0,~15,~8,~7] then ()
2.202 +else error "[8,~7,0,1] %-% 8 =[0,~15,~8,~7] changed";
2.203 +if ([8,7,6,5,4] %- ~2 )=[10,9,8,7,6] then ()
2.204 +else error "[8,7,6,5,4] %-% ~2 =[10,9,8,7,6] changed";
2.205 +
2.206 +"----------- fun %+ ------------------------------";
2.207 +"----------- fun %+ ------------------------------";
2.208 +"----------- fun %+ ------------------------------";
2.209 +if ([8,~7,0,1] %+ [~2,2,3,0])=[6,~5,3,1] then ()
2.210 +else error "[8,~7,0,1] %+ [~2,2,3,0]=[6,~5,3,1] changed";
2.211 +if ([8,7,6,5,4] %+ [2,2,3,1,1])=[10,9,9,6,5] then ()
2.212 +else error "[8,7,6,5,4] %+ [2,2,3,1,1]=[10,9,9,6,5] changed";
2.213 +
2.214 +"----------- fun CRA_2_poly -----------------------------";
2.215 +"----------- fun CRA_2_poly -----------------------------";
2.216 +"----------- fun CRA_2_poly -----------------------------";
2.217 +if (CRA_2_poly (5, 7)([2,2,4,3],[3,2,3,5]))=[17,2,24,33] then ()
2.218 +else error "CRA_2_poly (5, 7)([2,2,4,3],[3,2,3,5])=[17,2,24,33] changed";
2.219 +
2.220 +"----------- fun mod_div --------------------------------";
2.221 +"----------- fun mod_div --------------------------------";
2.222 +"----------- fun mod_div --------------------------------";
2.223 +if mod_div 21 4 5 = 4 then () else error "mod_div 21 4 5 = 4 changed";
2.224 +if mod_div 1 4 5 = 4 then () else error "mod_div 1 4 5 = 4 changed";
2.225 +if mod_div 0 4 5 = 0 then () else error "mod_div 0 4 5 = 0 changed";
2.226 +
2.227 +"----------- fun poly_mod_div ---------------------------";
2.228 +"----------- fun poly_mod_div ---------------------------";
2.229 +"----------- fun poly_mod_div ---------------------------";
2.230 +if poly_mod_div [21, 1, 0] 4 5 = [4, 4, 0] then ()
2.231 + else error "poly_mod_div [21, 1, 0]=[4,4,0] changed";
2.232 +
2.233 +"----------- fun mod_mult -------------------------------";
2.234 +"----------- fun mod_mult -------------------------------";
2.235 +"----------- fun mod_mult -------------------------------";
2.236 +if mod_mult 4 4 5 = 1 then () else error "mod_mult 4 4 5 = 1 changed";
2.237 +
2.238 +"----------- fun mod_minus ------------------------------";
2.239 +"----------- fun mod_minus ------------------------------";
2.240 +"----------- fun mod_minus ------------------------------";
2.241 +if mod_minus 9 3 5 = 1 then () else error "mod_minus 9 3 5 = 1 changed";
2.242 +
2.243 +"----------- fun mod_plus -------------------------------";
2.244 +"----------- fun mod_plus -------------------------------";
2.245 +"----------- fun mod_plus -------------------------------";
2.246 +if mod_plus 9 3 5 = 2 then () else error "mod_plus 9 3 5 =2 changed";
2.247 +
2.248 +"----------- fun lc_of_list_not_0 -----------------------";
2.249 +"----------- fun lc_of_list_not_0 -----------------------";
2.250 +"----------- fun lc_of_list_not_0 -----------------------";
2.251 +if lc_of_list_not_0 [0,1,2,3,4,5,0,0]=[1, 2, 3, 4, 5] then ()
2.252 +else error "lc_of_list_not_0 [0,1,2,3,4,5,0,0]=[1, 2, 3, 4, 5] changed";
2.253 +if lc_of_list_not_0 [0,1,2,3,4,5]=[1, 2, 3, 4, 5] then ()
2.254 +else error "lc_of_list_not_0 [0,1,2,3,4,5]=[ 1, 2, 3, 4, 5] changed";
2.255 +
2.256 +"----------- fun %/% -------------------------------";
2.257 +"----------- fun %/% -------------------------------";
2.258 +"----------- fun %/% -------------------------------";
2.259 +"~~~~~ fun %/% , args:"; val (poly1, poly2) = ([9,12,4],[3,2]);
2.260 +val b = (replicate (length poly1 - length poly2) 0 ) @ poly2 ;
2.261 +val c = lc poly1 div lc b;
2.262 +val rest = lc_of_list_not_0 ( poly1 %-% (b %* c ));
2.263 +rest = [];(*=false*)
2.264 +c<>0 andalso length rest >= length poly2; (*=true*)
2.265 +"~~~~~ fun %/% , args:"; val (poly1, poly2) = (rest, poly2);
2.266 +val b = (replicate (length poly1 - length poly2) 0 ) @ poly2 ;
2.267 +val c = lc poly1 div lc b;
2.268 +val rest = lc_of_list_not_0 ( poly1 %-% (b %* c ));
2.269 +rest = [];(*=true*)
2.270 +"~~~~~ to return val:"; val (divide) = (true);
2.271 +
2.272 +if ([6] %/% [4]) = false then () else error "[6] %/% [4] = false changed";
2.273 +if ([16,0] %/% [8]) = true then () else error "[16,0] %/% [8] = true changed";
2.274 +if ([0,0,9,12,4] %/% [3,2]) = true then ()
2.275 +else error "[0,0,9,12,4] %/% [3,2] = true changed";
2.276 +if ([16] %/% [8,0]) = true then () else error "[16] %/% [8,0] = true changed";
2.277 +
2.278 +"----------- fun mod_polynom_gcd ------------------------";
2.279 +"----------- fun mod_polynom_gcd ------------------------";
2.280 +"----------- fun mod_polynom_gcd ------------------------";
2.281 +"~~~~~ fun mod_polynom_gcd , args:";
2.282 +val (a,b,m) = ([ ~13, 2,12],[ ~14, 2],5);
2.283 +val moda = a poly_mod m;
2.284 +val modb = (replicate (length a - length(lc_of_list_not_0 b)) 0) @ (lc_of_list_not_0 b poly_mod m) ;
2.285 +val c = mod_div (lc moda) (lc modb) m;
2.286 +val rest = lc_of_list_not_0 (moda %-% (modb %* c) poly_mod m) poly_mod m;
2.287 +rest = []; (*=false*)
2.288 +length rest < length b; (*=false*)
2.289 +"~~~~~ fun mod_poly_gcd , args:";
2.290 +val (a,b,m,d) = (rest, b , m , [c]);
2.291 +val moda = a poly_mod m
2.292 +val modb = (replicate (length a - length(lc_of_list_not_0 b)) 0) @ (lc_of_list_not_0 b poly_mod m) ;
2.293 +val c = mod_div (lc moda) (lc modb) m
2.294 +val rest = lc_of_list_not_0 ((moda %-% (modb %* c)) poly_mod m);
2.295 +rest = [];(*=flase*)
2.296 +length rest < length b; (*=true*)
2.297 +"~~~~~ fun mod_poly_gcd , args:";
2.298 +val (a,b,m,d) = (b, rest, m, [c] @ d);
2.299 +val moda = a poly_mod m
2.300 +val modb = (replicate (length a - length(lc_of_list_not_0 b)) 0) @ (lc_of_list_not_0 b poly_mod m) ;
2.301 +val c = mod_div (lc moda) (lc modb) m
2.302 +val rest = lc_of_list_not_0 ((moda %-% (modb %* c)) poly_mod m);
2.303 +rest = [];(*=flase*)
2.304 +length rest < length b; (*=false*)
2.305 +"~~~~~ fun mod_poly_gcd , args:";
2.306 +val (a,b,m,d) = (b, rest, m, [c] @ d);
2.307 +val moda = a poly_mod m
2.308 +val modb = (replicate (length a - length(lc_of_list_not_0 b)) 0) @ (lc_of_list_not_0 b poly_mod m) ;
2.309 +val c = mod_div (lc moda) (lc modb) m
2.310 +val rest = lc_of_list_not_0 ((moda %-% (modb %* c)) poly_mod m);
2.311 +rest = [];(*=true*)
2.312 +"~~~~~ to return val:"; val ([b, rest, lsit])=[b, [0], [c] @ d];
2.313 +
2.314 +if ( mod_polynom_gcd [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 7)=[[2, 6, 0, 2, 6]]
2.315 +then () else error "( mod_polynom_gcd [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 7 = [[2, 6, 0, 2, 6] changed";
2.316 +if ( mod_polynom_gcd [8, 28, 22, ~11, ~14, 1, 2] [2, 6, 0, 2, 6] 7)=[[1, 3, 0, 1, 3], [0]] then ()
2.317 +else error "mod_polynom_gcd [8, 28, 22, ~11, ~14, 1, 2] [2, 6, 0, 2, 6] 7=[[1, 3, 0, 1, 3], [0]] changed";
2.318 +
2.319 +"----------- fun normirt_polynom ------------------------";
2.320 +"----------- fun normirt_polynom ------------------------";
2.321 +"----------- fun normirt_polynom ------------------------";
2.322 +if (normirt_polynom [~18, ~15, ~20, 12, 20, ~13, 2] 5) =[1, 0, 0, 1, 0, 1, 1 ] then ()
2.323 +else error "(normirt_polynom [~18, ~15, ~20, 12, 20, ~13, 2] 5) =[1, 0, 0, 1, 0, 1, 1 ] changed"
2.324 +
2.325 +"----------- fun poly_centr -----------------------------";
2.326 +"----------- fun poly_centr -----------------------------";
2.327 +"----------- fun poly_centr -----------------------------";
2.328 +if (poly_centr [7,3,5,8,1,3] 10) = [~3, 3, 5, ~2, 1, 3] then ()
2.329 +else error "poly_centr [7,3,5,8,1,3] 10 = [~3, 3, 5, ~2, 1, 3] cahnged";
2.330 +
2.331 +"----------- fun mod_gcd_p -------------------------------";
2.332 +"----------- fun mod_gcd_p -------------------------------";
2.333 +"----------- fun mod_gcd_p -------------------------------";
2.334 +if (mod_gcd_p [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 5) = [1, 1, 1, 1] then ()
2.335 +else error "mod_gcd_p [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 5 = [1, 1, 1, 1] changed";
2.336 +if (mod_gcd_p [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 7) = [5, 1, 0, 5, 1] then ()
2.337 +else error "mod_gcd_p [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 7 = [5, 1, 0, 5, 1] changed";
2.338 +
2.339 +"----------- fun gcd_n ----------------------------";
2.340 +"----------- fun gcd_n ----------------------------";
2.341 +"----------- fun gcd_n ----------------------------";
2.342 +if gcd_n [3,9,12,21] = 3 then () else error "gcd_n [3,9,12,21] = 3 changed";
2.343 +if gcd_n [12,16,32,44] = 4 then () else error "gcd_n [12,16,32,44] = 4 changed";
2.344 +if gcd_n [4,5,12] = 1 then () else error "gcd_n [4,5,12] = 1 changed";
2.345 +
2.346 +"----------- fun pp -------------------------------------";
2.347 +"----------- fun pp -------------------------------------";
2.348 +"----------- fun pp -------------------------------------";
2.349 +if pp [12,16,32,44] = [3, 4, 8, 11] then () else error "pp [12,16,32,44] = [3, 4, 8, 11] changed";
2.350 +if pp [4,5,12] = [4, 5, 12] then () else error "pp [4,5,12] = [4, 5, 12] changed";
2.351 +
2.352 +"----------- fun GCD_MOD -------------------------------";
2.353 +"----------- fun GCD_MOD -------------------------------";
2.354 +"----------- fun GCD_MOD -------------------------------";
2.355 +
2.356 +"~~~~~ fun GCD_MOD a b , args:"; val (a, b) = ([~18, ~15, ~20, 12, 20, ~13, 2], [8, 28, 22, ~11, ~14, 1, 2]);
2.357 +val d = INT_GCD (lc a) (lc b);
2.358 +val M =2*d* landau_mignotte_bound a b;
2.359 +(*val g = GOTO2 a b d (*p*)1 M*)
2.360 + "~~~~~ fun GOTO2 a b d M p , args:"; val (a, b, d, p, M) = (a,b,d,3,M);
2.361 + val p = find_next_prime_not_divide d p
2.362 + val c_p = normirt_polynom ( mod_gcd_p a b p) p
2.363 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
2.364 + (*val (p, g) = GOTO3 a b d M p g_p (* << 1.Mal ~~~~~~~~~~*)*)
2.365 + "~~~~~ fun GOTO3 a b d M p g_p , args:"; val (a, b, d, M, p, g_p) = (a,b,d,M,p,g_p);
2.366 + (deg g_p) = 0; (*=false*)
2.367 + val P = p
2.368 + val g = g_p;
2.369 + (*(p, g) = WHILE a b d M P g (* << 1.Mal -----*)*)
2.370 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,p,g);
2.371 + P > M ;(*false*)
2.372 + val p = find_next_prime_not_divide d p
2.373 + val c_p = normirt_polynom ( mod_gcd_p a b p) p
2.374 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
2.375 + deg g_p < deg g; (*=fasle*)
2.376 + deg g_p = deg g; (*false*);
2.377 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
2.378 + P > M ;(*false*)
2.379 + val p = find_next_prime_not_divide d p
2.380 + val c_p = normirt_polynom ( mod_gcd_p a b p) p
2.381 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
2.382 + deg g_p < deg g; (*=false*)
2.383 + deg g_p = deg g; (*=true*)
2.384 + val g = CRA_2_poly (P,p)(g,g_p)
2.385 + val P = P*p;
2.386 + val g = poly_centr(g poly_mod P)P;
2.387 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
2.388 + P > M ;(*false*)
2.389 + val p = find_next_prime_not_divide d p
2.390 + val c_p = normirt_polynom (mod_gcd_p a b p) p
2.391 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
2.392 + deg g_p < deg g; (*=true*)
2.393 + "~~~~~ fun GOTO3 a b d M p g_p , args:"; val (a, b, d, M, p, g_p) = (a,b,d,M,p,g_p);
2.394 + (deg g_p) = 0; (*false*)
2.395 + val P = p
2.396 + val g = g_p;
2.397 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
2.398 + P > M ;(*false*)
2.399 + val p = find_next_prime_not_divide d p
2.400 + val c_p = normirt_polynom (mod_gcd_p a b p) p
2.401 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
2.402 + deg g_p < deg g; (*=fasle*)
2.403 + deg g_p = deg g;(*true*)
2.404 + val g = CRA_2_poly (P,p)(g,g_p)
2.405 + val P = P*p;
2.406 + val g = poly_centr(g poly_mod P)P;
2.407 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
2.408 + P > M ;(*false*)
2.409 + val p = find_next_prime_not_divide d p
2.410 + val c_p = normirt_polynom (mod_gcd_p a b p) p
2.411 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
2.412 + deg g_p < deg g; (*=false*)
2.413 + deg g_p = deg g;(*true*)
2.414 + val g = CRA_2_poly (P,p)(g,g_p)
2.415 + val P = P*p;
2.416 + val g = poly_centr(g poly_mod P)P;
2.417 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
2.418 + P > M ;(*false*)
2.419 + val p = find_next_prime_not_divide d p
2.420 + val c_p = normirt_polynom (mod_gcd_p a b p) p
2.421 + val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;
2.422 + deg g_p < deg g; (*=false*)
2.423 + deg g_p = deg g;(*true*)
2.424 + val g = CRA_2_poly (P,p)(g,g_p)
2.425 + val P = P*p;
2.426 + val g = poly_centr(g poly_mod P)P;
2.427 + "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);
2.428 + P > M ;(*true*)
2.429 + "~~~~~ to return WHILE val:"; val (g,p) = (g,p);
2.430 + "~~~~~ to return GOTO3 val:"; val (g,p) = (g,p);
2.431 + val g = pp g;
2.432 + (a %/% g) andalso (b %/% g);(*=true*)
2.433 +"~~~~~ to return GOTO2 val:"; val (g) = (g);
2.434 +"~~~~~ to return GCD_MOD val:"; val (g) = (g);
2.435 +
2.436 +"----------- fun GCD_MOD --------------------------------";
2.437 +"----------- fun GCD_MOD --------------------------------";
2.438 +"----------- fun GCD_MOD --------------------------------";
2.439 +val a = [~18, ~15, ~20, 12, 20, ~13, 2];
2.440 +val b = [8, 28, 22, ~11, ~14, 1, 2];
2.441 +if GCD_MOD a b = [~2, ~1, 1] then ()
2.442 +else error "GCD_MOD a b = [~2, ~1, 1]";
2.443 +
2.444 +"----------- fun %/% -------------------------------";
2.445 +"----------- fun %/% -------------------------------";
2.446 +"----------- fun %/% -------------------------------";
2.447 +
2.448 +"----------- fun %/% -------------------------------";
2.449 +"----------- fun %/% -------------------------------";
2.450 +"----------- fun %/% -------------------------------";
2.451 +
2.452 +"----------- fun %/% -------------------------------";
2.453 +"----------- fun %/% -------------------------------";
2.454 +"----------- fun %/% -------------------------------";
2.455 +
2.456 +"----------- fun %/% -------------------------------";
2.457 +"----------- fun %/% -------------------------------";
2.458 +"----------- fun %/% -------------------------------";
2.459 +
2.460 +