(* 1st approach to Isabelle's new Polynomials

   continued in TPPL/material/AD_Polynomal.thy *)

theory GCD_Poly_OLD 

  imports "HOL-Computational_Algebra.Polynomial"

        (*"~~/src/HOL/Number_Theory/Primes" ??? rm ???, thys for code_gen missing*)

      (*"~~/src/HOL/Library/Code_Target_Numeral"  ...note (4)*)

begin

value "xxx" \<comment> \<open>TODO\<close>

subsection \<open>Euclidean Algorithm with generalised division\<close>

text \<open>

  "generalised division" addresses the generation of "polynomial remainder sequences"

  according to [1].p.83. These sequences allow tho generalise the Euclidean Algorithm

  to univariate polynomials in Z[x], i.e. keeping integer coefficients; see [1].p.84.

  References: 

  [1] Franz Winkler, Polynomial Algorithms in Computer Algebra. Springer 1996\<close>

subsubsection \<open>Auxiliary functions\<close>

text \<open>

  Some of these functions might be inappropriate due to lack of knowledge 

  about Isabelle's Polynomials at this first approach.

\<close>

fun gcds :: "int list \<Rightarrow> int" where "gcds ns = List.fold gcd ns (List.hd ns)"

(*

lift_definition sdiv :: "'a::comm_semiring_1 \<Rightarrow> 'a poly \<Rightarrow> 'a poly"

  is "\<lambda>a p n. (coeff p n) div a"

proof sorry

Lifting failed for the following term:

Term:  \<lambda>a p n. coeff p n div a

Type:  ?'b \<Rightarrow> ?'b poly \<Rightarrow> nat \<Rightarrow> ?'b

Reason:  The type of the term cannot be instancied to "'a \<Rightarrow> 'a poly \<Rightarrow> nat \<Rightarrow> 'a".

*)

(*/------------------------------------------------------------------------------\

  this is an odd bypass for the above problem: 

  divide a polynomial with (a factor of) the gcd; i.e. remainder is expected 0 *)

fun pConss :: "int list \<Rightarrow> int poly \<Rightarrow> int poly" where

"pConss [] p = p" |

"pConss (x # xs) p = pConss xs (pCons x p)"

value "pConss (rev \<comment> \<open>!!!\<close> [1,2,3]) [:0:]" \<comment> \<open>Poly [1, 2, 3]\<close>

value "pConss [] [:1,2,3:]"                \<comment> \<open>Poly [1, 2, 3]\<close>

value "pConss [0,0] [:1,2,3:]"             \<comment> \<open>Poly [0, 0, 1, 2, 3]\<close>

fun sdiv :: "int poly \<Rightarrow> int \<Rightarrow> int poly" where

"sdiv p n =

(let cs = map (\<lambda>i. i div n) (coeffs p)

 in pConss (rev cs) [:0:])"

value "sdiv [:2,4,6::int:] 2" \<comment> \<open>Poly [1, 2, 3]\<close>

(*\------------------------------------------------------------------------------/*)

(* find the primitive part of a polynomial, see [1].p.83

   primitive p = p'

    assumes p = [:c\<^isub>1, .., c\<^isub>n:]

    yields p' = [:c\<^isub>1, .., c\<^isub>n:] = [:c\<^isub>1', .., c\<^isub>n':] * d

      \<and>  d = gcds c\<^isub>1, .., c\<^isub>n *)

definition primitive :: "int poly \<Rightarrow> int poly" where

"primitive p = sdiv p (gcds (coeffs p))"

value "primitive [:1::int, 2, 3:]"             \<comment> \<open>Poly [1, 2, 3]\<close>

value "primitive [:2::int, 4, 6:]"             \<comment> \<open>Poly [1, 2, 3]\<close>

value "primitive [:2000::int, 4000, 6000:]"    \<comment> \<open>Poly [1, 2, 3]\<close>

(* arguments from the example from [1].p.94 *)

value "       primitive [:-5000595353600, -2500297676800, 2500297676800::int:]" 

  \<comment> \<open>= Poly [-2, -1, 1]  OK\<close>

value (*[simp]*) "primitive [:-5000595353600, -2500297676800, 2500297676800::int:]" 

  \<comment> \<open>= [:-2, -1, 1, 0:]  incorrect ???????????????????????????????????????????????????????????\<close>

(* goes forever:

value (*[nbe]*)  "primitive [:-5000595353600, -2500297676800, 2500297676800::int:]" 

*)

value (*[code]*) "primitive [:-5000595353600, -2500297676800, 2500297676800::int:]" 

  \<comment> \<open>= Poly [-2, -1, 1] ...ok\<close>

(*where does [simp] fail:*)

value (*[simp]*) "let p = [:-5000595353600, -2500297676800, 2500297676800::int:] 

in (coeffs p)"      \<comment> \<open>= "[-5000595353600, -2500297676800, 2500297676800]" :: "int list" OK\<close>

value (*[simp]*) "let p = [:-5000595353600, -2500297676800, 2500297676800::int:] 

in gcds (coeffs p)"  \<comment> \<open>= "2500297676800" :: "int" OK\<close>

(* declare [[simp_trace = true]]

value [simp] "let p = [:-5000595353600, -2500297676800, 2500297676800::int:] 

in sdiv p (gcds (coeffs p))"  --{* = "[:-2, -1, 1, 0:]" :: "int poly" ???*}

... gives an unreadable trace, unfortunately.*)

(*where does  [nbe] fail:*)

value (*[nbe]*)  "let p = [:-5000595353600, -2500297676800, 2500297676800::int:] 

in (coeffs p)" \<comment> \<open>= "coeffs (Poly (-5000595353600 # coeffs (Poly (-2500297676800 # coef... ???\<close>

(* value [nbe]  "let p = [:-5000595353600, -2500297676800, 2500297676800::int:] 

in gcds (coeffs p)"

... goes forever, clear due to the above line, which is unclear.*)

(* regard zeros occurring in the quotient, also in the last step of division *)

fun for_quot_regarding :: "int poly \<Rightarrow> int poly \<Rightarrow> int poly \<Rightarrow> int poly \<Rightarrow> int poly \<Rightarrow> nat" where

"for_quot_regarding p1 p2 p1' quot remd =

(let

  zeros = degree p1' - degree remd - 1\<comment> \<open>length [q]\<close>;

  max_qot_deg = degree p1 - degree p2

in 

  if zeros + 1 \<comment> \<open>length [q]\<close> + degree quot > max_qot_deg \<comment> \<open>possible in last step of div\<close>

  then max_qot_deg - degree quot - 1 \<comment> \<open>length [q]\<close>

  else zeros)"

(* multiply polynomial p such that it has degree d with d = deg p*)

fun mult_to_deg :: "nat \<Rightarrow> int poly \<Rightarrow> int poly" where

"mult_to_deg d p =

(let

  d' = degree p

in 

  if d \<ge> d' then pConss (replicate (d - d') 0) p else p)"

value "mult_to_deg 4 [:1, 2, 3::int:]" \<comment> \<open>Poly [0, 0, 1, 2, 3]\<close>

(* find a factor for the dividend, such the divisor divides to an integer quotient.

value "xxx" 

   fact_quot c1 c2 = (n, q)

    assumes c2 \<noteq> 0

    yields c1 * n = q * c2 *)

fun fact_quot :: "int \<Rightarrow> int \<Rightarrow> int \<times> int" where

"fact_quot c1 c2 =

(let

  d = gcd c1 c2;

  n = c2 div d;

  q = c1 div d

 in if n > 0 then (n, q) else (-1 * n, -1 * q))"

value "fact_quot 6 4 = (2, 3)"

definition leading_coeff :: "'a poly \<Rightarrow> 'a::zero" where

"leading_coeff p = coeff p (degree p)"

value "leading_coeff [:1, 2, 3::int:] = 3"

subsubsection \<open>Generalised division\<close>

text \<open>

  The "generalised division" is used to create "polynomial remainder sequences"

  which in turn allow to generalise the Euclidean Algorithm to Z[x].

\<close>

(* "generalized division" for "generalised Euclidean algorithms:

   divide in Z[x] while multiplying the dividend such that the quotient is again in Z[x], 

   i.e. has integer coefficients.

   fact_div :: unipoly \<Rightarrow> unipoly \<Rightarrow> (unipoly, unipoly, nat)

   fact_div p\<^isub>1 p\<^isub>2 = (n, q, r)

     assumes p\<^isub>2 \<noteq> [0]  \<and>  deg p\<^isub>1 \<ge> deg p\<^isub>2

     yields p\<^isub>1 *\<^isub>s n = q *\<^isub>p p\<^isub>2   +\<^isub>p   r *)

(* TODO: generalise to 'a poly *)

function div_int_up :: "int poly \<Rightarrow> int poly \<Rightarrow> int poly \<Rightarrow> int poly \<Rightarrow> int \<Rightarrow> int \<times> int poly \<times> int poly" where

"div_int_up p1 p1' p2 quot ns = 

(let

  (n, q) = fact_quot (leading_coeff p1') (leading_coeff p2);

  remd = smult n p1' - smult q (mult_to_deg (degree p1') p2);

  zeros = for_quot_regarding p1 p2 p1' quot remd;

  quot = pConss (replicate zeros 0) (pCons q (smult n quot));

  ns = n * ns

in 

  if degree remd \<ge> degree p2 

  then div_int_up p1 remd p2 quot ns

  else (ns, quot, remd))"

by pat_completeness (auto simp del: mult_to_deg.simps)

termination sorry

(*TODO how can we retain the operator with nice subscripts ?*)

definition mult_div :: "int poly \<Rightarrow> int poly \<Rightarrow> int \<times> int poly \<times> int poly" (infixl "div\<^sub>*\<^sub>P" 70) where

"p1 div\<^sub>*\<^sub>P p2 = (div_int_up p1 p1 p2 0 1)"

subsubsection \<open>Euclidean Algorithm applied to univariate polynomials\<close>

text \<open>

  See [1].p.84.

\<close>

function euclidean_algorithm :: "int poly \<Rightarrow> int poly \<Rightarrow> int poly" where

"euclidean_algorithm a b =

  (if b  = [:0:] then primitive a else

     if degree a < degree b then euclidean_algorithm b a else

     let (n, q, r) = a div\<^sub>*\<^sub>P b

     in euclidean_algorithm b r)"

by pat_completeness auto

termination sorry

(*The example used in the inital mail to Tobias Nipkow made univariate: y -> 1*)

value "euclidean_algorithm [:0::int, -1, 1:]  [:-1, 0, 1:]"    \<comment> \<open>Poly [1, -1]\<close>

(*                                -x + x^2    -1 + x^2                

                     gcd     -x * (1 - x  )  (-1-x)*(1 - x)  =       (1 - x) *)

(* example from [1].p.94: DOESN'T WORK SINCE REMOVAL OF BUG IN function div_int_up 

  in changeset 044419cc7bf8 this example works (with [code]), if in

  https://intra.ist.tugraz.at/hg/isa/file/044419cc7bf8/src/Tools/isac/Knowledge/GCD_Poly.thy

  in line #30 "2" is replaced by "n".

*)

(* these go forever (on Intel\<registered> Core™ i5-2410M CPU @ 2.30GHz \<times> 4, 3.8 GiB)

value 

  "euclidean_algorithm [:-18, -15, -20, 12, 20, -13, 2::int:]  [:8, 28, 22, -11, -14, 1, 2:]"

value [simp]

  "euclidean_algorithm [:-18, -15, -20, 12, 20, -13, 2::int:]  [:8, 28, 22, -11, -14, 1, 2:]"*)

value (*[nbe]*)

  "euclidean_algorithm [:-18, -15, -20, 12, 20, -13, 2::int:]  [:8, 28, 22, -11, -14, 1, 2:]"

(* = "euclidean_algorithm (Poly (-18 # coeffs (Poly (-15 # coeffs ...why not simplified ???*)

value (*[code]*)

  "euclidean_algorithm [:-18, -15, -20, 12, 20, -13, 2::int:]  [:8, 28, 22, -11, -14, 1, 2:]"

(* = "Poly [-2, -1, 1]"*)

subsubsection \<open>Observations and conclusions\<close>

text \<open>Which part of the code crucially slows down with large numerals ?

  For answering this questions we investigate the last iteration in the

  example from [1].p.94. The test-data are obtained from the respective test

  "--- last step in EUCLID -\<comment> \<open> in gcd_poly.sml.\<close>

\<close>

(* we step through the function:

"div_int_up p1 p1' p2 quot ns = 

(let

##1##   (n, q) = fact_quot (leading_coeff p1') (leading_coeff p2);

##2##   remd = smult (int n) p1' - smult q (mult_to_deg (degree p1') p2);

##3##   zeros = for_quot_regarding p1 p2 p1' quot remd;

##4##   quot = pConss (replicate zeros 0) (pCons q (smult (int n) quot));

##5##   ns = n * ns

in 

##6##   if degree remd \<ge> degree p2 

##7##   then div_int_up p1 remd p2 quot ns

##8##   else (ns, quot, remd))"

*)

(*

##1##   (n, q) = fact_quot (leading_coeff p1') (leading_coeff p2);

*)

value "let 

  p1' = pConss (rev [-1316434, -3708859, -867104, 1525321]) [:0:] 

in leading_coeff p1'"

value "let 

  p1' = pConss (rev [-1316434, -3708859, -867104, 1525321]) [:0:];

  p2 = pConss (rev [-5000595353600, -2500297676800, 2500297676800]) [:0:] 

in fact_quot (leading_coeff p1') (leading_coeff p2)"

(* note: this example shows that "value" selects the appropriate evaluation strategy

   very reasonably:

HANGS: gcd, div

value "let c1 = 1525321; c2 = 2500297676800 in gcd c1 c2" \<comment> \<open> = 343\<close>

                                               =========

value "let c2 = 2500297676800; d = 343 in c2 div d"       \<comment> \<open> = 7289497600\<close>

                                          ========                        *)

(*

##2##   remd = smult (int n) p1' - smult q (mult_to_deg (degree p1') p2);

*)

value "let 

  n = 7289497600::int \<comment> \<open>nat\<close>;

  p1' = pConss (rev [-1316434, -3708859, -867104, 1525321]) [:0:] 

in smult (\<comment> \<open>int\<close> n) p1'" \<comment> \<open>the type conversion hangs !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\<close>

value "let 

  q = 111::int;

  p1' = pConss (rev [-1316434, -3708859, -867104, 1525321]) [:0:];

  p2 = pConss (rev [-5000595353600, -2500297676800, 2500297676800]) [:0:] 

in smult q (mult_to_deg (degree p1') p2)"

(* 

##3##   zeros = for_quot_regarding p1 p2 p1' quot remd;

*)

value "let

  p1 = pConss (rev [-1316434, -3708859, -867104, 1525321]) [:0:];

  p2 = pConss (rev [-5000595353600, -2500297676800, 2500297676800]) [:0:];

  p1' = pConss (rev [-1316434, -3708859, -867104, 1525321]) [:0:];

  quot = [:0::int:];

  remd = pConss (rev [-9596142483558400, -4798071241779200, 4798071241779200]) [:0:]

in for_quot_regarding p1 p2 p1' quot remd"

(* 

##4##   quot = pConss (replicate zeros 0) (pCons q (smult (int n) quot));

*)

value "let

  zeros = 0::nat;

  q = 4447::int;

  n = 7289497600::\<comment> \<open>nat\<close> int;

  quot = [:0::int:]

in pConss (replicate zeros 0) (pCons q (smult (\<comment> \<open>int\<close> n) quot))"

(* 

##5##   ns = n * ns

*)

value "let 

  n = 7289497600::int\<comment> \<open>nat\<close>; \<comment> \<open>nat even cannot be initialised !!!!!!!!!!!!!!!!!!!!!!!!!!!!!\<close>

  ns = 1::int

in n * ns" \<comment> \<open>\<close>

\<comment> \<open>-------------------------------------------------------------------------------------------\<close>

value "xxx" \<comment> \<open>"\<close>

end