wneuper/isa: test/Tools/isac/Knowledge/gcd_poly

     1 (* Title: test/../gcd_poly_winkler

     2    Author: Diana Meindl

     3    Copyright (c) Diana Meindl 2011

     4 12345678901234567890123456789012345678901234567890123456789012345678901234567890

     5         10        20        30        40        50        60        70        80

     7 Programcode according to [1] for later lookup for mathematicians.

     8 The tests below remain according to [1],

     9 while "~~/test/Tools/isac/Knowledge/gcd_poly.sml" start exactly from the same state

    10 and in time follows development in "~~/src/Tools/isac/Knowledge/GCD_Poly.thy".

    12 [1] Franz Winkler. Polynomial Algorithms in Computer Algebra.

    13 Springer-Verlag/Wien 1996.

    14 *)

    16 (*fun nth _ []      = error "nth _ []" (*Isabelle2002, still saved the hours of update*)

    17     | nth 1 (x::_) = x

    18     | nth n (_::xs) = nth (n-1) xs;*)

    19   fun nth xs i = List.nth (xs, i);     (*recent Isabelle: TODO update all isac code   *)

    20 "--------------------------------------------------------";

    21 "table of contents --------------------------------------";

    22 "--------------------------------------------------------";

    23 "----------- auxiliary functions univariate -------------";

    24 "=========== code for [1] p.93: univariate ==============";

    25 "----------- auxiliary functions multivariate -----------";

    26 "=========== code for [1] p.95: multivariate ============";

    27 "--------------------------------------------------------";

    28 "--- begin of univariate polynomials --------------------";

    29 "----------- fun mod_inv --------------------------------";

    30 "----------- fun CRA_2 ----------------------------------";

    31 "----------- fun lc -------------------------------------";

    32 "----------- fun deg ------------------------------------";

    33 "----------- fun primeList ------------------------------";

    34 "----------- fun find_next_prime_not_divide -------------";

    35 "----------- fun poly_mod -------------------------------";

    36 "----------- fun %* -------------------------------------";

    37 "----------- fun %/ -------------------------------------";

    38 "----------- fun %-% ------------------------------------";

    39 "----------- fun %+% ------------------------------------";

    40 "----------- fun CRA_2_poly -----------------------------";

    41 "----------- fun mod_div --------------------------------";

    42 "----------- fun lc_of_unipoly_not_0 --------------------";

    43 "----------- fun %|% ------------------------------------";

    44 "----------- fun mod_poly_gcd ---------------------------";

    45 "----------- fun normirt_polynom ------------------------";

    46 "----------- fun poly_centr -----------------------------";

    47 "----------- fun mod_gcd_p ------------------------------";

    48 "----------- fun gcd_n ----------------------------------";

    49 "----------- fun pp -------------------------------------";

    50 "----------- fun GCD_MOD---------------------------------";

    51 "--- begin of multivariate polynomials ------------------";

    52 "----------- fun get_coef --------------------------- ---";

    53 "----------- fun get_varlist ----------------------------";

    54 "----------- fun get_coeflist ---------------------------";

    55 "----------- fun add_var_to_multipoly -------------------";

    56 "----------- fun cero_multipoly -------------------------";

    57 "----------- fun short_mv -------------------------------";

    58 "----------- fun greater_var ----------------------------";

    59 "----------- fun order_multipoly ------------------------";

    60 "----------- fun lc_multipoly ---------------------------";

    61 "----------- fun deg_multipoly --------------------------";

    62 "----------- fun max_deg_var ----------------------------";

    63 "----------- infix %%/ ----------------------------------";

    64 "----------- fun cont -----------------------------------";

    65 "----------- fun cont_multipoly -------------------------";

    66 "----------- fun evaluate_mv ----------------------------";

    67 "----------- fun mutli_to_uni ---------------------------";

    68 "----------- fun find_new_r  ----------------------------";

    69 "----------- fun newton_mv  -----------------------------";

    70 "----------- fun mult_with_new_var  ---------------------";

    71 "----------- fun greater_deg_multipoly  -----------------";

    72 "----------- infix %%+%%  -------------------------------";

    73 "----------- infix %%-%%  -------------------------------";

    74 "----------- infix %%*%%  -------------------------------";

    75 "----------- fun newton_first  --------------------------";

    76 "----------- fun newton_mv  -----------------------------";

    77 "----------- fun listgreater  ---------------------------";

    78 "----------- finfix  %%|%%  -----------------------------";

    79 "----------- fun GCD_MODm  ------------------------------";

    81 (*text {*

    82   Below we reference

    83   F. Winkler, Polynomial Algorithms. ...

    84   by page 93 and 95. The identifiers used in this book are re-used in this thesis

    85   in order to support reference (although some of these identifiers

    86   doe not conform with the Isabelle coding standards)

    87 *}*)

    89  print_depth 3; (*20*)

    90  type unipoly = int list;

    92 "----------- auxiliary functions univariate -------------";

    93 "----------- auxiliary functions univariate -------------";

    94 "----------- auxiliary functions univariate -------------";

    96 (*subsection {* calculations for integers *}*)

    98   infix div'

    99   fun a div' b =

   100     if a < 0 then abs a div b * ~1

   101     else a div b;

   103 (*subsection {* modulo calculations for integers *}*)

   104   fun mod_inv r m =

   105     let

   106       fun modi (r, rold, m, rinv) =

   107         if rinv < m

   108         then

   109           if r mod m = 1

   110           then rinv

   111           else modi (rold * (rinv + 1), rold, m, rinv + 1)

   112         else 0

   113      in modi (r, r, m, 1) end;

   115   fun mod_div a b m =

   116     a * (mod_inv b m) mod m;

   118   fun aprox_root a =

   119     let fun root' a n =

   120       if n*n < a then root' a (n+1) else n

   121     in root' a 1

   122   end

   124   (*  r = r1 mod m1 and r= r2 mod m2 *)

   125   fun CRA_2 (r1: int, m1: int, r2: int, m2: int ) =

   126    (r1 mod m1) + ((r2 - (r1 mod m1)) * (mod_inv m1 m2) mod m2) * m1

   128 (*subsection {* prime opertations *}*)

   130   fun is_prime primelist number =

   131     if length primelist >0

   132     then

   133       if (number mod (nth primelist 0))=0

   134       then false

   135       else is_prime (nth_drop 0 primelist) number

   136     else true

   138   fun primeList number =

   139     let

   140     fun make_primelist list last number =

   141       if (nth list (length list -1)) < number

   142       then

   143         if ( is_prime list (last + 2))

   144         then make_primelist (list @ [(last + 2)]) (last + 2) number

   145         else make_primelist  list (last + 2) number

   146       else list

   147     in

   148       if number < 3

   149       then [2]

   150       else make_primelist [2,3] 3 number end

   152   (*find a prime greater p not dividing the number a*)

   153   fun find_next_prime_not_divide a p  =

   154     let

   155       val next = nth (primeList (p + 1)) (length (primeList (p + 1)) - 1) ;

   156     in

   157       if a mod next  <> 0

   158         then next

   159       else find_next_prime_not_divide a next

   160   end

   162 (*subsection {* calculations for univariate polynomials *}*)

   164   fun lc (uvp: unipoly) =

   165     if nth uvp (length uvp-1) <>0

   166     then nth uvp (length uvp-1)

   167     else lc (nth_drop (length uvp-1) uvp);

   169   fun deg (uvp: unipoly) =

   170     if nth uvp (length uvp-1) <>0

   171     then length uvp - 1

   172     else deg (nth_drop (length uvp-1) uvp)

   174   fun lc_of_unipoly_not_0 [] = [](* and delete lc=0*)

   175     | lc_of_unipoly_not_0 (uvp: unipoly) =

   176       if nth uvp (length uvp -1) =0

   177       then lc_of_unipoly_not_0 (nth_drop (length uvp-1) uvp)

   178       else  uvp;

   180   fun normirt_polynom (poly1: unipoly) (m: int) =

   181     let

   182       val poly1 = lc_of_unipoly_not_0 poly1

   183       val lc_a=lc poly1;

   184       fun normirt poly1 b m lc_a i =

   185         if i=0

   186         then [mod_div (nth poly1 i) lc_a m]@b

   187         else normirt poly1 ( [mod_div(nth poly1 i) lc_a m]@b) m lc_a (i-1) ;

   188     in

   189       if length(poly1)=0

   190       then poly1

   191       else normirt poly1  [] m lc_a (length(poly1)-1)

   192   end

   194    infix %*

   195   fun (a: unipoly) %* b =  map2 Integer.mult a (replicate (length (a)) b)

   197   infix %/

   198   fun (poly: unipoly) %/ b = (* =quotient*)

   199     let fun division poly b = poly div' b;

   200     in map2 division poly (replicate (length (poly)) b) end

   202   infix %-

   203   fun (poly: unipoly) %- b =

   204     let fun minus poly b = poly - b;

   205   in map2 minus poly (replicate (length (poly)) b) end

   207   infix %+%

   208   fun (a: unipoly) %+% (b: unipoly) =  map2 Integer.add a b

   210   infix %-%

   211   fun (a: unipoly) %-% (b: unipoly) =

   212     let fun minus a b = a-b;

   213     in  map2 minus a b end

   215   (* if poly2|poly1 *)

   216   infix %|%

   217   fun [b: int] %|% [a: int] =

   218     if abs b<= abs a andalso a mod b = 0

   219     then true

   220     else false

   221   | (poly2: unipoly) %|% (poly1: unipoly) =

   222     let

   223     val b = (replicate (length poly1 - length(lc_of_unipoly_not_0 poly2)) 0) @ lc_of_unipoly_not_0 poly2 ;

   224     val c = lc poly1 div' lc b;

   225     val rest = lc_of_unipoly_not_0 (poly1 %-% (b %* c));

   226     in

   227       if rest = []

   228       then true

   229       else

   230         if c<>0 andalso length rest >= length poly2

   231         then poly2 %|% rest

   232         else false

   233     end

   235   fun poly_centr (poly: unipoly) (m: int) =

   236     let

   237       fun midle m =

   238         let val mid = m div' 2

   239         in

   240           if m mod mid = 0

   241           then  mid

   242           else  mid+1

   243         end

   244       fun centr a m mid =

   245         if a > mid

   246         then a - m

   247         else a

   248       fun polyCentr poly poly' m mid counter =

   249         if length(poly) > counter

   250         then polyCentr poly (poly' @ [centr (nth poly counter) m mid]) m mid (counter+1)

   251         else (poly': unipoly)

   252     in polyCentr poly [] m (midle m) 0

   253   end

   255   fun sum_lmb (a: unipoly) exp =

   256     let fun  sum' [a] _ = a

   257            | sum' (a: int list) exp =

   258       sum' ([((nth a 0)) + ( Integer.pow exp (nth a 1))] @ nth_drop 0 (nth_drop 0 a)) exp

   259     in sum' ([Integer.pow exp (nth a 0)] @ nth_drop 0 a) exp

   260   end;

   261 Integer.min ;

   262   fun landau_mignotte_bound a b =

   263   (Integer.pow (Integer.min (deg a) (deg b)) 2) * (abs (Integer.gcd (lc a) (lc b)))*

   264   (Int.min (abs((aprox_root (sum_lmb a 2)) div ~(lc a)), abs(( aprox_root (sum_lmb b 2)) div ~(lc b))));

   266 (*subsection {* modulo calculations for polynomials *}*)

   268   fun CRA_2_poly (ma, mb) (r1, r2) =

   269   let

   270     fun CRA_2' (ma, mb) (r1, r2) = CRA_2 (r1, ma,r2, mb)

   271   in  map (CRA_2' (ma, mb)) (r1 ~~ r2) end

   273   infix poly_mod

   274   fun uvp poly_mod m =

   275     let fun poly_mod' uvp m n =

   276       if n < length (uvp)

   277       then poly_mod' ((nth_drop 0 uvp)@[(nth uvp 0) mod m]) m ( n + 1 )

   278       else uvp (*end of poly_mod'*)

   279     in

   280       poly_mod' uvp m 0 end

   282   fun mod_poly_gcd (upoly1: unipoly) (upoly2: unipoly) (m: int) =

   283     let

   284       val moda = upoly1 poly_mod m;

   285       val modb = (replicate (length upoly1 - length(lc_of_unipoly_not_0 (upoly2 poly_mod m))) 0)

   286                   @ (lc_of_unipoly_not_0 (upoly2 poly_mod m)) ;

   287       val c =  mod_div (lc moda) (lc modb) m

   288       val rest = lc_of_unipoly_not_0 ((moda %-% (modb %*  c)) poly_mod m)

   289     in if rest = []

   290         then [upoly2, [0]]

   291         else

   292           if length rest < length upoly2

   293           then mod_poly_gcd upoly2 rest m

   294           else mod_poly_gcd rest upoly2 m

   295     end

   297   fun mod_gcd_p (poly1: unipoly) (poly2: unipoly) (p: int) =

   298     let val gcd_p = mod_poly_gcd poly1 poly2 p

   299     in normirt_polynom (nth gcd_p 0) p

   300     end

   302   fun gcd_n (up: unipoly) =

   303     if length up = 2

   304     then abs (Integer.gcd (nth up 0)(nth up 1))

   305     else gcd_n ([abs (Integer.gcd (nth up 0)(nth up 1))]@(nth_drop 0 (nth_drop 0 up)))

   307 fun gcd_n up = abs (Integer.gcds up);

   309   fun pp [_: int] = [1]

   310     | pp (poly1: unipoly) =

   311     if (poly1 poly_mod (gcd_n poly1)) = (replicate (length (poly1)) 0)

   312     then  poly1 %/ (gcd_n poly1)

   313     else poly1;

   315 "=========== code for [1] p.93: univariate ==============";

   316 "=========== code for [1] p.93: univariate ==============";

   317 "=========== code for [1] p.93: univariate ==============";

   318 (*subsection {* GCD_MOD Algorgithmus, code for [1] p.93: univariate *}*)

   320  fun GCD_MOD (a: unipoly) (b: unipoly) =

   321    if a = [0] orelse  b = [0] then [1] else

   322    if a = b then a else

   323    let

   324 (*1*)val d = abs (Integer.gcd (lc a) (lc b));

   325      val M  = 2*d*landau_mignotte_bound a b;

   326     fun GOTO2 a b d M p =   (*==============================*)

   327        let

   328 (*2*)    val p = find_next_prime_not_divide d p

   329          val c_p = normirt_polynom ( mod_gcd_p a b p) p

   330          val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p

   331          fun GOTO3 a b d M p g_p =  (*~~~~~~~~~~~~~~~~~~~~~~*)

   332 (*3*)      if (deg g_p) = 0

   333            then  [1]

   334            else

   335              let

   336                val P = p

   337                val g = g_p

   338                fun WHILE a b d M P g p = (*------------------*)

   339 (*4*)            if P > M

   340                  then g

   341                  else

   342                    let

   343                      val p = find_next_prime_not_divide d p

   344                      val c_p = normirt_polynom ( mod_gcd_p a b p) p

   345                      val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p

   346                    in

   347                      if deg g_p < deg g

   348                      then GOTO3 a b d M p g_p

   349                      else

   350                        if deg g_p = deg g

   351                        then

   352                          let

   353                            val g = CRA_2_poly (P,p)(g,g_p)

   354                            val P = P*p

   355                            val g = poly_centr(g poly_mod P)P

   356                          in WHILE a b d M P g p end

   357                        else WHILE a b d M P g p

   358                end (*----------------------------------*)

   361                val g = WHILE a b d M P g p (* << 1.Mal -----*)

   362              in g end (*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)

   364          val g = GOTO3 a b d M p g_p (* << 1.Mal ~~~~~~~~~~*)

   365 (*5*)    val g = pp g

   366         in

   367           if (g %|% a) andalso (g %|% b)

   368           then  g

   369           else GOTO2 a b d M p

   370         end (*==============================================*)

   373      val  g = GOTO2 a b d  M (*p*)1(* << 1. Mal  =============*)

   374    in g end;

   376 "----------- auxiliary functions multivariate -----------";

   377 "----------- auxiliary functions multivariate -----------";

   378 "----------- auxiliary functions multivariate -----------";

   379 (*subsection {* GCD_MODm Algorgithmus, auxiliary functions multivariate *}*)

   381 type monom = (int * int list);

   382 type multipoly = monom list;

   384 (*subsection {* calculations for multivariate polynomials *}*)

   386   fun listgreater a b =

   387     let fun gr a b = a>=b;

   388     fun all [] = true

   389       | all list =

   390       if nth list 0 then all (nth_drop 0 list)  else false

   391     in all (map2 gr a b)

   392   end

   394   fun get_coef (mvp: multipoly) (place: int) =

   395     let fun coef ((coef,_): monom) = coef;

   396     in coef (nth mvp place) end

   398   fun get_varlist (mvp: multipoly) (place: int) =

   399     let fun varlist ((_,var): monom) = var;

   400     in varlist (nth mvp place) end

   402   fun get_coeflist (poly: multipoly) =

   403     let

   404       fun get_coeflist' (poly: multipoly) list =

   405         if poly = []

   406          then

   407           list

   408          else

   409            get_coeflist' (nth_drop 0 poly) (list @ [get_coef poly 0])

   410     in

   411       get_coeflist' poly []

   412   end

   414   fun add_var_to_multipoly (mvp: multipoly) (order: int) =

   415     let fun add ([]: multipoly) (_: int) (new: multipoly) = new

   416           | add (mvp: multipoly) (order: int) (new: multipoly) =

   417           let val (first, last) = chop order (get_varlist mvp 0)

   418           in add (nth_drop 0 mvp) (order) (new @ [((get_coef mvp 0),(first @ [0] @ last))]) end

   419     in add mvp order [] end

   421   fun cero_multipoly 1 = [(0,[0])]

   422     | cero_multipoly (diferent_var: int) =

   423       add_var_to_multipoly (cero_multipoly (diferent_var-1)) 0;

   425   fun short_mv (mvp: multipoly) =

   426     let fun short (mvp: multipoly) (new: multipoly) =

   427               if length mvp =1

   428               then if (get_coef mvp 0) = 0

   429                    then if new = [] then cero_multipoly (length (get_varlist mvp 0)) else new

   430                    else  new @ mvp

   431               else if get_varlist  mvp 0 = get_varlist mvp 1

   432                 then short ( [(get_coef mvp 0 + get_coef mvp 1,get_varlist  mvp 0)]

   433                              @ nth_drop 0 (nth_drop 0 mvp))   new

   434                 else if (get_coef mvp 0) = 0 then short (nth_drop 0 mvp) new

   435                   else short (nth_drop 0 mvp) (new @ [nth mvp 0])

   436     in short mvp [] end

   438   (* if a is greater than b *)

   439   fun greater_var [a: int] [b: int] =

   440     a > b

   441     | greater_var (a: int list) (b: int list) =

   442     if  (nth a (length a -1))= (nth b (length b -1))

   443       then greater_var (nth_drop (length a-1) a) (nth_drop (length b-1) b)

   444       else (nth a (length a -1)) > (nth b (length b -1))

   446   fun order_multipoly (a: multipoly)=

   447     let

   448       val ordered = [nth a 0];

   449       fun order_mp [] [] = cero_multipoly (length (get_varlist a 0))

   450         | order_mp [] (ordered: multipoly) = short_mv ordered

   451         | order_mp a (ordered: multipoly) =

   452         if  greater_var (get_varlist a 0) (get_varlist ordered (length ordered -1))

   453           then order_mp (nth_drop 0 a)(ordered @ [nth a 0])

   454           else let

   455             val rest = [nth ordered (length ordered -1)];

   456             fun order_mp' [] (new: multipoly) (rest: multipoly) = new @ rest

   457               | order_mp' (ordered': multipoly) (new: multipoly) (rest: multipoly) =

   458                 if greater_var (get_varlist new 0) (get_varlist ordered' (length ordered' -1))

   459                   then ordered' @ new @ rest

   460                   else order_mp' (nth_drop (length ordered' -1) ordered') new

   461                                  ([nth ordered' (length ordered' -1)]@ rest)

   462             in order_mp (nth_drop 0 a)

   463                         (order_mp' (nth_drop (length ordered -1) ordered) [nth a 0] rest)

   464           end

   465     in order_mp (nth_drop 0 a) ordered

   466     end

   468   fun lc_multipoly (mpoly: multipoly) =

   469    get_coef (order_multipoly mpoly) (length (order_multipoly mpoly) - 1);

   471  (* greatest variablegroup  *)

   472   fun deg_multipoly (mvp: multipoly) =

   473     get_varlist (order_multipoly mvp) (length (order_multipoly mvp) -1)

   475   fun  max_deg_var [m: monom] (x: int) = nth (get_varlist [m] 0) x |

   476    max_deg_var (mvp: multipoly) (x: int) =

   477     let

   478       fun max_monom (m1: monom) (m2: monom) (x: int)=

   479         if nth (get_varlist [m1] 0) x < nth (get_varlist [m2] 0) x  then m2  else m1;

   480     in

   481     max_deg_var ([max_monom (nth(mvp) 0)( nth(mvp) 1) x] @ nth_drop 0 (nth_drop 0 mvp)) x

   482     end

   484   fun greater_deg (monom1: monom) (monom2: monom)=

   485     if (max_deg_var [monom1] (length(get_varlist [monom1] 0) -1)) >

   486        (max_deg_var [monom2] (length (get_varlist [monom2] 0) -1))

   487       then 1

   488       else if (max_deg_var [monom1] (length(get_varlist [monom1] 0) -1)) <

   489               (max_deg_var [monom2] (length (get_varlist [monom2] 0) -1))

   490         then 2

   491         else if length (get_varlist [monom1] 0) >1

   492           then

   493             greater_deg (1,(nth_drop 0 (get_varlist [monom1] 0))) (1,(nth_drop 0 (get_varlist [monom2] 0)))

   494           else 0

   496   infix %%/

   497   fun (poly: multipoly) %%/ b = (* =quotient*)

   498     let fun division monom b = (((get_coef [monom] 0) div' b),get_varlist [monom] 0);

   499     in order_multipoly(map2 division poly (replicate (length(poly)) b)) end

   501   fun cont [a: int] = a

   502     | cont (poly1: unipoly) = gcd_n poly1

   503   fun cont_multipoly [(a,_): monom] = a

   504     | cont_multipoly (poly: multipoly) = gcd_n (get_coeflist poly)

   506   fun evaluate_mv (mvp: multipoly) (var: int) (value: int) =

   507     let fun eval ([]: multipoly) (_: int) (_: int) (new: multipoly) = order_multipoly new |

   508       eval (mvp: multipoly) (var: int) (value: int) (new: multipoly) =

   509       eval (nth_drop 0 mvp) var value

   510            (new @  [((Integer.pow  (nth (get_varlist mvp 0) var)value) * get_coef mvp 0,

   511                     nth_drop var (get_varlist mvp 0))]);

   512     in eval mvp var value [] end

   514   fun multi_to_uni (mvp: multipoly) =

   515     if length (get_varlist mvp 0) = 1

   516     then let fun mtu ([]: multipoly) (uvp: unipoly) = uvp |

   517                  mtu (mvp: multipoly) (uvp: unipoly) =

   518                if length uvp = (nth(get_varlist mvp 0) 0)

   519                then mtu (nth_drop 0 mvp) (uvp @ [get_coef mvp 0])

   520                else mtu mvp (uvp @ [0])

   521          in mtu (order_multipoly mvp) [] end

   522     else error "Polynom has more than one variable!";

   524   fun uni_to_multi (uvp: unipoly) =

   525     let fun utm ([]: unipoly) (mvp: multipoly) (_: int)=  short_mv mvp

   526           | utm (uvp: unipoly) (mvp: multipoly) (counter: int) =

   527           utm (nth_drop 0 uvp) (mvp @ [((nth uvp 0),[counter])]) (counter+1)

   528     in  utm uvp [] 0 end

   530   fun find_new_r (mvp1: multipoly) (mvp2: multipoly) (old: int) =

   531     let val poly1 = evaluate_mv mvp1 (length (get_varlist mvp1 0) - 2) (old + 1)

   532         val poly2 = evaluate_mv mvp2 (length (get_varlist mvp2 0) - 2) (old + 1);

   533     in

   534       if max_deg_var poly1 (length (get_varlist poly1 0) - 1)= max_deg_var mvp1 (length (get_varlist mvp1 0) - 1) orelse

   535          max_deg_var poly2 (length (get_varlist poly2 0) - 1)= max_deg_var mvp2 (length (get_varlist mvp2 0) - 1)

   536       then old + 1

   537       else find_new_r (mvp1) (mvp2) (old + 1)

   538     end

   540   fun mult_with_new_var ([]: multipoly) (_: unipoly) (_: int) = []

   541     | mult_with_new_var (mvp: multipoly) (uvp: unipoly) (order: int) =

   542     let fun mult ([]: multipoly) (_: unipoly) (_: int) (new: multipoly) (_: int) = short_mv new

   543           | mult (mvp: multipoly) (uvp: unipoly) (order: int) (new: multipoly) (_: int)  =

   544           let fun mult' (_: multipoly) ([]: unipoly) (_: int) (new': multipoly) (_) = order_multipoly new'

   545                 | mult' (mvp': multipoly) (uvp': unipoly) (order: int) (new': multipoly) (c2': int) =

   546                 let val (first, last) = chop order (get_varlist mvp' 0)

   547                 in mult' mvp' (nth_drop 0 uvp') order

   548                         (new' @ [(((get_coef mvp' 0) * (nth uvp' 0)),(first @ [c2'] @ last))]) (c2'+1) end

   549           in mult (nth_drop 0 mvp) uvp order (new @ (mult' mvp uvp order [] 0)) (0) end

   550     in mult mvp uvp order [] 0 end

   552   fun greater_deg_multipoly (var1: int list) (var2: int list) =

   553     if var1 = [] andalso var2 =[] then 0

   554     else if (nth var1 (length var1 -1)) = (nth var2 (length var1 -1) )

   555          then greater_deg_multipoly (nth_drop (length var1 -1) var1) (nth_drop (length var1 -1) var2)

   556          else if (nth var1 (length var1 -1)) > (nth var2 (length var1 -1))

   557               then 1 else 2 ;

   559   infix %%+%%

   560   fun ([]: multipoly) %%+%% (mvp2: multipoly) = mvp2

   561     | (mvp1: multipoly) %%+%% ([]: multipoly) = mvp1

   562     | (mvp1: multipoly) %%+%% (mvp2: multipoly) =

   563     let fun plus ([]: multipoly) (mvp2: multipoly) (new: multipoly) = order_multipoly (new @ mvp2)

   564           | plus (mvp1: multipoly) ([]: multipoly) (new: multipoly) = order_multipoly (new @ mvp1)

   565           | plus (mvp1: multipoly) (mvp2: multipoly) (new: multipoly) =

   566           if greater_deg_multipoly (get_varlist mvp1 0) (get_varlist mvp2 0) = 0

   567           then plus (nth_drop 0 mvp1) (nth_drop 0 mvp2)

   568                     (new @ [(((get_coef mvp1 0) + (get_coef mvp2 0)), get_varlist mvp1 0)])

   569           else if greater_deg_multipoly (get_varlist mvp1 0) (get_varlist mvp2 0) = 1

   570                then plus mvp1 (nth_drop 0 mvp2) (new @ [nth mvp2 0])

   571                else plus (nth_drop 0 mvp1) mvp2 (new @ [nth mvp1 0])

   572     in plus mvp1 mvp2 [] end

   574   infix %%-%%

   575   fun (mvp1: multipoly) %%-%% (mvp2: multipoly) =

   576     let fun neg ([]: multipoly) (new: multipoly) = order_multipoly new

   577           | neg (mvp: multipoly) (new: multipoly) =

   578           neg (nth_drop 0 mvp) (new @ [(((get_coef mvp 0) * ~1), get_varlist mvp 0)])

   579         val neg_mvp2 = neg mvp2

   580     in mvp1 %%+%% (neg_mvp2 [])  end

   582   infix %%*%%

   583   fun (mvp1: multipoly) %%*%% (mvp2: multipoly) =

   584     let fun mult ([]: multipoly) (_: multipoly) (_: multipoly) (new: multipoly) = order_multipoly new

   585           | mult (mvp1: multipoly) ([]: multipoly) (regular_mvp2: multipoly) (new: multipoly) = mult (nth_drop 0 mvp1) regular_mvp2 regular_mvp2 new

   586           | mult (mvp1: multipoly) (mvp2: multipoly) (regular_mvp2: multipoly) (new: multipoly) =

   587           mult mvp1 (nth_drop 0 mvp2) regular_mvp2 (new @ [(((get_coef mvp1 0) * (get_coef mvp2 0)), ((get_varlist mvp1 0) %+% (get_varlist mvp2 0)))])

   588   in if (length mvp1) > (length mvp2) then mult mvp1 mvp2 mvp2 [] else mult mvp2 mvp1 mvp1 [] end

   590  (* if poly1|poly2 *)

   591   infix %%|%%

   592   fun [(coef2, var2): monom] %%|%%  [(coef1, var1): monom]  =

   593     ( (listgreater var1 var2)

   594         andalso (coef1 mod coef2) = 0)

   595     | (_: multipoly) %%|%% [(0,_)]=  true

   596     | (poly2: multipoly)  %%|%% (poly1: multipoly) =

   597     if [nth poly2 (length poly2 -1)]   %%|%% [nth poly1 (length poly1 -1)]

   598     then poly2 %%|%% (poly1 %%-%%

   599       ([((get_coef poly1 (length poly1 -1)) div (get_coef poly2 (length poly2 -1)),

   600       (get_varlist poly1 (length poly1 -1)) %-%(get_varlist poly2 (length poly2 -1)))] %%*%%

   601       poly2))

   602     else false

   604 (*subsection {* Newtoninterpolation *}*)

   605  (* first step *)

   606   fun newton_first (x: int list) (f: multipoly list) (order: int) =

   607     let val polynom =(add_var_to_multipoly (nth f 0) order) %%+%%

   608             ((mult_with_new_var (((nth f 1)%%-%% (nth f 0))%%/

   609                (nth x 1) - (nth x 0))) [(nth x 0) * ~1, 1] order);

   610         val new_value_poly =   [(nth x 0) * ~1, 1];

   611         val steps = [((nth f 1) %%-%%(nth f 0))%%/ ((nth x 1) - (nth x 0))]

   612     in (polynom, new_value_poly, steps) end

   614   fun newton_mv (x: int list) (f: multipoly list) (steps: multipoly list) (t: unipoly) (p: multipoly) (order: int) =

   615     if length x = 2

   616     then let val (polynom, new_value_poly, steps) =  newton_first x [(nth f 0), (nth f 1)] order

   617          in (polynom, new_value_poly, steps) end

   618     else let val new_value_poly = multi_to_uni((uni_to_multi t)  %%*%% (uni_to_multi [(nth x (length x -2) )* ~1, 1]));

   619              val new_steps = [((nth f (length f -1))  %%-%% (nth f (length f -2)))  %%/ ((nth x (length x - 1)) - (nth x (length x - 2)))];

   620              fun next_step ([]: multipoly list) (new_steps: multipoly list) (_: int list) = new_steps

   621                | next_step (steps: multipoly list) (new_steps: multipoly list) (x': int list) =

   622                  next_step (nth_drop 0 steps)

   623                            (new_steps @ [(((nth new_steps (length new_steps -1)) %%-%%(nth steps 0)) %%/

   624                                     ((nth x' (length x' - 1)) - (nth x' (length x' - 3))))])

   625                            ( nth_drop (length x' -2) x')

   626              val steps = next_step steps new_steps x;

   627              val polynom' =  (p %%+%% (mult_with_new_var (nth steps (length steps -1)) new_value_poly order));

   628        in (order_multipoly(polynom'), new_value_poly, steps) end;

   630 "=========== code for [1] p.95: multivariate ============";

   631 "=========== code for [1] p.95: multivariate ============";

   632 "=========== code for [1] p.95: multivariate ============";

   633 (*subsection {* GCD_MODm algorithm, code for [1] p.95: multivariate *}*)

   635 fun GCD_MODm a b n s r=

   636   if greater_var (deg_multipoly b) (deg_multipoly a) then GCD_MODm b a n s r else

   637 (*0*)  if s = 0

   638     then  uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %*

   639           (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))))

   640       else

   641       let

   642 (*1*)   val M = 1 + Int.min (max_deg_var a (length(get_varlist a 0)-2), max_deg_var b (length(get_varlist b 0)-2));

   643 (*2*)     fun GOTO2 a b n s M r_list steps =  (*==============================*)

   644               let

   645                 val r = find_new_r a b r;

   646                 val r_list = r_list @ [r];

   647                 val g_r = GCD_MODm (order_multipoly (evaluate_mv a (s-1) r)) ( order_multipoly (evaluate_mv b (s-1) r)) (n-1) (s-1) 0

   648                 val steps = steps @ [g_r];

   649 (*3*)           fun GOTO3 a b n s M g_r r r_list steps = (* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ *)

   650                   let

   651                     val m = 1;

   652                     fun WHILE a b n s M m r r_list steps g g_n mult =   (* ----------------------- *)

   653                       if m > M

   654                         then

   655                          if g_n %%|%% a andalso g_n %%|%% b

   656                           then  g_n

   657                           else GOTO2 a b n s M r_list steps

   658                         else

   659                         let

   660                         val r = find_new_r a b r;

   661                         val r_list = r_list @ [r];

   662                         val g_r = GCD_MODm (order_multipoly  (evaluate_mv a (s-1) r))

   663                                             (order_multipoly  (evaluate_mv b (s-1) r)) (n-1) (s-1) 0

   664                         in

   665                           if greater_var  (deg_multipoly g) (deg_multipoly g_r)

   666                             then GOTO3 a b n s M g_r r r_list steps

   667                             else if (deg_multipoly g)= (deg_multipoly g_r)

   668                               then

   669                                 let

   670                                   val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s-1)

   671                                 in if (nth steps (length steps -1)) = (cero_multipoly s)

   672                                    then WHILE a b n s M (M+1) r r_list steps g_r g_n new

   673                                    else WHILE a b n s M (m+1) r r_list steps g_r g_n new

   674                                 end

   675                               else WHILE a b n s M (m+1) r r_list steps g  g_n mult

   676                         end (* WHILE *)

   677                   in (* GOTO3*) (* ----------------------- *)

   678                     WHILE a b n s M m r r_list steps g_r ( cero_multipoly (s+1)) [1]

   679                   end (*GOTO3*)

   680               in (* GOTO2*)  (* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ *)

   681                 GOTO3 a b n s M g_r r r_list steps

   682               end (*GOTO2*)

   683       in  (*==============================*)

   684       GOTO2 a b n s M [] []

   685     end (*end 0*);

   687 "--- begin of univariate polynomials --------------------";

   688 "--- begin of univariate polynomials --------------------";

   689 "--- begin of univariate polynomials --------------------";

   691 "----------- fun mod_inv --------------------------------";

   692 "----------- fun mod_inv --------------------------------";

   693 "----------- fun mod_inv--- -----------------------------";

   694 "~~~~~ fun mod_inv, args:"; val (r, m) = (5,7);

   695 "~~~~~ fun modi, args:"; val (r, rold, m, rinv) = (r, r, m, 1);

   696 rinv < m; (*=true*)

   697 r mod m = 1; (*=false*)

   698 "~~~~~ fun modi, args:"; val (r, rold, m, rinv) = (rold * (rinv + 1), rold, m, rinv + 1);

   699 rinv < m;(*=true*)

   700 r mod m = 1; (*=false*)

   701 "~~~~~ fun modi, args:"; val (r, rold, m, rinv) = (rold * (rinv + 1), rold, m, rinv + 1);

   702 rinv < m;(*=true*)

   703 r mod m = 1; (*=true*)

   704 "~~~~~ to mod_inv return val:"; val (rinv) = (rinv);

   706 if mod_inv 5 7 = 3 then () else error "mod_inv 5 7 = 3 changed";

   707 if mod_inv 3 7 = 5 then () else error "mod_inv 3 7 = 5 changed";

   708 if mod_inv 4 339 = 85 then () else error "mod_inv 4 339 = 85 changed";

   710 "----------- fun CRA_2 ----------------------------------";

   711 "----------- fun CRA_2 ----------------------------------";

   712 "----------- fun CRA_2 ----------------------------------";

   714 if  CRA_2(17,9,3,4) = 35 then () else error "CRA_2(17,9,3,4)=35 changed";

   715 if  CRA_2(7,2,6,11) = 17  then () else error "CRA_2(7,2,6,11) = 17 changed";

   717 "----------- fun lc -------------------------------------";

   718 "----------- fun lc -------------------------------------";

   719 "----------- fun lc -------------------------------------";

   720 if lc [3,4,5,6] = 6 then () else error "lc (3,4,5,6) = 6 changed" ;

   721 if lc [3,4,5,6,0] = 6 then () else error "lc (3,4,5,6,0) = 6 changed"  ;

   723 "----------- fun deg ------------------------------------";

   724 "----------- fun deg ------------------------------------";

   725 "----------- fun deg ------------------------------------";

   726 if deg [3,4,5,6] = 3 then () else error "lc (3,4,5,6) = 3 changed" ;

   727 if deg [3,4,5,6,0] = 3 then () else error "lc (3,4,5,6) = 6 changed";

   729 "----------- fun primeList ------------------------------";

   730 "----------- fun primeList ------------------------------";

   731 "----------- fun primeList ------------------------------";

   732 if primeList 1 = [2] then () else error " primeList 1 = [2] changed";

   733 if primeList 3 = [2,3] then () else error " primeList 3 = [2,3] changed";

   734 if primeList 6 = [2,3,5,7] then () else error " primeList 6 = [2,3,5,7] changed";

   735 if primeList 7 = [2,3,5,7] then () else error " primeList 7 = [2,3,5,7] changed";

   737 "----------- fun find_next_prime_not_divide -----------------";

   738 "----------- fun find_next_prime_not_divide -----------------";

   739 "----------- fun find_next_prime_not_divide -----------------";

   740 if find_next_prime_not_divide 15 2 = 7 then () else error "findNextPrimeNotdivide 15 2 = 7 changed";

   741 if find_next_prime_not_divide 90 1 = 7 then () else error "findNextPrimeNotdivide 90 1 = 7 changed";

   743 "----------- fun poly_mod -------------------------------";

   744 "----------- fun poly_mod -------------------------------";

   745 "----------- fun poly_mod -------------------------------";

   746 if ([5,4,7,8,1] poly_mod 5) = [0, 4, 2, 3, 1] then ()

   747 else error "[5,4,7,8,1] poly_mod 5 = [0, 4, 2, 3, 1] changed" ;

   748 if ([5,4,~7,8,~1] poly_mod 5) = [0, 4, 3, 3, 4] then ()

   749 else error "[5,4,~7,8,~1] poly_mod 5  = [0, 4, 3, 3, 4] changed" ;

   751 "----------- fun %* ------------------------------";

   752 "----------- fun %* ------------------------------";

   753 "----------- fun %* ------------------------------";

   754 if ([5,4,7,8,1] %* 5) = [25, 20, 35, 40, 5] then ()

   755 else error "[5,4,7,8,1] %* 5 = [25, 20, 35, 40, 5] changed";

   756 if ([5,4,~7,8,~1] %* 5) = [25, 20, ~35, 40, ~5] then ()

   757 else error "[5,4,~7,8,~1] %* 5 = [25, 20, ~35, 40, ~5] changed";

   759 "----------- fun %/ -------------------------------";

   760 "----------- fun %/ -------------------------------";

   761 "----------- fun %/ -------------------------------";

   762 if ([4,3,2,5,6] %/ 3) =[1, 1, 0, 1, 2] then ()

   763 else error "%/ [4,3,2,5,6] 3 =[1, 1, 0, 1, 2] changed";

   764 if ([4,3,2,0] %/ 3) =[1, 1, 0, 0] then ()

   765 else error "%/ [4,3,2,0] 3 =[1, 1, 0, 0] changed";

   767 "----------- fun %-% ------------------------";

   768 "----------- fun %-% ------------------------";

   769 "----------- fun %-% ------------------------";

   770 if ([8,~7,0,1] %-% [~2,2,3,0])=[10,~9,~3,1] then ()

   771 else error "[8,~7,0,1] %-% [~2,2,3,0]=[10,~9,~3,1] changed";

   772 if ([8,7,6,5,4] %-% [2,2,3,1,1])=[6,5,3,4,3] then ()

   773 else error "[8,7,6,5,4] %-% [2,2,3,1,1]=[6,5,3,4,3] changed";

   775 "----------- fun %+% ------------------------------";

   776 "----------- fun %+% ------------------------------";

   777 "----------- fun %+% ------------------------------";

   778 if ([8,~7,0,1] %+% [~2,2,3,0])=[6,~5,3,1] then ()

   779 else error "[8,~7,0,1] %+% [~2,2,3,0]=[6,~5,3,1] changed";

   780 if ([8,7,6,5,4] %+% [2,2,3,1,1])=[10,9,9,6,5] then ()

   781 else error "[8,7,6,5,4] %+% [2,2,3,1,1]=[10,9,9,6,5] changed";

   783 "----------- fun CRA_2_poly -----------------------------";

   784 "----------- fun CRA_2_poly -----------------------------";

   785 "----------- fun CRA_2_poly -----------------------------";

   786 if (CRA_2_poly (5, 7)([2,2,4,3],[3,2,3,5]))=[17,2,24,33] then ()

   787 else error "CRA_2_poly (5, 7)([2,2,4,3],[3,2,3,5])=[17,2,24,33] changed";

   789 "----------- fun mod_div --------------------------------";

   790 "----------- fun mod_div --------------------------------";

   791 "----------- fun mod_div --------------------------------";

   792 if mod_div 21 4 5 = 4 then () else error "mod_div 21 4 5 = 4 changed";

   793 if mod_div 1 4 5 = 4 then () else error "mod_div 1 4 5 = 4 changed";

   794 if mod_div 0 4 5 = 0 then () else error "mod_div 0 4 5 = 0 changed";

   796 "----------- fun lc_of_unipoly_not_0 --------------------";

   797 "----------- fun lc_of_unipoly_not_0 --------------------";

   798 "----------- fun lc_of_unipoly_not_0 --------------------";

   799 if lc_of_unipoly_not_0 [0,1,2,3,4,5,0,0]=[0,1, 2, 3, 4, 5] then ()

   800 else error "lc_of_unipoly_not_0 [0,1,2,3,4,5,0,0]=[0,1, 2, 3, 4, 5] changed";

   801 if lc_of_unipoly_not_0 [0,1,2,3,4,5]=[0,1, 2, 3, 4, 5] then ()

   802 else error "lc_of_unipoly_not_0 [0,1,2,3,4,5]=[0, 1, 2, 3, 4, 5] changed";

   804 "----------- fun %|% -------------------------------";

   805 "----------- fun %|% -------------------------------";

   806 "----------- fun %|% -------------------------------";

   807 "~~~~~ fun %|% , args:"; val (poly1, poly2) = ([9,12,4],[3,2]);

   808 val b = (replicate (length poly1 - length poly2) 0 ) @ poly2 ;

   809 val c = lc poly1 div lc b;

   810 val rest = lc_of_unipoly_not_0 ( poly1 %-% (b %* c ));

   811 rest = [];(*=false*)

   812 c<>0 andalso length rest >= length poly2; (*=true*)

   813 "~~~~~ fun %|% , args:"; val (poly1, poly2) = (rest, poly2);

   814 val b = (replicate (length poly1 - length poly2) 0 ) @ poly2 ;

   815 val c = lc poly1 div lc b;

   816 val rest = lc_of_unipoly_not_0 ( poly1 %-% (b %* c ));

   817 rest = [];(*=true*)

   818 "~~~~~ to  return val:"; val (divide) = (true);

   820 if ([4] %|% [6]) = false then () else error "[4] %|% [6] = false changed";

   821 if ( [8] %|%[16,0]) = true then () else error "[8] %|%[16,0] = true changed";

   822 if ([3,2] %|%[0,0,9,12,4] ) = true then ()

   823 else error "[3,2] %|%[0,0,9,12,4]  = true changed";

   824 if ([8,0] %|% [16]) = true then () else error "[8,0] %|% [16] = true changed";

   826 "----------- fun mod_poly_gcd ------------------------";

   827 "----------- fun mod_poly_gcd ------------------------";

   828 "----------- fun mod_poly_gcd ------------------------";

   830 if ( mod_poly_gcd [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 7)=[[2, 6, 0, 2, 6], [0]]

   831 then () else error "( mod_poly_gcd [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 7 = [ [5, 1, 0, 5, 1], [0]] changed";

   832 if ( mod_poly_gcd [8, 28, 22, ~11, ~14, 1, 2] [2, 6, 0, 2, 6] 7)=[[2, 6, 0, 2, 6], [0]] then ()

   833 else error "mod_poly_gcd [8, 28, 22, ~11, ~14, 1, 2] [2, 6, 0, 2, 6] 7=[[1, 3, 0, 1, 3], [0]] changed";

   834 if mod_poly_gcd [20,15,8,6] [8,~2,~2,3] 2 = [[0, 1], [0]] then ()

   835 else error " mod_poly_gcd [20,15,8,6] [8,~2,~2,3] 2 = [[0, 1], [0]] changed";

   837 "~~~~~ fun mod_poly_gcd , args:";

   838 val (a,b,m) = ([ ~13, 2,12],[ ~14, 2],5);

   839 val moda = a poly_mod m;

   840 val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b poly_mod m) ;

   841 val c =  mod_div (lc moda) (lc modb) m;

   842 val rest = lc_of_unipoly_not_0 (moda %-% (modb %* c) poly_mod m) poly_mod m;

   843 rest = []; (*=false*)

   844 length rest < length b; (*=false*)

   845 "~~~~~ fun  mod_poly_gcd , args:";

   846 val (a,b,m,d) = (rest, b ,  m , [c]);

   847 val moda = a poly_mod m

   848 val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b poly_mod m) ;

   849 val c =  mod_div (lc moda) (lc modb) m

   850 val rest = lc_of_unipoly_not_0 ((moda %-% (modb %*  c)) poly_mod m);

   851 rest = [];(*=flase*)

   852 length rest < length b; (*=true*)

   853 "~~~~~ fun  mod_poly_gcd , args:";

   854 val (a,b,m,d) = (b, rest, m, [c] @ d);

   855 val moda = a poly_mod m

   856 val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b poly_mod m) ;

   857 val c =  mod_div (lc moda) (lc modb) m

   858 val rest = lc_of_unipoly_not_0 ((moda %-% (modb %*  c)) poly_mod m);

   859 rest = [];(*=flase*)

   860 length rest < length b; (*=false*)

   861 "~~~~~ fun  mod_poly_gcd , args:";

   862 val (a,b,m,d) = (b, rest, m, [c] @ d);

   863 val moda = a poly_mod m

   864 val modb = (replicate (length a - length(lc_of_unipoly_not_0 b)) 0) @ (lc_of_unipoly_not_0 b poly_mod m) ;

   865 val c =  mod_div (lc moda) (lc modb) m

   866 val rest = lc_of_unipoly_not_0 ((moda %-% (modb %*  c)) poly_mod m);

   867 rest = [];(*=true*)

   868 "~~~~~ to  return val:"; val ([b, rest, lsit])=[b, [0], [c] @ d];

   871 "----------- fun normirt_polynom ------------------------";

   872 "----------- fun normirt_polynom ------------------------";

   873 "----------- fun normirt_polynom ------------------------";

   874 if (normirt_polynom [~18, ~15, ~20, 12, 20, ~13, 2] 5) =[1, 0, 0, 1, 0, 1, 1 ] then ()

   875 else error "(normirt_polynom [~18, ~15, ~20, 12, 20, ~13, 2] 5) =[1, 0, 0, 1, 0, 1, 1 ] changed"

   877 "----------- fun poly_centr -----------------------------";

   878 "----------- fun poly_centr -----------------------------";

   879 "----------- fun poly_centr -----------------------------";

   880 if (poly_centr [7,3,5,8,1,3] 10) = [~3, 3, 5, ~2, 1, 3] then ()

   881 else error "poly_centr [7,3,5,8,1,3] 10 = [~3, 3, 5, ~2, 1, 3] cahnged";

   883 "----------- fun mod_gcd_p -------------------------------";

   884 "----------- fun mod_gcd_p -------------------------------";

   885 "----------- fun mod_gcd_p -------------------------------";

   886 if  (mod_gcd_p [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] 5) = [1, 1, 1, 1] then ()

   887 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";

   888 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 ()

   889 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";

   892 "----------- fun gcd_n ----------------------------";

   893 "----------- fun gcd_n ----------------------------";

   894 "----------- fun gcd_n ----------------------------";

   895 if gcd_n [3,9,12,21] = 3 then () else error "gcd_n [3,9,12,21] = 3 changed";

   896 if gcd_n [12,16,32,44] = 4 then () else error "gcd_n [12,16,32,44] = 4 changed";

   897 if gcd_n [4,5,12] = 1 then () else error "gcd_n [4,5,12] = 1 changed";

   899 "----------- fun pp -------------------------------------";

   900 "----------- fun pp -------------------------------------";

   901 "----------- fun pp -------------------------------------";

   902 if pp [12,16,32,44] = [3, 4, 8, 11] then () else error "pp [12,16,32,44] = [3, 4, 8, 11] changed";

   903 if pp [4,5,12] =  [4, 5, 12] then () else error "pp [4,5,12] =  [4, 5, 12] changed";

   905 "----------- fun GCD_MOD --------------------------------";

   906 "----------- fun GCD_MOD --------------------------------";

   907 "----------- fun GCD_MOD --------------------------------";

   909 "~~~~~ fun GCD_MOD a b , args:"; val (a, b) = ([~18, ~15, ~20, 12, 20, ~13, 2], [8, 28, 22, ~11, ~14, 1, 2]);

   910 val d = abs (Integer.gcd (lc a) (lc b));

   911 val M  =2*d* landau_mignotte_bound a b;

   912 (*val  g = GOTO2 a b d (*p*)1 M*)

   913    "~~~~~ fun GOTO2 a b d M p , args:"; val (a, b, d, p, M) = (a,b,d,3,M);

   914   val p = find_next_prime_not_divide d p

   915   val c_p = normirt_polynom ( mod_gcd_p a b p) p

   916   val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;

   917   (*val (p, g) = GOTO3 a b d M p g_p (* << 1.Mal ~~~~~~~~~~*)*)

   918      "~~~~~ 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);

   919     (deg g_p) = 0; (*=false*)

   920     val P = p

   921     val g = g_p;

   922     (*(p, g) = WHILE a b d M P g (* << 1.Mal -----*)*)

   923        "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,p,g);

   924        P > M ;(*false*)

   925       val p = find_next_prime_not_divide d p

   926       val c_p = normirt_polynom ( mod_gcd_p a b p) p

   927       val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;

   928       deg g_p < deg g;   (*=fasle*)

   929       deg g_p = deg g; (*false*);

   930       "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);

   931        P > M ;(*false*)

   932       val p = find_next_prime_not_divide d p

   933       val c_p = normirt_polynom ( mod_gcd_p a b p) p

   934       val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;

   935       deg g_p < deg g;   (*=false*)

   936       deg g_p = deg g; (*=true*)

   937       val g = CRA_2_poly (P,p)(g,g_p)

   938       val P = P*p;

   939       val g = poly_centr(g poly_mod P)P;

   940       "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);

   941        P > M ;(*false*)

   942       val p = find_next_prime_not_divide d p

   943       val c_p = normirt_polynom (mod_gcd_p a b p) p

   944       val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;

   945       deg g_p < deg g;   (*=true*)

   946        "~~~~~ 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);

   947       (deg g_p) = 0; (*false*)

   948       val P = p

   949       val g = g_p;

   950       "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);

   951        P > M ;(*false*)

   952       val p = find_next_prime_not_divide d p

   953       val c_p = normirt_polynom (mod_gcd_p a b p) p

   954       val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;

   955       deg g_p < deg g;   (*=fasle*)

   956       deg g_p = deg g;(*true*)

   957       val g = CRA_2_poly (P,p)(g,g_p)

   958       val P = P*p;

   959       val g = poly_centr(g poly_mod P)P;

   960       "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);

   961        P > M ;(*false*)

   962       val p = find_next_prime_not_divide d p

   963       val c_p = normirt_polynom (mod_gcd_p a b p) p

   964       val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;

   965       deg g_p < deg g;   (*=false*)

   966       deg g_p = deg g;(*true*)

   967       val g = CRA_2_poly (P,p)(g,g_p)

   968       val P = P*p;

   969       val g = poly_centr(g poly_mod P)P;

   970       "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);

   971        P > M ;(*false*)

   972       val p = find_next_prime_not_divide d p

   973       val c_p = normirt_polynom (mod_gcd_p a b p) p

   974       val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p;

   975       deg g_p < deg g;   (*=false*)

   976       deg g_p = deg g;(*true*)

   977       val g = CRA_2_poly (P,p)(g,g_p)

   978       val P = P*p;

   979       val g = poly_centr(g poly_mod P)P;

   980       "~~~~~ fun WHILE a b d M p g , args:"; val (a, b, d, M, P, g) = (a,b,d,M,P,g);

   981        P > M ;(*true*)

   982     "~~~~~ to  return WHILE val:"; val (g,p) = (g,p);

   983   "~~~~~ to  return GOTO3 val:"; val (g,p) = (g,p);

   984    val g = pp g;

   985    (g %|% a) andalso (g %|% b);(*=true*)

   986 "~~~~~ to  return GOTO2 val:"; val (g) = (g);

   987 "~~~~~ to  return GCD_MOD val:"; val (g) = (g);

   989 "----------- fun GCD_MOD---------------------------------";

   990 "----------- fun GCD_MOD --------------------------------";

   991 "----------- fun GCD_MOD --------------------------------";

   992 if GCD_MOD [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] = [~2, ~1, 1] then ()

   993 else error "GCD_MOD [~18, ~15, ~20, 12, 20, ~13, 2] [8, 28, 22, ~11, ~14, 1, 2] = [~2, ~1, 1] changed";

   995 "------------------------------------------------------------";

   996 "--------------------- Multivariate Case --------------------";

   997 "------------------------------------------------------------";

   999 "----------- fun get_coef --------------------------- ---";

  1000 "----------- fun get_coef -------------------------------";

  1001 "----------- fun get_coef -------------------------------";

  1002 if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 0 = 1 then () else

  1003 error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 0 = 1 changed";

  1004 if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 1 = 2 then () else

  1005 error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 1 = 2 changed";

  1006 if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 2 = 3 then () else

  1007 error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 2 = 3 changed";

  1008 if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 3 = 4 then () else

  1009 error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 3 = 4 changed";

  1010 if get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 4 = 5 then () else

  1011 error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 4 = 5 changed";

  1014 "----------- fun get_varlist ----------------------------";

  1015 "----------- fun get_varlist ----------------------------";

  1016 "----------- fun get_varlist ----------------------------";

  1017 if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 0 = [1,1] then () else

  1018 error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 0 = [1,1] changed";

  1019 if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 1 = [1,2] then () else

  1020 error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 1 = [1,2] changed";

  1021 if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 2 = [1,3] then () else

  1022 error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 2 = [1,3] changed";

  1023 if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 3 = [1,4] then () else

  1024 error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 3 = [1,4] changed";

  1025 if get_varlist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 4 = [1,5] then () else

  1026 error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] 4 = [1,5] changed";

  1028 "----------- fun get_coeflist ---------------------------";

  1029 "----------- fun get_coeflist ---------------------------";

  1030 "----------- fun get_coeflist ---------------------------";

  1031 if get_coeflist [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])]  = [1,2,3,4,5] then () else

  1032 error "get_coef [(1,[1,1]),(2,[1,2]),(3,[1,3]),(4,[1,4]),(5,[1,5])] = [1,2,3,4,5] changed";

  1034 "----------- fun add_var_to_multipoly -------------------";

  1035 "----------- fun add_var_to_multipoly -------------------";

  1036 "----------- fun add_var_to_multipoly -------------------";

  1037 if add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 0 = [(1, [0, 1, 2]), (2, [0, 2, 3])] then () else

  1038 error "add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 0 = [(1, [0, 1, 2]), (2, [0, 2, 3])] changed";

  1039 if add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 1 = [(1, [1, 0,  2]), (2, [2, 0, 3])] then () else

  1040 error "add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 1 = [(1, [1, 0,  2]), (2, [2, 0, 3]) changed";

  1041 if add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 2 = [(1, [1, 2, 0]), (2, [2, 3, 0])] then () else

  1042 error "add_var_to_multipoly [(1,[1,2]),(2,[2,3])] 2 = [(1, [1, 2, 0]), (2, [2, 3, 0])] changed";

  1044 "----------- fun cero_multipoly -------------------------";

  1045 "----------- fun cero_multipoly -------------------------";

  1046 "----------- fun cero_multipoly -------------------------";

  1047 if cero_multipoly 1 = [(0, [0])] then () else

  1048 error "cero_multipoly 1 = [(0, [0])] changed";

  1049 if cero_multipoly 5 = [(0, [0, 0, 0, 0, 0])] then () else

  1050 error "cero_multipoly 1 = [(0, [0, 0, 0, 0, 0])] changed";

  1052 "----------- fun short_mv -------------------------------";

  1053 "----------- fun short_mv -------------------------------";

  1054 "----------- fun short_mv -------------------------------";

  1055 "~~~~~ fun short_mv, args:"; val (mvp) = ([(~12, [1]), (2, [1]), (4, [0])]);

  1056 "~~~~~ fun short , args:"; val (mvp, new) = (mvp,([]:multipoly));

  1057 length mvp =1(*false*);

  1058  get_varlist  mvp 0 = get_varlist  mvp 1; (* true*)

  1059 "~~~~~ fun short , args:"; val (mvp, new) =

  1060 ([(get_coef  mvp 0 + get_coef  mvp 1,get_varlist mvp 0)] @ nth_drop 0 (nth_drop 0 mvp), new );

  1061 length mvp =1(*false*);

  1062 get_varlist  mvp 0 = get_varlist mvp 1;(*false*)

  1063 "~~~~~ fun short , args:"; val (mvp, new) = ((nth_drop 0 mvp),new @ [nth mvp 0]);

  1064 length mvp =1(*true*);

  1065 "~~~~~ to  return val:"; val (shortform) = (new @ mvp);

  1067 if short_mv [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = [(1, [0, 0]), (~1, [1, 2])]

  1068 then () else error "short_mv [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = [(1, [0, 0]), (~1, [1, 2])] changed"

  1070 "----------- fun greater_var ----------------------------";

  1071 "----------- fun greater_var ----------------------------";

  1072 "----------- fun greater_var ----------------------------";

  1073 if greater_var [3] [4]

  1074 then error " greater_var [3] [4] changed = false"  else ();

  1075 if greater_var [1,2] [0,3]

  1076 then error " greater_var [1,2] [0,3] changed = false"  else ();

  1077 if greater_var  [1,2] [3,0]

  1078 then ()  else error " greater_var  [1,2] [3,0] = true changed";

  1079 if greater_var [1,2] [1,2]

  1080 then error " greater_var [1,2] [1,2] changed = false"  else ();

  1082 "----------- fun order_multipoly ------------------------";

  1083 "----------- fun order_multipoly ------------------------";

  1084 "----------- fun order_multipoly ------------------------";

  1085 if order_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0]),(~3,[1,1]),(~3,[1,2]),(4,[0,0])] = [(1, [0, 0]), (~1, [1, 2])]

  1086 then () else error "order_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0]),(~3,[1,1]),(~3,[1,2]),(4,[0,0])] = [(1, [0, 0]), (~1, [1, 2])] changed"

  1088 "----------- fun lc_multipoly ---------------------------";

  1089 "----------- fun lc_multipoly ---------------------------";

  1090 "----------- fun lc_multipoly ---------------------------";

  1091 if lc_multipoly [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = ~1

  1092 then () else error "lc_multipoly [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = ~1  changed";

  1093 if lc_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0])] = 2

  1094 then () else error "lc_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0])] = 2   changed";

  1096 "----------- fun deg_multipoly -------------------------";

  1097 "----------- fun deg_multipoly -------------------------";

  1098 "----------- fun deg_multipoly -------------------------";

  1099 if deg_multipoly [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = [1,2]

  1100 then () else error "lc_multipoly [(~3,[0,0]),(4,[0,0]),(3,[1,1]),(~3,[1,1]),(2,[1,2]),(~3,[1,2])] = ~1  changed";

  1101 if deg_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0])] = [1,2]

  1102 then () else error "lc_multipoly [(3,[1,1]),(2,[1,2]),(~3,[0,0])] = 2   changed";

  1106 "----------- fun max_deg_var ----------------------------";

  1107 "----------- fun max_deg_var ----------------------------";

  1108 "----------- fun max_deg_var ----------------------------";

  1109 if max_deg_var [(1,[1,2]),(2,[2,3]),(3,[5,4]),(1,[0,0])] 0 = 5 then ()

  1110 else error " max_deg_var [(1,[1,2]),(2,[2,3]),(3,[5,4]),(1,[0,0])] 0 = 5 changed";

  1111 if max_deg_var [(1,[1,2]),(2,[2,3]),(3,[5,4]),(1,[0,0])] 1 = 4 then ()

  1112 else error " max_deg_var [(1,[1,2]),(2,[2,3]),(3,[5,4]),(1,[0,0])] 1 = 4 changed";

  1114 "----------- infix %%/ ----------------------------------";

  1115 "----------- infix %%/ ----------------------------------";

  1116 "----------- infix %%/ ----------------------------------";

  1117 if ([(~3, [2, 0]), (~6, [1, 1]), (~1, [3, 1]),(1, [5, 0]), (2, [4, 1])] %%/ 2) = [(~1, [2, 0]), (~3, [1, 1]), (1, [4, 1])] then ()

  1118 else error "[(~3, [2, 0]), (~6, [1, 1]), (~1, [3, 1]),(1, [5, 0]), (2, [4, 1])] %%/ 2  =  [(~1, [2, 0]), (~3, [1, 1]), (1, [4, 1])] changed";

  1120 "----------- fun cont -----------------------------------";

  1121 "----------- fun cont -----------------------------------";

  1122 "----------- fun cont -----------------------------------";

  1123 if cont [4,8,12,~2,0,~4] = 2 then () else error " cont [4,8,12,~2,0,~4] = 2 changed";

  1124 if cont [3,6,9,19] = 1 then () else error " cont [3,6,9,19] = 1 changed";

  1126 "----------- fun cont_multipoly -------------------------";

  1127 "----------- fun cont_multipoly -------------------------";

  1128 "----------- fun cont_multipoly -------------------------";

  1129 if cont_multipoly [(3,[1,2,3,1])] = 3 then () else error " cont_multipoly [(3,[1,2,3,1])] = 3 changed";

  1130 if cont_multipoly [(~2, [2, 0]), (4, [5, 0]),(~4, [3, 2]), (2, [2, 3])] = 2 then ()

  1131 else error "cont_multipoly [(~2, [2, 0]), (4, [5, 0]),(~4, [3, 2]), (2, [2, 3])] = 2 changed";

  1132 if cont_multipoly [(~2, [2, 0]), (5, [5, 0]),(~4, [3, 2]), (2, [2, 3])] = 1 then ()

  1133 else error "cont_multipoly [(~2, [2, 0]), (5, [5, 0]),(~4, [3, 2]), (2, [2, 3])] = 1 changed";

  1135 "----------- fun evaluate_mv ----------------------------";

  1136 "----------- fun evaluate_mv ----------------------------";

  1137 "----------- fun evaluate_mv ----------------------------";

  1138 if evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 0 1 = [(1, [0]), (~1, [1])]

  1139 then () else error " evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 0 1 = [(1, [0]), (~1, [1])] changed";

  1140 if evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 1 1 = [(2, [0]), (~2, [2])]

  1141 then () else error " evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 1 1 =[(2, [0]), (~2, [2])] changed";

  1142 if evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 0 3 = [(9, [0]), (~25, [1])]

  1143 then () else error " evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 0 3 = [(9, [0]), (~25, [1])] changed";

  1144 if evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 1 3 = [ (6, [0]), (~8, [2])]

  1145 then () else error " evaluate_mv [(~3, [2, 1]), (2, [0, 1]),(1,[2,0])] 1 3 = [ (6, [0]), (~8, [2])] changed";

  1147 "----------- fun mutli_to_uni ---------------------------";

  1148 "----------- fun mutli_to_uni ---------------------------";

  1149 "----------- fun mutli_to_uni ---------------------------";

  1150 if multi_to_uni [(~3, [1]), (2, [1]), (1, [0])] = [1, ~1]

  1151 then () else error " multi_to_uni [(~3, [1]), (2, [1]), (1, [0])] = [1, ~1] changed";

  1153 "~~~~~ fun multi_to_uni , args:"; val (mvp) = ([(~3, [0]), (1, [0]), (3, [1]), (~6, [1]),(2,[3])]);

  1154 val mvp = order_multipoly mvp;

  1155 "~~~~~ fun short, args:"; val (mvp, uvp) = ((short_mv mvp), ([]:unipoly));

  1156 length uvp = (nth(get_varlist  mvp 0) 0);(*true*)

  1157 "~~~~~ fun short, args:"; val (mvp, uvp) = ((nth_drop 0 mvp), (uvp @ [get_coef mvp 0]));

  1158 length uvp = (nth(get_varlist  mvp 0) 0);(*true*)

  1159 "~~~~~ fun short, args:"; val (mvp, uvp) = ((nth_drop 0 mvp), (uvp @ [get_coef mvp 0]));

  1160 length uvp = (nth(get_varlist  mvp 0) 0);(*false*)

  1161 "~~~~~ fun short, args:"; val (mvp, uvp) = (mvp, (uvp @ [0]));

  1162 length uvp = (nth(get_varlist  mvp 0) 0);(*true*)

  1163 "~~~~~ fun short, args:"; val (mvp, uvp) = ((nth_drop 0 mvp), (uvp @ [get_coef mvp 0]));

  1164 "~~~~~ to  return val:"; val (unipoly) = (uvp);

  1168 "----------- fun find_new_r  ----------------------------";

  1169 "----------- fun find_new_r  ----------------------------";

  1170 "----------- fun find_new_r  ----------------------------";

  1171 if find_new_r [(2,[0,0]),(1,[1,2])] [(4,[2,2]),(1,[1,2])] ~1 = 1 then ()

  1172 else error " find_new_r [(2,[0,0]),(1,[1,2])] [(4,[2,2]),(1,[1,2])] ~1 = 1 changed";

  1173 if find_new_r [(2,[0,0]),(1,[1,2])] [(4,[2,2]),(1,[1,2])] 1 = 2 then ()

  1174 else error "find_new_r [(2,[0,0]),(1,[1,2])] [(4,[2,2]),(1,[1,2])] 1 = 2 changed";

  1175 if find_new_r [(2,[1,2]),(~4,[0,2]),(2,[1,1])] [(1,[1,3]),(~2,[0,3])] 1 = 3 then ()

  1176 else error "find_new_r [(2,[1,2]),(~4,[0,2]),(2,[1,1])] [(1,[1,3]),(~2,[0,3])] 1 = 3 changed";

  1178 "----------- fun newton_mv  -----------------------------";

  1179 "----------- fun newton_mv  -----------------------------";

  1180 "----------- fun newton_mv  -----------------------------";

  1181 if uni_to_multi [1,2,1,1,3,4] = [(1, [0]), (2, [1]), (1, [2]), (1, [3]), (3, [4]), (4, [5])] then ()

  1182 else error "uni_to_multi [1,2,1,1,3,4] = [(1, [0]), (2, [1]), (1, [2]), (1, [3]), (3, [4]), (4, [5])] changed";

  1183 if uni_to_multi [1,2,0,0,5,6] = [(1, [0]), (2, [1]), (5, [4]), (6, [5])] then ()

  1184 else error "uni_to_multi [1,2,0,0,5,6] = [(1, [0]), (2, [1]), (5, [4]), (6, [5])] changed";

  1186 "----------- fun mult_with_new_var  ---------------------";

  1187 "----------- fun mult_with_new_var  ---------------------";

  1188 "----------- fun mult_with_new_var  ---------------------";

  1189 if mult_with_new_var [(3,[0]),(2,[1])] [~1,1] 0 = [(~3, [0, 0]), (3, [1, 0]), (~2, [0, 1]), (2, [1, 1])] then ()

  1190 else error "mult_with_new_var [(3,[0]),(2,[1])] [~1,1] 0 = [(~3, [0, 0]), (3, [1, 0]), (~2, [0, 1]), (2, [1, 1])] changed";

  1192 "----------- fun greater_deg_multipoly  -----------------";

  1193 "----------- fun greater_deg_multipoly  -----------------";

  1194 "----------- fun greater_deg_multipoly  -----------------";

  1195 greater_deg_multipoly [1,2,3] [1,2,4] =2;

  1197 "----------- infix %%+%%  -------------------------------";

  1198 "----------- infix %%+%%  -------------------------------";

  1199 "----------- infix %%+%%  -------------------------------";

  1200 if ([(3,[0,0]),(2,[3,2]),(~3,[1,3])] %%+%% [(3,[3,2]),(~1,[0,0]),(2,[1,1])]) = [(2, [0, 0]), (2, [1, 1]), (5, [3, 2]), (~3, [1, 3])]

  1201 then () else error "([(3,[0,0]),(2,[3,2]),(~3,[1,3])] %%+%% [(3,[3,2]),(~1,[0,0]),(2,[1,1])]) = [(2, [0, 0]), (2, [1, 1]), (5, [3, 2]), (~3, [1, 3])] changed";

  1203 "----------- infix %%-%%  -------------------------------";

  1204 "----------- infix %%-%%  -------------------------------";

  1205 "----------- infix %%-%%  -------------------------------";

  1206 if ([(2, [0, 0]), (2, [1, 1]), (5, [3, 2]), (~3, [1, 3])] %%-%% [(3,[3,2]),(~1,[0,0]),(2,[1,1])]) = [(3, [0, 0]), (2, [3, 2]), (~3, [1, 3])]

  1207 then () else error "([(2, [0, 0]), (2, [1, 1]), (5, [3, 2]), (~3, [1, 3])] %%-%% [(3,[3,2]),(~1,[0,0]),(2,[1,1])]) = [(3, [0, 0]), (2, [3, 2]), (~3, [1, 3])] changed";

  1209 "----------- infix %%*%%  -------------------------------";

  1210 "----------- infix %%*%%  -------------------------------";

  1211 "----------- infix %%*%%  -------------------------------";

  1212 if ([(~1,[0]),(1,[1])] %%*%% [(~2,[0]),(1,[1])]) = [(2, [0]), (~3, [1]), (1, [2])]

  1213 then () else error "([(~1,[0]),(1,[1])] %%*%% [(~2,[0]),(1,[1])]) = [(2, [0]), (~3, [1]), (1, [2])] changed";

  1214 if ([(~2,[0]),(1,[1])] %%*%% [(1,[1]),(2,[0])]) = [(~4, [0]), (1, [2])] then ()

  1215 else error "([(~2,[0]),(1,[1])] %%*%% [(1,[1]),(2,[0])]) = [(~4, [0]), (1, [2])] changed";

  1216 if ([(3,[0,0]),(2,[0,1])] %%*%% [(1,[1,0]),(~1,[0,0])]) = mult_with_new_var [(3,[0]),(2,[1])] [~1,1] 0

  1217 then () else error

  1218 "([(3,[0,0]),(2,[0,1])] %%*%% [(1,[1,0]),(~1,[0,0])]) = mult_with_new_var [(3,[0]),(2,[1])] [~1,1] 0 changed";

  1220 "----------- fun newton_first  --------------------------";

  1221 "----------- fun newton_first  --------------------------";

  1222 "----------- fun newton_first  --------------------------";

  1223 if newton_first [1, 2] [[(1,[1,1])], [(4,[1,1])]] 1 =

  1224 ([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]])

  1225 then () else error

  1226 "newton_first [1, 2] [[(1,[1,1])], [(4,[1,1])]] 1 = "

  1227 "([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]]) changed";

  1229 "----------- fun newton_mv  -----------------------------";

  1230 "----------- fun newton_mv  -----------------------------";

  1231 "----------- fun newton_mv  -----------------------------";

  1232 if newton_mv [1, 2 ] [[(1,[1,1])], [(4,[1,1])]] [] [1] (cero_multipoly 2) 1 =

  1233 ([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]])

  1234 then () else error

  1235 "newton_mv [1, 2 ] [[(1,[1,1])], [(4,[1,1])]] [] [1] (cero_multipoly 2) 1 = "

  1236 "([(~2, [1, 0, 1]), (3, [1, 1, 1])], [~1, 1], [[(3, [1, 1])]]) changed";

  1237 if newton_mv [1, 2, 3 ] [[(4,[1,1])], [(9,[1,1])]] [[(3, [1, 1])]] [~1, 1] [(~2, [1, 0, 1]), (3, [1, 1, 1])] 1 =

  1238 ([(1, [1, 2, 1])], [2, ~3, 1], [[(5, [1, 1])], [(1, [1, 1])]])

  1239 then () else error

  1240 "newton_mv [1, 2, 3 ] [[(4,[1,1])], [(9,[1,1])]] [[(3, [1, 1])]] [~1, 1] [(~2, [1, 0, 1]), (3, [1, 1, 1])] 1 ="

  1241 " ([(1, [1, 2, 1])], [2, ~3, 1], [[(5, [1, 1])], [(1, [1, 1])]]) changed";

  1244 "~~~~~ fun newton_mv, args:"; val (x,f,steps,t,p,order) = ([1, 2, 3, 4],[[(9, [0]), (5, [1])], [(16, [0]), (7, [1])]], [[(5, [0]), (2, [1])], [(1, [0])]], [2, ~3, 1], [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])], 0  );

  1245 length x = 2; (* false *)

  1246 val new_value_poly = multi_to_uni((uni_to_multi t)  %%*%% (uni_to_multi [(nth x (length x -2) )* ~1, 1]));

  1247 val new_steps = [((nth f (length f -1)) %%/ ((nth x (length x - 1)) - (nth x (length x - 2)))) %%-%% ((nth f (length f -2)))];

  1249 "~~~~~ fun next_step, args:"; val (steps,new_steps, x') = (steps, new_steps, x);

  1250 steps = []; (*false*)

  1252 "~~~~~ fun next_step, args:"; val (steps,new_steps, x') = (nth_drop 0 steps, new_steps @ [(((nth new_steps (length new_steps -1)) %%-%%(nth steps 0))) %%/

  1253                                   ((nth x' (length x' - 1)) - (nth x' (length x' - 3)))], nth_drop (length x' -2) x');steps = []; (*false*)

  1255 "~~~~~ fun next_step, args:"; val (steps,new_steps, x') = (nth_drop 0 steps, new_steps @ [(((nth new_steps (length new_steps -1)) %%-%%(nth steps 0))) %%/

  1256                                   ((nth x' (length x' - 1)) - (nth x' (length x' - 3)))], nth_drop (length x' -2) x');

  1257 steps = []; (*true*)

  1258 val steps = new_steps;

  1259 val polynom' =  p %%+%% (mult_with_new_var (nth steps (length steps -1)) new_value_poly order);

  1261 "----------- fun listgreater  ---------------------------";

  1262 "----------- fun listgreater  ---------------------------";

  1263 "----------- fun listgreater  ---------------------------";

  1264 if listgreater  [1,2,3,4,5] [1,2,3,4,5] = true then ()

  1265 else error " listgreater  [1,2,3,4,5] [1,2,3,4,5] = true changed";

  1266 if listgreater [1,2,3,4] [1,2,3,5] = false then ()

  1267 else error "listgreater [1,2,3,4] [1,2,3,5] = false changed"  ;

  1268 if listgreater [1,4,5,4] [0,3,4,5] = false then ()

  1269 else error "listgreater [1,2,3,4] [0,3,4,5] = false changed ";

  1271 "----------- fun greater_deg_multipoly  -----------------";

  1272 "----------- fun greater_deg_multipoly  -----------------";

  1273 "----------- fun greater_deg_multipoly  -----------------";

  1274 if greater_deg_multipoly  [1,2,3,4,5] [1,2,3,4,5] = 0 then ()

  1275 else error " greater_deg_multipoly  [1,2,3,4,5] [1,2,3,4,5] = 0 changed";

  1276 if greater_deg_multipoly [1,2,3,4] [5,2,8,1] = 1 then ()

  1277 else error "greater_deg_multipoly [1,2,3,4] [1,2,3,5] = 1 changed"  ;

  1278 if greater_deg_multipoly [1,4,5,4] [0,3,4,5] = 2 then ()

  1279 else error "greater_deg_multipoly [1,2,3,4] [0,3,4,5] = 2 changed ";

  1281 "----------- finfix  %%|%%  ------------------------------";

  1282 "----------- finfix  %%|%%  ------------------------------";

  1283 "----------- finfix  %%|%%  ------------------------------";

  1284 if [(1, [0, 0]), (~1, [0, 1])] %%|%% [(1, [0, 0]), (~1, [0, 1])] then ()

  1285 else error "[(1, [0, 0]), (~1, [0, 1])] %%|%% [(1, [0, 0]), (~1, [0, 1])] = true changed";

  1286 if [(3,[1,0])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])] then ()

  1287 else error "[(3,[1,0])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])] = true changed";

  1288 if [(3,[2,1])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])]

  1289 then error "[(3,[2,1])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])] = false changed" else ();

  1291 "----------- fun GCD_MODm  ------------------------------";

  1292 "----------- fun GCD_MODm  ------------------------------";

  1293 "----------- fun GCD_MODm  ------------------------------";

  1295 if GCD_MODm

  1296   [(~3,[2,0]),(1,[5,0]),(3,[0,1]),(~6,[1,1]),(~1,[3,1]),(2,[4,1]),(1,[3,2]),(~1,[1,3]),(2,[2,3])]

  1297   [(2,[2,0]),(~2,[0,1]),(4,[1,1]),(~1,[3,1]),(1,[1,2]),(~1,[2,2]),(~1,[0,3]),(2,[1,3])] 2 1 0

  1298   = [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])] then () else error

  1299 "GCD_MODm [(~3,[2,0]),(1,[5,0]),(3,[0,1]),(~6,[1,1]),(~1,[3,1]),(2,[4,1]),(1,[3,2]),(~1,[1,3]),(2,[2,3])]"

  1300 "[(2,[2,0]),(~2,[0,1]),(4,[1,1]),(~1,[3,1]),(1,[1,2]),(~1,[2,2]),(~1,[0,3]),(2,[1,3])] 2 1 0 "

  1301 "= [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])] changed";

  1302 (* -xy +xy^2z+yz - 1*)(* xy +1*) (*=*) (*xy -1*)

  1303 if GCD_MODm [(~1,[0,0,0]),(1,[0,1,1]),(1,[1,2,1]),(~1,[1,1,0])] [(1,[0,0,0]),(1,[1,1,0])] 3 2 0

  1304 = [(1, [0, 0, 0]), (1, [1, 1, 0])] then () else error

  1305 "GCD_MODm [(~1,[0,0,0]),(1,[0,1,1]),(1,[1,2,1]),(~1,[1,1,0])] [(1,[0,0,0]),(1,[1,1,0])] 3 2 0 "

  1306 "= [(1, [0, 0, 0]), (1, [1, 1, 0])] changed";

  1308 " ~~~ 1 ~~~";

  1309 "~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ([(1,[1,2,1])], [(1,[1,2,1])], 3, 2, 0);

  1310 s = 0; (*false*)

  1311 val M = 1 +  Int.min (max_deg_var a (length(get_varlist a 0)-2), max_deg_var b (length(get_varlist b 0)-2));

  1312  " ~~~ 1 ~~~";

  1313 "~~~~~ fun GOTO2 , args:"; val(a, b, n, s, M, r_list, steps ) =(a, b, n, s, M, [], []);

  1314 val r = find_new_r a b r;

  1315 val r_list = r_list @ [r];

  1316 val (a',b',n',s',r',r_list',steps',M') = ( a,b,n,s,r,r_list,steps,M);

  1317 (*g_1=*)

  1318  " ~~~ 1.1 ~~~";

  1319 "~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ((order_multipoly (evaluate_mv a (s-1) r)), ( order_multipoly (evaluate_mv b (s-1) r)), (n-1), (s-1), 0);

  1320 s = 0; (*false*)

  1322 val M = 1 +  Int.min (max_deg_var a (length(get_varlist a 0)-2), max_deg_var b (length(get_varlist b 0)-2));

  1323  " ~~~ 1.1 ~~~";

  1324 "~~~~~ fun GOTO2 , args:"; val(a, b, n, s, M, r_list, steps ) =(a, b, n, s, M, [], []);

  1325 val r = find_new_r a b r;

  1326 val r_list = r_list @ [r];

  1327 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1328 (*g_2=*)

  1329  " ~~~ 1.1.1 ~~~";

  1330 "~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ((order_multipoly (evaluate_mv a (s-1) r)), ( order_multipoly (evaluate_mv b (s-1) r)), (n-1), (s-1), 0);

  1331 s = 0; (*true*)

  1332 val g_2 =uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))))

  1333 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1334 val g_r = g_2;

  1335 val steps = steps @ [g_r];

  1336  " ~~~ 1.1 ~~~";

  1337 "~~~~~ fun GOTO3 , args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =

  1338 (a, b, n, s, M, g_r, r, r_list, steps);

  1339 val m = 1;

  1340  " ~~~ 1.1.W1 ~~~";

  1341 "~~~~~ fun WHILE , args:"; val (a ,b ,n ,s, M, m, r, r_list, steps, g, g_n ,mult) =

  1342 (a, b, n, s, M, m,r, r_list, steps, g_r, ( cero_multipoly (s+1)), [1]);

  1343 ; (*false*)

  1344 val r = find_new_r a b r;

  1345 val r_list = r_list @ [r];

  1346 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1347 (*g_3=*)

  1348  " ~~~ 1.1.W1.1 ~~~";

  1349 "~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ((order_multipoly (evaluate_mv a (s-1) r)), ( order_multipoly (evaluate_mv b (s-1) r)), (n-1), (s-1), 0);

  1350 s = 0; (*true*)

  1351 val g_3 =uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));

  1352  " ~~~ 1.1.W1 ~~~";

  1353 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1354 val g_r = g_3;

  1355 greater_var  (deg_multipoly g) (deg_multipoly g_r) (*false*);

  1356 (deg_multipoly g)= (deg_multipoly g_r) (*true*);

  1357 val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s-1);

  1358 if g_n = [(1,[1,1])] then () else error "g_n is false" ;

  1359 (nth steps (length steps -1)) = (cero_multipoly s) (* false*);

  1360  " ~~~ 1.1.W2 ~~~";

  1361 "~~~~~ fun WHILE , args:"; val (a ,b ,n ,s, M, m, r, r_list, steps, g, g_n ,mult) =

  1362 (a, b, n, s, M, (m+1), r, r_list, steps, g_r, g_n, new);

  1363 m > M; (*false*)

  1364 val r = find_new_r a b r;

  1365 val r_list = r_list @ [r];

  1366 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1367 (*g_4=*)

  1368 " ~~~ 1.1.W2.1 ~~~";

  1369 "~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ((order_multipoly (evaluate_mv a (s-1) r)), ( order_multipoly (evaluate_mv b (s-1) r)), (n-1), (s-1), 0);

  1370 s = 0; (*true*)

  1371 val g_r =uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));

  1372 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1373 " ~~~ 1.1.W2 ~~~";

  1374 greater_var  (deg_multipoly g) (deg_multipoly g_r) (*false*);

  1375 (deg_multipoly g)= (deg_multipoly g_r) (*true*);

  1376 val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s-1);

  1377 if g_n = [(1,[1,1])] then () else error "g_n is false" ;

  1378 (nth steps (length steps -1)) = (cero_multipoly s) (* true*);

  1379  " ~~~ 1.1.W3 ~~~";

  1380 "~~~~~ fun WHILE , args:"; val (a ,b ,n ,s, M, m, r, r_list, steps, g, g_n ,mult) =

  1381 (a, b, n, s, M, (M+1), r, r_list, steps, g_r, g_n, new);

  1382 m > M; (*true*)

  1383 g_n %%|%% a andalso g_n %%|%% b(*true*);

  1384  " ~~~ 1 ~~~";

  1385 val ( a,b,n,s,r,r_list,steps,M) = (a',b',n',s',r',r_list',steps',M');

  1386 val g_1=g_n;

  1387 val g_r = g_1;

  1388 val steps = steps @ [g_r];

  1389  " ~~~ 1 ~~~";

  1390 "~~~~~ fun GOTO3, args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =

  1391 (a, b, n, s, M, g_r, r, r_list, steps);

  1392 val m = 1;

  1393  " ~~~ 1.W1 ~~~";

  1394 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1395 (a, b, n, s, M, m, r, r_list, steps ,g_r, ( cero_multipoly (s+1)), [1]);

  1396 val g' = g;

  1397 m > M; (*false*)

  1398 val r = find_new_r a b r;

  1399 val r_list = r_list @ [r];

  1400 val (a',b',n',s',r',r_list',steps',m') = ( a,b,n,s,r,r_list,steps,m);

  1401 (* g_2 = GCD_MODm *)

  1402  " ~~~ 1.W1.1 ~~~";

  1403 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1404 ((order_multipoly  (evaluate_mv a (s-1) r)),(order_multipoly  (evaluate_mv b (s-1) r)), (n-1), (s-1), 0);

  1405 s = 0; (*false*)

  1406 val M = 1 +  Int.min(max_deg_var a (length(get_varlist a 0)-2), max_deg_var b (length(get_varlist b 0)-2));

  1407  " ~~~ 1.W1.1 ~~~";

  1408 "~~~~~ fun GOTO2, args:"; val (a,b,n,s,M,r_list,steps) =

  1409 (a, b, n, s, M, [], []);

  1410 val r = find_new_r a b r;

  1411 val r_list = r_list @ [r];

  1412 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1413  " ~~~ 1.W1.1.1 ~~~";

  1414  (*g_r=*)

  1415 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1416 ((order_multipoly  (evaluate_mv a (s-1) r)),(order_multipoly  (evaluate_mv b (s-1) r)), (n-1), (s-1), 0);

  1417 s = 0; (*true*)

  1418 val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));

  1419  " ~~~ 1.W1.1 ~~~";

  1420 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1421 val steps = steps @ [g_r];

  1422  " ~~~ 1.W1.1 ~~~";

  1423 "~~~~~ fun GOTO3, args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =

  1424 (a, b, n, s, M, g_r, r, r_list, steps);

  1425 val m = 1;

  1426 val (g'',g_n'',g_r'', m'',new'') = (g,g_n,g_r,m,mult);

  1427  " ~~~ 1.W1.1.W1 ~~~";

  1428 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1429 (a, b, n, s, M, m, r, r_list, steps ,g_r, ( cero_multipoly (s+1)), [1]);

  1431 m > M; (*false*)

  1432 val r = find_new_r a b r;

  1433 val r_list = r_list @ [r];

  1434 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1435 (* g_r = GCD_MODm *)

  1436  " ~~~ 1.W1.1.W1.1 ~~~";

  1437  (*g_r=*)

  1438 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1439 ((order_multipoly  (evaluate_mv a (s-1) r)),(order_multipoly  (evaluate_mv b (s-1) r)), (n-1), (s-1), 0);

  1440 s = 0; (*true*)

  1441 val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));

  1442  " ~~~ 1.W1.1.W1 ~~~";

  1443 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1444 greater_var  (deg_multipoly g) (deg_multipoly g_r); (*false*)

  1445 (deg_multipoly g)= (deg_multipoly g_r); (*true*)

  1446 val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s-1);

  1447 (nth steps (length steps -1)) = (cero_multipoly s); (*false*)

  1448  " ~~~ 1.W1.1.W2 ~~~";

  1449 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1450 (a, b, n, s, M, m+1, r, r_list, steps ,g_r, g_n,new);

  1451 m > M; (*false*)

  1452 val r = find_new_r a b r;

  1453 val r_list = r_list @ [r];

  1454 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1455 (* g_r = GCD_MODm *)

  1456  " ~~~ 1.W1.1.W2.1 ~~~";

  1457  (*g_r=*)

  1458 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1459 ((order_multipoly  (evaluate_mv a (s-1) r)),(order_multipoly  (evaluate_mv b (s-1) r)), (n-1), (s-1), 0);

  1460 s = 0; (*true*)

  1461 val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));

  1462  " ~~~ 1.W1.1.W2 ~~~";

  1463 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1464 greater_var  (deg_multipoly g) (deg_multipoly g_r); (*false*)

  1465 (deg_multipoly g)= (deg_multipoly g_r); (*true*)

  1466 val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s-1);

  1467 (nth steps (length steps -1)) = (cero_multipoly s); (*true*)

  1468  " ~~~ 1.W1.1.W3 ~~~";

  1469 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1470 (a, b, n, s, M, M+1, r, r_list, steps ,g_r, g_n,new);

  1471 m>M; (*true*)

  1472 g_n %%|%% a andalso g_n %%|%% b; (*true*)

  1473  " ~~~ 1.W1.1 ~~~";

  1474 val ( a,b,n,s,r,r_list,g,M) = (a',b',n',s',r',r_list',g',M');

  1475 val g_r = g_n;

  1476 greater_var  (deg_multipoly g) (deg_multipoly g_r); (*false*)

  1477 (deg_multipoly g)= (deg_multipoly g_r); (*true*)

  1478 val (g_n,steps,m,new) = (g_n'',steps',m'',new'');

  1479  " ~~~ 1.W2 ~~~";

  1480 val (g_n, new, steps) = newton_mv r_list [g, g_r] steps new g_n (s-1);

  1481 (nth steps (length steps -1)) = (cero_multipoly s); (*false*)

  1482 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1483 (a, b, n, s, M, m+1, r, r_list, steps ,g_r, g_n,new);

  1484 m > M; (*false*)

  1485 val r = find_new_r a b r;

  1486 val r_list = r_list @ [r];

  1487 val (a',b',n',s',r',r_list',steps',m',M',g') = ( a,b,n,s,r,r_list,steps,m,M,g);

  1488 (* g_2 = GCD_MODm *)

  1489  " ~~~ 1.W2.1 ~~~";

  1490 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1491 ((order_multipoly  (evaluate_mv a (s-1) r)),(order_multipoly  (evaluate_mv b (s-1) r)), (n-1), (s-1), 0);

  1492 s = 0; (*false*)

  1493 val M = 1 +  Int.min(max_deg_var a (length(get_varlist a 0)-2), max_deg_var b (length(get_varlist b 0)-2));

  1494  " ~~~ 1.W2.1 ~~~";

  1495 "~~~~~ fun GOTO2, args:"; val (a,b,n,s,M,r_list,steps) =

  1496 (a, b, n, s, M, [], []);

  1497 val r = find_new_r a b r;

  1498 val r_list = r_list @ [r];

  1499 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1500  " ~~~ 1.W2.1.1 ~~~";

  1501  (*g_r=*)

  1502 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1503 ((order_multipoly  (evaluate_mv a (s-1) r)),(order_multipoly  (evaluate_mv b (s-1) r)), (n-1), (s-1), 0);

  1504 s = 0; (*true*)

  1505 val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));

  1506  " ~~~ 1.W2.1 ~~~";

  1507 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1508 val steps = steps @ [g_r];

  1509  " ~~~ 1.W2.1 ~~~";

  1510 "~~~~~ fun GOTO3, args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =

  1511 (a, b, n, s, M, g_r, r, r_list, steps);

  1512 val m = 1;

  1513 val (g'',g_n'',g_r'', m'',new'') = (g,g_n,g_r,m,mult);

  1514  " ~~~ 1.W2.1.W1 ~~~";

  1515 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1516 (a, b, n, s, M, m, r, r_list, steps ,g_r, ( cero_multipoly (s+1)), [1]);

  1517 m > M; (*false*)

  1518 val r = find_new_r a b r;

  1519 val r_list = r_list @ [r];

  1520 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1521 (* g_r = GCD_MODm *)

  1522  " ~~~ 1.W2.1.W1.1 ~~~";

  1523  (*g_r=*)

  1524 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1525 ((order_multipoly  (evaluate_mv a (s-1) r)),(order_multipoly  (evaluate_mv b (s-1) r)), (n-1), (s-1), 0);

  1526 s = 0; (*true*)

  1527 val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));

  1528  " ~~~ 1.W2.1.W1 ~~~";

  1529 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1530 greater_var  (deg_multipoly g) (deg_multipoly g_r); (*false*)

  1531 (deg_multipoly g)= (deg_multipoly g_r); (*true*)

  1532 val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s-1);

  1533 (nth steps (length steps -1)) = (cero_multipoly s); (*false*)

  1534  " ~~~ 1.W2.1.W2 ~~~";

  1535 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1536 (a, b, n, s, M, m+1, r, r_list, steps ,g_r, g_n,new);

  1537 m > M; (*false*)

  1538 val r = find_new_r a b r;

  1539 val r_list = r_list @ [r];

  1540 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1541 (* g_r = GCD_MODm *)

  1542  " ~~~ 1.W2.1.W2.1 ~~~";

  1543  (*g_r=*)

  1544 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1545 ((order_multipoly  (evaluate_mv a (s-1) r)),(order_multipoly  (evaluate_mv b (s-1) r)), (n-1), (s-1), 0);

  1546 s = 0; (*true*)

  1547 val g_r = uni_to_multi((GCD_MOD (pp(multi_to_uni a)) (pp(multi_to_uni b))) %* (abs (Integer.gcd (cont (multi_to_uni a)) (cont( multi_to_uni b)))));

  1548  " ~~~ 1.W2.1.W2 ~~~";

  1549 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1550 greater_var  (deg_multipoly g) (deg_multipoly g_r); (*false*)

  1551 (deg_multipoly g)= (deg_multipoly g_r); (*true*)

  1552 val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s-1);

  1553 (nth steps (length steps -1)) = (cero_multipoly s); (*true*)

  1554  " ~~~ 1.W2.1.W3 ~~~";

  1555 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1556 (a, b, n, s, M, M+1, r, r_list, steps ,g_r, g_n,new);

  1557 m>M; (*true*)

  1558 g_n %%|%% a andalso g_n %%|%% b; (*true*)

  1559  " ~~~ 1.W2.1 ~~~";

  1560 val ( a,b,n,s,r,r_list,g) = (a',b',n',s',r',r_list',g');

  1561 val g_r = g_n;

  1562 greater_var  (deg_multipoly g) (deg_multipoly g_r); (*false*)

  1563 (deg_multipoly g)= (deg_multipoly g_r); (*true*)

  1564 val (g_n,steps,m,mult) = (g_n'',steps',m'',new'');

  1565  " ~~~ 1.W2 ~~~";

  1566  val (g_n, new, steps) = newton_mv r_list [g, g_r] steps mult g_n (s-1);

  1567 "~~~~~ to  return val:"; val (g) = (g_n);

  1571 if g =  [(1, [1, 2, 1])] then () else error "GCD_MODm [(1, [1, 2, 1])] [(1, [1, 2, 1])] has changes.";

  1573 " ==========================  END  ========================== ";

  1574   fun nth _ []      = error "nth _ []" (*Isabelle2002, still saved the hours of update*)

  1575     | nth 1 (x::_) = x

  1576     | nth n (_::xs) = nth (n-1) xs;

  1577 (*fun nth xs i = List.nth (xs, i);       recent Isabelle: TODO update all isac code   *)

author	Walther Neuper <neuper@ist.tugraz.at>
	Mon, 02 Sep 2013 16:16:08 +0200
changeset 52101	c3f399ce32af
parent 52074	044419cc7bf8
child 59111	c730b643bc0e
permissions	-rw-r--r--