(* rationals, i.e. fractions of multivariate polynomials over the real field

   author: (c) Diana Meindl

   Use is subject to license terms. 

1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890

        10        20        30        40        50        60        70        80         90      100

*)

theory Rational2 imports Complex_Main begin 

chapter {* The Isabelle Theory *}

text {*

  Below we reference \cite{winkler96} by page. From this book we adopt the 

  algorithms as closely as possible. Isabelle coding standards \cite{TODO isarimpl}

  require to change the identifiers; further changes concern the transition from

  Winkler's imperative pseudocode to efficient functional code.

  However, in test/../reational.sml there is an almost original version of the

  algorithms, including the identifiers; this version was the first one

  developed in this thesis, the version is kept for simple reference to 

  \cite{winkler96}.

*}

section {* Auxiliary Functions *}

ML {*

  type unipoly = int list;

*}

subsection {* modulo calculations etc. *}

ML {*

  (*The "inverse modulo" functions: "mod_inv 5 7 = 3" 

  because 5*3 mod 7 = 1, or more precisely 5^-1 mod 7 = 3.*)

  fun mod_inv r m =

    let

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

        if rinv < m

        then

          if r mod m = 1

          then rinv

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

        else 0

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

  fun leading_coeff (uvp: unipoly) =

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

    then nth uvp (length uvp - 1)

    else leading_coeff (nth_drop (length uvp - 1) uvp);

  fun deg (uvp:unipoly) = 

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

    then length uvp - 1

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

  (* gcd is: greatest common divisor *)

  fun gcd a b = if b = 0 then a else gcd b (a mod b);

  (* Chinese remainder: given r1, r2, m1, m2 

     find r such that r = r1 mod m1 and r = r2 mod m2 *)

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

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

*}

subsection {* primes *}

ML{*

  (* is n divisble by some number in primelist ? *)

  fun is_prime primelist n = 

    if length primelist > 0

    then 

      if (n mod (nth primelist 0)) = 0

      then false

      else is_prime (nth_drop 0 primelist) n

    else true

  (* sieve of erathostenes *)

  fun primeList p =

    let 

    fun make_primelist list last p =

      if (nth list (length list - 1)) < p 

      then

        if ( is_prime list (last + 2)) 

        then make_primelist (list @ [(last + 2)]) (last + 2) p

        else make_primelist  list (last + 2) p

      else list

    in

      if p < 3

      then [2]

      else make_primelist [2,3] 3 p

    end

*}

text {*

fun primeList' 2 =[2,3]|

 primeList' x= if x<2 then [2] else

  let 

  fun make_primelist list last x =

    if (nth list (length list -1)) < x 

    then

      if ( is_prime list (last + 2)) 

      then make_primelist (list @ [(last + 2)]) (last + 2) x

      else make_primelist  list (last + 2) x

    else list

  in

     make_primelist [2,3] 3 x end

*}

text {*

	To find now a new prime number not dividing a number d, we find the next prime number higher than the old one with (primeList old+1) and test again if it divide our d. If not we found a new prime not dividing d, else we repeat it with the next prime.\\

*}

ML {*

fun find_next_prime_not_divide a p  = 

  let

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

  in

    if a mod next  <> 0

      then next

    else find_next_prime_not_divide a next

end

*}

text{*

find_next_prime_not_divide 12 2;

find_next_prime_not_divide 15 2;

find_next_prime_not_divide 90 7;

find_next_prime_not_divide 90 1;

find_next_prime_not_divide 2 7;

primeList 15;

primeList 1;

primeList 8;

primeList 2;

primeList 0;

primeList' 15;

primeList' 1;

primeList' 8;

primeList' 2;

*}

ML {*

is_prime (primeList 7) 9;

*}

ML {*

fun landau_mignotte_bound a b = 2603;

(* 2^min(m,n)gcd(leading_coeff(a),leading_coeff(b)) min(1/(abs leading_coeff(a)) sqrt(sum(a_^2), 1/(abs leading_coeff(b)) sqrt(sum(b_i^2))  *)

*}

ML {*

infix poly_mod

fun uvp poly_mod m =

  let fun poly_mod' uvp m n =

    if n<length (uvp)

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

    else uvp (*end of poly_mod'*)

  in

    poly_mod' uvp m 0 end (*poly_mod*)

*}                 

ML {*

infix %*

fun (a:unipoly) %* (b:int) =

  let fun mul a b = a*b;

  in map2 mul a (replicate (length(a)) b)end

infix %/

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

  let fun division poly b = poly div b;

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

*}                 

ML {*

infix %-%

fun a %-% b =

  let fun minus a b = a-b;

  in  map2 minus a b end

infix %-

fun p %- b =

  let fun minus p b = p-b;

in map2 minus p (replicate (length(p)) b) end

infix %+

fun a %+ b =

  let fun plus a b = a+b;

  in  map2 plus a b end

*}

ML{*

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

let 

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

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

*}

ML {*(*a/b mod m*)

fun mod_div (a:int) (b:int) (m:int) =

  a*(mod_inv b m) mod m

fun  poly_mod_div (a:unipoly) (b:int)  m =

  let fun lm_div a b c m i =

    if length c = length a

    then c

    else lm_div a b (c @ [ mod_div (nth a i) b m ]) m (i+1)

  in  lm_div a b [] m 0 end

*}

ML {*

fun mod_mult a b m = (a * b) mod m;

fun mod_minus a b m = (a - b) mod m;

fun mod_plus a b m = (a + b) mod m;

*}

ML{*

fun lc_of_list_not_0 [] = [](* and delete leading_coeff=0*)

| lc_of_list_not_0 (uvp:unipoly) = 

  if nth uvp (length uvp -1) =0

  then lc_of_list_not_0 (nth_drop (length uvp-1) uvp)

  else 

    if nth uvp 0 = 0

    then  lc_of_list_not_0 (nth_drop 0 uvp)

    else uvp;

*}

ML {* (* if poly2|poly1 *)

infix %/%

fun [a:int] %/% [b:int] = 

  if b<a andalso a mod b = 0

  then true

  else false

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

  let 

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

  val c = leading_coeff poly1 div leading_coeff b;

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

  in 

    if rest = []

    then true

    else

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

      then rest %/% poly2 

      else false

  end

*}

ML {*

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

  let 

    val moda = upoly1 poly_mod m

    val modb = (replicate (length upoly1 - length(lc_of_list_not_0 upoly2)) 0)

                @ (lc_of_list_not_0 upoly2 poly_mod m) ;

    val c =  mod_div (leading_coeff moda) (leading_coeff modb) m

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

  in 

    if rest = []

    then [upoly1, [0] ]

    else

      if length rest < length upoly2

      then mod_poly_gcd (upoly2, rest, m) (*[ rest,[c]@d]*)

      else mod_poly_gcd (rest, upoly2, m)

  end

*}

ML {*

fun  mod_polynom_gcd (poly1:unipoly) (poly2:unipoly) (m:int) =

  let 

    val moda = poly1 poly_mod m;

    val modb = (replicate (length poly1 - length poly2) 0) @ (poly2 poly_mod m) ;

    val c =  mod_div (leading_coeff moda) (leading_coeff modb) m;

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

  in (* nth_drop ( length ( upoly1 ) - 1 )*)

    if rest = []

    then [poly2, [0]]

    else

      if length rest < length poly2

      then [rest]

      else mod_poly_gcd (rest, poly2 ,  m)

  end

*}

ML{*

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

  let

    val poly1 = lc_of_list_not_0 poly1

    val lc_a=leading_coeff poly1;

    fun normirt poly1 b m lc_a i =

      if i=0

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

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

  in 

    if length(poly1)=0

    then poly1

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

end

*}

ML{*

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

  let

    val rest = mod_polynom_gcd poly1 poly2 p;

  in

    if length rest = 2

    then normirt_polynom (nth rest 0) p

    else   mod_gcd_p  poly2 (nth rest 0) p

end

*}

ML{*

fun gcd_n (list:int list) =

  if length list = 2

  then gcd (nth list 0)(nth list 1)

  else gcd_n ([gcd (nth list 0)(nth list 1)]@(nth_drop 0 (nth_drop 0 list)))

*}

ML {*

fun pp (poly1:unipoly) = 

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

  then  poly1 %/ (gcd_n poly1)

  else poly1

*}

ML {* 

pp [4,8,2,6];

*}

ML {* 

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

  let

    fun midle m =

      let val mid = m div 2

      in 

        if m mod mid = 0

        then  mid

        else  mid+1

      end

    fun centr a m mid =

      if a > mid

      then a - m

      else a

    fun polyCentr poly poly' m mid counter =

      if length(poly) > counter

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

      else (poly':unipoly)

  in polyCentr poly [] m (midle m) 0

end

*}

subsection {* GCD_MOD Algorgithmus *}

ML {*

 fun GCD_MOD a b =

   let

(*1*)val d = gcd (leading_coeff a) (leading_coeff b);

     val M  = 2*d*landau_mignotte_bound a b;

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

       let 

(*2*)    val p = find_next_prime_not_divide d p

         val c_p = normirt_polynom ( mod_gcd_p a b p) p

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

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

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

           then  [1]

           else

             let

               val P = p

               val g = g_p

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

(*4*)            if P > M 

                 then g

                 else

                   let 

                     val p = find_next_prime_not_divide d p

                     val c_p = normirt_polynom ( mod_gcd_p a b p) p

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

                   in

                     if deg g_p < deg g

                     then GOTO3 a b d M p g_p

                     else

                       if deg g_p = deg g

                       then 

                         let 

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

                           val P = P*p

                           val g = poly_centr(g poly_mod P)P

                         in WHILE a b d M P g end

                       else WHILE a b d M P g

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

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

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

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

(*5*)     val g = pp g

        in

          if (a %/% g) andalso (b %/% g)

          then  g

          else GOTO2 a b d M p

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

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

   in g end;

*}

ML{*

 fun GCD_MOD' a b =

   let

(*1*)val d = gcd (leading_coeff a) (leading_coeff b);

     val M  = 2*d*landau_mignotte_bound a b;

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

       let 

(*2*)    val p = find_next_prime_not_divide d p

         val c_p = normirt_polynom ( mod_gcd_p a b p) p

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

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

(*4*)      if P > M 

           then 

             if a %/% (pp g) andalso b %/% (pp g)

              then pp g

              else GOTO2 a b d M p

           else

             let 

               val p = find_next_prime_not_divide d p

               val c_p = normirt_polynom (mod_gcd_p a b p) p

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

             in

              if deg g_p < deg g

               then 

(*3*)            if deg g_p = 0

                 then [1]

                 else WHILE a b d M p p g_p

               else

                 if deg g_p = deg g

                 then 

                   WHILE a b d M (P*p) p ( poly_centr(CRA_2_poly (P,p)(g,g_p))(P*p)) 

                 else WHILE a b d M P p g

         end (*--------------WHILE--------------------*)

       in 

(*3*)      if deg g_p = 0 then [1]

           else WHILE a b d M p p g_p

    end (*========================GOTO2======================*)

in  GOTO2 a b d (*p*)1 M end

*}

subsection {* GCD_MODm Algorgithmus *}

ML {*

type monom = (int * int list);

type multipoly = monom list;

val mpoly1 = [(~3,[0,2]),(3,[1,0]),(1,[0,5]),(~6,[1,1]),(~1,[1,3]),(2,[1,4]),(1,[2,3]),(~1,[3,1]),(2,[3,2])];

val mpoly2 = [(2,[3,1]),(~1,[3,0]),(~1,[2,2]),(1,[2,1]),(~1,[1,3]),(4,[1,1]),(~2,[1,0]),(2,[0,2])];

val (test:monom) =(22,[5,6,7]);

*}

ML {*

fun get_coef ((coef,var):monom) = coef;

fun get_varlist ((coef,var):monom) = var;

fun get_coeflist (poly:multipoly) = 

  let 

    fun get_coeflist' (poly:multipoly) list = 

      if poly = []

       then 

        list

       else          

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

  in 

    get_coeflist' poly []

end

(* Eine Funktion für verschiedene Typen? *)

fun lc_mv (poly:multipoly) = 

  if get_coef (nth mpoly1 (length mpoly1 - 1)) <> 0

  then get_coef (nth poly (length poly-1))

  else lc_mv (nth_drop (length poly-1) poly);

fun deg_mv (mvp:multipoly) = 

  if get_coef (nth mvp (length mvp-1)) <> 0

  then nth (get_varlist (nth mvp (length mvp-1))) 0

  else  deg_mv (nth_drop (length mvp-1) mvp)

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

  if  deg_mv [monom1] >  deg_mv [monom2]

    then 1

    else if  deg_mv [monom1] <  deg_mv [monom2]

      then 2

      else if length (get_varlist monom1) >0

        then

          greater_deg (1,(nth_drop 0 (get_varlist monom1))) (1,(nth_drop 0 (get_varlist monom2)))

        else 0

*}

ML {*

get_coef test;

get_varlist test;

val t = [[1,2,4,5],[4,5,6,7]];

nth (nth t 1) 1;

lc_mv mpoly1;

 deg_mv mpoly1;

 deg_mv [(5,[4,5,6])];

get_coeflist mpoly1;

val monom =(5,[4,5,6]) ;

 deg_mv [monom];

get_varlist monom;

 [(1,(nth_drop 0 (get_varlist monom)))];

greater_deg (5,[4,5,6]) (5,[4,5,5]);

*}

ML {*

infix %%/

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

  let fun division monom b = ((get_coef monom div b),get_varlist monom);

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

*}

ML {*

mpoly1;

mpoly1 %%/ 2;

~11 div 2;

11 div 2;

[3,4,~4,~3,~1] %/ 2;

*}

ML {*

fun cont (poly1:unipoly) = gcd_n poly1

fun cont' (poly:multipoly) = gcd_n (get_coeflist poly)

*}

ML {*

[3,4,5] = [3,4,5];

[3,4,5]=[5,4,3];

[3,4,5] @ [get_coef (nth mpoly1 (length mpoly1 - 1))

                                - get_coef (nth mpoly2 (length mpoly2 - 1))];

mpoly1 = [] orelse mpoly2 =[];

get_varlist (nth mpoly1 (length mpoly1 - 1)) = get_varlist (nth mpoly2 (length mpoly2 - 1))

*}

ML {*

infix %%-

fun (mpoly1:multipoly) %%- (mpoly2:multipoly) = 

  let 

    fun min' (mpoly1:multipoly) (mpoly2:multipoly) (result:multipoly) = 

      if mpoly1 = [] andalso mpoly2 =[] 

        then result

        else if mpoly1 = []

          then let

            val result = [(0 - get_coef (nth mpoly2 (length mpoly2 - 1)),

                               get_varlist (nth mpoly2 (length mpoly2 - 1)))] @ result

            in 

               min'  mpoly1 (nth_drop (length mpoly2 - 1) mpoly2) result

            end

          else if mpoly2 = []

            then 

               mpoly1 @ result

       else 

         if get_varlist (nth mpoly1 (length mpoly1 - 1)) = get_varlist (nth mpoly2 (length mpoly2 - 1))

           then

            let

             val result =  [(get_coef (nth mpoly1 (length mpoly1 - 1))

                             - get_coef (nth mpoly2 (length mpoly2 - 1)),

                           get_varlist (nth mpoly1 (length mpoly1 - 1)))] @ result

            in 

              min' (nth_drop (length mpoly1 - 1) mpoly1) (nth_drop (length mpoly2 - 1) mpoly2) result

            end

         else

           if greater_deg (nth mpoly1 (length mpoly1 - 1)) (nth mpoly2 (length mpoly2 - 1)) = 1 

             then 

               min' (nth_drop (length mpoly1 - 1) mpoly1) mpoly2 ( [nth mpoly1 (length mpoly1 - 1)] @ result)

             else 

               let

                 val result =  [(0 - get_coef (nth mpoly2 (length mpoly2 - 1)),

                               get_varlist (nth mpoly2 (length mpoly2 - 1)))]@ result

               in 

                  min'  mpoly1 (nth_drop (length mpoly2 - 1) mpoly2) result

            end

    in 

      min' mpoly1 mpoly2 []

end

*}

ML {*

[(~6,[1,1]),(~1,[1,3])]%%-[(~4,[1,1]),(3,[1,3])];

[(~6,[1,1]),(~1,[1,5])]%%-[(~4,[1,1]),(3,[2,3])];

*}

ML {*

*}

ML {*

*}

ML {*

*}

ML {*

*}

ML {*

*}

ML {*

*}

ML {*

*}

text {*

fun GCD_MODm a b n s =

(*0*)  if s = 0

    then g = (gcd (cont a) (cont b)) %* (GCD_MOD pp(a) pp(b))

    else 

      let

(*1*)   val M = 1 + Int.min (DEG_X a, DEG_X b);

(*2*)     fun GOTO2 a b n s x M =

              let

                val r = AN_Integer_DEG_H;

                val g_r = GCD_MODm (REDUCED a)(REDUCED b) n (s-1)

(*3*)           fun GOTO3 a b n s x M g_r r=

                  let

                    val m = 1;

                    val g = g_r

                    fun WHILE a b n s x M m r =

                      if m > M

                        then 

                         if g TEILT a andalso g TEILT a

                          then  g

                          else GOTO2 a b n s x M

                        else

                        let

                        val r = NEW_INTEGER;

                        val g_r = GCD_MODm (REDUCED a)(REDUCED b) n (s-1);

                        in

                          if DEG_XN g_r < DEG_XN g

                            then GOTO3 a b n s x M g_r r

                            else if DEG_XN g_r < DEG_XN g

                              then INCORPORATE g_r

                              else WHILE a b n s x M (m+1) r

                        end (* WHILE *)

                  in (* GOTO3*)

                    WHILE a b n s x M m r

                  end (*GOTO3*)

              in (* GOTO2*)

                GOTO3 a b n s x M g_r r

              end (*GOTO2*)

            in

             GOTO2 a b n s x M 

          end (*end 0*)

*}

ML {*

  val a = [~18, ~15, ~20, 12, 20, ~13, 2];

  val b = [8, 28, 22, ~11, ~14, 1, 2];

  if GCD_MOD a b = [~2, ~1, 1] then ()

  else error "GCD_MOD a b = [~2, ~1, 1]";

*}

end