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

     1 (* Title: src/../gcd_poly_old.sml

     2    Author: Diana Meindl

     3    Copyright (c) Diana Meindl 2011

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

     6 The tests below remain according to [1],

     7 while "$ISABELLE_ISAC_TEST/Tools/isac/Knowledge/gcd_poly_ml.sml" start exactly from the same state

     8 and in time follows development in "$ISABELLE_ISAC/Knowledge/GCD_Poly_ML.thy".

    10 This file includes *** tests ***.

    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);

    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  (*default_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*);

   688 (******************************************** tests ********************************************)

   689 (******************************************** tests ********************************************)

   690 (******************************************** tests ********************************************)

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

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

   694 "----------- fun mod_inv--- -----------------------------";

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

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

   697 rinv < m; (*=true*)

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

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

   700 rinv < m;(*=true*)

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

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

   703 rinv < m;(*=true*)

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

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

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

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

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

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

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

   713 "----------- fun CRA_2 ----------------------------------";

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

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

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

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

   720 "----------- fun lc -------------------------------------";

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

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

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

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

   726 "----------- fun deg ------------------------------------";

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

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

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

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

   732 "----------- fun primeList ------------------------------";

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

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

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

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

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

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

   740 "----------- fun find_next_prime_not_divide -----------------";

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

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

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

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

   746 "----------- fun poly_mod -------------------------------";

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

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

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

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

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

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

   754 "----------- fun %* ------------------------------";

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

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

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

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

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

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

   762 "----------- fun %/ -------------------------------";

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

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

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

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

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

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

   770 "----------- fun %-% ------------------------";

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

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

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

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

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

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

   778 "----------- fun %+% ------------------------------";

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

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

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

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

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

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

   786 "----------- fun CRA_2_poly -----------------------------";

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

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

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

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

   792 "----------- fun mod_div --------------------------------";

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

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

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

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

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

   799 "----------- fun lc_of_unipoly_not_0 --------------------";

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

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

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

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

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

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

   807 "----------- fun %|% -------------------------------";

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

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

   810 val c = lc poly1 div lc b;

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

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

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

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

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

   816 val c = lc poly1 div lc b;

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

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

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

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

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

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

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

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

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

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

   829 "----------- fun mod_poly_gcd ------------------------";

   831 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]]

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

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

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

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

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

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

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

   840 val moda = a poly_mod m;

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

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

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

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

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

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

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

   848 val moda = a poly_mod m

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

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

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

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

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

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

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

   856 val moda = a poly_mod m

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

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

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

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

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

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

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

   864 val moda = a poly_mod m

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

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

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

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

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

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

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

   874 "----------- fun normirt_polynom ------------------------";

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

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

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

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

   880 "----------- fun poly_centr -----------------------------";

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

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

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

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

   886 "----------- fun mod_gcd_p -------------------------------";

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

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

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

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

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

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

   895 "----------- fun gcd_n ----------------------------";

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

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

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

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

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

   902 "----------- fun pp -------------------------------------";

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

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

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

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

   908 "----------- fun GCD_MOD --------------------------------";

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

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

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

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

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

   915   val p = find_next_prime_not_divide d p

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

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

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

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

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

   921     val P = p

   922     val g = g_p;

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

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

   925        P > M ;(*false*)

   926       val p = find_next_prime_not_divide d p

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

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

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

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

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

   932        P > M ;(*false*)

   933       val p = find_next_prime_not_divide d p

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

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

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

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

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

   939       val P = P*p;

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

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

   942        P > M ;(*false*)

   943       val p = find_next_prime_not_divide d p

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

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

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

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

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

   949       val P = p

   950       val g = g_p;

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

   952        P > M ;(*false*)

   953       val p = find_next_prime_not_divide d p

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

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

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

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

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

   959       val P = P*p;

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

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

   962        P > M ;(*false*)

   963       val p = find_next_prime_not_divide d p

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

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

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

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

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

   969       val P = P*p;

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

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

   972        P > M ;(*false*)

   973       val p = find_next_prime_not_divide d p

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

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

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

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

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

   979       val P = P*p;

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

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

   982        P > M ;(*true*)

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

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

   985    val g = pp g;

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

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

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

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

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

   992 "----------- fun GCD_MOD --------------------------------";

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

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

   996 "------------------------------------------------------------";

   997 "--------------------- Multivariate Case --------------------";

   998 "------------------------------------------------------------";

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

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

  1002 "----------- fun get_coef -------------------------------";

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

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

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

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

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

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

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

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

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

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

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

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

  1017 "----------- fun get_varlist ----------------------------";

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

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

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

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

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

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

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

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

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

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

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

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

  1031 "----------- fun get_coeflist ---------------------------";

  1032 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

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

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

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

  1037 "----------- fun add_var_to_multipoly -------------------";

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

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

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

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

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

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

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

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

  1047 "----------- fun cero_multipoly -------------------------";

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

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

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

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

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

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

  1055 "----------- fun short_mv -------------------------------";

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

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

  1058 length mvp =1(*false*);

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

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

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

  1062 length mvp =1(*false*);

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

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

  1065 length mvp =1(*true*);

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

  1068 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])]

  1069 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"

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

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

  1073 "----------- fun greater_var ----------------------------";

  1074 if greater_var [3] [4]

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

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

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

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

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

  1080 if greater_var [1,2] [1,2]

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

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

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

  1085 "----------- fun order_multipoly ------------------------";

  1086 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])]

  1087 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"

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

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

  1091 "----------- fun lc_multipoly ---------------------------";

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

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

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

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

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

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

  1099 "----------- fun deg_multipoly -------------------------";

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

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

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

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

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

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

  1109 "----------- fun max_deg_var ----------------------------";

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

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

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

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

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

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

  1117 "----------- infix %%/ ----------------------------------";

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

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

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

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

  1123 "----------- fun cont -----------------------------------";

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

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

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

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

  1129 "----------- fun cont_multipoly -------------------------";

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

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

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

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

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

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

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

  1138 "----------- fun evaluate_mv ----------------------------";

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

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

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

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

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

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

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

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

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

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

  1150 "----------- fun mutli_to_uni ---------------------------";

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

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

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

  1155 val mvp = order_multipoly mvp;

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

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

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

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

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

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

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

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

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

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

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

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

  1171 "----------- fun find_new_r  ----------------------------";

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

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

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

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

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

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

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

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

  1181 "----------- fun newton_mv  -----------------------------";

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

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

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

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

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

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

  1189 "----------- fun mult_with_new_var  ---------------------";

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

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

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

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

  1195 "----------- fun greater_deg_multipoly  -----------------";

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

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

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

  1200 "----------- infix %%+%%  -------------------------------";

  1201 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])]

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

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

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

  1206 "----------- infix %%-%%  -------------------------------";

  1207 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])]

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

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

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

  1212 "----------- infix %%*%%  -------------------------------";

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

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

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

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

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

  1218 then () else error

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

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

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

  1223 "----------- fun newton_first  --------------------------";

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

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

  1226 then () else error

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

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

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

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

  1232 "----------- fun newton_mv  -----------------------------";

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

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

  1235 then () else error

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

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

  1238 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 =

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

  1240 then () else error

  1241 "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 ="

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

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

  1246 length x = 2; (* false *)

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

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

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

  1251 steps = []; (*false*)

  1253 "~~~~~ 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))) %%/

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

  1256 "~~~~~ 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))) %%/

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

  1258 steps = []; (*true*)

  1259 val steps = new_steps;

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

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

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

  1264 "----------- fun listgreater  ---------------------------";

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

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

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

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

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

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

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

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

  1274 "----------- fun greater_deg_multipoly  -----------------";

  1275 if greater_deg_multipoly  [1,2,3,4,5] [1,2,3,4,5] = 0 then ()

  1276 else error " greater_deg_multipoly  [1,2,3,4,5] [1,2,3,4,5] = 0 changed";

  1277 if greater_deg_multipoly [1,2,3,4] [5,2,8,1] = 1 then ()

  1278 else error "greater_deg_multipoly [1,2,3,4] [1,2,3,5] = 1 changed"  ;

  1279 if greater_deg_multipoly [1,4,5,4] [0,3,4,5] = 2 then ()

  1280 else error "greater_deg_multipoly [1,2,3,4] [0,3,4,5] = 2 changed ";

  1282 "----------- finfix  %%|%%  ------------------------------";

  1283 "----------- finfix  %%|%%  ------------------------------";

  1284 "----------- finfix  %%|%%  ------------------------------";

  1285 if [(1, [0, 0]), (~1, [0, 1])] %%|%% [(1, [0, 0]), (~1, [0, 1])] then ()

  1286 else error "[(1, [0, 0]), (~1, [0, 1])] %%|%% [(1, [0, 0]), (~1, [0, 1])] = true changed";

  1287 if [(3,[1,0])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])] then ()

  1288 else error "[(3,[1,0])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])] = true changed";

  1289 if [(3,[2,1])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])]

  1290 then error "[(3,[2,1])] %%|%% [(9,[1,1]),(12,[2,1]),(~3,[1,2])] = false changed" else ();

  1292 "----------- fun GCD_MODm  ------------------------------";

  1293 "----------- fun GCD_MODm  ------------------------------";

  1294 "----------- fun GCD_MODm  ------------------------------";

  1296 if GCD_MODm

  1297   [(~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])]

  1298   [(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

  1299   = [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])] then () else error

  1300 "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])]"

  1301 "[(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 "

  1302 "= [(1, [2, 0]), (~1, [0, 1]), (2, [1, 1])] changed";

  1303 (* -xy +xy^2z+yz - 1*)(* xy +1*) (*=*) (*xy - 1*)

  1304 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

  1305 = [(1, [0, 0, 0]), (1, [1, 1, 0])] then () else error

  1306 "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 "

  1307 "= [(1, [0, 0, 0]), (1, [1, 1, 0])] changed";

  1309 " ~~~ 1 ~~~";

  1310 "~~~~~ fun GCD_MODm , args:"; val (a,b,n,s,r) = ([(1,[1,2,1])], [(1,[1,2,1])], 3, 2, 0);

  1311 s = 0; (*false*)

  1312 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));

  1313  " ~~~ 1 ~~~";

  1314 "~~~~~ fun GOTO2 , args:"; val(a, b, n, s, M, r_list, steps ) =(a, b, n, s, M, [], []);

  1315 val r = find_new_r a b r;

  1316 val r_list = r_list @ [r];

  1317 val (a',b',n',s',r',r_list',steps',M') = ( a,b,n,s,r,r_list,steps,M);

  1318 (*g_1=*)

  1319  " ~~~ 1.1 ~~~";

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

  1321 s = 0; (*false*)

  1323 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));

  1324  " ~~~ 1.1 ~~~";

  1325 "~~~~~ fun GOTO2 , args:"; val(a, b, n, s, M, r_list, steps ) =(a, b, n, s, M, [], []);

  1326 val r = find_new_r a b r;

  1327 val r_list = r_list @ [r];

  1328 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1329 (*g_2=*)

  1330  " ~~~ 1.1.1 ~~~";

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

  1332 s = 0; (*true*)

  1333 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)))))

  1334 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1335 val g_r = g_2;

  1336 val steps = steps @ [g_r];

  1337  " ~~~ 1.1 ~~~";

  1338 "~~~~~ fun GOTO3 , args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =

  1339 (a, b, n, s, M, g_r, r, r_list, steps);

  1340 val m = 1;

  1341  " ~~~ 1.1.W1 ~~~";

  1342 "~~~~~ fun WHILE , args:"; val (a ,b ,n ,s, M, m, r, r_list, steps, g, g_n ,mult) =

  1343 (a, b, n, s, M, m,r, r_list, steps, g_r, ( cero_multipoly (s+1)), [1]);

  1344 ; (*false*)

  1345 val r = find_new_r a b r;

  1346 val r_list = r_list @ [r];

  1347 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1348 (*g_3=*)

  1349  " ~~~ 1.1.W1.1 ~~~";

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

  1351 s = 0; (*true*)

  1352 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)))));

  1353  " ~~~ 1.1.W1 ~~~";

  1354 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1355 val g_r = g_3;

  1356 greater_var  (deg_multipoly g) (deg_multipoly g_r) (*false*);

  1357 (deg_multipoly g)= (deg_multipoly g_r) (*true*);

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

  1359 if g_n = [(1,[1,1])] then () else error "g_n is false" ;

  1360 (nth steps (length steps - 1)) = (cero_multipoly s) (* false*);

  1361  " ~~~ 1.1.W2 ~~~";

  1362 "~~~~~ fun WHILE , args:"; val (a ,b ,n ,s, M, m, r, r_list, steps, g, g_n ,mult) =

  1363 (a, b, n, s, M, (m+1), r, r_list, steps, g_r, g_n, new);

  1364 m > M; (*false*)

  1365 val r = find_new_r a b r;

  1366 val r_list = r_list @ [r];

  1367 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1368 (*g_4=*)

  1369 " ~~~ 1.1.W2.1 ~~~";

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

  1371 s = 0; (*true*)

  1372 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)))));

  1373 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1374 " ~~~ 1.1.W2 ~~~";

  1375 greater_var  (deg_multipoly g) (deg_multipoly g_r) (*false*);

  1376 (deg_multipoly g)= (deg_multipoly g_r) (*true*);

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

  1378 if g_n = [(1,[1,1])] then () else error "g_n is false" ;

  1379 (nth steps (length steps - 1)) = (cero_multipoly s) (* true*);

  1380  " ~~~ 1.1.W3 ~~~";

  1381 "~~~~~ fun WHILE , args:"; val (a ,b ,n ,s, M, m, r, r_list, steps, g, g_n ,mult) =

  1382 (a, b, n, s, M, (M+1), r, r_list, steps, g_r, g_n, new);

  1383 m > M; (*true*)

  1384 g_n %%|%% a andalso g_n %%|%% b(*true*);

  1385  " ~~~ 1 ~~~";

  1386 val ( a,b,n,s,r,r_list,steps,M) = (a',b',n',s',r',r_list',steps',M');

  1387 val g_1=g_n;

  1388 val g_r = g_1;

  1389 val steps = steps @ [g_r];

  1390  " ~~~ 1 ~~~";

  1391 "~~~~~ fun GOTO3, args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =

  1392 (a, b, n, s, M, g_r, r, r_list, steps);

  1393 val m = 1;

  1394  " ~~~ 1.W1 ~~~";

  1395 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1396 (a, b, n, s, M, m, r, r_list, steps ,g_r, ( cero_multipoly (s+1)), [1]);

  1397 val g' = g;

  1398 m > M; (*false*)

  1399 val r = find_new_r a b r;

  1400 val r_list = r_list @ [r];

  1401 val (a',b',n',s',r',r_list',steps',m') = ( a,b,n,s,r,r_list,steps,m);

  1402 (* g_2 = GCD_MODm *)

  1403  " ~~~ 1.W1.1 ~~~";

  1404 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1405 ((order_multipoly  (evaluate_mv a (s- 1) r)),(order_multipoly  (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);

  1406 s = 0; (*false*)

  1407 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));

  1408  " ~~~ 1.W1.1 ~~~";

  1409 "~~~~~ fun GOTO2, args:"; val (a,b,n,s,M,r_list,steps) =

  1410 (a, b, n, s, M, [], []);

  1411 val r = find_new_r a b r;

  1412 val r_list = r_list @ [r];

  1413 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1414  " ~~~ 1.W1.1.1 ~~~";

  1415  (*g_r=*)

  1416 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1417 ((order_multipoly  (evaluate_mv a (s- 1) r)),(order_multipoly  (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);

  1418 s = 0; (*true*)

  1419 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)))));

  1420  " ~~~ 1.W1.1 ~~~";

  1421 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1422 val steps = steps @ [g_r];

  1423  " ~~~ 1.W1.1 ~~~";

  1424 "~~~~~ fun GOTO3, args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =

  1425 (a, b, n, s, M, g_r, r, r_list, steps);

  1426 val m = 1;

  1427 val (g'',g_n'',g_r'', m'',new'') = (g,g_n,g_r,m,mult);

  1428  " ~~~ 1.W1.1.W1 ~~~";

  1429 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1430 (a, b, n, s, M, m, r, r_list, steps ,g_r, ( cero_multipoly (s+1)), [1]);

  1432 m > M; (*false*)

  1433 val r = find_new_r a b r;

  1434 val r_list = r_list @ [r];

  1435 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1436 (* g_r = GCD_MODm *)

  1437  " ~~~ 1.W1.1.W1.1 ~~~";

  1438  (*g_r=*)

  1439 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1440 ((order_multipoly  (evaluate_mv a (s- 1) r)),(order_multipoly  (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);

  1441 s = 0; (*true*)

  1442 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)))));

  1443  " ~~~ 1.W1.1.W1 ~~~";

  1444 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1445 greater_var  (deg_multipoly g) (deg_multipoly g_r); (*false*)

  1446 (deg_multipoly g)= (deg_multipoly g_r); (*true*)

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

  1448 (nth steps (length steps - 1)) = (cero_multipoly s); (*false*)

  1449  " ~~~ 1.W1.1.W2 ~~~";

  1450 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1451 (a, b, n, s, M, m+1, r, r_list, steps ,g_r, g_n,new);

  1452 m > M; (*false*)

  1453 val r = find_new_r a b r;

  1454 val r_list = r_list @ [r];

  1455 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1456 (* g_r = GCD_MODm *)

  1457  " ~~~ 1.W1.1.W2.1 ~~~";

  1458  (*g_r=*)

  1459 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1460 ((order_multipoly  (evaluate_mv a (s- 1) r)),(order_multipoly  (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);

  1461 s = 0; (*true*)

  1462 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)))));

  1463  " ~~~ 1.W1.1.W2 ~~~";

  1464 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1465 greater_var  (deg_multipoly g) (deg_multipoly g_r); (*false*)

  1466 (deg_multipoly g)= (deg_multipoly g_r); (*true*)

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

  1468 (nth steps (length steps - 1)) = (cero_multipoly s); (*true*)

  1469  " ~~~ 1.W1.1.W3 ~~~";

  1470 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1471 (a, b, n, s, M, M+1, r, r_list, steps ,g_r, g_n,new);

  1472 m>M; (*true*)

  1473 g_n %%|%% a andalso g_n %%|%% b; (*true*)

  1474  " ~~~ 1.W1.1 ~~~";

  1475 val ( a,b,n,s,r,r_list,g,M) = (a',b',n',s',r',r_list',g',M');

  1476 val g_r = g_n;

  1477 greater_var  (deg_multipoly g) (deg_multipoly g_r); (*false*)

  1478 (deg_multipoly g)= (deg_multipoly g_r); (*true*)

  1479 val (g_n,steps,m,new) = (g_n'',steps',m'',new'');

  1480  " ~~~ 1.W2 ~~~";

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

  1482 (nth steps (length steps - 1)) = (cero_multipoly s); (*false*)

  1483 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1484 (a, b, n, s, M, m+1, r, r_list, steps ,g_r, g_n,new);

  1485 m > M; (*false*)

  1486 val r = find_new_r a b r;

  1487 val r_list = r_list @ [r];

  1488 val (a',b',n',s',r',r_list',steps',m',M',g') = ( a,b,n,s,r,r_list,steps,m,M,g);

  1489 (* g_2 = GCD_MODm *)

  1490  " ~~~ 1.W2.1 ~~~";

  1491 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1492 ((order_multipoly  (evaluate_mv a (s- 1) r)),(order_multipoly  (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);

  1493 s = 0; (*false*)

  1494 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));

  1495  " ~~~ 1.W2.1 ~~~";

  1496 "~~~~~ fun GOTO2, args:"; val (a,b,n,s,M,r_list,steps) =

  1497 (a, b, n, s, M, [], []);

  1498 val r = find_new_r a b r;

  1499 val r_list = r_list @ [r];

  1500 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1501  " ~~~ 1.W2.1.1 ~~~";

  1502  (*g_r=*)

  1503 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1504 ((order_multipoly  (evaluate_mv a (s- 1) r)),(order_multipoly  (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);

  1505 s = 0; (*true*)

  1506 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)))));

  1507  " ~~~ 1.W2.1 ~~~";

  1508 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1509 val steps = steps @ [g_r];

  1510  " ~~~ 1.W2.1 ~~~";

  1511 "~~~~~ fun GOTO3, args:"; val (a, b, n, s, M, g_r, r, r_list, steps) =

  1512 (a, b, n, s, M, g_r, r, r_list, steps);

  1513 val m = 1;

  1514 val (g'',g_n'',g_r'', m'',new'') = (g,g_n,g_r,m,mult);

  1515  " ~~~ 1.W2.1.W1 ~~~";

  1516 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1517 (a, b, n, s, M, m, r, r_list, steps ,g_r, ( cero_multipoly (s+1)), [1]);

  1518 m > M; (*false*)

  1519 val r = find_new_r a b r;

  1520 val r_list = r_list @ [r];

  1521 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1522 (* g_r = GCD_MODm *)

  1523  " ~~~ 1.W2.1.W1.1 ~~~";

  1524  (*g_r=*)

  1525 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1526 ((order_multipoly  (evaluate_mv a (s- 1) r)),(order_multipoly  (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);

  1527 s = 0; (*true*)

  1528 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)))));

  1529  " ~~~ 1.W2.1.W1 ~~~";

  1530 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1531 greater_var  (deg_multipoly g) (deg_multipoly g_r); (*false*)

  1532 (deg_multipoly g)= (deg_multipoly g_r); (*true*)

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

  1534 (nth steps (length steps - 1)) = (cero_multipoly s); (*false*)

  1535  " ~~~ 1.W2.1.W2 ~~~";

  1536 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1537 (a, b, n, s, M, m+1, r, r_list, steps ,g_r, g_n,new);

  1538 m > M; (*false*)

  1539 val r = find_new_r a b r;

  1540 val r_list = r_list @ [r];

  1541 val (a'',b'',n'',s'',r'',r_list'') = ( a,b,n,s,r,r_list);

  1542 (* g_r = GCD_MODm *)

  1543  " ~~~ 1.W2.1.W2.1 ~~~";

  1544  (*g_r=*)

  1545 "~~~~~ fun GCD_MODm, args:"; val  (a,b,n,s,r) =

  1546 ((order_multipoly  (evaluate_mv a (s- 1) r)),(order_multipoly  (evaluate_mv b (s- 1) r)), (n- 1), (s- 1), 0);

  1547 s = 0; (*true*)

  1548 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)))));

  1549  " ~~~ 1.W2.1.W2 ~~~";

  1550 val ( a,b,n,s,r,r_list) = (a'',b'',n'',s'',r'',r_list'');

  1551 greater_var  (deg_multipoly g) (deg_multipoly g_r); (*false*)

  1552 (deg_multipoly g)= (deg_multipoly g_r); (*true*)

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

  1554 (nth steps (length steps - 1)) = (cero_multipoly s); (*true*)

  1555  " ~~~ 1.W2.1.W3 ~~~";

  1556 "~~~~~ fun WHILE, args:"; val (a, b, n, s, M, m, r, r_list, steps ,g, g_n, mult) =

  1557 (a, b, n, s, M, M+1, r, r_list, steps ,g_r, g_n,new);

  1558 m>M; (*true*)

  1559 g_n %%|%% a andalso g_n %%|%% b; (*true*)

  1560  " ~~~ 1.W2.1 ~~~";

  1561 val ( a,b,n,s,r,r_list,g) = (a',b',n',s',r',r_list',g');

  1562 val g_r = g_n;

  1563 greater_var  (deg_multipoly g) (deg_multipoly g_r); (*false*)

  1564 (deg_multipoly g)= (deg_multipoly g_r); (*true*)

  1565 val (g_n,steps,m,mult) = (g_n'',steps',m'',new'');

  1566  " ~~~ 1.W2 ~~~";

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

  1568 "~~~~~ to  return val:"; val (g) = (g_n);

  1572 if g =  [(1, [1, 2, 1])] then () else error "GCD_MODm [(1, [1, 2, 1])] [(1, [1, 2, 1])] has changes.";

  1574 " ==========================  END  ========================== ";

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

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

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

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

author	wneuper <Walther.Neuper@jku.at>
	Thu, 04 Aug 2022 12:48:37 +0200
changeset 60509	2e0b7ca391dc
parent 60354	716dd4a05cc8
child 60771	1b072aab8f4e
permissions	-rw-r--r--