src/HOL/Library/Polynomial.thy
changeset 65417 fc41a5650fb1
parent 65416 f707dbcf11e3
child 65418 c821f1f3d92d
child 65419 457e4fbed731
equal deleted inserted replaced
65416:f707dbcf11e3 65417:fc41a5650fb1
     1 (*  Title:      HOL/Library/Polynomial.thy
       
     2     Author:     Brian Huffman
       
     3     Author:     Clemens Ballarin
       
     4     Author:     Amine Chaieb
       
     5     Author:     Florian Haftmann
       
     6 *)
       
     7 
       
     8 section \<open>Polynomials as type over a ring structure\<close>
       
     9 
       
    10 theory Polynomial
       
    11   imports
       
    12     "~~/src/HOL/Deriv"
       
    13     More_List
       
    14     Infinite_Set
       
    15 begin
       
    16 
       
    17 subsection \<open>Auxiliary: operations for lists (later) representing coefficients\<close>
       
    18 
       
    19 definition cCons :: "'a::zero \<Rightarrow> 'a list \<Rightarrow> 'a list"  (infixr "##" 65)
       
    20   where "x ## xs = (if xs = [] \<and> x = 0 then [] else x # xs)"
       
    21 
       
    22 lemma cCons_0_Nil_eq [simp]: "0 ## [] = []"
       
    23   by (simp add: cCons_def)
       
    24 
       
    25 lemma cCons_Cons_eq [simp]: "x ## y # ys = x # y # ys"
       
    26   by (simp add: cCons_def)
       
    27 
       
    28 lemma cCons_append_Cons_eq [simp]: "x ## xs @ y # ys = x # xs @ y # ys"
       
    29   by (simp add: cCons_def)
       
    30 
       
    31 lemma cCons_not_0_eq [simp]: "x \<noteq> 0 \<Longrightarrow> x ## xs = x # xs"
       
    32   by (simp add: cCons_def)
       
    33 
       
    34 lemma strip_while_not_0_Cons_eq [simp]:
       
    35   "strip_while (\<lambda>x. x = 0) (x # xs) = x ## strip_while (\<lambda>x. x = 0) xs"
       
    36 proof (cases "x = 0")
       
    37   case False
       
    38   then show ?thesis by simp
       
    39 next
       
    40   case True
       
    41   show ?thesis
       
    42   proof (induct xs rule: rev_induct)
       
    43     case Nil
       
    44     with True show ?case by simp
       
    45   next
       
    46     case (snoc y ys)
       
    47     then show ?case
       
    48       by (cases "y = 0") (simp_all add: append_Cons [symmetric] del: append_Cons)
       
    49   qed
       
    50 qed
       
    51 
       
    52 lemma tl_cCons [simp]: "tl (x ## xs) = xs"
       
    53   by (simp add: cCons_def)
       
    54 
       
    55 
       
    56 subsection \<open>Definition of type \<open>poly\<close>\<close>
       
    57 
       
    58 typedef (overloaded) 'a poly = "{f :: nat \<Rightarrow> 'a::zero. \<forall>\<^sub>\<infinity> n. f n = 0}"
       
    59   morphisms coeff Abs_poly
       
    60   by (auto intro!: ALL_MOST)
       
    61 
       
    62 setup_lifting type_definition_poly
       
    63 
       
    64 lemma poly_eq_iff: "p = q \<longleftrightarrow> (\<forall>n. coeff p n = coeff q n)"
       
    65   by (simp add: coeff_inject [symmetric] fun_eq_iff)
       
    66 
       
    67 lemma poly_eqI: "(\<And>n. coeff p n = coeff q n) \<Longrightarrow> p = q"
       
    68   by (simp add: poly_eq_iff)
       
    69 
       
    70 lemma MOST_coeff_eq_0: "\<forall>\<^sub>\<infinity> n. coeff p n = 0"
       
    71   using coeff [of p] by simp
       
    72 
       
    73 
       
    74 subsection \<open>Degree of a polynomial\<close>
       
    75 
       
    76 definition degree :: "'a::zero poly \<Rightarrow> nat"
       
    77   where "degree p = (LEAST n. \<forall>i>n. coeff p i = 0)"
       
    78 
       
    79 lemma coeff_eq_0:
       
    80   assumes "degree p < n"
       
    81   shows "coeff p n = 0"
       
    82 proof -
       
    83   have "\<exists>n. \<forall>i>n. coeff p i = 0"
       
    84     using MOST_coeff_eq_0 by (simp add: MOST_nat)
       
    85   then have "\<forall>i>degree p. coeff p i = 0"
       
    86     unfolding degree_def by (rule LeastI_ex)
       
    87   with assms show ?thesis by simp
       
    88 qed
       
    89 
       
    90 lemma le_degree: "coeff p n \<noteq> 0 \<Longrightarrow> n \<le> degree p"
       
    91   by (erule contrapos_np, rule coeff_eq_0, simp)
       
    92 
       
    93 lemma degree_le: "\<forall>i>n. coeff p i = 0 \<Longrightarrow> degree p \<le> n"
       
    94   unfolding degree_def by (erule Least_le)
       
    95 
       
    96 lemma less_degree_imp: "n < degree p \<Longrightarrow> \<exists>i>n. coeff p i \<noteq> 0"
       
    97   unfolding degree_def by (drule not_less_Least, simp)
       
    98 
       
    99 
       
   100 subsection \<open>The zero polynomial\<close>
       
   101 
       
   102 instantiation poly :: (zero) zero
       
   103 begin
       
   104 
       
   105 lift_definition zero_poly :: "'a poly"
       
   106   is "\<lambda>_. 0"
       
   107   by (rule MOST_I) simp
       
   108 
       
   109 instance ..
       
   110 
       
   111 end
       
   112 
       
   113 lemma coeff_0 [simp]: "coeff 0 n = 0"
       
   114   by transfer rule
       
   115 
       
   116 lemma degree_0 [simp]: "degree 0 = 0"
       
   117   by (rule order_antisym [OF degree_le le0]) simp
       
   118 
       
   119 lemma leading_coeff_neq_0:
       
   120   assumes "p \<noteq> 0"
       
   121   shows "coeff p (degree p) \<noteq> 0"
       
   122 proof (cases "degree p")
       
   123   case 0
       
   124   from \<open>p \<noteq> 0\<close> obtain n where "coeff p n \<noteq> 0"
       
   125     by (auto simp add: poly_eq_iff)
       
   126   then have "n \<le> degree p"
       
   127     by (rule le_degree)
       
   128   with \<open>coeff p n \<noteq> 0\<close> and \<open>degree p = 0\<close> show "coeff p (degree p) \<noteq> 0"
       
   129     by simp
       
   130 next
       
   131   case (Suc n)
       
   132   from \<open>degree p = Suc n\<close> have "n < degree p"
       
   133     by simp
       
   134   then have "\<exists>i>n. coeff p i \<noteq> 0"
       
   135     by (rule less_degree_imp)
       
   136   then obtain i where "n < i" and "coeff p i \<noteq> 0"
       
   137     by blast
       
   138   from \<open>degree p = Suc n\<close> and \<open>n < i\<close> have "degree p \<le> i"
       
   139     by simp
       
   140   also from \<open>coeff p i \<noteq> 0\<close> have "i \<le> degree p"
       
   141     by (rule le_degree)
       
   142   finally have "degree p = i" .
       
   143   with \<open>coeff p i \<noteq> 0\<close> show "coeff p (degree p) \<noteq> 0" by simp
       
   144 qed
       
   145 
       
   146 lemma leading_coeff_0_iff [simp]: "coeff p (degree p) = 0 \<longleftrightarrow> p = 0"
       
   147   by (cases "p = 0") (simp_all add: leading_coeff_neq_0)
       
   148 
       
   149 lemma eq_zero_or_degree_less:
       
   150   assumes "degree p \<le> n" and "coeff p n = 0"
       
   151   shows "p = 0 \<or> degree p < n"
       
   152 proof (cases n)
       
   153   case 0
       
   154   with \<open>degree p \<le> n\<close> and \<open>coeff p n = 0\<close> have "coeff p (degree p) = 0"
       
   155     by simp
       
   156   then have "p = 0" by simp
       
   157   then show ?thesis ..
       
   158 next
       
   159   case (Suc m)
       
   160   from \<open>degree p \<le> n\<close> have "\<forall>i>n. coeff p i = 0"
       
   161     by (simp add: coeff_eq_0)
       
   162   with \<open>coeff p n = 0\<close> have "\<forall>i\<ge>n. coeff p i = 0"
       
   163     by (simp add: le_less)
       
   164   with \<open>n = Suc m\<close> have "\<forall>i>m. coeff p i = 0"
       
   165     by (simp add: less_eq_Suc_le)
       
   166   then have "degree p \<le> m"
       
   167     by (rule degree_le)
       
   168   with \<open>n = Suc m\<close> have "degree p < n"
       
   169     by (simp add: less_Suc_eq_le)
       
   170   then show ?thesis ..
       
   171 qed
       
   172 
       
   173 lemma coeff_0_degree_minus_1: "coeff rrr dr = 0 \<Longrightarrow> degree rrr \<le> dr \<Longrightarrow> degree rrr \<le> dr - 1"
       
   174   using eq_zero_or_degree_less by fastforce
       
   175 
       
   176 
       
   177 subsection \<open>List-style constructor for polynomials\<close>
       
   178 
       
   179 lift_definition pCons :: "'a::zero \<Rightarrow> 'a poly \<Rightarrow> 'a poly"
       
   180   is "\<lambda>a p. case_nat a (coeff p)"
       
   181   by (rule MOST_SucD) (simp add: MOST_coeff_eq_0)
       
   182 
       
   183 lemmas coeff_pCons = pCons.rep_eq
       
   184 
       
   185 lemma coeff_pCons_0 [simp]: "coeff (pCons a p) 0 = a"
       
   186   by transfer simp
       
   187 
       
   188 lemma coeff_pCons_Suc [simp]: "coeff (pCons a p) (Suc n) = coeff p n"
       
   189   by (simp add: coeff_pCons)
       
   190 
       
   191 lemma degree_pCons_le: "degree (pCons a p) \<le> Suc (degree p)"
       
   192   by (rule degree_le) (simp add: coeff_eq_0 coeff_pCons split: nat.split)
       
   193 
       
   194 lemma degree_pCons_eq: "p \<noteq> 0 \<Longrightarrow> degree (pCons a p) = Suc (degree p)"
       
   195   apply (rule order_antisym [OF degree_pCons_le])
       
   196   apply (rule le_degree, simp)
       
   197   done
       
   198 
       
   199 lemma degree_pCons_0: "degree (pCons a 0) = 0"
       
   200   apply (rule order_antisym [OF _ le0])
       
   201   apply (rule degree_le, simp add: coeff_pCons split: nat.split)
       
   202   done
       
   203 
       
   204 lemma degree_pCons_eq_if [simp]: "degree (pCons a p) = (if p = 0 then 0 else Suc (degree p))"
       
   205   apply (cases "p = 0", simp_all)
       
   206   apply (rule order_antisym [OF _ le0])
       
   207   apply (rule degree_le, simp add: coeff_pCons split: nat.split)
       
   208   apply (rule order_antisym [OF degree_pCons_le])
       
   209   apply (rule le_degree, simp)
       
   210   done
       
   211 
       
   212 lemma pCons_0_0 [simp]: "pCons 0 0 = 0"
       
   213   by (rule poly_eqI) (simp add: coeff_pCons split: nat.split)
       
   214 
       
   215 lemma pCons_eq_iff [simp]: "pCons a p = pCons b q \<longleftrightarrow> a = b \<and> p = q"
       
   216 proof safe
       
   217   assume "pCons a p = pCons b q"
       
   218   then have "coeff (pCons a p) 0 = coeff (pCons b q) 0"
       
   219     by simp
       
   220   then show "a = b"
       
   221     by simp
       
   222 next
       
   223   assume "pCons a p = pCons b q"
       
   224   then have "coeff (pCons a p) (Suc n) = coeff (pCons b q) (Suc n)" for n
       
   225     by simp
       
   226   then show "p = q"
       
   227     by (simp add: poly_eq_iff)
       
   228 qed
       
   229 
       
   230 lemma pCons_eq_0_iff [simp]: "pCons a p = 0 \<longleftrightarrow> a = 0 \<and> p = 0"
       
   231   using pCons_eq_iff [of a p 0 0] by simp
       
   232 
       
   233 lemma pCons_cases [cases type: poly]:
       
   234   obtains (pCons) a q where "p = pCons a q"
       
   235 proof
       
   236   show "p = pCons (coeff p 0) (Abs_poly (\<lambda>n. coeff p (Suc n)))"
       
   237     by transfer
       
   238       (simp_all add: MOST_inj[where f=Suc and P="\<lambda>n. p n = 0" for p] fun_eq_iff Abs_poly_inverse
       
   239         split: nat.split)
       
   240 qed
       
   241 
       
   242 lemma pCons_induct [case_names 0 pCons, induct type: poly]:
       
   243   assumes zero: "P 0"
       
   244   assumes pCons: "\<And>a p. a \<noteq> 0 \<or> p \<noteq> 0 \<Longrightarrow> P p \<Longrightarrow> P (pCons a p)"
       
   245   shows "P p"
       
   246 proof (induct p rule: measure_induct_rule [where f=degree])
       
   247   case (less p)
       
   248   obtain a q where "p = pCons a q" by (rule pCons_cases)
       
   249   have "P q"
       
   250   proof (cases "q = 0")
       
   251     case True
       
   252     then show "P q" by (simp add: zero)
       
   253   next
       
   254     case False
       
   255     then have "degree (pCons a q) = Suc (degree q)"
       
   256       by (rule degree_pCons_eq)
       
   257     with \<open>p = pCons a q\<close> have "degree q < degree p"
       
   258       by simp
       
   259     then show "P q"
       
   260       by (rule less.hyps)
       
   261   qed
       
   262   have "P (pCons a q)"
       
   263   proof (cases "a \<noteq> 0 \<or> q \<noteq> 0")
       
   264     case True
       
   265     with \<open>P q\<close> show ?thesis by (auto intro: pCons)
       
   266   next
       
   267     case False
       
   268     with zero show ?thesis by simp
       
   269   qed
       
   270   with \<open>p = pCons a q\<close> show ?case
       
   271     by simp
       
   272 qed
       
   273 
       
   274 lemma degree_eq_zeroE:
       
   275   fixes p :: "'a::zero poly"
       
   276   assumes "degree p = 0"
       
   277   obtains a where "p = pCons a 0"
       
   278 proof -
       
   279   obtain a q where p: "p = pCons a q"
       
   280     by (cases p)
       
   281   with assms have "q = 0"
       
   282     by (cases "q = 0") simp_all
       
   283   with p have "p = pCons a 0"
       
   284     by simp
       
   285   then show thesis ..
       
   286 qed
       
   287 
       
   288 
       
   289 subsection \<open>Quickcheck generator for polynomials\<close>
       
   290 
       
   291 quickcheck_generator poly constructors: "0 :: _ poly", pCons
       
   292 
       
   293 
       
   294 subsection \<open>List-style syntax for polynomials\<close>
       
   295 
       
   296 syntax "_poly" :: "args \<Rightarrow> 'a poly"  ("[:(_):]")
       
   297 translations
       
   298   "[:x, xs:]" \<rightleftharpoons> "CONST pCons x [:xs:]"
       
   299   "[:x:]" \<rightleftharpoons> "CONST pCons x 0"
       
   300   "[:x:]" \<leftharpoondown> "CONST pCons x (_constrain 0 t)"
       
   301 
       
   302 
       
   303 subsection \<open>Representation of polynomials by lists of coefficients\<close>
       
   304 
       
   305 primrec Poly :: "'a::zero list \<Rightarrow> 'a poly"
       
   306   where
       
   307     [code_post]: "Poly [] = 0"
       
   308   | [code_post]: "Poly (a # as) = pCons a (Poly as)"
       
   309 
       
   310 lemma Poly_replicate_0 [simp]: "Poly (replicate n 0) = 0"
       
   311   by (induct n) simp_all
       
   312 
       
   313 lemma Poly_eq_0: "Poly as = 0 \<longleftrightarrow> (\<exists>n. as = replicate n 0)"
       
   314   by (induct as) (auto simp add: Cons_replicate_eq)
       
   315 
       
   316 lemma Poly_append_replicate_zero [simp]: "Poly (as @ replicate n 0) = Poly as"
       
   317   by (induct as) simp_all
       
   318 
       
   319 lemma Poly_snoc_zero [simp]: "Poly (as @ [0]) = Poly as"
       
   320   using Poly_append_replicate_zero [of as 1] by simp
       
   321 
       
   322 lemma Poly_cCons_eq_pCons_Poly [simp]: "Poly (a ## p) = pCons a (Poly p)"
       
   323   by (simp add: cCons_def)
       
   324 
       
   325 lemma Poly_on_rev_starting_with_0 [simp]: "hd as = 0 \<Longrightarrow> Poly (rev (tl as)) = Poly (rev as)"
       
   326   by (cases as) simp_all
       
   327 
       
   328 lemma degree_Poly: "degree (Poly xs) \<le> length xs"
       
   329   by (induct xs) simp_all
       
   330 
       
   331 lemma coeff_Poly_eq [simp]: "coeff (Poly xs) = nth_default 0 xs"
       
   332   by (induct xs) (simp_all add: fun_eq_iff coeff_pCons split: nat.splits)
       
   333 
       
   334 definition coeffs :: "'a poly \<Rightarrow> 'a::zero list"
       
   335   where "coeffs p = (if p = 0 then [] else map (\<lambda>i. coeff p i) [0 ..< Suc (degree p)])"
       
   336 
       
   337 lemma coeffs_eq_Nil [simp]: "coeffs p = [] \<longleftrightarrow> p = 0"
       
   338   by (simp add: coeffs_def)
       
   339 
       
   340 lemma not_0_coeffs_not_Nil: "p \<noteq> 0 \<Longrightarrow> coeffs p \<noteq> []"
       
   341   by simp
       
   342 
       
   343 lemma coeffs_0_eq_Nil [simp]: "coeffs 0 = []"
       
   344   by simp
       
   345 
       
   346 lemma coeffs_pCons_eq_cCons [simp]: "coeffs (pCons a p) = a ## coeffs p"
       
   347 proof -
       
   348   have *: "\<forall>m\<in>set ms. m > 0 \<Longrightarrow> map (case_nat x f) ms = map f (map (\<lambda>n. n - 1) ms)"
       
   349     for ms :: "nat list" and f :: "nat \<Rightarrow> 'a" and x :: "'a"
       
   350     by (induct ms) (auto split: nat.split)
       
   351   show ?thesis
       
   352     by (simp add: * coeffs_def upt_conv_Cons coeff_pCons map_decr_upt del: upt_Suc)
       
   353 qed
       
   354 
       
   355 lemma length_coeffs: "p \<noteq> 0 \<Longrightarrow> length (coeffs p) = degree p + 1"
       
   356   by (simp add: coeffs_def)
       
   357 
       
   358 lemma coeffs_nth: "p \<noteq> 0 \<Longrightarrow> n \<le> degree p \<Longrightarrow> coeffs p ! n = coeff p n"
       
   359   by (auto simp: coeffs_def simp del: upt_Suc)
       
   360 
       
   361 lemma coeff_in_coeffs: "p \<noteq> 0 \<Longrightarrow> n \<le> degree p \<Longrightarrow> coeff p n \<in> set (coeffs p)"
       
   362   using coeffs_nth [of p n, symmetric] by (simp add: length_coeffs)
       
   363 
       
   364 lemma not_0_cCons_eq [simp]: "p \<noteq> 0 \<Longrightarrow> a ## coeffs p = a # coeffs p"
       
   365   by (simp add: cCons_def)
       
   366 
       
   367 lemma Poly_coeffs [simp, code abstype]: "Poly (coeffs p) = p"
       
   368   by (induct p) auto
       
   369 
       
   370 lemma coeffs_Poly [simp]: "coeffs (Poly as) = strip_while (HOL.eq 0) as"
       
   371 proof (induct as)
       
   372   case Nil
       
   373   then show ?case by simp
       
   374 next
       
   375   case (Cons a as)
       
   376   from replicate_length_same [of as 0] have "(\<forall>n. as \<noteq> replicate n 0) \<longleftrightarrow> (\<exists>a\<in>set as. a \<noteq> 0)"
       
   377     by (auto dest: sym [of _ as])
       
   378   with Cons show ?case by auto
       
   379 qed
       
   380 
       
   381 lemma no_trailing_coeffs [simp]:
       
   382   "no_trailing (HOL.eq 0) (coeffs p)"
       
   383   by (induct p)  auto
       
   384 
       
   385 lemma strip_while_coeffs [simp]:
       
   386   "strip_while (HOL.eq 0) (coeffs p) = coeffs p"
       
   387   by simp
       
   388 
       
   389 lemma coeffs_eq_iff: "p = q \<longleftrightarrow> coeffs p = coeffs q"
       
   390   (is "?P \<longleftrightarrow> ?Q")
       
   391 proof
       
   392   assume ?P
       
   393   then show ?Q by simp
       
   394 next
       
   395   assume ?Q
       
   396   then have "Poly (coeffs p) = Poly (coeffs q)" by simp
       
   397   then show ?P by simp
       
   398 qed
       
   399 
       
   400 lemma nth_default_coeffs_eq: "nth_default 0 (coeffs p) = coeff p"
       
   401   by (simp add: fun_eq_iff coeff_Poly_eq [symmetric])
       
   402 
       
   403 lemma [code]: "coeff p = nth_default 0 (coeffs p)"
       
   404   by (simp add: nth_default_coeffs_eq)
       
   405 
       
   406 lemma coeffs_eqI:
       
   407   assumes coeff: "\<And>n. coeff p n = nth_default 0 xs n"
       
   408   assumes zero: "no_trailing (HOL.eq 0) xs"
       
   409   shows "coeffs p = xs"
       
   410 proof -
       
   411   from coeff have "p = Poly xs"
       
   412     by (simp add: poly_eq_iff)
       
   413   with zero show ?thesis by simp
       
   414 qed
       
   415 
       
   416 lemma degree_eq_length_coeffs [code]: "degree p = length (coeffs p) - 1"
       
   417   by (simp add: coeffs_def)
       
   418 
       
   419 lemma length_coeffs_degree: "p \<noteq> 0 \<Longrightarrow> length (coeffs p) = Suc (degree p)"
       
   420   by (induct p) (auto simp: cCons_def)
       
   421 
       
   422 lemma [code abstract]: "coeffs 0 = []"
       
   423   by (fact coeffs_0_eq_Nil)
       
   424 
       
   425 lemma [code abstract]: "coeffs (pCons a p) = a ## coeffs p"
       
   426   by (fact coeffs_pCons_eq_cCons)
       
   427 
       
   428 instantiation poly :: ("{zero, equal}") equal
       
   429 begin
       
   430 
       
   431 definition [code]: "HOL.equal (p::'a poly) q \<longleftrightarrow> HOL.equal (coeffs p) (coeffs q)"
       
   432 
       
   433 instance
       
   434   by standard (simp add: equal equal_poly_def coeffs_eq_iff)
       
   435 
       
   436 end
       
   437 
       
   438 lemma [code nbe]: "HOL.equal (p :: _ poly) p \<longleftrightarrow> True"
       
   439   by (fact equal_refl)
       
   440 
       
   441 definition is_zero :: "'a::zero poly \<Rightarrow> bool"
       
   442   where [code]: "is_zero p \<longleftrightarrow> List.null (coeffs p)"
       
   443 
       
   444 lemma is_zero_null [code_abbrev]: "is_zero p \<longleftrightarrow> p = 0"
       
   445   by (simp add: is_zero_def null_def)
       
   446 
       
   447 
       
   448 subsubsection \<open>Reconstructing the polynomial from the list\<close>
       
   449   \<comment> \<open>contributed by Sebastiaan J.C. Joosten and René Thiemann\<close>
       
   450 
       
   451 definition poly_of_list :: "'a::comm_monoid_add list \<Rightarrow> 'a poly"
       
   452   where [simp]: "poly_of_list = Poly"
       
   453 
       
   454 lemma poly_of_list_impl [code abstract]: "coeffs (poly_of_list as) = strip_while (HOL.eq 0) as"
       
   455   by simp
       
   456 
       
   457 
       
   458 subsection \<open>Fold combinator for polynomials\<close>
       
   459 
       
   460 definition fold_coeffs :: "('a::zero \<Rightarrow> 'b \<Rightarrow> 'b) \<Rightarrow> 'a poly \<Rightarrow> 'b \<Rightarrow> 'b"
       
   461   where "fold_coeffs f p = foldr f (coeffs p)"
       
   462 
       
   463 lemma fold_coeffs_0_eq [simp]: "fold_coeffs f 0 = id"
       
   464   by (simp add: fold_coeffs_def)
       
   465 
       
   466 lemma fold_coeffs_pCons_eq [simp]: "f 0 = id \<Longrightarrow> fold_coeffs f (pCons a p) = f a \<circ> fold_coeffs f p"
       
   467   by (simp add: fold_coeffs_def cCons_def fun_eq_iff)
       
   468 
       
   469 lemma fold_coeffs_pCons_0_0_eq [simp]: "fold_coeffs f (pCons 0 0) = id"
       
   470   by (simp add: fold_coeffs_def)
       
   471 
       
   472 lemma fold_coeffs_pCons_coeff_not_0_eq [simp]:
       
   473   "a \<noteq> 0 \<Longrightarrow> fold_coeffs f (pCons a p) = f a \<circ> fold_coeffs f p"
       
   474   by (simp add: fold_coeffs_def)
       
   475 
       
   476 lemma fold_coeffs_pCons_not_0_0_eq [simp]:
       
   477   "p \<noteq> 0 \<Longrightarrow> fold_coeffs f (pCons a p) = f a \<circ> fold_coeffs f p"
       
   478   by (simp add: fold_coeffs_def)
       
   479 
       
   480 
       
   481 subsection \<open>Canonical morphism on polynomials -- evaluation\<close>
       
   482 
       
   483 definition poly :: "'a::comm_semiring_0 poly \<Rightarrow> 'a \<Rightarrow> 'a"
       
   484   where "poly p = fold_coeffs (\<lambda>a f x. a + x * f x) p (\<lambda>x. 0)" \<comment> \<open>The Horner Schema\<close>
       
   485 
       
   486 lemma poly_0 [simp]: "poly 0 x = 0"
       
   487   by (simp add: poly_def)
       
   488 
       
   489 lemma poly_pCons [simp]: "poly (pCons a p) x = a + x * poly p x"
       
   490   by (cases "p = 0 \<and> a = 0") (auto simp add: poly_def)
       
   491 
       
   492 lemma poly_altdef: "poly p x = (\<Sum>i\<le>degree p. coeff p i * x ^ i)"
       
   493   for x :: "'a::{comm_semiring_0,semiring_1}"
       
   494 proof (induction p rule: pCons_induct)
       
   495   case 0
       
   496   then show ?case
       
   497     by simp
       
   498 next
       
   499   case (pCons a p)
       
   500   show ?case
       
   501   proof (cases "p = 0")
       
   502     case True
       
   503     then show ?thesis by simp
       
   504   next
       
   505     case False
       
   506     let ?p' = "pCons a p"
       
   507     note poly_pCons[of a p x]
       
   508     also note pCons.IH
       
   509     also have "a + x * (\<Sum>i\<le>degree p. coeff p i * x ^ i) =
       
   510         coeff ?p' 0 * x^0 + (\<Sum>i\<le>degree p. coeff ?p' (Suc i) * x^Suc i)"
       
   511       by (simp add: field_simps sum_distrib_left coeff_pCons)
       
   512     also note sum_atMost_Suc_shift[symmetric]
       
   513     also note degree_pCons_eq[OF \<open>p \<noteq> 0\<close>, of a, symmetric]
       
   514     finally show ?thesis .
       
   515   qed
       
   516 qed
       
   517 
       
   518 lemma poly_0_coeff_0: "poly p 0 = coeff p 0"
       
   519   by (cases p) (auto simp: poly_altdef)
       
   520 
       
   521 
       
   522 subsection \<open>Monomials\<close>
       
   523 
       
   524 lift_definition monom :: "'a \<Rightarrow> nat \<Rightarrow> 'a::zero poly"
       
   525   is "\<lambda>a m n. if m = n then a else 0"
       
   526   by (simp add: MOST_iff_cofinite)
       
   527 
       
   528 lemma coeff_monom [simp]: "coeff (monom a m) n = (if m = n then a else 0)"
       
   529   by transfer rule
       
   530 
       
   531 lemma monom_0: "monom a 0 = pCons a 0"
       
   532   by (rule poly_eqI) (simp add: coeff_pCons split: nat.split)
       
   533 
       
   534 lemma monom_Suc: "monom a (Suc n) = pCons 0 (monom a n)"
       
   535   by (rule poly_eqI) (simp add: coeff_pCons split: nat.split)
       
   536 
       
   537 lemma monom_eq_0 [simp]: "monom 0 n = 0"
       
   538   by (rule poly_eqI) simp
       
   539 
       
   540 lemma monom_eq_0_iff [simp]: "monom a n = 0 \<longleftrightarrow> a = 0"
       
   541   by (simp add: poly_eq_iff)
       
   542 
       
   543 lemma monom_eq_iff [simp]: "monom a n = monom b n \<longleftrightarrow> a = b"
       
   544   by (simp add: poly_eq_iff)
       
   545 
       
   546 lemma degree_monom_le: "degree (monom a n) \<le> n"
       
   547   by (rule degree_le, simp)
       
   548 
       
   549 lemma degree_monom_eq: "a \<noteq> 0 \<Longrightarrow> degree (monom a n) = n"
       
   550   apply (rule order_antisym [OF degree_monom_le])
       
   551   apply (rule le_degree)
       
   552   apply simp
       
   553   done
       
   554 
       
   555 lemma coeffs_monom [code abstract]:
       
   556   "coeffs (monom a n) = (if a = 0 then [] else replicate n 0 @ [a])"
       
   557   by (induct n) (simp_all add: monom_0 monom_Suc)
       
   558 
       
   559 lemma fold_coeffs_monom [simp]: "a \<noteq> 0 \<Longrightarrow> fold_coeffs f (monom a n) = f 0 ^^ n \<circ> f a"
       
   560   by (simp add: fold_coeffs_def coeffs_monom fun_eq_iff)
       
   561 
       
   562 lemma poly_monom: "poly (monom a n) x = a * x ^ n"
       
   563   for a x :: "'a::comm_semiring_1"
       
   564   by (cases "a = 0", simp_all) (induct n, simp_all add: mult.left_commute poly_def)
       
   565 
       
   566 lemma monom_eq_iff': "monom c n = monom d m \<longleftrightarrow>  c = d \<and> (c = 0 \<or> n = m)"
       
   567   by (auto simp: poly_eq_iff)
       
   568 
       
   569 lemma monom_eq_const_iff: "monom c n = [:d:] \<longleftrightarrow> c = d \<and> (c = 0 \<or> n = 0)"
       
   570   using monom_eq_iff'[of c n d 0] by (simp add: monom_0)
       
   571 
       
   572 
       
   573 subsection \<open>Leading coefficient\<close>
       
   574 
       
   575 abbreviation lead_coeff:: "'a::zero poly \<Rightarrow> 'a"
       
   576   where "lead_coeff p \<equiv> coeff p (degree p)"
       
   577 
       
   578 lemma lead_coeff_pCons[simp]:
       
   579   "p \<noteq> 0 \<Longrightarrow> lead_coeff (pCons a p) = lead_coeff p"
       
   580   "p = 0 \<Longrightarrow> lead_coeff (pCons a p) = a"
       
   581   by auto
       
   582 
       
   583 lemma lead_coeff_monom [simp]: "lead_coeff (monom c n) = c"
       
   584   by (cases "c = 0") (simp_all add: degree_monom_eq)
       
   585 
       
   586 
       
   587 subsection \<open>Addition and subtraction\<close>
       
   588 
       
   589 instantiation poly :: (comm_monoid_add) comm_monoid_add
       
   590 begin
       
   591 
       
   592 lift_definition plus_poly :: "'a poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly"
       
   593   is "\<lambda>p q n. coeff p n + coeff q n"
       
   594 proof -
       
   595   fix q p :: "'a poly"
       
   596   show "\<forall>\<^sub>\<infinity>n. coeff p n + coeff q n = 0"
       
   597     using MOST_coeff_eq_0[of p] MOST_coeff_eq_0[of q] by eventually_elim simp
       
   598 qed
       
   599 
       
   600 lemma coeff_add [simp]: "coeff (p + q) n = coeff p n + coeff q n"
       
   601   by (simp add: plus_poly.rep_eq)
       
   602 
       
   603 instance
       
   604 proof
       
   605   fix p q r :: "'a poly"
       
   606   show "(p + q) + r = p + (q + r)"
       
   607     by (simp add: poly_eq_iff add.assoc)
       
   608   show "p + q = q + p"
       
   609     by (simp add: poly_eq_iff add.commute)
       
   610   show "0 + p = p"
       
   611     by (simp add: poly_eq_iff)
       
   612 qed
       
   613 
       
   614 end
       
   615 
       
   616 instantiation poly :: (cancel_comm_monoid_add) cancel_comm_monoid_add
       
   617 begin
       
   618 
       
   619 lift_definition minus_poly :: "'a poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly"
       
   620   is "\<lambda>p q n. coeff p n - coeff q n"
       
   621 proof -
       
   622   fix q p :: "'a poly"
       
   623   show "\<forall>\<^sub>\<infinity>n. coeff p n - coeff q n = 0"
       
   624     using MOST_coeff_eq_0[of p] MOST_coeff_eq_0[of q] by eventually_elim simp
       
   625 qed
       
   626 
       
   627 lemma coeff_diff [simp]: "coeff (p - q) n = coeff p n - coeff q n"
       
   628   by (simp add: minus_poly.rep_eq)
       
   629 
       
   630 instance
       
   631 proof
       
   632   fix p q r :: "'a poly"
       
   633   show "p + q - p = q"
       
   634     by (simp add: poly_eq_iff)
       
   635   show "p - q - r = p - (q + r)"
       
   636     by (simp add: poly_eq_iff diff_diff_eq)
       
   637 qed
       
   638 
       
   639 end
       
   640 
       
   641 instantiation poly :: (ab_group_add) ab_group_add
       
   642 begin
       
   643 
       
   644 lift_definition uminus_poly :: "'a poly \<Rightarrow> 'a poly"
       
   645   is "\<lambda>p n. - coeff p n"
       
   646 proof -
       
   647   fix p :: "'a poly"
       
   648   show "\<forall>\<^sub>\<infinity>n. - coeff p n = 0"
       
   649     using MOST_coeff_eq_0 by simp
       
   650 qed
       
   651 
       
   652 lemma coeff_minus [simp]: "coeff (- p) n = - coeff p n"
       
   653   by (simp add: uminus_poly.rep_eq)
       
   654 
       
   655 instance
       
   656 proof
       
   657   fix p q :: "'a poly"
       
   658   show "- p + p = 0"
       
   659     by (simp add: poly_eq_iff)
       
   660   show "p - q = p + - q"
       
   661     by (simp add: poly_eq_iff)
       
   662 qed
       
   663 
       
   664 end
       
   665 
       
   666 lemma add_pCons [simp]: "pCons a p + pCons b q = pCons (a + b) (p + q)"
       
   667   by (rule poly_eqI) (simp add: coeff_pCons split: nat.split)
       
   668 
       
   669 lemma minus_pCons [simp]: "- pCons a p = pCons (- a) (- p)"
       
   670   by (rule poly_eqI) (simp add: coeff_pCons split: nat.split)
       
   671 
       
   672 lemma diff_pCons [simp]: "pCons a p - pCons b q = pCons (a - b) (p - q)"
       
   673   by (rule poly_eqI) (simp add: coeff_pCons split: nat.split)
       
   674 
       
   675 lemma degree_add_le_max: "degree (p + q) \<le> max (degree p) (degree q)"
       
   676   by (rule degree_le) (auto simp add: coeff_eq_0)
       
   677 
       
   678 lemma degree_add_le: "degree p \<le> n \<Longrightarrow> degree q \<le> n \<Longrightarrow> degree (p + q) \<le> n"
       
   679   by (auto intro: order_trans degree_add_le_max)
       
   680 
       
   681 lemma degree_add_less: "degree p < n \<Longrightarrow> degree q < n \<Longrightarrow> degree (p + q) < n"
       
   682   by (auto intro: le_less_trans degree_add_le_max)
       
   683 
       
   684 lemma degree_add_eq_right: "degree p < degree q \<Longrightarrow> degree (p + q) = degree q"
       
   685   apply (cases "q = 0")
       
   686    apply simp
       
   687   apply (rule order_antisym)
       
   688    apply (simp add: degree_add_le)
       
   689   apply (rule le_degree)
       
   690   apply (simp add: coeff_eq_0)
       
   691   done
       
   692 
       
   693 lemma degree_add_eq_left: "degree q < degree p \<Longrightarrow> degree (p + q) = degree p"
       
   694   using degree_add_eq_right [of q p] by (simp add: add.commute)
       
   695 
       
   696 lemma degree_minus [simp]: "degree (- p) = degree p"
       
   697   by (simp add: degree_def)
       
   698 
       
   699 lemma lead_coeff_add_le: "degree p < degree q \<Longrightarrow> lead_coeff (p + q) = lead_coeff q"
       
   700   by (metis coeff_add coeff_eq_0 monoid_add_class.add.left_neutral degree_add_eq_right)
       
   701 
       
   702 lemma lead_coeff_minus: "lead_coeff (- p) = - lead_coeff p"
       
   703   by (metis coeff_minus degree_minus)
       
   704 
       
   705 lemma degree_diff_le_max: "degree (p - q) \<le> max (degree p) (degree q)"
       
   706   for p q :: "'a::ab_group_add poly"
       
   707   using degree_add_le [where p=p and q="-q"] by simp
       
   708 
       
   709 lemma degree_diff_le: "degree p \<le> n \<Longrightarrow> degree q \<le> n \<Longrightarrow> degree (p - q) \<le> n"
       
   710   for p q :: "'a::ab_group_add poly"
       
   711   using degree_add_le [of p n "- q"] by simp
       
   712 
       
   713 lemma degree_diff_less: "degree p < n \<Longrightarrow> degree q < n \<Longrightarrow> degree (p - q) < n"
       
   714   for p q :: "'a::ab_group_add poly"
       
   715   using degree_add_less [of p n "- q"] by simp
       
   716 
       
   717 lemma add_monom: "monom a n + monom b n = monom (a + b) n"
       
   718   by (rule poly_eqI) simp
       
   719 
       
   720 lemma diff_monom: "monom a n - monom b n = monom (a - b) n"
       
   721   by (rule poly_eqI) simp
       
   722 
       
   723 lemma minus_monom: "- monom a n = monom (- a) n"
       
   724   by (rule poly_eqI) simp
       
   725 
       
   726 lemma coeff_sum: "coeff (\<Sum>x\<in>A. p x) i = (\<Sum>x\<in>A. coeff (p x) i)"
       
   727   by (induct A rule: infinite_finite_induct) simp_all
       
   728 
       
   729 lemma monom_sum: "monom (\<Sum>x\<in>A. a x) n = (\<Sum>x\<in>A. monom (a x) n)"
       
   730   by (rule poly_eqI) (simp add: coeff_sum)
       
   731 
       
   732 fun plus_coeffs :: "'a::comm_monoid_add list \<Rightarrow> 'a list \<Rightarrow> 'a list"
       
   733   where
       
   734     "plus_coeffs xs [] = xs"
       
   735   | "plus_coeffs [] ys = ys"
       
   736   | "plus_coeffs (x # xs) (y # ys) = (x + y) ## plus_coeffs xs ys"
       
   737 
       
   738 lemma coeffs_plus_eq_plus_coeffs [code abstract]:
       
   739   "coeffs (p + q) = plus_coeffs (coeffs p) (coeffs q)"
       
   740 proof -
       
   741   have *: "nth_default 0 (plus_coeffs xs ys) n = nth_default 0 xs n + nth_default 0 ys n"
       
   742     for xs ys :: "'a list" and n
       
   743   proof (induct xs ys arbitrary: n rule: plus_coeffs.induct)
       
   744     case (3 x xs y ys n)
       
   745     then show ?case
       
   746       by (cases n) (auto simp add: cCons_def)
       
   747   qed simp_all
       
   748   have **: "no_trailing (HOL.eq 0) (plus_coeffs xs ys)"
       
   749     if "no_trailing (HOL.eq 0) xs" and "no_trailing (HOL.eq 0) ys"
       
   750     for xs ys :: "'a list"
       
   751     using that by (induct xs ys rule: plus_coeffs.induct) (simp_all add: cCons_def)
       
   752   show ?thesis
       
   753     by (rule coeffs_eqI) (auto simp add: * nth_default_coeffs_eq intro: **)
       
   754 qed
       
   755 
       
   756 lemma coeffs_uminus [code abstract]:
       
   757   "coeffs (- p) = map uminus (coeffs p)"
       
   758 proof -
       
   759   have eq_0: "HOL.eq 0 \<circ> uminus = HOL.eq (0::'a)"
       
   760     by (simp add: fun_eq_iff)
       
   761   show ?thesis
       
   762     by (rule coeffs_eqI) (simp_all add: nth_default_map_eq nth_default_coeffs_eq no_trailing_map eq_0)
       
   763 qed
       
   764 
       
   765 lemma [code]: "p - q = p + - q"
       
   766   for p q :: "'a::ab_group_add poly"
       
   767   by (fact diff_conv_add_uminus)
       
   768 
       
   769 lemma poly_add [simp]: "poly (p + q) x = poly p x + poly q x"
       
   770   apply (induct p arbitrary: q)
       
   771    apply simp
       
   772   apply (case_tac q, simp, simp add: algebra_simps)
       
   773   done
       
   774 
       
   775 lemma poly_minus [simp]: "poly (- p) x = - poly p x"
       
   776   for x :: "'a::comm_ring"
       
   777   by (induct p) simp_all
       
   778 
       
   779 lemma poly_diff [simp]: "poly (p - q) x = poly p x - poly q x"
       
   780   for x :: "'a::comm_ring"
       
   781   using poly_add [of p "- q" x] by simp
       
   782 
       
   783 lemma poly_sum: "poly (\<Sum>k\<in>A. p k) x = (\<Sum>k\<in>A. poly (p k) x)"
       
   784   by (induct A rule: infinite_finite_induct) simp_all
       
   785 
       
   786 lemma degree_sum_le: "finite S \<Longrightarrow> (\<And>p. p \<in> S \<Longrightarrow> degree (f p) \<le> n) \<Longrightarrow> degree (sum f S) \<le> n"
       
   787 proof (induct S rule: finite_induct)
       
   788   case empty
       
   789   then show ?case by simp
       
   790 next
       
   791   case (insert p S)
       
   792   then have "degree (sum f S) \<le> n" "degree (f p) \<le> n"
       
   793     by auto
       
   794   then show ?case
       
   795     unfolding sum.insert[OF insert(1-2)] by (metis degree_add_le)
       
   796 qed
       
   797 
       
   798 lemma poly_as_sum_of_monoms':
       
   799   assumes "degree p \<le> n"
       
   800   shows "(\<Sum>i\<le>n. monom (coeff p i) i) = p"
       
   801 proof -
       
   802   have eq: "\<And>i. {..n} \<inter> {i} = (if i \<le> n then {i} else {})"
       
   803     by auto
       
   804   from assms show ?thesis
       
   805     by (simp add: poly_eq_iff coeff_sum coeff_eq_0 sum.If_cases eq
       
   806         if_distrib[where f="\<lambda>x. x * a" for a])
       
   807 qed
       
   808 
       
   809 lemma poly_as_sum_of_monoms: "(\<Sum>i\<le>degree p. monom (coeff p i) i) = p"
       
   810   by (intro poly_as_sum_of_monoms' order_refl)
       
   811 
       
   812 lemma Poly_snoc: "Poly (xs @ [x]) = Poly xs + monom x (length xs)"
       
   813   by (induct xs) (simp_all add: monom_0 monom_Suc)
       
   814 
       
   815 
       
   816 subsection \<open>Multiplication by a constant, polynomial multiplication and the unit polynomial\<close>
       
   817 
       
   818 lift_definition smult :: "'a::comm_semiring_0 \<Rightarrow> 'a poly \<Rightarrow> 'a poly"
       
   819   is "\<lambda>a p n. a * coeff p n"
       
   820 proof -
       
   821   fix a :: 'a and p :: "'a poly"
       
   822   show "\<forall>\<^sub>\<infinity> i. a * coeff p i = 0"
       
   823     using MOST_coeff_eq_0[of p] by eventually_elim simp
       
   824 qed
       
   825 
       
   826 lemma coeff_smult [simp]: "coeff (smult a p) n = a * coeff p n"
       
   827   by (simp add: smult.rep_eq)
       
   828 
       
   829 lemma degree_smult_le: "degree (smult a p) \<le> degree p"
       
   830   by (rule degree_le) (simp add: coeff_eq_0)
       
   831 
       
   832 lemma smult_smult [simp]: "smult a (smult b p) = smult (a * b) p"
       
   833   by (rule poly_eqI) (simp add: mult.assoc)
       
   834 
       
   835 lemma smult_0_right [simp]: "smult a 0 = 0"
       
   836   by (rule poly_eqI) simp
       
   837 
       
   838 lemma smult_0_left [simp]: "smult 0 p = 0"
       
   839   by (rule poly_eqI) simp
       
   840 
       
   841 lemma smult_1_left [simp]: "smult (1::'a::comm_semiring_1) p = p"
       
   842   by (rule poly_eqI) simp
       
   843 
       
   844 lemma smult_add_right: "smult a (p + q) = smult a p + smult a q"
       
   845   by (rule poly_eqI) (simp add: algebra_simps)
       
   846 
       
   847 lemma smult_add_left: "smult (a + b) p = smult a p + smult b p"
       
   848   by (rule poly_eqI) (simp add: algebra_simps)
       
   849 
       
   850 lemma smult_minus_right [simp]: "smult a (- p) = - smult a p"
       
   851   for a :: "'a::comm_ring"
       
   852   by (rule poly_eqI) simp
       
   853 
       
   854 lemma smult_minus_left [simp]: "smult (- a) p = - smult a p"
       
   855   for a :: "'a::comm_ring"
       
   856   by (rule poly_eqI) simp
       
   857 
       
   858 lemma smult_diff_right: "smult a (p - q) = smult a p - smult a q"
       
   859   for a :: "'a::comm_ring"
       
   860   by (rule poly_eqI) (simp add: algebra_simps)
       
   861 
       
   862 lemma smult_diff_left: "smult (a - b) p = smult a p - smult b p"
       
   863   for a b :: "'a::comm_ring"
       
   864   by (rule poly_eqI) (simp add: algebra_simps)
       
   865 
       
   866 lemmas smult_distribs =
       
   867   smult_add_left smult_add_right
       
   868   smult_diff_left smult_diff_right
       
   869 
       
   870 lemma smult_pCons [simp]: "smult a (pCons b p) = pCons (a * b) (smult a p)"
       
   871   by (rule poly_eqI) (simp add: coeff_pCons split: nat.split)
       
   872 
       
   873 lemma smult_monom: "smult a (monom b n) = monom (a * b) n"
       
   874   by (induct n) (simp_all add: monom_0 monom_Suc)
       
   875 
       
   876 lemma smult_Poly: "smult c (Poly xs) = Poly (map (op * c) xs)"
       
   877   by (auto simp: poly_eq_iff nth_default_def)
       
   878 
       
   879 lemma degree_smult_eq [simp]: "degree (smult a p) = (if a = 0 then 0 else degree p)"
       
   880   for a :: "'a::{comm_semiring_0,semiring_no_zero_divisors}"
       
   881   by (cases "a = 0") (simp_all add: degree_def)
       
   882 
       
   883 lemma smult_eq_0_iff [simp]: "smult a p = 0 \<longleftrightarrow> a = 0 \<or> p = 0"
       
   884   for a :: "'a::{comm_semiring_0,semiring_no_zero_divisors}"
       
   885   by (simp add: poly_eq_iff)
       
   886 
       
   887 lemma coeffs_smult [code abstract]:
       
   888   "coeffs (smult a p) = (if a = 0 then [] else map (Groups.times a) (coeffs p))"
       
   889   for p :: "'a::{comm_semiring_0,semiring_no_zero_divisors} poly"
       
   890 proof -
       
   891   have eq_0: "HOL.eq 0 \<circ> times a = HOL.eq (0::'a)" if "a \<noteq> 0"
       
   892     using that by (simp add: fun_eq_iff)
       
   893   show ?thesis
       
   894     by (rule coeffs_eqI) (auto simp add: no_trailing_map nth_default_map_eq nth_default_coeffs_eq eq_0)
       
   895 qed  
       
   896 
       
   897 lemma smult_eq_iff:
       
   898   fixes b :: "'a :: field"
       
   899   assumes "b \<noteq> 0"
       
   900   shows "smult a p = smult b q \<longleftrightarrow> smult (a / b) p = q"
       
   901     (is "?lhs \<longleftrightarrow> ?rhs")
       
   902 proof
       
   903   assume ?lhs
       
   904   also from assms have "smult (inverse b) \<dots> = q"
       
   905     by simp
       
   906   finally show ?rhs
       
   907     by (simp add: field_simps)
       
   908 next
       
   909   assume ?rhs
       
   910   with assms show ?lhs by auto
       
   911 qed
       
   912 
       
   913 instantiation poly :: (comm_semiring_0) comm_semiring_0
       
   914 begin
       
   915 
       
   916 definition "p * q = fold_coeffs (\<lambda>a p. smult a q + pCons 0 p) p 0"
       
   917 
       
   918 lemma mult_poly_0_left: "(0::'a poly) * q = 0"
       
   919   by (simp add: times_poly_def)
       
   920 
       
   921 lemma mult_pCons_left [simp]: "pCons a p * q = smult a q + pCons 0 (p * q)"
       
   922   by (cases "p = 0 \<and> a = 0") (auto simp add: times_poly_def)
       
   923 
       
   924 lemma mult_poly_0_right: "p * (0::'a poly) = 0"
       
   925   by (induct p) (simp_all add: mult_poly_0_left)
       
   926 
       
   927 lemma mult_pCons_right [simp]: "p * pCons a q = smult a p + pCons 0 (p * q)"
       
   928   by (induct p) (simp_all add: mult_poly_0_left algebra_simps)
       
   929 
       
   930 lemmas mult_poly_0 = mult_poly_0_left mult_poly_0_right
       
   931 
       
   932 lemma mult_smult_left [simp]: "smult a p * q = smult a (p * q)"
       
   933   by (induct p) (simp_all add: mult_poly_0 smult_add_right)
       
   934 
       
   935 lemma mult_smult_right [simp]: "p * smult a q = smult a (p * q)"
       
   936   by (induct q) (simp_all add: mult_poly_0 smult_add_right)
       
   937 
       
   938 lemma mult_poly_add_left: "(p + q) * r = p * r + q * r"
       
   939   for p q r :: "'a poly"
       
   940   by (induct r) (simp_all add: mult_poly_0 smult_distribs algebra_simps)
       
   941 
       
   942 instance
       
   943 proof
       
   944   fix p q r :: "'a poly"
       
   945   show 0: "0 * p = 0"
       
   946     by (rule mult_poly_0_left)
       
   947   show "p * 0 = 0"
       
   948     by (rule mult_poly_0_right)
       
   949   show "(p + q) * r = p * r + q * r"
       
   950     by (rule mult_poly_add_left)
       
   951   show "(p * q) * r = p * (q * r)"
       
   952     by (induct p) (simp_all add: mult_poly_0 mult_poly_add_left)
       
   953   show "p * q = q * p"
       
   954     by (induct p) (simp_all add: mult_poly_0)
       
   955 qed
       
   956 
       
   957 end
       
   958 
       
   959 lemma coeff_mult_degree_sum:
       
   960   "coeff (p * q) (degree p + degree q) = coeff p (degree p) * coeff q (degree q)"
       
   961   by (induct p) (simp_all add: coeff_eq_0)
       
   962 
       
   963 instance poly :: ("{comm_semiring_0,semiring_no_zero_divisors}") semiring_no_zero_divisors
       
   964 proof
       
   965   fix p q :: "'a poly"
       
   966   assume "p \<noteq> 0" and "q \<noteq> 0"
       
   967   have "coeff (p * q) (degree p + degree q) = coeff p (degree p) * coeff q (degree q)"
       
   968     by (rule coeff_mult_degree_sum)
       
   969   also from \<open>p \<noteq> 0\<close> \<open>q \<noteq> 0\<close> have "coeff p (degree p) * coeff q (degree q) \<noteq> 0"
       
   970     by simp
       
   971   finally have "\<exists>n. coeff (p * q) n \<noteq> 0" ..
       
   972   then show "p * q \<noteq> 0"
       
   973     by (simp add: poly_eq_iff)
       
   974 qed
       
   975 
       
   976 instance poly :: (comm_semiring_0_cancel) comm_semiring_0_cancel ..
       
   977 
       
   978 lemma coeff_mult: "coeff (p * q) n = (\<Sum>i\<le>n. coeff p i * coeff q (n-i))"
       
   979 proof (induct p arbitrary: n)
       
   980   case 0
       
   981   show ?case by simp
       
   982 next
       
   983   case (pCons a p n)
       
   984   then show ?case
       
   985     by (cases n) (simp_all add: sum_atMost_Suc_shift del: sum_atMost_Suc)
       
   986 qed
       
   987 
       
   988 lemma degree_mult_le: "degree (p * q) \<le> degree p + degree q"
       
   989   apply (rule degree_le)
       
   990   apply (induct p)
       
   991    apply simp
       
   992   apply (simp add: coeff_eq_0 coeff_pCons split: nat.split)
       
   993   done
       
   994 
       
   995 lemma mult_monom: "monom a m * monom b n = monom (a * b) (m + n)"
       
   996   by (induct m) (simp add: monom_0 smult_monom, simp add: monom_Suc)
       
   997 
       
   998 instantiation poly :: (comm_semiring_1) comm_semiring_1
       
   999 begin
       
  1000 
       
  1001 definition one_poly_def: "1 = pCons 1 0"
       
  1002 
       
  1003 instance
       
  1004 proof
       
  1005   show "1 * p = p" for p :: "'a poly"
       
  1006     by (simp add: one_poly_def)
       
  1007   show "0 \<noteq> (1::'a poly)"
       
  1008     by (simp add: one_poly_def)
       
  1009 qed
       
  1010 
       
  1011 end
       
  1012 
       
  1013 instance poly :: ("{comm_semiring_1,semiring_1_no_zero_divisors}") semiring_1_no_zero_divisors ..
       
  1014 instance poly :: (comm_ring) comm_ring ..
       
  1015 instance poly :: (comm_ring_1) comm_ring_1 ..
       
  1016 instance poly :: (comm_ring_1) comm_semiring_1_cancel ..
       
  1017 
       
  1018 lemma coeff_1 [simp]: "coeff 1 n = (if n = 0 then 1 else 0)"
       
  1019   by (simp add: one_poly_def coeff_pCons split: nat.split)
       
  1020 
       
  1021 lemma monom_eq_1 [simp]: "monom 1 0 = 1"
       
  1022   by (simp add: monom_0 one_poly_def)
       
  1023 
       
  1024 lemma monom_eq_1_iff: "monom c n = 1 \<longleftrightarrow> c = 1 \<and> n = 0"
       
  1025   using monom_eq_const_iff[of c n 1] by (auto simp: one_poly_def)
       
  1026 
       
  1027 lemma monom_altdef: "monom c n = smult c ([:0, 1:]^n)"
       
  1028   by (induct n) (simp_all add: monom_0 monom_Suc one_poly_def)
       
  1029 
       
  1030 lemma degree_1 [simp]: "degree 1 = 0"
       
  1031   unfolding one_poly_def by (rule degree_pCons_0)
       
  1032 
       
  1033 lemma coeffs_1_eq [simp, code abstract]: "coeffs 1 = [1]"
       
  1034   by (simp add: one_poly_def)
       
  1035 
       
  1036 lemma degree_power_le: "degree (p ^ n) \<le> degree p * n"
       
  1037   by (induct n) (auto intro: order_trans degree_mult_le)
       
  1038 
       
  1039 lemma coeff_0_power: "coeff (p ^ n) 0 = coeff p 0 ^ n"
       
  1040   by (induct n) (simp_all add: coeff_mult)
       
  1041 
       
  1042 lemma poly_smult [simp]: "poly (smult a p) x = a * poly p x"
       
  1043   by (induct p) (simp_all add: algebra_simps)
       
  1044 
       
  1045 lemma poly_mult [simp]: "poly (p * q) x = poly p x * poly q x"
       
  1046   by (induct p) (simp_all add: algebra_simps)
       
  1047 
       
  1048 lemma poly_1 [simp]: "poly 1 x = 1"
       
  1049   by (simp add: one_poly_def)
       
  1050 
       
  1051 lemma poly_power [simp]: "poly (p ^ n) x = poly p x ^ n"
       
  1052   for p :: "'a::comm_semiring_1 poly"
       
  1053   by (induct n) simp_all
       
  1054 
       
  1055 lemma poly_prod: "poly (\<Prod>k\<in>A. p k) x = (\<Prod>k\<in>A. poly (p k) x)"
       
  1056   by (induct A rule: infinite_finite_induct) simp_all
       
  1057 
       
  1058 lemma degree_prod_sum_le: "finite S \<Longrightarrow> degree (prod f S) \<le> sum (degree o f) S"
       
  1059 proof (induct S rule: finite_induct)
       
  1060   case empty
       
  1061   then show ?case by simp
       
  1062 next
       
  1063   case (insert a S)
       
  1064   show ?case
       
  1065     unfolding prod.insert[OF insert(1-2)] sum.insert[OF insert(1-2)]
       
  1066     by (rule le_trans[OF degree_mult_le]) (use insert in auto)
       
  1067 qed
       
  1068 
       
  1069 lemma coeff_0_prod_list: "coeff (prod_list xs) 0 = prod_list (map (\<lambda>p. coeff p 0) xs)"
       
  1070   by (induct xs) (simp_all add: coeff_mult)
       
  1071 
       
  1072 lemma coeff_monom_mult: "coeff (monom c n * p) k = (if k < n then 0 else c * coeff p (k - n))"
       
  1073 proof -
       
  1074   have "coeff (monom c n * p) k = (\<Sum>i\<le>k. (if n = i then c else 0) * coeff p (k - i))"
       
  1075     by (simp add: coeff_mult)
       
  1076   also have "\<dots> = (\<Sum>i\<le>k. (if n = i then c * coeff p (k - i) else 0))"
       
  1077     by (intro sum.cong) simp_all
       
  1078   also have "\<dots> = (if k < n then 0 else c * coeff p (k - n))"
       
  1079     by (simp add: sum.delta')
       
  1080   finally show ?thesis .
       
  1081 qed
       
  1082 
       
  1083 lemma monom_1_dvd_iff': "monom 1 n dvd p \<longleftrightarrow> (\<forall>k<n. coeff p k = 0)"
       
  1084 proof
       
  1085   assume "monom 1 n dvd p"
       
  1086   then obtain r where "p = monom 1 n * r"
       
  1087     by (rule dvdE)
       
  1088   then show "\<forall>k<n. coeff p k = 0"
       
  1089     by (simp add: coeff_mult)
       
  1090 next
       
  1091   assume zero: "(\<forall>k<n. coeff p k = 0)"
       
  1092   define r where "r = Abs_poly (\<lambda>k. coeff p (k + n))"
       
  1093   have "\<forall>\<^sub>\<infinity>k. coeff p (k + n) = 0"
       
  1094     by (subst cofinite_eq_sequentially, subst eventually_sequentially_seg,
       
  1095         subst cofinite_eq_sequentially [symmetric]) transfer
       
  1096   then have coeff_r [simp]: "coeff r k = coeff p (k + n)" for k
       
  1097     unfolding r_def by (subst poly.Abs_poly_inverse) simp_all
       
  1098   have "p = monom 1 n * r"
       
  1099     by (rule poly_eqI, subst coeff_monom_mult) (simp_all add: zero)
       
  1100   then show "monom 1 n dvd p" by simp
       
  1101 qed
       
  1102 
       
  1103 
       
  1104 subsection \<open>Mapping polynomials\<close>
       
  1105 
       
  1106 definition map_poly :: "('a :: zero \<Rightarrow> 'b :: zero) \<Rightarrow> 'a poly \<Rightarrow> 'b poly"
       
  1107   where "map_poly f p = Poly (map f (coeffs p))"
       
  1108 
       
  1109 lemma map_poly_0 [simp]: "map_poly f 0 = 0"
       
  1110   by (simp add: map_poly_def)
       
  1111 
       
  1112 lemma map_poly_1: "map_poly f 1 = [:f 1:]"
       
  1113   by (simp add: map_poly_def)
       
  1114 
       
  1115 lemma map_poly_1' [simp]: "f 1 = 1 \<Longrightarrow> map_poly f 1 = 1"
       
  1116   by (simp add: map_poly_def one_poly_def)
       
  1117 
       
  1118 lemma coeff_map_poly:
       
  1119   assumes "f 0 = 0"
       
  1120   shows "coeff (map_poly f p) n = f (coeff p n)"
       
  1121   by (auto simp: assms map_poly_def nth_default_def coeffs_def not_less Suc_le_eq coeff_eq_0
       
  1122       simp del: upt_Suc)
       
  1123 
       
  1124 lemma coeffs_map_poly [code abstract]:
       
  1125   "coeffs (map_poly f p) = strip_while (op = 0) (map f (coeffs p))"
       
  1126   by (simp add: map_poly_def)
       
  1127 
       
  1128 lemma coeffs_map_poly':
       
  1129   assumes "\<And>x. x \<noteq> 0 \<Longrightarrow> f x \<noteq> 0"
       
  1130   shows "coeffs (map_poly f p) = map f (coeffs p)"
       
  1131   using assms by (simp add: coeffs_map_poly no_trailing_map strip_while_idem_iff)
       
  1132     (metis comp_def no_leading_def no_trailing_coeffs)
       
  1133 
       
  1134 lemma set_coeffs_map_poly:
       
  1135   "(\<And>x. f x = 0 \<longleftrightarrow> x = 0) \<Longrightarrow> set (coeffs (map_poly f p)) = f ` set (coeffs p)"
       
  1136   by (simp add: coeffs_map_poly')
       
  1137 
       
  1138 lemma degree_map_poly:
       
  1139   assumes "\<And>x. x \<noteq> 0 \<Longrightarrow> f x \<noteq> 0"
       
  1140   shows "degree (map_poly f p) = degree p"
       
  1141   by (simp add: degree_eq_length_coeffs coeffs_map_poly' assms)
       
  1142 
       
  1143 lemma map_poly_eq_0_iff:
       
  1144   assumes "f 0 = 0" "\<And>x. x \<in> set (coeffs p) \<Longrightarrow> x \<noteq> 0 \<Longrightarrow> f x \<noteq> 0"
       
  1145   shows "map_poly f p = 0 \<longleftrightarrow> p = 0"
       
  1146 proof -
       
  1147   have "(coeff (map_poly f p) n = 0) = (coeff p n = 0)" for n
       
  1148   proof -
       
  1149     have "coeff (map_poly f p) n = f (coeff p n)"
       
  1150       by (simp add: coeff_map_poly assms)
       
  1151     also have "\<dots> = 0 \<longleftrightarrow> coeff p n = 0"
       
  1152     proof (cases "n < length (coeffs p)")
       
  1153       case True
       
  1154       then have "coeff p n \<in> set (coeffs p)"
       
  1155         by (auto simp: coeffs_def simp del: upt_Suc)
       
  1156       with assms show "f (coeff p n) = 0 \<longleftrightarrow> coeff p n = 0"
       
  1157         by auto
       
  1158     next
       
  1159       case False
       
  1160       then show ?thesis
       
  1161         by (auto simp: assms length_coeffs nth_default_coeffs_eq [symmetric] nth_default_def)
       
  1162     qed
       
  1163     finally show ?thesis .
       
  1164   qed
       
  1165   then show ?thesis by (auto simp: poly_eq_iff)
       
  1166 qed
       
  1167 
       
  1168 lemma map_poly_smult:
       
  1169   assumes "f 0 = 0""\<And>c x. f (c * x) = f c * f x"
       
  1170   shows "map_poly f (smult c p) = smult (f c) (map_poly f p)"
       
  1171   by (intro poly_eqI) (simp_all add: assms coeff_map_poly)
       
  1172 
       
  1173 lemma map_poly_pCons:
       
  1174   assumes "f 0 = 0"
       
  1175   shows "map_poly f (pCons c p) = pCons (f c) (map_poly f p)"
       
  1176   by (intro poly_eqI) (simp_all add: assms coeff_map_poly coeff_pCons split: nat.splits)
       
  1177 
       
  1178 lemma map_poly_map_poly:
       
  1179   assumes "f 0 = 0" "g 0 = 0"
       
  1180   shows "map_poly f (map_poly g p) = map_poly (f \<circ> g) p"
       
  1181   by (intro poly_eqI) (simp add: coeff_map_poly assms)
       
  1182 
       
  1183 lemma map_poly_id [simp]: "map_poly id p = p"
       
  1184   by (simp add: map_poly_def)
       
  1185 
       
  1186 lemma map_poly_id' [simp]: "map_poly (\<lambda>x. x) p = p"
       
  1187   by (simp add: map_poly_def)
       
  1188 
       
  1189 lemma map_poly_cong:
       
  1190   assumes "(\<And>x. x \<in> set (coeffs p) \<Longrightarrow> f x = g x)"
       
  1191   shows "map_poly f p = map_poly g p"
       
  1192 proof -
       
  1193   from assms have "map f (coeffs p) = map g (coeffs p)"
       
  1194     by (intro map_cong) simp_all
       
  1195   then show ?thesis
       
  1196     by (simp only: coeffs_eq_iff coeffs_map_poly)
       
  1197 qed
       
  1198 
       
  1199 lemma map_poly_monom: "f 0 = 0 \<Longrightarrow> map_poly f (monom c n) = monom (f c) n"
       
  1200   by (intro poly_eqI) (simp_all add: coeff_map_poly)
       
  1201 
       
  1202 lemma map_poly_idI:
       
  1203   assumes "\<And>x. x \<in> set (coeffs p) \<Longrightarrow> f x = x"
       
  1204   shows "map_poly f p = p"
       
  1205   using map_poly_cong[OF assms, of _ id] by simp
       
  1206 
       
  1207 lemma map_poly_idI':
       
  1208   assumes "\<And>x. x \<in> set (coeffs p) \<Longrightarrow> f x = x"
       
  1209   shows "p = map_poly f p"
       
  1210   using map_poly_cong[OF assms, of _ id] by simp
       
  1211 
       
  1212 lemma smult_conv_map_poly: "smult c p = map_poly (\<lambda>x. c * x) p"
       
  1213   by (intro poly_eqI) (simp_all add: coeff_map_poly)
       
  1214 
       
  1215 
       
  1216 subsection \<open>Conversions from @{typ nat} and @{typ int} numbers\<close>
       
  1217 
       
  1218 lemma of_nat_poly: "of_nat n = [:of_nat n :: 'a :: comm_semiring_1:]"
       
  1219 proof (induct n)
       
  1220   case 0
       
  1221   then show ?case by simp
       
  1222 next
       
  1223   case (Suc n)
       
  1224   then have "of_nat (Suc n) = 1 + (of_nat n :: 'a poly)"
       
  1225     by simp
       
  1226   also have "(of_nat n :: 'a poly) = [: of_nat n :]"
       
  1227     by (subst Suc) (rule refl)
       
  1228   also have "1 = [:1:]"
       
  1229     by (simp add: one_poly_def)
       
  1230   finally show ?case
       
  1231     by (subst (asm) add_pCons) simp
       
  1232 qed
       
  1233 
       
  1234 lemma degree_of_nat [simp]: "degree (of_nat n) = 0"
       
  1235   by (simp add: of_nat_poly)
       
  1236 
       
  1237 lemma lead_coeff_of_nat [simp]:
       
  1238   "lead_coeff (of_nat n) = (of_nat n :: 'a::{comm_semiring_1,semiring_char_0})"
       
  1239   by (simp add: of_nat_poly)
       
  1240 
       
  1241 lemma of_int_poly: "of_int k = [:of_int k :: 'a :: comm_ring_1:]"
       
  1242   by (simp only: of_int_of_nat of_nat_poly) simp
       
  1243 
       
  1244 lemma degree_of_int [simp]: "degree (of_int k) = 0"
       
  1245   by (simp add: of_int_poly)
       
  1246 
       
  1247 lemma lead_coeff_of_int [simp]:
       
  1248   "lead_coeff (of_int k) = (of_int k :: 'a::{comm_ring_1,ring_char_0})"
       
  1249   by (simp add: of_int_poly)
       
  1250 
       
  1251 lemma numeral_poly: "numeral n = [:numeral n:]"
       
  1252   by (subst of_nat_numeral [symmetric], subst of_nat_poly) simp
       
  1253 
       
  1254 lemma degree_numeral [simp]: "degree (numeral n) = 0"
       
  1255   by (subst of_nat_numeral [symmetric], subst of_nat_poly) simp
       
  1256 
       
  1257 lemma lead_coeff_numeral [simp]:
       
  1258   "lead_coeff (numeral n) = numeral n"
       
  1259   by (simp add: numeral_poly)
       
  1260 
       
  1261 
       
  1262 subsection \<open>Lemmas about divisibility\<close>
       
  1263 
       
  1264 lemma dvd_smult:
       
  1265   assumes "p dvd q"
       
  1266   shows "p dvd smult a q"
       
  1267 proof -
       
  1268   from assms obtain k where "q = p * k" ..
       
  1269   then have "smult a q = p * smult a k" by simp
       
  1270   then show "p dvd smult a q" ..
       
  1271 qed
       
  1272 
       
  1273 lemma dvd_smult_cancel: "p dvd smult a q \<Longrightarrow> a \<noteq> 0 \<Longrightarrow> p dvd q"
       
  1274   for a :: "'a::field"
       
  1275   by (drule dvd_smult [where a="inverse a"]) simp
       
  1276 
       
  1277 lemma dvd_smult_iff: "a \<noteq> 0 \<Longrightarrow> p dvd smult a q \<longleftrightarrow> p dvd q"
       
  1278   for a :: "'a::field"
       
  1279   by (safe elim!: dvd_smult dvd_smult_cancel)
       
  1280 
       
  1281 lemma smult_dvd_cancel:
       
  1282   assumes "smult a p dvd q"
       
  1283   shows "p dvd q"
       
  1284 proof -
       
  1285   from assms obtain k where "q = smult a p * k" ..
       
  1286   then have "q = p * smult a k" by simp
       
  1287   then show "p dvd q" ..
       
  1288 qed
       
  1289 
       
  1290 lemma smult_dvd: "p dvd q \<Longrightarrow> a \<noteq> 0 \<Longrightarrow> smult a p dvd q"
       
  1291   for a :: "'a::field"
       
  1292   by (rule smult_dvd_cancel [where a="inverse a"]) simp
       
  1293 
       
  1294 lemma smult_dvd_iff: "smult a p dvd q \<longleftrightarrow> (if a = 0 then q = 0 else p dvd q)"
       
  1295   for a :: "'a::field"
       
  1296   by (auto elim: smult_dvd smult_dvd_cancel)
       
  1297 
       
  1298 lemma is_unit_smult_iff: "smult c p dvd 1 \<longleftrightarrow> c dvd 1 \<and> p dvd 1"
       
  1299 proof -
       
  1300   have "smult c p = [:c:] * p" by simp
       
  1301   also have "\<dots> dvd 1 \<longleftrightarrow> c dvd 1 \<and> p dvd 1"
       
  1302   proof safe
       
  1303     assume *: "[:c:] * p dvd 1"
       
  1304     then show "p dvd 1"
       
  1305       by (rule dvd_mult_right)
       
  1306     from * obtain q where q: "1 = [:c:] * p * q"
       
  1307       by (rule dvdE)
       
  1308     have "c dvd c * (coeff p 0 * coeff q 0)"
       
  1309       by simp
       
  1310     also have "\<dots> = coeff ([:c:] * p * q) 0"
       
  1311       by (simp add: mult.assoc coeff_mult)
       
  1312     also note q [symmetric]
       
  1313     finally have "c dvd coeff 1 0" .
       
  1314     then show "c dvd 1" by simp
       
  1315   next
       
  1316     assume "c dvd 1" "p dvd 1"
       
  1317     from this(1) obtain d where "1 = c * d"
       
  1318       by (rule dvdE)
       
  1319     then have "1 = [:c:] * [:d:]"
       
  1320       by (simp add: one_poly_def mult_ac)
       
  1321     then have "[:c:] dvd 1"
       
  1322       by (rule dvdI)
       
  1323     from mult_dvd_mono[OF this \<open>p dvd 1\<close>] show "[:c:] * p dvd 1"
       
  1324       by simp
       
  1325   qed
       
  1326   finally show ?thesis .
       
  1327 qed
       
  1328 
       
  1329 
       
  1330 subsection \<open>Polynomials form an integral domain\<close>
       
  1331 
       
  1332 instance poly :: (idom) idom ..
       
  1333 
       
  1334 lemma degree_mult_eq: "p \<noteq> 0 \<Longrightarrow> q \<noteq> 0 \<Longrightarrow> degree (p * q) = degree p + degree q"
       
  1335   for p q :: "'a::{comm_semiring_0,semiring_no_zero_divisors} poly"
       
  1336   by (rule order_antisym [OF degree_mult_le le_degree]) (simp add: coeff_mult_degree_sum)
       
  1337 
       
  1338 lemma degree_mult_eq_0:
       
  1339   "degree (p * q) = 0 \<longleftrightarrow> p = 0 \<or> q = 0 \<or> (p \<noteq> 0 \<and> q \<noteq> 0 \<and> degree p = 0 \<and> degree q = 0)"
       
  1340   for p q :: "'a::{comm_semiring_0,semiring_no_zero_divisors} poly"
       
  1341   by (auto simp: degree_mult_eq)
       
  1342 
       
  1343 lemma degree_mult_right_le:
       
  1344   fixes p q :: "'a::{comm_semiring_0,semiring_no_zero_divisors} poly"
       
  1345   assumes "q \<noteq> 0"
       
  1346   shows "degree p \<le> degree (p * q)"
       
  1347   using assms by (cases "p = 0") (simp_all add: degree_mult_eq)
       
  1348 
       
  1349 lemma coeff_degree_mult: "coeff (p * q) (degree (p * q)) = coeff q (degree q) * coeff p (degree p)"
       
  1350   for p q :: "'a::{comm_semiring_0,semiring_no_zero_divisors} poly"
       
  1351   by (cases "p = 0 \<or> q = 0") (auto simp: degree_mult_eq coeff_mult_degree_sum mult_ac)
       
  1352 
       
  1353 lemma dvd_imp_degree_le: "p dvd q \<Longrightarrow> q \<noteq> 0 \<Longrightarrow> degree p \<le> degree q"
       
  1354   for p q :: "'a::{comm_semiring_1,semiring_no_zero_divisors} poly"
       
  1355   by (erule dvdE, hypsubst, subst degree_mult_eq) auto
       
  1356 
       
  1357 lemma divides_degree:
       
  1358   fixes p q :: "'a ::{comm_semiring_1,semiring_no_zero_divisors} poly"
       
  1359   assumes "p dvd q"
       
  1360   shows "degree p \<le> degree q \<or> q = 0"
       
  1361   by (metis dvd_imp_degree_le assms)
       
  1362 
       
  1363 lemma const_poly_dvd_iff:
       
  1364   fixes c :: "'a::{comm_semiring_1,semiring_no_zero_divisors}"
       
  1365   shows "[:c:] dvd p \<longleftrightarrow> (\<forall>n. c dvd coeff p n)"
       
  1366 proof (cases "c = 0 \<or> p = 0")
       
  1367   case True
       
  1368   then show ?thesis
       
  1369     by (auto intro!: poly_eqI)
       
  1370 next
       
  1371   case False
       
  1372   show ?thesis
       
  1373   proof
       
  1374     assume "[:c:] dvd p"
       
  1375     then show "\<forall>n. c dvd coeff p n"
       
  1376       by (auto elim!: dvdE simp: coeffs_def)
       
  1377   next
       
  1378     assume *: "\<forall>n. c dvd coeff p n"
       
  1379     define mydiv where "mydiv x y = (SOME z. x = y * z)" for x y :: 'a
       
  1380     have mydiv: "x = y * mydiv x y" if "y dvd x" for x y
       
  1381       using that unfolding mydiv_def dvd_def by (rule someI_ex)
       
  1382     define q where "q = Poly (map (\<lambda>a. mydiv a c) (coeffs p))"
       
  1383     from False * have "p = q * [:c:]"
       
  1384       by (intro poly_eqI)
       
  1385         (auto simp: q_def nth_default_def not_less length_coeffs_degree coeffs_nth
       
  1386           intro!: coeff_eq_0 mydiv)
       
  1387     then show "[:c:] dvd p"
       
  1388       by (simp only: dvd_triv_right)
       
  1389   qed
       
  1390 qed
       
  1391 
       
  1392 lemma const_poly_dvd_const_poly_iff [simp]: "[:a:] dvd [:b:] \<longleftrightarrow> a dvd b"
       
  1393   for a b :: "'a::{comm_semiring_1,semiring_no_zero_divisors}"
       
  1394   by (subst const_poly_dvd_iff) (auto simp: coeff_pCons split: nat.splits)
       
  1395 
       
  1396 lemma lead_coeff_mult: "lead_coeff (p * q) = lead_coeff p * lead_coeff q"
       
  1397   for p q :: "'a::{comm_semiring_0, semiring_no_zero_divisors} poly"
       
  1398   by (cases "p = 0 \<or> q = 0") (auto simp: coeff_mult_degree_sum degree_mult_eq)
       
  1399 
       
  1400 lemma lead_coeff_smult: "lead_coeff (smult c p) = c * lead_coeff p"
       
  1401   for p :: "'a::{comm_semiring_0,semiring_no_zero_divisors} poly"
       
  1402 proof -
       
  1403   have "smult c p = [:c:] * p" by simp
       
  1404   also have "lead_coeff \<dots> = c * lead_coeff p"
       
  1405     by (subst lead_coeff_mult) simp_all
       
  1406   finally show ?thesis .
       
  1407 qed
       
  1408 
       
  1409 lemma lead_coeff_1 [simp]: "lead_coeff 1 = 1"
       
  1410   by simp
       
  1411 
       
  1412 lemma lead_coeff_power: "lead_coeff (p ^ n) = lead_coeff p ^ n"
       
  1413   for p :: "'a::{comm_semiring_1,semiring_no_zero_divisors} poly"
       
  1414   by (induct n) (simp_all add: lead_coeff_mult)
       
  1415 
       
  1416 
       
  1417 subsection \<open>Polynomials form an ordered integral domain\<close>
       
  1418 
       
  1419 definition pos_poly :: "'a::linordered_semidom poly \<Rightarrow> bool"
       
  1420   where "pos_poly p \<longleftrightarrow> 0 < coeff p (degree p)"
       
  1421 
       
  1422 lemma pos_poly_pCons: "pos_poly (pCons a p) \<longleftrightarrow> pos_poly p \<or> (p = 0 \<and> 0 < a)"
       
  1423   by (simp add: pos_poly_def)
       
  1424 
       
  1425 lemma not_pos_poly_0 [simp]: "\<not> pos_poly 0"
       
  1426   by (simp add: pos_poly_def)
       
  1427 
       
  1428 lemma pos_poly_add: "pos_poly p \<Longrightarrow> pos_poly q \<Longrightarrow> pos_poly (p + q)"
       
  1429   apply (induct p arbitrary: q)
       
  1430    apply simp
       
  1431   apply (case_tac q)
       
  1432   apply (force simp add: pos_poly_pCons add_pos_pos)
       
  1433   done
       
  1434 
       
  1435 lemma pos_poly_mult: "pos_poly p \<Longrightarrow> pos_poly q \<Longrightarrow> pos_poly (p * q)"
       
  1436   unfolding pos_poly_def
       
  1437   apply (subgoal_tac "p \<noteq> 0 \<and> q \<noteq> 0")
       
  1438    apply (simp add: degree_mult_eq coeff_mult_degree_sum)
       
  1439   apply auto
       
  1440   done
       
  1441 
       
  1442 lemma pos_poly_total: "p = 0 \<or> pos_poly p \<or> pos_poly (- p)"
       
  1443   for p :: "'a::linordered_idom poly"
       
  1444   by (induct p) (auto simp: pos_poly_pCons)
       
  1445 
       
  1446 lemma last_coeffs_eq_coeff_degree: "p \<noteq> 0 \<Longrightarrow> last (coeffs p) = coeff p (degree p)"
       
  1447   by (simp add: coeffs_def)
       
  1448 
       
  1449 lemma pos_poly_coeffs [code]: "pos_poly p \<longleftrightarrow> (let as = coeffs p in as \<noteq> [] \<and> last as > 0)"
       
  1450   (is "?lhs \<longleftrightarrow> ?rhs")
       
  1451 proof
       
  1452   assume ?rhs
       
  1453   then show ?lhs
       
  1454     by (auto simp add: pos_poly_def last_coeffs_eq_coeff_degree)
       
  1455 next
       
  1456   assume ?lhs
       
  1457   then have *: "0 < coeff p (degree p)"
       
  1458     by (simp add: pos_poly_def)
       
  1459   then have "p \<noteq> 0"
       
  1460     by auto
       
  1461   with * show ?rhs
       
  1462     by (simp add: last_coeffs_eq_coeff_degree)
       
  1463 qed
       
  1464 
       
  1465 instantiation poly :: (linordered_idom) linordered_idom
       
  1466 begin
       
  1467 
       
  1468 definition "x < y \<longleftrightarrow> pos_poly (y - x)"
       
  1469 
       
  1470 definition "x \<le> y \<longleftrightarrow> x = y \<or> pos_poly (y - x)"
       
  1471 
       
  1472 definition "\<bar>x::'a poly\<bar> = (if x < 0 then - x else x)"
       
  1473 
       
  1474 definition "sgn (x::'a poly) = (if x = 0 then 0 else if 0 < x then 1 else - 1)"
       
  1475 
       
  1476 instance
       
  1477 proof
       
  1478   fix x y z :: "'a poly"
       
  1479   show "x < y \<longleftrightarrow> x \<le> y \<and> \<not> y \<le> x"
       
  1480     unfolding less_eq_poly_def less_poly_def
       
  1481     apply safe
       
  1482      apply simp
       
  1483     apply (drule (1) pos_poly_add)
       
  1484     apply simp
       
  1485     done
       
  1486   show "x \<le> x"
       
  1487     by (simp add: less_eq_poly_def)
       
  1488   show "x \<le> y \<Longrightarrow> y \<le> z \<Longrightarrow> x \<le> z"
       
  1489     unfolding less_eq_poly_def
       
  1490     apply safe
       
  1491     apply (drule (1) pos_poly_add)
       
  1492     apply (simp add: algebra_simps)
       
  1493     done
       
  1494   show "x \<le> y \<Longrightarrow> y \<le> x \<Longrightarrow> x = y"
       
  1495     unfolding less_eq_poly_def
       
  1496     apply safe
       
  1497     apply (drule (1) pos_poly_add)
       
  1498     apply simp
       
  1499     done
       
  1500   show "x \<le> y \<Longrightarrow> z + x \<le> z + y"
       
  1501     unfolding less_eq_poly_def
       
  1502     apply safe
       
  1503     apply (simp add: algebra_simps)
       
  1504     done
       
  1505   show "x \<le> y \<or> y \<le> x"
       
  1506     unfolding less_eq_poly_def
       
  1507     using pos_poly_total [of "x - y"]
       
  1508     by auto
       
  1509   show "x < y \<Longrightarrow> 0 < z \<Longrightarrow> z * x < z * y"
       
  1510     by (simp add: less_poly_def right_diff_distrib [symmetric] pos_poly_mult)
       
  1511   show "\<bar>x\<bar> = (if x < 0 then - x else x)"
       
  1512     by (rule abs_poly_def)
       
  1513   show "sgn x = (if x = 0 then 0 else if 0 < x then 1 else - 1)"
       
  1514     by (rule sgn_poly_def)
       
  1515 qed
       
  1516 
       
  1517 end
       
  1518 
       
  1519 text \<open>TODO: Simplification rules for comparisons\<close>
       
  1520 
       
  1521 
       
  1522 subsection \<open>Synthetic division and polynomial roots\<close>
       
  1523 
       
  1524 subsubsection \<open>Synthetic division\<close>
       
  1525 
       
  1526 text \<open>Synthetic division is simply division by the linear polynomial @{term "x - c"}.\<close>
       
  1527 
       
  1528 definition synthetic_divmod :: "'a::comm_semiring_0 poly \<Rightarrow> 'a \<Rightarrow> 'a poly \<times> 'a"
       
  1529   where "synthetic_divmod p c = fold_coeffs (\<lambda>a (q, r). (pCons r q, a + c * r)) p (0, 0)"
       
  1530 
       
  1531 definition synthetic_div :: "'a::comm_semiring_0 poly \<Rightarrow> 'a \<Rightarrow> 'a poly"
       
  1532   where "synthetic_div p c = fst (synthetic_divmod p c)"
       
  1533 
       
  1534 lemma synthetic_divmod_0 [simp]: "synthetic_divmod 0 c = (0, 0)"
       
  1535   by (simp add: synthetic_divmod_def)
       
  1536 
       
  1537 lemma synthetic_divmod_pCons [simp]:
       
  1538   "synthetic_divmod (pCons a p) c = (\<lambda>(q, r). (pCons r q, a + c * r)) (synthetic_divmod p c)"
       
  1539   by (cases "p = 0 \<and> a = 0") (auto simp add: synthetic_divmod_def)
       
  1540 
       
  1541 lemma synthetic_div_0 [simp]: "synthetic_div 0 c = 0"
       
  1542   by (simp add: synthetic_div_def)
       
  1543 
       
  1544 lemma synthetic_div_unique_lemma: "smult c p = pCons a p \<Longrightarrow> p = 0"
       
  1545   by (induct p arbitrary: a) simp_all
       
  1546 
       
  1547 lemma snd_synthetic_divmod: "snd (synthetic_divmod p c) = poly p c"
       
  1548   by (induct p) (simp_all add: split_def)
       
  1549 
       
  1550 lemma synthetic_div_pCons [simp]:
       
  1551   "synthetic_div (pCons a p) c = pCons (poly p c) (synthetic_div p c)"
       
  1552   by (simp add: synthetic_div_def split_def snd_synthetic_divmod)
       
  1553 
       
  1554 lemma synthetic_div_eq_0_iff: "synthetic_div p c = 0 \<longleftrightarrow> degree p = 0"
       
  1555 proof (induct p)
       
  1556   case 0
       
  1557   then show ?case by simp
       
  1558 next
       
  1559   case (pCons a p)
       
  1560   then show ?case by (cases p) simp
       
  1561 qed
       
  1562 
       
  1563 lemma degree_synthetic_div: "degree (synthetic_div p c) = degree p - 1"
       
  1564   by (induct p) (simp_all add: synthetic_div_eq_0_iff)
       
  1565 
       
  1566 lemma synthetic_div_correct:
       
  1567   "p + smult c (synthetic_div p c) = pCons (poly p c) (synthetic_div p c)"
       
  1568   by (induct p) simp_all
       
  1569 
       
  1570 lemma synthetic_div_unique: "p + smult c q = pCons r q \<Longrightarrow> r = poly p c \<and> q = synthetic_div p c"
       
  1571   apply (induct p arbitrary: q r)
       
  1572    apply simp
       
  1573    apply (frule synthetic_div_unique_lemma)
       
  1574    apply simp
       
  1575   apply (case_tac q, force)
       
  1576   done
       
  1577 
       
  1578 lemma synthetic_div_correct': "[:-c, 1:] * synthetic_div p c + [:poly p c:] = p"
       
  1579   for c :: "'a::comm_ring_1"
       
  1580   using synthetic_div_correct [of p c] by (simp add: algebra_simps)
       
  1581 
       
  1582 
       
  1583 subsubsection \<open>Polynomial roots\<close>
       
  1584 
       
  1585 lemma poly_eq_0_iff_dvd: "poly p c = 0 \<longleftrightarrow> [:- c, 1:] dvd p"
       
  1586   (is "?lhs \<longleftrightarrow> ?rhs")
       
  1587   for c :: "'a::comm_ring_1"
       
  1588 proof
       
  1589   assume ?lhs
       
  1590   with synthetic_div_correct' [of c p] have "p = [:-c, 1:] * synthetic_div p c" by simp
       
  1591   then show ?rhs ..
       
  1592 next
       
  1593   assume ?rhs
       
  1594   then obtain k where "p = [:-c, 1:] * k" by (rule dvdE)
       
  1595   then show ?lhs by simp
       
  1596 qed
       
  1597 
       
  1598 lemma dvd_iff_poly_eq_0: "[:c, 1:] dvd p \<longleftrightarrow> poly p (- c) = 0"
       
  1599   for c :: "'a::comm_ring_1"
       
  1600   by (simp add: poly_eq_0_iff_dvd)
       
  1601 
       
  1602 lemma poly_roots_finite: "p \<noteq> 0 \<Longrightarrow> finite {x. poly p x = 0}"
       
  1603   for p :: "'a::{comm_ring_1,ring_no_zero_divisors} poly"
       
  1604 proof (induct n \<equiv> "degree p" arbitrary: p)
       
  1605   case 0
       
  1606   then obtain a where "a \<noteq> 0" and "p = [:a:]"
       
  1607     by (cases p) (simp split: if_splits)
       
  1608   then show "finite {x. poly p x = 0}"
       
  1609     by simp
       
  1610 next
       
  1611   case (Suc n)
       
  1612   show "finite {x. poly p x = 0}"
       
  1613   proof (cases "\<exists>x. poly p x = 0")
       
  1614     case False
       
  1615     then show "finite {x. poly p x = 0}" by simp
       
  1616   next
       
  1617     case True
       
  1618     then obtain a where "poly p a = 0" ..
       
  1619     then have "[:-a, 1:] dvd p"
       
  1620       by (simp only: poly_eq_0_iff_dvd)
       
  1621     then obtain k where k: "p = [:-a, 1:] * k" ..
       
  1622     with \<open>p \<noteq> 0\<close> have "k \<noteq> 0"
       
  1623       by auto
       
  1624     with k have "degree p = Suc (degree k)"
       
  1625       by (simp add: degree_mult_eq del: mult_pCons_left)
       
  1626     with \<open>Suc n = degree p\<close> have "n = degree k"
       
  1627       by simp
       
  1628     from this \<open>k \<noteq> 0\<close> have "finite {x. poly k x = 0}"
       
  1629       by (rule Suc.hyps)
       
  1630     then have "finite (insert a {x. poly k x = 0})"
       
  1631       by simp
       
  1632     then show "finite {x. poly p x = 0}"
       
  1633       by (simp add: k Collect_disj_eq del: mult_pCons_left)
       
  1634   qed
       
  1635 qed
       
  1636 
       
  1637 lemma poly_eq_poly_eq_iff: "poly p = poly q \<longleftrightarrow> p = q"
       
  1638   (is "?lhs \<longleftrightarrow> ?rhs")
       
  1639   for p q :: "'a::{comm_ring_1,ring_no_zero_divisors,ring_char_0} poly"
       
  1640 proof
       
  1641   assume ?rhs
       
  1642   then show ?lhs by simp
       
  1643 next
       
  1644   assume ?lhs
       
  1645   have "poly p = poly 0 \<longleftrightarrow> p = 0" for p :: "'a poly"
       
  1646     apply (cases "p = 0")
       
  1647      apply simp_all
       
  1648     apply (drule poly_roots_finite)
       
  1649     apply (auto simp add: infinite_UNIV_char_0)
       
  1650     done
       
  1651   from \<open>?lhs\<close> and this [of "p - q"] show ?rhs
       
  1652     by auto
       
  1653 qed
       
  1654 
       
  1655 lemma poly_all_0_iff_0: "(\<forall>x. poly p x = 0) \<longleftrightarrow> p = 0"
       
  1656   for p :: "'a::{ring_char_0,comm_ring_1,ring_no_zero_divisors} poly"
       
  1657   by (auto simp add: poly_eq_poly_eq_iff [symmetric])
       
  1658 
       
  1659 
       
  1660 subsubsection \<open>Order of polynomial roots\<close>
       
  1661 
       
  1662 definition order :: "'a::idom \<Rightarrow> 'a poly \<Rightarrow> nat"
       
  1663   where "order a p = (LEAST n. \<not> [:-a, 1:] ^ Suc n dvd p)"
       
  1664 
       
  1665 lemma coeff_linear_power: "coeff ([:a, 1:] ^ n) n = 1"
       
  1666   for a :: "'a::comm_semiring_1"
       
  1667   apply (induct n)
       
  1668    apply simp_all
       
  1669   apply (subst coeff_eq_0)
       
  1670    apply (auto intro: le_less_trans degree_power_le)
       
  1671   done
       
  1672 
       
  1673 lemma degree_linear_power: "degree ([:a, 1:] ^ n) = n"
       
  1674   for a :: "'a::comm_semiring_1"
       
  1675   apply (rule order_antisym)
       
  1676    apply (rule ord_le_eq_trans [OF degree_power_le])
       
  1677    apply simp
       
  1678   apply (rule le_degree)
       
  1679   apply (simp add: coeff_linear_power)
       
  1680   done
       
  1681 
       
  1682 lemma order_1: "[:-a, 1:] ^ order a p dvd p"
       
  1683   apply (cases "p = 0")
       
  1684    apply simp
       
  1685   apply (cases "order a p")
       
  1686    apply simp
       
  1687   apply (subgoal_tac "nat < (LEAST n. \<not> [:-a, 1:] ^ Suc n dvd p)")
       
  1688    apply (drule not_less_Least)
       
  1689    apply simp
       
  1690   apply (fold order_def)
       
  1691   apply simp
       
  1692   done
       
  1693 
       
  1694 lemma order_2: "p \<noteq> 0 \<Longrightarrow> \<not> [:-a, 1:] ^ Suc (order a p) dvd p"
       
  1695   unfolding order_def
       
  1696   apply (rule LeastI_ex)
       
  1697   apply (rule_tac x="degree p" in exI)
       
  1698   apply (rule notI)
       
  1699   apply (drule (1) dvd_imp_degree_le)
       
  1700   apply (simp only: degree_linear_power)
       
  1701   done
       
  1702 
       
  1703 lemma order: "p \<noteq> 0 \<Longrightarrow> [:-a, 1:] ^ order a p dvd p \<and> \<not> [:-a, 1:] ^ Suc (order a p) dvd p"
       
  1704   by (rule conjI [OF order_1 order_2])
       
  1705 
       
  1706 lemma order_degree:
       
  1707   assumes p: "p \<noteq> 0"
       
  1708   shows "order a p \<le> degree p"
       
  1709 proof -
       
  1710   have "order a p = degree ([:-a, 1:] ^ order a p)"
       
  1711     by (simp only: degree_linear_power)
       
  1712   also from order_1 p have "\<dots> \<le> degree p"
       
  1713     by (rule dvd_imp_degree_le)
       
  1714   finally show ?thesis .
       
  1715 qed
       
  1716 
       
  1717 lemma order_root: "poly p a = 0 \<longleftrightarrow> p = 0 \<or> order a p \<noteq> 0"
       
  1718   apply (cases "p = 0")
       
  1719    apply simp_all
       
  1720   apply (rule iffI)
       
  1721    apply (metis order_2 not_gr0 poly_eq_0_iff_dvd power_0 power_Suc_0 power_one_right)
       
  1722   unfolding poly_eq_0_iff_dvd
       
  1723   apply (metis dvd_power dvd_trans order_1)
       
  1724   done
       
  1725 
       
  1726 lemma order_0I: "poly p a \<noteq> 0 \<Longrightarrow> order a p = 0"
       
  1727   by (subst (asm) order_root) auto
       
  1728 
       
  1729 lemma order_unique_lemma:
       
  1730   fixes p :: "'a::idom poly"
       
  1731   assumes "[:-a, 1:] ^ n dvd p" "\<not> [:-a, 1:] ^ Suc n dvd p"
       
  1732   shows "n = order a p"
       
  1733   unfolding Polynomial.order_def
       
  1734   apply (rule Least_equality [symmetric])
       
  1735    apply (fact assms)
       
  1736   apply (rule classical)
       
  1737   apply (erule notE)
       
  1738   unfolding not_less_eq_eq
       
  1739   using assms(1)
       
  1740   apply (rule power_le_dvd)
       
  1741   apply assumption
       
  1742   done
       
  1743 
       
  1744 lemma order_mult: "p * q \<noteq> 0 \<Longrightarrow> order a (p * q) = order a p + order a q"
       
  1745 proof -
       
  1746   define i where "i = order a p"
       
  1747   define j where "j = order a q"
       
  1748   define t where "t = [:-a, 1:]"
       
  1749   have t_dvd_iff: "\<And>u. t dvd u \<longleftrightarrow> poly u a = 0"
       
  1750     by (simp add: t_def dvd_iff_poly_eq_0)
       
  1751   assume "p * q \<noteq> 0"
       
  1752   then show "order a (p * q) = i + j"
       
  1753     apply clarsimp
       
  1754     apply (drule order [where a=a and p=p, folded i_def t_def])
       
  1755     apply (drule order [where a=a and p=q, folded j_def t_def])
       
  1756     apply clarify
       
  1757     apply (erule dvdE)+
       
  1758     apply (rule order_unique_lemma [symmetric], fold t_def)
       
  1759      apply (simp_all add: power_add t_dvd_iff)
       
  1760     done
       
  1761 qed
       
  1762 
       
  1763 lemma order_smult:
       
  1764   assumes "c \<noteq> 0"
       
  1765   shows "order x (smult c p) = order x p"
       
  1766 proof (cases "p = 0")
       
  1767   case True
       
  1768   then show ?thesis
       
  1769     by simp
       
  1770 next
       
  1771   case False
       
  1772   have "smult c p = [:c:] * p" by simp
       
  1773   also from assms False have "order x \<dots> = order x [:c:] + order x p"
       
  1774     by (subst order_mult) simp_all
       
  1775   also have "order x [:c:] = 0"
       
  1776     by (rule order_0I) (use assms in auto)
       
  1777   finally show ?thesis
       
  1778     by simp
       
  1779 qed
       
  1780 
       
  1781 (* Next two lemmas contributed by Wenda Li *)
       
  1782 lemma order_1_eq_0 [simp]:"order x 1 = 0"
       
  1783   by (metis order_root poly_1 zero_neq_one)
       
  1784 
       
  1785 lemma order_power_n_n: "order a ([:-a,1:]^n)=n"
       
  1786 proof (induct n) (*might be proved more concisely using nat_less_induct*)
       
  1787   case 0
       
  1788   then show ?case
       
  1789     by (metis order_root poly_1 power_0 zero_neq_one)
       
  1790 next
       
  1791   case (Suc n)
       
  1792   have "order a ([:- a, 1:] ^ Suc n) = order a ([:- a, 1:] ^ n) + order a [:-a,1:]"
       
  1793     by (metis (no_types, hide_lams) One_nat_def add_Suc_right monoid_add_class.add.right_neutral
       
  1794       one_neq_zero order_mult pCons_eq_0_iff power_add power_eq_0_iff power_one_right)
       
  1795   moreover have "order a [:-a,1:] = 1"
       
  1796     unfolding order_def
       
  1797   proof (rule Least_equality, rule notI)
       
  1798     assume "[:- a, 1:] ^ Suc 1 dvd [:- a, 1:]"
       
  1799     then have "degree ([:- a, 1:] ^ Suc 1) \<le> degree ([:- a, 1:])"
       
  1800       by (rule dvd_imp_degree_le) auto
       
  1801     then show False
       
  1802       by auto
       
  1803   next
       
  1804     fix y
       
  1805     assume *: "\<not> [:- a, 1:] ^ Suc y dvd [:- a, 1:]"
       
  1806     show "1 \<le> y"
       
  1807     proof (rule ccontr)
       
  1808       assume "\<not> 1 \<le> y"
       
  1809       then have "y = 0" by auto
       
  1810       then have "[:- a, 1:] ^ Suc y dvd [:- a, 1:]" by auto
       
  1811       with * show False by auto
       
  1812     qed
       
  1813   qed
       
  1814   ultimately show ?case
       
  1815     using Suc by auto
       
  1816 qed
       
  1817 
       
  1818 lemma order_0_monom [simp]: "c \<noteq> 0 \<Longrightarrow> order 0 (monom c n) = n"
       
  1819   using order_power_n_n[of 0 n] by (simp add: monom_altdef order_smult)
       
  1820 
       
  1821 lemma dvd_imp_order_le: "q \<noteq> 0 \<Longrightarrow> p dvd q \<Longrightarrow> Polynomial.order a p \<le> Polynomial.order a q"
       
  1822   by (auto simp: order_mult elim: dvdE)
       
  1823 
       
  1824 text \<open>Now justify the standard squarefree decomposition, i.e. \<open>f / gcd f f'\<close>.\<close>
       
  1825 
       
  1826 lemma order_divides: "[:-a, 1:] ^ n dvd p \<longleftrightarrow> p = 0 \<or> n \<le> order a p"
       
  1827   apply (cases "p = 0")
       
  1828   apply auto
       
  1829    apply (drule order_2 [where a=a and p=p])
       
  1830    apply (metis not_less_eq_eq power_le_dvd)
       
  1831   apply (erule power_le_dvd [OF order_1])
       
  1832   done
       
  1833 
       
  1834 lemma order_decomp:
       
  1835   assumes "p \<noteq> 0"
       
  1836   shows "\<exists>q. p = [:- a, 1:] ^ order a p * q \<and> \<not> [:- a, 1:] dvd q"
       
  1837 proof -
       
  1838   from assms have *: "[:- a, 1:] ^ order a p dvd p"
       
  1839     and **: "\<not> [:- a, 1:] ^ Suc (order a p) dvd p"
       
  1840     by (auto dest: order)
       
  1841   from * obtain q where q: "p = [:- a, 1:] ^ order a p * q" ..
       
  1842   with ** have "\<not> [:- a, 1:] ^ Suc (order a p) dvd [:- a, 1:] ^ order a p * q"
       
  1843     by simp
       
  1844   then have "\<not> [:- a, 1:] ^ order a p * [:- a, 1:] dvd [:- a, 1:] ^ order a p * q"
       
  1845     by simp
       
  1846   with idom_class.dvd_mult_cancel_left [of "[:- a, 1:] ^ order a p" "[:- a, 1:]" q]
       
  1847   have "\<not> [:- a, 1:] dvd q" by auto
       
  1848   with q show ?thesis by blast
       
  1849 qed
       
  1850 
       
  1851 lemma monom_1_dvd_iff: "p \<noteq> 0 \<Longrightarrow> monom 1 n dvd p \<longleftrightarrow> n \<le> order 0 p"
       
  1852   using order_divides[of 0 n p] by (simp add: monom_altdef)
       
  1853 
       
  1854 
       
  1855 subsection \<open>Additional induction rules on polynomials\<close>
       
  1856 
       
  1857 text \<open>
       
  1858   An induction rule for induction over the roots of a polynomial with a certain property.
       
  1859   (e.g. all positive roots)
       
  1860 \<close>
       
  1861 lemma poly_root_induct [case_names 0 no_roots root]:
       
  1862   fixes p :: "'a :: idom poly"
       
  1863   assumes "Q 0"
       
  1864     and "\<And>p. (\<And>a. P a \<Longrightarrow> poly p a \<noteq> 0) \<Longrightarrow> Q p"
       
  1865     and "\<And>a p. P a \<Longrightarrow> Q p \<Longrightarrow> Q ([:a, -1:] * p)"
       
  1866   shows "Q p"
       
  1867 proof (induction "degree p" arbitrary: p rule: less_induct)
       
  1868   case (less p)
       
  1869   show ?case
       
  1870   proof (cases "p = 0")
       
  1871     case True
       
  1872     with assms(1) show ?thesis by simp
       
  1873   next
       
  1874     case False
       
  1875     show ?thesis
       
  1876     proof (cases "\<exists>a. P a \<and> poly p a = 0")
       
  1877       case False
       
  1878       then show ?thesis by (intro assms(2)) blast
       
  1879     next
       
  1880       case True
       
  1881       then obtain a where a: "P a" "poly p a = 0"
       
  1882         by blast
       
  1883       then have "-[:-a, 1:] dvd p"
       
  1884         by (subst minus_dvd_iff) (simp add: poly_eq_0_iff_dvd)
       
  1885       then obtain q where q: "p = [:a, -1:] * q" by (elim dvdE) simp
       
  1886       with False have "q \<noteq> 0" by auto
       
  1887       have "degree p = Suc (degree q)"
       
  1888         by (subst q, subst degree_mult_eq) (simp_all add: \<open>q \<noteq> 0\<close>)
       
  1889       then have "Q q" by (intro less) simp
       
  1890       with a(1) have "Q ([:a, -1:] * q)"
       
  1891         by (rule assms(3))
       
  1892       with q show ?thesis by simp
       
  1893     qed
       
  1894   qed
       
  1895 qed
       
  1896 
       
  1897 lemma dropWhile_replicate_append:
       
  1898   "dropWhile (op = a) (replicate n a @ ys) = dropWhile (op = a) ys"
       
  1899   by (induct n) simp_all
       
  1900 
       
  1901 lemma Poly_append_replicate_0: "Poly (xs @ replicate n 0) = Poly xs"
       
  1902   by (subst coeffs_eq_iff) (simp_all add: strip_while_def dropWhile_replicate_append)
       
  1903 
       
  1904 text \<open>
       
  1905   An induction rule for simultaneous induction over two polynomials,
       
  1906   prepending one coefficient in each step.
       
  1907 \<close>
       
  1908 lemma poly_induct2 [case_names 0 pCons]:
       
  1909   assumes "P 0 0" "\<And>a p b q. P p q \<Longrightarrow> P (pCons a p) (pCons b q)"
       
  1910   shows "P p q"
       
  1911 proof -
       
  1912   define n where "n = max (length (coeffs p)) (length (coeffs q))"
       
  1913   define xs where "xs = coeffs p @ (replicate (n - length (coeffs p)) 0)"
       
  1914   define ys where "ys = coeffs q @ (replicate (n - length (coeffs q)) 0)"
       
  1915   have "length xs = length ys"
       
  1916     by (simp add: xs_def ys_def n_def)
       
  1917   then have "P (Poly xs) (Poly ys)"
       
  1918     by (induct rule: list_induct2) (simp_all add: assms)
       
  1919   also have "Poly xs = p"
       
  1920     by (simp add: xs_def Poly_append_replicate_0)
       
  1921   also have "Poly ys = q"
       
  1922     by (simp add: ys_def Poly_append_replicate_0)
       
  1923   finally show ?thesis .
       
  1924 qed
       
  1925 
       
  1926 
       
  1927 subsection \<open>Composition of polynomials\<close>
       
  1928 
       
  1929 (* Several lemmas contributed by René Thiemann and Akihisa Yamada *)
       
  1930 
       
  1931 definition pcompose :: "'a::comm_semiring_0 poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly"
       
  1932   where "pcompose p q = fold_coeffs (\<lambda>a c. [:a:] + q * c) p 0"
       
  1933 
       
  1934 notation pcompose (infixl "\<circ>\<^sub>p" 71)
       
  1935 
       
  1936 lemma pcompose_0 [simp]: "pcompose 0 q = 0"
       
  1937   by (simp add: pcompose_def)
       
  1938 
       
  1939 lemma pcompose_pCons: "pcompose (pCons a p) q = [:a:] + q * pcompose p q"
       
  1940   by (cases "p = 0 \<and> a = 0") (auto simp add: pcompose_def)
       
  1941 
       
  1942 lemma pcompose_1: "pcompose 1 p = 1"
       
  1943   for p :: "'a::comm_semiring_1 poly"
       
  1944   by (auto simp: one_poly_def pcompose_pCons)
       
  1945 
       
  1946 lemma poly_pcompose: "poly (pcompose p q) x = poly p (poly q x)"
       
  1947   by (induct p) (simp_all add: pcompose_pCons)
       
  1948 
       
  1949 lemma degree_pcompose_le: "degree (pcompose p q) \<le> degree p * degree q"
       
  1950   apply (induct p)
       
  1951    apply simp
       
  1952   apply (simp add: pcompose_pCons)
       
  1953   apply clarify
       
  1954   apply (rule degree_add_le)
       
  1955    apply simp
       
  1956   apply (rule order_trans [OF degree_mult_le])
       
  1957   apply simp
       
  1958   done
       
  1959 
       
  1960 lemma pcompose_add: "pcompose (p + q) r = pcompose p r + pcompose q r"
       
  1961   for p q r :: "'a::{comm_semiring_0, ab_semigroup_add} poly"
       
  1962 proof (induction p q rule: poly_induct2)
       
  1963   case 0
       
  1964   then show ?case by simp
       
  1965 next
       
  1966   case (pCons a p b q)
       
  1967   have "pcompose (pCons a p + pCons b q) r = [:a + b:] + r * pcompose p r + r * pcompose q r"
       
  1968     by (simp_all add: pcompose_pCons pCons.IH algebra_simps)
       
  1969   also have "[:a + b:] = [:a:] + [:b:]" by simp
       
  1970   also have "\<dots> + r * pcompose p r + r * pcompose q r =
       
  1971     pcompose (pCons a p) r + pcompose (pCons b q) r"
       
  1972     by (simp only: pcompose_pCons add_ac)
       
  1973   finally show ?case .
       
  1974 qed
       
  1975 
       
  1976 lemma pcompose_uminus: "pcompose (-p) r = -pcompose p r"
       
  1977   for p r :: "'a::comm_ring poly"
       
  1978   by (induct p) (simp_all add: pcompose_pCons)
       
  1979 
       
  1980 lemma pcompose_diff: "pcompose (p - q) r = pcompose p r - pcompose q r"
       
  1981   for p q r :: "'a::comm_ring poly"
       
  1982   using pcompose_add[of p "-q"] by (simp add: pcompose_uminus)
       
  1983 
       
  1984 lemma pcompose_smult: "pcompose (smult a p) r = smult a (pcompose p r)"
       
  1985   for p r :: "'a::comm_semiring_0 poly"
       
  1986   by (induct p) (simp_all add: pcompose_pCons pcompose_add smult_add_right)
       
  1987 
       
  1988 lemma pcompose_mult: "pcompose (p * q) r = pcompose p r * pcompose q r"
       
  1989   for p q r :: "'a::comm_semiring_0 poly"
       
  1990   by (induct p arbitrary: q) (simp_all add: pcompose_add pcompose_smult pcompose_pCons algebra_simps)
       
  1991 
       
  1992 lemma pcompose_assoc: "pcompose p (pcompose q r) = pcompose (pcompose p q) r"
       
  1993   for p q r :: "'a::comm_semiring_0 poly"
       
  1994   by (induct p arbitrary: q) (simp_all add: pcompose_pCons pcompose_add pcompose_mult)
       
  1995 
       
  1996 lemma pcompose_idR[simp]: "pcompose p [: 0, 1 :] = p"
       
  1997   for p :: "'a::comm_semiring_1 poly"
       
  1998   by (induct p) (simp_all add: pcompose_pCons)
       
  1999 
       
  2000 lemma pcompose_sum: "pcompose (sum f A) p = sum (\<lambda>i. pcompose (f i) p) A"
       
  2001   by (induct A rule: infinite_finite_induct) (simp_all add: pcompose_1 pcompose_add)
       
  2002 
       
  2003 lemma pcompose_prod: "pcompose (prod f A) p = prod (\<lambda>i. pcompose (f i) p) A"
       
  2004   by (induct A rule: infinite_finite_induct) (simp_all add: pcompose_1 pcompose_mult)
       
  2005 
       
  2006 lemma pcompose_const [simp]: "pcompose [:a:] q = [:a:]"
       
  2007   by (subst pcompose_pCons) simp
       
  2008 
       
  2009 lemma pcompose_0': "pcompose p 0 = [:coeff p 0:]"
       
  2010   by (induct p) (auto simp add: pcompose_pCons)
       
  2011 
       
  2012 lemma degree_pcompose: "degree (pcompose p q) = degree p * degree q"
       
  2013   for p q :: "'a::{comm_semiring_0,semiring_no_zero_divisors} poly"
       
  2014 proof (induct p)
       
  2015   case 0
       
  2016   then show ?case by auto
       
  2017 next
       
  2018   case (pCons a p)
       
  2019   consider "degree (q * pcompose p q) = 0" | "degree (q * pcompose p q) > 0"
       
  2020     by blast
       
  2021   then show ?case
       
  2022   proof cases
       
  2023     case prems: 1
       
  2024     show ?thesis
       
  2025     proof (cases "p = 0")
       
  2026       case True
       
  2027       then show ?thesis by auto
       
  2028     next
       
  2029       case False
       
  2030       from prems have "degree q = 0 \<or> pcompose p q = 0"
       
  2031         by (auto simp add: degree_mult_eq_0)
       
  2032       moreover have False if "pcompose p q = 0" "degree q \<noteq> 0"
       
  2033       proof -
       
  2034         from pCons.hyps(2) that have "degree p = 0"
       
  2035           by auto
       
  2036         then obtain a1 where "p = [:a1:]"
       
  2037           by (metis degree_pCons_eq_if old.nat.distinct(2) pCons_cases)
       
  2038         with \<open>pcompose p q = 0\<close> \<open>p \<noteq> 0\<close> show False
       
  2039           by auto
       
  2040       qed
       
  2041       ultimately have "degree (pCons a p) * degree q = 0"
       
  2042         by auto
       
  2043       moreover have "degree (pcompose (pCons a p) q) = 0"
       
  2044       proof -
       
  2045         from prems have "0 = max (degree [:a:]) (degree (q * pcompose p q))"
       
  2046           by simp
       
  2047         also have "\<dots> \<ge> degree ([:a:] + q * pcompose p q)"
       
  2048           by (rule degree_add_le_max)
       
  2049         finally show ?thesis
       
  2050           by (auto simp add: pcompose_pCons)
       
  2051       qed
       
  2052       ultimately show ?thesis by simp
       
  2053     qed
       
  2054   next
       
  2055     case prems: 2
       
  2056     then have "p \<noteq> 0" "q \<noteq> 0" "pcompose p q \<noteq> 0"
       
  2057       by auto
       
  2058     from prems degree_add_eq_right [of "[:a:]"]
       
  2059     have "degree (pcompose (pCons a p) q) = degree (q * pcompose p q)"
       
  2060       by (auto simp: pcompose_pCons)
       
  2061     with pCons.hyps(2) degree_mult_eq[OF \<open>q\<noteq>0\<close> \<open>pcompose p q\<noteq>0\<close>] show ?thesis
       
  2062       by auto
       
  2063   qed
       
  2064 qed
       
  2065 
       
  2066 lemma pcompose_eq_0:
       
  2067   fixes p q :: "'a::{comm_semiring_0,semiring_no_zero_divisors} poly"
       
  2068   assumes "pcompose p q = 0" "degree q > 0"
       
  2069   shows "p = 0"
       
  2070 proof -
       
  2071   from assms degree_pcompose [of p q] have "degree p = 0"
       
  2072     by auto
       
  2073   then obtain a where "p = [:a:]"
       
  2074     by (metis degree_pCons_eq_if gr0_conv_Suc neq0_conv pCons_cases)
       
  2075   with assms(1) have "a = 0"
       
  2076     by auto
       
  2077   with \<open>p = [:a:]\<close> show ?thesis
       
  2078     by simp
       
  2079 qed
       
  2080 
       
  2081 lemma lead_coeff_comp:
       
  2082   fixes p q :: "'a::{comm_semiring_1,semiring_no_zero_divisors} poly"
       
  2083   assumes "degree q > 0"
       
  2084   shows "lead_coeff (pcompose p q) = lead_coeff p * lead_coeff q ^ (degree p)"
       
  2085 proof (induct p)
       
  2086   case 0
       
  2087   then show ?case by auto
       
  2088 next
       
  2089   case (pCons a p)
       
  2090   consider "degree (q * pcompose p q) = 0" | "degree (q * pcompose p q) > 0"
       
  2091     by blast
       
  2092   then show ?case
       
  2093   proof cases
       
  2094     case prems: 1
       
  2095     then have "pcompose p q = 0"
       
  2096       by (metis assms degree_0 degree_mult_eq_0 neq0_conv)
       
  2097     with pcompose_eq_0[OF _ \<open>degree q > 0\<close>] have "p = 0"
       
  2098       by simp
       
  2099     then show ?thesis
       
  2100       by auto
       
  2101   next
       
  2102     case prems: 2
       
  2103     then have "degree [:a:] < degree (q * pcompose p q)"
       
  2104       by simp
       
  2105     then have "lead_coeff ([:a:] + q * p \<circ>\<^sub>p q) = lead_coeff (q * p \<circ>\<^sub>p q)"
       
  2106       by (rule lead_coeff_add_le)
       
  2107     then have "lead_coeff (pcompose (pCons a p) q) = lead_coeff (q * pcompose p q)"
       
  2108       by (simp add: pcompose_pCons)
       
  2109     also have "\<dots> = lead_coeff q * (lead_coeff p * lead_coeff q ^ degree p)"
       
  2110       using pCons.hyps(2) lead_coeff_mult[of q "pcompose p q"] by simp
       
  2111     also have "\<dots> = lead_coeff p * lead_coeff q ^ (degree p + 1)"
       
  2112       by (auto simp: mult_ac)
       
  2113     finally show ?thesis by auto
       
  2114   qed
       
  2115 qed
       
  2116 
       
  2117 
       
  2118 subsection \<open>Shifting polynomials\<close>
       
  2119 
       
  2120 definition poly_shift :: "nat \<Rightarrow> 'a::zero poly \<Rightarrow> 'a poly"
       
  2121   where "poly_shift n p = Abs_poly (\<lambda>i. coeff p (i + n))"
       
  2122 
       
  2123 lemma nth_default_drop: "nth_default x (drop n xs) m = nth_default x xs (m + n)"
       
  2124   by (auto simp add: nth_default_def add_ac)
       
  2125 
       
  2126 lemma nth_default_take: "nth_default x (take n xs) m = (if m < n then nth_default x xs m else x)"
       
  2127   by (auto simp add: nth_default_def add_ac)
       
  2128 
       
  2129 lemma coeff_poly_shift: "coeff (poly_shift n p) i = coeff p (i + n)"
       
  2130 proof -
       
  2131   from MOST_coeff_eq_0[of p] obtain m where "\<forall>k>m. coeff p k = 0"
       
  2132     by (auto simp: MOST_nat)
       
  2133   then have "\<forall>k>m. coeff p (k + n) = 0"
       
  2134     by auto
       
  2135   then have "\<forall>\<^sub>\<infinity>k. coeff p (k + n) = 0"
       
  2136     by (auto simp: MOST_nat)
       
  2137   then show ?thesis
       
  2138     by (simp add: poly_shift_def poly.Abs_poly_inverse)
       
  2139 qed
       
  2140 
       
  2141 lemma poly_shift_id [simp]: "poly_shift 0 = (\<lambda>x. x)"
       
  2142   by (simp add: poly_eq_iff fun_eq_iff coeff_poly_shift)
       
  2143 
       
  2144 lemma poly_shift_0 [simp]: "poly_shift n 0 = 0"
       
  2145   by (simp add: poly_eq_iff coeff_poly_shift)
       
  2146 
       
  2147 lemma poly_shift_1: "poly_shift n 1 = (if n = 0 then 1 else 0)"
       
  2148   by (simp add: poly_eq_iff coeff_poly_shift)
       
  2149 
       
  2150 lemma poly_shift_monom: "poly_shift n (monom c m) = (if m \<ge> n then monom c (m - n) else 0)"
       
  2151   by (auto simp add: poly_eq_iff coeff_poly_shift)
       
  2152 
       
  2153 lemma coeffs_shift_poly [code abstract]:
       
  2154   "coeffs (poly_shift n p) = drop n (coeffs p)"
       
  2155 proof (cases "p = 0")
       
  2156   case True
       
  2157   then show ?thesis by simp
       
  2158 next
       
  2159   case False
       
  2160   then show ?thesis
       
  2161     by (intro coeffs_eqI)
       
  2162       (simp_all add: coeff_poly_shift nth_default_drop nth_default_coeffs_eq)
       
  2163 qed
       
  2164 
       
  2165 
       
  2166 subsection \<open>Truncating polynomials\<close>
       
  2167 
       
  2168 definition poly_cutoff
       
  2169   where "poly_cutoff n p = Abs_poly (\<lambda>k. if k < n then coeff p k else 0)"
       
  2170 
       
  2171 lemma coeff_poly_cutoff: "coeff (poly_cutoff n p) k = (if k < n then coeff p k else 0)"
       
  2172   unfolding poly_cutoff_def
       
  2173   by (subst poly.Abs_poly_inverse) (auto simp: MOST_nat intro: exI[of _ n])
       
  2174 
       
  2175 lemma poly_cutoff_0 [simp]: "poly_cutoff n 0 = 0"
       
  2176   by (simp add: poly_eq_iff coeff_poly_cutoff)
       
  2177 
       
  2178 lemma poly_cutoff_1 [simp]: "poly_cutoff n 1 = (if n = 0 then 0 else 1)"
       
  2179   by (simp add: poly_eq_iff coeff_poly_cutoff)
       
  2180 
       
  2181 lemma coeffs_poly_cutoff [code abstract]:
       
  2182   "coeffs (poly_cutoff n p) = strip_while (op = 0) (take n (coeffs p))"
       
  2183 proof (cases "strip_while (op = 0) (take n (coeffs p)) = []")
       
  2184   case True
       
  2185   then have "coeff (poly_cutoff n p) k = 0" for k
       
  2186     unfolding coeff_poly_cutoff
       
  2187     by (auto simp: nth_default_coeffs_eq [symmetric] nth_default_def set_conv_nth)
       
  2188   then have "poly_cutoff n p = 0"
       
  2189     by (simp add: poly_eq_iff)
       
  2190   then show ?thesis
       
  2191     by (subst True) simp_all
       
  2192 next
       
  2193   case False
       
  2194   have "no_trailing (op = 0) (strip_while (op = 0) (take n (coeffs p)))"
       
  2195     by simp
       
  2196   with False have "last (strip_while (op = 0) (take n (coeffs p))) \<noteq> 0"
       
  2197     unfolding no_trailing_unfold by auto
       
  2198   then show ?thesis
       
  2199     by (intro coeffs_eqI)
       
  2200       (simp_all add: coeff_poly_cutoff nth_default_take nth_default_coeffs_eq)
       
  2201 qed
       
  2202 
       
  2203 
       
  2204 subsection \<open>Reflecting polynomials\<close>
       
  2205 
       
  2206 definition reflect_poly :: "'a::zero poly \<Rightarrow> 'a poly"
       
  2207   where "reflect_poly p = Poly (rev (coeffs p))"
       
  2208 
       
  2209 lemma coeffs_reflect_poly [code abstract]:
       
  2210   "coeffs (reflect_poly p) = rev (dropWhile (op = 0) (coeffs p))"
       
  2211   by (simp add: reflect_poly_def)
       
  2212 
       
  2213 lemma reflect_poly_0 [simp]: "reflect_poly 0 = 0"
       
  2214   by (simp add: reflect_poly_def)
       
  2215 
       
  2216 lemma reflect_poly_1 [simp]: "reflect_poly 1 = 1"
       
  2217   by (simp add: reflect_poly_def one_poly_def)
       
  2218 
       
  2219 lemma coeff_reflect_poly:
       
  2220   "coeff (reflect_poly p) n = (if n > degree p then 0 else coeff p (degree p - n))"
       
  2221   by (cases "p = 0")
       
  2222     (auto simp add: reflect_poly_def nth_default_def
       
  2223       rev_nth degree_eq_length_coeffs coeffs_nth not_less
       
  2224       dest: le_imp_less_Suc)
       
  2225 
       
  2226 lemma coeff_0_reflect_poly_0_iff [simp]: "coeff (reflect_poly p) 0 = 0 \<longleftrightarrow> p = 0"
       
  2227   by (simp add: coeff_reflect_poly)
       
  2228 
       
  2229 lemma reflect_poly_at_0_eq_0_iff [simp]: "poly (reflect_poly p) 0 = 0 \<longleftrightarrow> p = 0"
       
  2230   by (simp add: coeff_reflect_poly poly_0_coeff_0)
       
  2231 
       
  2232 lemma reflect_poly_pCons':
       
  2233   "p \<noteq> 0 \<Longrightarrow> reflect_poly (pCons c p) = reflect_poly p + monom c (Suc (degree p))"
       
  2234   by (intro poly_eqI)
       
  2235     (auto simp: coeff_reflect_poly coeff_pCons not_less Suc_diff_le split: nat.split)
       
  2236 
       
  2237 lemma reflect_poly_const [simp]: "reflect_poly [:a:] = [:a:]"
       
  2238   by (cases "a = 0") (simp_all add: reflect_poly_def)
       
  2239 
       
  2240 lemma poly_reflect_poly_nz:
       
  2241   "x \<noteq> 0 \<Longrightarrow> poly (reflect_poly p) x = x ^ degree p * poly p (inverse x)"
       
  2242   for x :: "'a::field"
       
  2243   by (induct rule: pCons_induct) (simp_all add: field_simps reflect_poly_pCons' poly_monom)
       
  2244 
       
  2245 lemma coeff_0_reflect_poly [simp]: "coeff (reflect_poly p) 0 = lead_coeff p"
       
  2246   by (simp add: coeff_reflect_poly)
       
  2247 
       
  2248 lemma poly_reflect_poly_0 [simp]: "poly (reflect_poly p) 0 = lead_coeff p"
       
  2249   by (simp add: poly_0_coeff_0)
       
  2250 
       
  2251 lemma reflect_poly_reflect_poly [simp]: "coeff p 0 \<noteq> 0 \<Longrightarrow> reflect_poly (reflect_poly p) = p"
       
  2252   by (cases p rule: pCons_cases) (simp add: reflect_poly_def )
       
  2253 
       
  2254 lemma degree_reflect_poly_le: "degree (reflect_poly p) \<le> degree p"
       
  2255   by (simp add: degree_eq_length_coeffs coeffs_reflect_poly length_dropWhile_le diff_le_mono)
       
  2256 
       
  2257 lemma reflect_poly_pCons: "a \<noteq> 0 \<Longrightarrow> reflect_poly (pCons a p) = Poly (rev (a # coeffs p))"
       
  2258   by (subst coeffs_eq_iff) (simp add: coeffs_reflect_poly)
       
  2259 
       
  2260 lemma degree_reflect_poly_eq [simp]: "coeff p 0 \<noteq> 0 \<Longrightarrow> degree (reflect_poly p) = degree p"
       
  2261   by (cases p rule: pCons_cases) (simp add: reflect_poly_pCons degree_eq_length_coeffs)
       
  2262 
       
  2263 (* TODO: does this work with zero divisors as well? Probably not. *)
       
  2264 lemma reflect_poly_mult: "reflect_poly (p * q) = reflect_poly p * reflect_poly q"
       
  2265   for p q :: "'a::{comm_semiring_0,semiring_no_zero_divisors} poly"
       
  2266 proof (cases "p = 0 \<or> q = 0")
       
  2267   case False
       
  2268   then have [simp]: "p \<noteq> 0" "q \<noteq> 0" by auto
       
  2269   show ?thesis
       
  2270   proof (rule poly_eqI)
       
  2271     show "coeff (reflect_poly (p * q)) i = coeff (reflect_poly p * reflect_poly q) i" for i
       
  2272     proof (cases "i \<le> degree (p * q)")
       
  2273       case True
       
  2274       define A where "A = {..i} \<inter> {i - degree q..degree p}"
       
  2275       define B where "B = {..degree p} \<inter> {degree p - i..degree (p*q) - i}"
       
  2276       let ?f = "\<lambda>j. degree p - j"
       
  2277 
       
  2278       from True have "coeff (reflect_poly (p * q)) i = coeff (p * q) (degree (p * q) - i)"
       
  2279         by (simp add: coeff_reflect_poly)
       
  2280       also have "\<dots> = (\<Sum>j\<le>degree (p * q) - i. coeff p j * coeff q (degree (p * q) - i - j))"
       
  2281         by (simp add: coeff_mult)
       
  2282       also have "\<dots> = (\<Sum>j\<in>B. coeff p j * coeff q (degree (p * q) - i - j))"
       
  2283         by (intro sum.mono_neutral_right) (auto simp: B_def degree_mult_eq not_le coeff_eq_0)
       
  2284       also from True have "\<dots> = (\<Sum>j\<in>A. coeff p (degree p - j) * coeff q (degree q - (i - j)))"
       
  2285         by (intro sum.reindex_bij_witness[of _ ?f ?f])
       
  2286           (auto simp: A_def B_def degree_mult_eq add_ac)
       
  2287       also have "\<dots> =
       
  2288         (\<Sum>j\<le>i.
       
  2289           if j \<in> {i - degree q..degree p}
       
  2290           then coeff p (degree p - j) * coeff q (degree q - (i - j))
       
  2291           else 0)"
       
  2292         by (subst sum.inter_restrict [symmetric]) (simp_all add: A_def)
       
  2293       also have "\<dots> = coeff (reflect_poly p * reflect_poly q) i"
       
  2294         by (fastforce simp: coeff_mult coeff_reflect_poly intro!: sum.cong)
       
  2295       finally show ?thesis .
       
  2296     qed (auto simp: coeff_mult coeff_reflect_poly coeff_eq_0 degree_mult_eq intro!: sum.neutral)
       
  2297   qed
       
  2298 qed auto
       
  2299 
       
  2300 lemma reflect_poly_smult: "reflect_poly (smult c p) = smult c (reflect_poly p)"
       
  2301   for p :: "'a::{comm_semiring_0,semiring_no_zero_divisors} poly"
       
  2302   using reflect_poly_mult[of "[:c:]" p] by simp
       
  2303 
       
  2304 lemma reflect_poly_power: "reflect_poly (p ^ n) = reflect_poly p ^ n"
       
  2305   for p :: "'a::{comm_semiring_1,semiring_no_zero_divisors} poly"
       
  2306   by (induct n) (simp_all add: reflect_poly_mult)
       
  2307 
       
  2308 lemma reflect_poly_prod: "reflect_poly (prod f A) = prod (\<lambda>x. reflect_poly (f x)) A"
       
  2309   for f :: "_ \<Rightarrow> _::{comm_semiring_0,semiring_no_zero_divisors} poly"
       
  2310   by (induct A rule: infinite_finite_induct) (simp_all add: reflect_poly_mult)
       
  2311 
       
  2312 lemma reflect_poly_prod_list: "reflect_poly (prod_list xs) = prod_list (map reflect_poly xs)"
       
  2313   for xs :: "_::{comm_semiring_0,semiring_no_zero_divisors} poly list"
       
  2314   by (induct xs) (simp_all add: reflect_poly_mult)
       
  2315 
       
  2316 lemma reflect_poly_Poly_nz:
       
  2317   "no_trailing (HOL.eq 0) xs \<Longrightarrow> reflect_poly (Poly xs) = Poly (rev xs)"
       
  2318   by (simp add: reflect_poly_def)
       
  2319 
       
  2320 lemmas reflect_poly_simps =
       
  2321   reflect_poly_0 reflect_poly_1 reflect_poly_const reflect_poly_smult reflect_poly_mult
       
  2322   reflect_poly_power reflect_poly_prod reflect_poly_prod_list
       
  2323 
       
  2324 
       
  2325 subsection \<open>Derivatives\<close>
       
  2326 
       
  2327 function pderiv :: "('a :: {comm_semiring_1,semiring_no_zero_divisors}) poly \<Rightarrow> 'a poly"
       
  2328   where "pderiv (pCons a p) = (if p = 0 then 0 else p + pCons 0 (pderiv p))"
       
  2329   by (auto intro: pCons_cases)
       
  2330 
       
  2331 termination pderiv
       
  2332   by (relation "measure degree") simp_all
       
  2333 
       
  2334 declare pderiv.simps[simp del]
       
  2335 
       
  2336 lemma pderiv_0 [simp]: "pderiv 0 = 0"
       
  2337   using pderiv.simps [of 0 0] by simp
       
  2338 
       
  2339 lemma pderiv_pCons: "pderiv (pCons a p) = p + pCons 0 (pderiv p)"
       
  2340   by (simp add: pderiv.simps)
       
  2341 
       
  2342 lemma pderiv_1 [simp]: "pderiv 1 = 0"
       
  2343   by (simp add: one_poly_def pderiv_pCons)
       
  2344 
       
  2345 lemma pderiv_of_nat [simp]: "pderiv (of_nat n) = 0"
       
  2346   and pderiv_numeral [simp]: "pderiv (numeral m) = 0"
       
  2347   by (simp_all add: of_nat_poly numeral_poly pderiv_pCons)
       
  2348 
       
  2349 lemma coeff_pderiv: "coeff (pderiv p) n = of_nat (Suc n) * coeff p (Suc n)"
       
  2350   by (induct p arbitrary: n)
       
  2351     (auto simp add: pderiv_pCons coeff_pCons algebra_simps split: nat.split)
       
  2352 
       
  2353 fun pderiv_coeffs_code :: "'a::{comm_semiring_1,semiring_no_zero_divisors} \<Rightarrow> 'a list \<Rightarrow> 'a list"
       
  2354   where
       
  2355     "pderiv_coeffs_code f (x # xs) = cCons (f * x) (pderiv_coeffs_code (f+1) xs)"
       
  2356   | "pderiv_coeffs_code f [] = []"
       
  2357 
       
  2358 definition pderiv_coeffs :: "'a::{comm_semiring_1,semiring_no_zero_divisors} list \<Rightarrow> 'a list"
       
  2359   where "pderiv_coeffs xs = pderiv_coeffs_code 1 (tl xs)"
       
  2360 
       
  2361 (* Efficient code for pderiv contributed by René Thiemann and Akihisa Yamada *)
       
  2362 lemma pderiv_coeffs_code:
       
  2363   "nth_default 0 (pderiv_coeffs_code f xs) n = (f + of_nat n) * nth_default 0 xs n"
       
  2364 proof (induct xs arbitrary: f n)
       
  2365   case Nil
       
  2366   then show ?case by simp
       
  2367 next
       
  2368   case (Cons x xs)
       
  2369   show ?case
       
  2370   proof (cases n)
       
  2371     case 0
       
  2372     then show ?thesis
       
  2373       by (cases "pderiv_coeffs_code (f + 1) xs = [] \<and> f * x = 0") (auto simp: cCons_def)
       
  2374   next
       
  2375     case n: (Suc m)
       
  2376     show ?thesis
       
  2377     proof (cases "pderiv_coeffs_code (f + 1) xs = [] \<and> f * x = 0")
       
  2378       case False
       
  2379       then have "nth_default 0 (pderiv_coeffs_code f (x # xs)) n =
       
  2380           nth_default 0 (pderiv_coeffs_code (f + 1) xs) m"
       
  2381         by (auto simp: cCons_def n)
       
  2382       also have "\<dots> = (f + of_nat n) * nth_default 0 xs m"
       
  2383         by (simp add: Cons n add_ac)
       
  2384       finally show ?thesis
       
  2385         by (simp add: n)
       
  2386     next
       
  2387       case True
       
  2388       have empty: "pderiv_coeffs_code g xs = [] \<Longrightarrow> g + of_nat m = 0 \<or> nth_default 0 xs m = 0" for g
       
  2389       proof (induct xs arbitrary: g m)
       
  2390         case Nil
       
  2391         then show ?case by simp
       
  2392       next
       
  2393         case (Cons x xs)
       
  2394         from Cons(2) have empty: "pderiv_coeffs_code (g + 1) xs = []" and g: "g = 0 \<or> x = 0"
       
  2395           by (auto simp: cCons_def split: if_splits)
       
  2396         note IH = Cons(1)[OF empty]
       
  2397         from IH[of m] IH[of "m - 1"] g show ?case
       
  2398           by (cases m) (auto simp: field_simps)
       
  2399       qed
       
  2400       from True have "nth_default 0 (pderiv_coeffs_code f (x # xs)) n = 0"
       
  2401         by (auto simp: cCons_def n)
       
  2402       moreover from True have "(f + of_nat n) * nth_default 0 (x # xs) n = 0"
       
  2403         by (simp add: n) (use empty[of "f+1"] in \<open>auto simp: field_simps\<close>)
       
  2404       ultimately show ?thesis by simp
       
  2405     qed
       
  2406   qed
       
  2407 qed
       
  2408 
       
  2409 lemma map_upt_Suc: "map f [0 ..< Suc n] = f 0 # map (\<lambda>i. f (Suc i)) [0 ..< n]"
       
  2410   by (induct n arbitrary: f) auto
       
  2411 
       
  2412 lemma coeffs_pderiv_code [code abstract]: "coeffs (pderiv p) = pderiv_coeffs (coeffs p)"
       
  2413   unfolding pderiv_coeffs_def
       
  2414 proof (rule coeffs_eqI, unfold pderiv_coeffs_code coeff_pderiv, goal_cases)
       
  2415   case (1 n)
       
  2416   have id: "coeff p (Suc n) = nth_default 0 (map (\<lambda>i. coeff p (Suc i)) [0..<degree p]) n"
       
  2417     by (cases "n < degree p") (auto simp: nth_default_def coeff_eq_0)
       
  2418   show ?case
       
  2419     unfolding coeffs_def map_upt_Suc by (auto simp: id)
       
  2420 next
       
  2421   case 2
       
  2422   obtain n :: 'a and xs where defs: "tl (coeffs p) = xs" "1 = n"
       
  2423     by simp
       
  2424   from 2 show ?case
       
  2425     unfolding defs by (induct xs arbitrary: n) (auto simp: cCons_def)
       
  2426 qed
       
  2427 
       
  2428 lemma pderiv_eq_0_iff: "pderiv p = 0 \<longleftrightarrow> degree p = 0"
       
  2429   for p :: "'a::{comm_semiring_1,semiring_no_zero_divisors,semiring_char_0} poly"
       
  2430   apply (rule iffI)
       
  2431    apply (cases p)
       
  2432    apply simp
       
  2433    apply (simp add: poly_eq_iff coeff_pderiv del: of_nat_Suc)
       
  2434   apply (simp add: poly_eq_iff coeff_pderiv coeff_eq_0)
       
  2435   done
       
  2436 
       
  2437 lemma degree_pderiv: "degree (pderiv p) = degree p - 1"
       
  2438   for p :: "'a::{comm_semiring_1,semiring_no_zero_divisors,semiring_char_0} poly"
       
  2439   apply (rule order_antisym [OF degree_le])
       
  2440    apply (simp add: coeff_pderiv coeff_eq_0)
       
  2441   apply (cases "degree p")
       
  2442    apply simp
       
  2443   apply (rule le_degree)
       
  2444   apply (simp add: coeff_pderiv del: of_nat_Suc)
       
  2445   apply (metis degree_0 leading_coeff_0_iff nat.distinct(1))
       
  2446   done
       
  2447 
       
  2448 lemma not_dvd_pderiv:
       
  2449   fixes p :: "'a::{comm_semiring_1,semiring_no_zero_divisors,semiring_char_0} poly"
       
  2450   assumes "degree p \<noteq> 0"
       
  2451   shows "\<not> p dvd pderiv p"
       
  2452 proof
       
  2453   assume dvd: "p dvd pderiv p"
       
  2454   then obtain q where p: "pderiv p = p * q"
       
  2455     unfolding dvd_def by auto
       
  2456   from dvd have le: "degree p \<le> degree (pderiv p)"
       
  2457     by (simp add: assms dvd_imp_degree_le pderiv_eq_0_iff)
       
  2458   from assms and this [unfolded degree_pderiv]
       
  2459     show False by auto
       
  2460 qed
       
  2461 
       
  2462 lemma dvd_pderiv_iff [simp]: "p dvd pderiv p \<longleftrightarrow> degree p = 0"
       
  2463   for p :: "'a::{comm_semiring_1,semiring_no_zero_divisors,semiring_char_0} poly"
       
  2464   using not_dvd_pderiv[of p] by (auto simp: pderiv_eq_0_iff [symmetric])
       
  2465 
       
  2466 lemma pderiv_singleton [simp]: "pderiv [:a:] = 0"
       
  2467   by (simp add: pderiv_pCons)
       
  2468 
       
  2469 lemma pderiv_add: "pderiv (p + q) = pderiv p + pderiv q"
       
  2470   by (rule poly_eqI) (simp add: coeff_pderiv algebra_simps)
       
  2471 
       
  2472 lemma pderiv_minus: "pderiv (- p :: 'a :: idom poly) = - pderiv p"
       
  2473   by (rule poly_eqI) (simp add: coeff_pderiv algebra_simps)
       
  2474 
       
  2475 lemma pderiv_diff: "pderiv ((p :: _ :: idom poly) - q) = pderiv p - pderiv q"
       
  2476   by (rule poly_eqI) (simp add: coeff_pderiv algebra_simps)
       
  2477 
       
  2478 lemma pderiv_smult: "pderiv (smult a p) = smult a (pderiv p)"
       
  2479   by (rule poly_eqI) (simp add: coeff_pderiv algebra_simps)
       
  2480 
       
  2481 lemma pderiv_mult: "pderiv (p * q) = p * pderiv q + q * pderiv p"
       
  2482   by (induct p) (auto simp: pderiv_add pderiv_smult pderiv_pCons algebra_simps)
       
  2483 
       
  2484 lemma pderiv_power_Suc: "pderiv (p ^ Suc n) = smult (of_nat (Suc n)) (p ^ n) * pderiv p"
       
  2485   apply (induct n)
       
  2486    apply simp
       
  2487   apply (subst power_Suc)
       
  2488   apply (subst pderiv_mult)
       
  2489   apply (erule ssubst)
       
  2490   apply (simp only: of_nat_Suc smult_add_left smult_1_left)
       
  2491   apply (simp add: algebra_simps)
       
  2492   done
       
  2493 
       
  2494 lemma pderiv_prod: "pderiv (prod f (as)) = (\<Sum>a\<in>as. prod f (as - {a}) * pderiv (f a))"
       
  2495 proof (induct as rule: infinite_finite_induct)
       
  2496   case (insert a as)
       
  2497   then have id: "prod f (insert a as) = f a * prod f as"
       
  2498     "\<And>g. sum g (insert a as) = g a + sum g as"
       
  2499     "insert a as - {a} = as"
       
  2500     by auto
       
  2501   have "prod f (insert a as - {b}) = f a * prod f (as - {b})" if "b \<in> as" for b
       
  2502   proof -
       
  2503     from \<open>a \<notin> as\<close> that have *: "insert a as - {b} = insert a (as - {b})"
       
  2504       by auto
       
  2505     show ?thesis
       
  2506       unfolding * by (subst prod.insert) (use insert in auto)
       
  2507   qed
       
  2508   then show ?case
       
  2509     unfolding id pderiv_mult insert(3) sum_distrib_left
       
  2510     by (auto simp add: ac_simps intro!: sum.cong)
       
  2511 qed auto
       
  2512 
       
  2513 lemma DERIV_pow2: "DERIV (\<lambda>x. x ^ Suc n) x :> real (Suc n) * (x ^ n)"
       
  2514   by (rule DERIV_cong, rule DERIV_pow) simp
       
  2515 declare DERIV_pow2 [simp] DERIV_pow [simp]
       
  2516 
       
  2517 lemma DERIV_add_const: "DERIV f x :> D \<Longrightarrow> DERIV (\<lambda>x. a + f x :: 'a::real_normed_field) x :> D"
       
  2518   by (rule DERIV_cong, rule DERIV_add) auto
       
  2519 
       
  2520 lemma poly_DERIV [simp]: "DERIV (\<lambda>x. poly p x) x :> poly (pderiv p) x"
       
  2521   by (induct p) (auto intro!: derivative_eq_intros simp add: pderiv_pCons)
       
  2522 
       
  2523 lemma continuous_on_poly [continuous_intros]:
       
  2524   fixes p :: "'a :: {real_normed_field} poly"
       
  2525   assumes "continuous_on A f"
       
  2526   shows "continuous_on A (\<lambda>x. poly p (f x))"
       
  2527 proof -
       
  2528   have "continuous_on A (\<lambda>x. (\<Sum>i\<le>degree p. (f x) ^ i * coeff p i))"
       
  2529     by (intro continuous_intros assms)
       
  2530   also have "\<dots> = (\<lambda>x. poly p (f x))"
       
  2531     by (rule ext) (simp add: poly_altdef mult_ac)
       
  2532   finally show ?thesis .
       
  2533 qed
       
  2534 
       
  2535 text \<open>Consequences of the derivative theorem above.\<close>
       
  2536 
       
  2537 lemma poly_differentiable[simp]: "(\<lambda>x. poly p x) differentiable (at x)"
       
  2538   for x :: real
       
  2539   by (simp add: real_differentiable_def) (blast intro: poly_DERIV)
       
  2540 
       
  2541 lemma poly_isCont[simp]: "isCont (\<lambda>x. poly p x) x"
       
  2542   for x :: real
       
  2543   by (rule poly_DERIV [THEN DERIV_isCont])
       
  2544 
       
  2545 lemma poly_IVT_pos: "a < b \<Longrightarrow> poly p a < 0 \<Longrightarrow> 0 < poly p b \<Longrightarrow> \<exists>x. a < x \<and> x < b \<and> poly p x = 0"
       
  2546   for a b :: real
       
  2547   using IVT_objl [of "poly p" a 0 b] by (auto simp add: order_le_less)
       
  2548 
       
  2549 lemma poly_IVT_neg: "a < b \<Longrightarrow> 0 < poly p a \<Longrightarrow> poly p b < 0 \<Longrightarrow> \<exists>x. a < x \<and> x < b \<and> poly p x = 0"
       
  2550   for a b :: real
       
  2551   using poly_IVT_pos [where p = "- p"] by simp
       
  2552 
       
  2553 lemma poly_IVT: "a < b \<Longrightarrow> poly p a * poly p b < 0 \<Longrightarrow> \<exists>x>a. x < b \<and> poly p x = 0"
       
  2554   for p :: "real poly"
       
  2555   by (metis less_not_sym mult_less_0_iff poly_IVT_neg poly_IVT_pos)
       
  2556 
       
  2557 lemma poly_MVT: "a < b \<Longrightarrow> \<exists>x. a < x \<and> x < b \<and> poly p b - poly p a = (b - a) * poly (pderiv p) x"
       
  2558   for a b :: real
       
  2559   using MVT [of a b "poly p"]
       
  2560   apply auto
       
  2561   apply (rule_tac x = z in exI)
       
  2562   apply (auto simp add: mult_left_cancel poly_DERIV [THEN DERIV_unique])
       
  2563   done
       
  2564 
       
  2565 lemma poly_MVT':
       
  2566   fixes a b :: real
       
  2567   assumes "{min a b..max a b} \<subseteq> A"
       
  2568   shows "\<exists>x\<in>A. poly p b - poly p a = (b - a) * poly (pderiv p) x"
       
  2569 proof (cases a b rule: linorder_cases)
       
  2570   case less
       
  2571   from poly_MVT[OF less, of p] guess x by (elim exE conjE)
       
  2572   then show ?thesis by (intro bexI[of _ x]) (auto intro!: subsetD[OF assms])
       
  2573 next
       
  2574   case greater
       
  2575   from poly_MVT[OF greater, of p] guess x by (elim exE conjE)
       
  2576   then show ?thesis by (intro bexI[of _ x]) (auto simp: algebra_simps intro!: subsetD[OF assms])
       
  2577 qed (use assms in auto)
       
  2578 
       
  2579 lemma poly_pinfty_gt_lc:
       
  2580   fixes p :: "real poly"
       
  2581   assumes "lead_coeff p > 0"
       
  2582   shows "\<exists>n. \<forall> x \<ge> n. poly p x \<ge> lead_coeff p"
       
  2583   using assms
       
  2584 proof (induct p)
       
  2585   case 0
       
  2586   then show ?case by auto
       
  2587 next
       
  2588   case (pCons a p)
       
  2589   from this(1) consider "a \<noteq> 0" "p = 0" | "p \<noteq> 0" by auto
       
  2590   then show ?case
       
  2591   proof cases
       
  2592     case 1
       
  2593     then show ?thesis by auto
       
  2594   next
       
  2595     case 2
       
  2596     with pCons obtain n1 where gte_lcoeff: "\<forall>x\<ge>n1. lead_coeff p \<le> poly p x"
       
  2597       by auto
       
  2598     from pCons(3) \<open>p \<noteq> 0\<close> have gt_0: "lead_coeff p > 0" by auto
       
  2599     define n where "n = max n1 (1 + \<bar>a\<bar> / lead_coeff p)"
       
  2600     have "lead_coeff (pCons a p) \<le> poly (pCons a p) x" if "n \<le> x" for x
       
  2601     proof -
       
  2602       from gte_lcoeff that have "lead_coeff p \<le> poly p x"
       
  2603         by (auto simp: n_def)
       
  2604       with gt_0 have "\<bar>a\<bar> / lead_coeff p \<ge> \<bar>a\<bar> / poly p x" and "poly p x > 0"
       
  2605         by (auto intro: frac_le)
       
  2606       with \<open>n \<le> x\<close>[unfolded n_def] have "x \<ge> 1 + \<bar>a\<bar> / poly p x"
       
  2607         by auto
       
  2608       with \<open>lead_coeff p \<le> poly p x\<close> \<open>poly p x > 0\<close> \<open>p \<noteq> 0\<close>
       
  2609       show "lead_coeff (pCons a p) \<le> poly (pCons a p) x"
       
  2610         by (auto simp: field_simps)
       
  2611     qed
       
  2612     then show ?thesis by blast
       
  2613   qed
       
  2614 qed
       
  2615 
       
  2616 lemma lemma_order_pderiv1:
       
  2617   "pderiv ([:- a, 1:] ^ Suc n * q) = [:- a, 1:] ^ Suc n * pderiv q +
       
  2618     smult (of_nat (Suc n)) (q * [:- a, 1:] ^ n)"
       
  2619   by (simp only: pderiv_mult pderiv_power_Suc) (simp del: power_Suc of_nat_Suc add: pderiv_pCons)
       
  2620 
       
  2621 lemma lemma_order_pderiv:
       
  2622   fixes p :: "'a :: field_char_0 poly"
       
  2623   assumes n: "0 < n"
       
  2624     and pd: "pderiv p \<noteq> 0"
       
  2625     and pe: "p = [:- a, 1:] ^ n * q"
       
  2626     and nd: "\<not> [:- a, 1:] dvd q"
       
  2627   shows "n = Suc (order a (pderiv p))"
       
  2628 proof -
       
  2629   from assms have "pderiv ([:- a, 1:] ^ n * q) \<noteq> 0"
       
  2630     by auto
       
  2631   from assms obtain n' where "n = Suc n'" "0 < Suc n'" "pderiv ([:- a, 1:] ^ Suc n' * q) \<noteq> 0"
       
  2632     by (cases n) auto
       
  2633   have *: "k dvd k * pderiv q + smult (of_nat (Suc n')) l \<Longrightarrow> k dvd l" for k l
       
  2634     by (auto simp del: of_nat_Suc simp: dvd_add_right_iff dvd_smult_iff)
       
  2635   have "n' = order a (pderiv ([:- a, 1:] ^ Suc n' * q))"
       
  2636   proof (rule order_unique_lemma)
       
  2637     show "[:- a, 1:] ^ n' dvd pderiv ([:- a, 1:] ^ Suc n' * q)"
       
  2638       apply (subst lemma_order_pderiv1)
       
  2639       apply (rule dvd_add)
       
  2640        apply (metis dvdI dvd_mult2 power_Suc2)
       
  2641       apply (metis dvd_smult dvd_triv_right)
       
  2642       done
       
  2643     show "\<not> [:- a, 1:] ^ Suc n' dvd pderiv ([:- a, 1:] ^ Suc n' * q)"
       
  2644       apply (subst lemma_order_pderiv1)
       
  2645       apply (metis * nd dvd_mult_cancel_right power_not_zero pCons_eq_0_iff power_Suc zero_neq_one)
       
  2646       done
       
  2647   qed
       
  2648   then show ?thesis
       
  2649     by (metis \<open>n = Suc n'\<close> pe)
       
  2650 qed
       
  2651 
       
  2652 lemma order_pderiv: "pderiv p \<noteq> 0 \<Longrightarrow> order a p \<noteq> 0 \<Longrightarrow> order a p = Suc (order a (pderiv p))"
       
  2653   for p :: "'a::field_char_0 poly"
       
  2654   apply (cases "p = 0")
       
  2655    apply simp
       
  2656   apply (drule_tac a = a and p = p in order_decomp)
       
  2657   using neq0_conv
       
  2658   apply (blast intro: lemma_order_pderiv)
       
  2659   done
       
  2660 
       
  2661 lemma poly_squarefree_decomp_order:
       
  2662   fixes p :: "'a::field_char_0 poly"
       
  2663   assumes "pderiv p \<noteq> 0"
       
  2664     and p: "p = q * d"
       
  2665     and p': "pderiv p = e * d"
       
  2666     and d: "d = r * p + s * pderiv p"
       
  2667   shows "order a q = (if order a p = 0 then 0 else 1)"
       
  2668 proof (rule classical)
       
  2669   assume 1: "\<not> ?thesis"
       
  2670   from \<open>pderiv p \<noteq> 0\<close> have "p \<noteq> 0" by auto
       
  2671   with p have "order a p = order a q + order a d"
       
  2672     by (simp add: order_mult)
       
  2673   with 1 have "order a p \<noteq> 0"
       
  2674     by (auto split: if_splits)
       
  2675   from \<open>pderiv p \<noteq> 0\<close> \<open>pderiv p = e * d\<close> have "order a (pderiv p) = order a e + order a d"
       
  2676     by (simp add: order_mult)
       
  2677   from \<open>pderiv p \<noteq> 0\<close> \<open>order a p \<noteq> 0\<close> have "order a p = Suc (order a (pderiv p))"
       
  2678     by (rule order_pderiv)
       
  2679   from \<open>p \<noteq> 0\<close> \<open>p = q * d\<close> have "d \<noteq> 0"
       
  2680     by simp
       
  2681   have "([:-a, 1:] ^ (order a (pderiv p))) dvd d"
       
  2682     apply (simp add: d)
       
  2683     apply (rule dvd_add)
       
  2684      apply (rule dvd_mult)
       
  2685      apply (simp add: order_divides \<open>p \<noteq> 0\<close> \<open>order a p = Suc (order a (pderiv p))\<close>)
       
  2686     apply (rule dvd_mult)
       
  2687     apply (simp add: order_divides)
       
  2688     done
       
  2689   with \<open>d \<noteq> 0\<close> have "order a (pderiv p) \<le> order a d"
       
  2690     by (simp add: order_divides)
       
  2691   show ?thesis
       
  2692     using \<open>order a p = order a q + order a d\<close>
       
  2693       and \<open>order a (pderiv p) = order a e + order a d\<close>
       
  2694       and \<open>order a p = Suc (order a (pderiv p))\<close>
       
  2695       and \<open>order a (pderiv p) \<le> order a d\<close>
       
  2696     by auto
       
  2697 qed
       
  2698 
       
  2699 lemma poly_squarefree_decomp_order2:
       
  2700   "pderiv p \<noteq> 0 \<Longrightarrow> p = q * d \<Longrightarrow> pderiv p = e * d \<Longrightarrow>
       
  2701     d = r * p + s * pderiv p \<Longrightarrow> \<forall>a. order a q = (if order a p = 0 then 0 else 1)"
       
  2702   for p :: "'a::field_char_0 poly"
       
  2703   by (blast intro: poly_squarefree_decomp_order)
       
  2704 
       
  2705 lemma order_pderiv2:
       
  2706   "pderiv p \<noteq> 0 \<Longrightarrow> order a p \<noteq> 0 \<Longrightarrow> order a (pderiv p) = n \<longleftrightarrow> order a p = Suc n"
       
  2707   for p :: "'a::field_char_0 poly"
       
  2708   by (auto dest: order_pderiv)
       
  2709 
       
  2710 definition rsquarefree :: "'a::idom poly \<Rightarrow> bool"
       
  2711   where "rsquarefree p \<longleftrightarrow> p \<noteq> 0 \<and> (\<forall>a. order a p = 0 \<or> order a p = 1)"
       
  2712 
       
  2713 lemma pderiv_iszero: "pderiv p = 0 \<Longrightarrow> \<exists>h. p = [:h:]"
       
  2714   for p :: "'a::{semidom,semiring_char_0} poly"
       
  2715   by (cases p) (auto simp: pderiv_eq_0_iff split: if_splits)
       
  2716 
       
  2717 lemma rsquarefree_roots: "rsquarefree p \<longleftrightarrow> (\<forall>a. \<not> (poly p a = 0 \<and> poly (pderiv p) a = 0))"
       
  2718   for p :: "'a::field_char_0 poly"
       
  2719   apply (simp add: rsquarefree_def)
       
  2720   apply (case_tac "p = 0")
       
  2721    apply simp
       
  2722   apply simp
       
  2723   apply (case_tac "pderiv p = 0")
       
  2724    apply simp
       
  2725    apply (drule pderiv_iszero, clarsimp)
       
  2726    apply (metis coeff_0 coeff_pCons_0 degree_pCons_0 le0 le_antisym order_degree)
       
  2727   apply (force simp add: order_root order_pderiv2)
       
  2728   done
       
  2729 
       
  2730 lemma poly_squarefree_decomp:
       
  2731   fixes p :: "'a::field_char_0 poly"
       
  2732   assumes "pderiv p \<noteq> 0"
       
  2733     and "p = q * d"
       
  2734     and "pderiv p = e * d"
       
  2735     and "d = r * p + s * pderiv p"
       
  2736   shows "rsquarefree q \<and> (\<forall>a. poly q a = 0 \<longleftrightarrow> poly p a = 0)"
       
  2737 proof -
       
  2738   from \<open>pderiv p \<noteq> 0\<close> have "p \<noteq> 0" by auto
       
  2739   with \<open>p = q * d\<close> have "q \<noteq> 0" by simp
       
  2740   from assms have "\<forall>a. order a q = (if order a p = 0 then 0 else 1)"
       
  2741     by (rule poly_squarefree_decomp_order2)
       
  2742   with \<open>p \<noteq> 0\<close> \<open>q \<noteq> 0\<close> show ?thesis
       
  2743     by (simp add: rsquarefree_def order_root)
       
  2744 qed
       
  2745 
       
  2746 
       
  2747 subsection \<open>Algebraic numbers\<close>
       
  2748 
       
  2749 text \<open>
       
  2750   Algebraic numbers can be defined in two equivalent ways: all real numbers that are
       
  2751   roots of rational polynomials or of integer polynomials. The Algebraic-Numbers AFP entry
       
  2752   uses the rational definition, but we need the integer definition.
       
  2753 
       
  2754   The equivalence is obvious since any rational polynomial can be multiplied with the
       
  2755   LCM of its coefficients, yielding an integer polynomial with the same roots.
       
  2756 \<close>
       
  2757 
       
  2758 definition algebraic :: "'a :: field_char_0 \<Rightarrow> bool"
       
  2759   where "algebraic x \<longleftrightarrow> (\<exists>p. (\<forall>i. coeff p i \<in> \<int>) \<and> p \<noteq> 0 \<and> poly p x = 0)"
       
  2760 
       
  2761 lemma algebraicI: "(\<And>i. coeff p i \<in> \<int>) \<Longrightarrow> p \<noteq> 0 \<Longrightarrow> poly p x = 0 \<Longrightarrow> algebraic x"
       
  2762   unfolding algebraic_def by blast
       
  2763 
       
  2764 lemma algebraicE:
       
  2765   assumes "algebraic x"
       
  2766   obtains p where "\<And>i. coeff p i \<in> \<int>" "p \<noteq> 0" "poly p x = 0"
       
  2767   using assms unfolding algebraic_def by blast
       
  2768 
       
  2769 lemma algebraic_altdef: "algebraic x \<longleftrightarrow> (\<exists>p. (\<forall>i. coeff p i \<in> \<rat>) \<and> p \<noteq> 0 \<and> poly p x = 0)"
       
  2770   for p :: "'a::field_char_0 poly"
       
  2771 proof safe
       
  2772   fix p
       
  2773   assume rat: "\<forall>i. coeff p i \<in> \<rat>" and root: "poly p x = 0" and nz: "p \<noteq> 0"
       
  2774   define cs where "cs = coeffs p"
       
  2775   from rat have "\<forall>c\<in>range (coeff p). \<exists>c'. c = of_rat c'"
       
  2776     unfolding Rats_def by blast
       
  2777   then obtain f where f: "coeff p i = of_rat (f (coeff p i))" for i
       
  2778     by (subst (asm) bchoice_iff) blast
       
  2779   define cs' where "cs' = map (quotient_of \<circ> f) (coeffs p)"
       
  2780   define d where "d = Lcm (set (map snd cs'))"
       
  2781   define p' where "p' = smult (of_int d) p"
       
  2782 
       
  2783   have "coeff p' n \<in> \<int>" for n
       
  2784   proof (cases "n \<le> degree p")
       
  2785     case True
       
  2786     define c where "c = coeff p n"
       
  2787     define a where "a = fst (quotient_of (f (coeff p n)))"
       
  2788     define b where "b = snd (quotient_of (f (coeff p n)))"
       
  2789     have b_pos: "b > 0"
       
  2790       unfolding b_def using quotient_of_denom_pos' by simp
       
  2791     have "coeff p' n = of_int d * coeff p n"
       
  2792       by (simp add: p'_def)
       
  2793     also have "coeff p n = of_rat (of_int a / of_int b)"
       
  2794       unfolding a_def b_def
       
  2795       by (subst quotient_of_div [of "f (coeff p n)", symmetric]) (simp_all add: f [symmetric])
       
  2796     also have "of_int d * \<dots> = of_rat (of_int (a*d) / of_int b)"
       
  2797       by (simp add: of_rat_mult of_rat_divide)
       
  2798     also from nz True have "b \<in> snd ` set cs'"
       
  2799       by (force simp: cs'_def o_def b_def coeffs_def simp del: upt_Suc)
       
  2800     then have "b dvd (a * d)"
       
  2801       by (simp add: d_def)
       
  2802     then have "of_int (a * d) / of_int b \<in> (\<int> :: rat set)"
       
  2803       by (rule of_int_divide_in_Ints)
       
  2804     then have "of_rat (of_int (a * d) / of_int b) \<in> \<int>" by (elim Ints_cases) auto
       
  2805     finally show ?thesis .
       
  2806   next
       
  2807     case False
       
  2808     then show ?thesis
       
  2809       by (auto simp: p'_def not_le coeff_eq_0)
       
  2810   qed
       
  2811   moreover have "set (map snd cs') \<subseteq> {0<..}"
       
  2812     unfolding cs'_def using quotient_of_denom_pos' by (auto simp: coeffs_def simp del: upt_Suc)
       
  2813   then have "d \<noteq> 0"
       
  2814     unfolding d_def by (induct cs') simp_all
       
  2815   with nz have "p' \<noteq> 0" by (simp add: p'_def)
       
  2816   moreover from root have "poly p' x = 0"
       
  2817     by (simp add: p'_def)
       
  2818   ultimately show "algebraic x"
       
  2819     unfolding algebraic_def by blast
       
  2820 next
       
  2821   assume "algebraic x"
       
  2822   then obtain p where p: "coeff p i \<in> \<int>" "poly p x = 0" "p \<noteq> 0" for i
       
  2823     by (force simp: algebraic_def)
       
  2824   moreover have "coeff p i \<in> \<int> \<Longrightarrow> coeff p i \<in> \<rat>" for i
       
  2825     by (elim Ints_cases) simp
       
  2826   ultimately show "\<exists>p. (\<forall>i. coeff p i \<in> \<rat>) \<and> p \<noteq> 0 \<and> poly p x = 0" by auto
       
  2827 qed
       
  2828 
       
  2829 
       
  2830 subsection \<open>Content and primitive part of a polynomial\<close>
       
  2831 
       
  2832 definition content :: "'a::semiring_gcd poly \<Rightarrow> 'a"
       
  2833   where "content p = gcd_list (coeffs p)"
       
  2834 
       
  2835 lemma content_eq_fold_coeffs [code]: "content p = fold_coeffs gcd p 0"
       
  2836   by (simp add: content_def Gcd_fin.set_eq_fold fold_coeffs_def foldr_fold fun_eq_iff ac_simps)
       
  2837 
       
  2838 lemma content_0 [simp]: "content 0 = 0"
       
  2839   by (simp add: content_def)
       
  2840 
       
  2841 lemma content_1 [simp]: "content 1 = 1"
       
  2842   by (simp add: content_def)
       
  2843 
       
  2844 lemma content_const [simp]: "content [:c:] = normalize c"
       
  2845   by (simp add: content_def cCons_def)
       
  2846 
       
  2847 lemma const_poly_dvd_iff_dvd_content: "[:c:] dvd p \<longleftrightarrow> c dvd content p"
       
  2848   for c :: "'a::semiring_gcd"
       
  2849 proof (cases "p = 0")
       
  2850   case True
       
  2851   then show ?thesis by simp
       
  2852 next
       
  2853   case False
       
  2854   have "[:c:] dvd p \<longleftrightarrow> (\<forall>n. c dvd coeff p n)"
       
  2855     by (rule const_poly_dvd_iff)
       
  2856   also have "\<dots> \<longleftrightarrow> (\<forall>a\<in>set (coeffs p). c dvd a)"
       
  2857   proof safe
       
  2858     fix n :: nat
       
  2859     assume "\<forall>a\<in>set (coeffs p). c dvd a"
       
  2860     then show "c dvd coeff p n"
       
  2861       by (cases "n \<le> degree p") (auto simp: coeff_eq_0 coeffs_def split: if_splits)
       
  2862   qed (auto simp: coeffs_def simp del: upt_Suc split: if_splits)
       
  2863   also have "\<dots> \<longleftrightarrow> c dvd content p"
       
  2864     by (simp add: content_def dvd_Gcd_fin_iff dvd_mult_unit_iff)
       
  2865   finally show ?thesis .
       
  2866 qed
       
  2867 
       
  2868 lemma content_dvd [simp]: "[:content p:] dvd p"
       
  2869   by (subst const_poly_dvd_iff_dvd_content) simp_all
       
  2870 
       
  2871 lemma content_dvd_coeff [simp]: "content p dvd coeff p n"
       
  2872 proof (cases "p = 0")
       
  2873   case True
       
  2874   then show ?thesis
       
  2875     by simp
       
  2876 next
       
  2877   case False
       
  2878   then show ?thesis
       
  2879     by (cases "n \<le> degree p")
       
  2880       (auto simp add: content_def not_le coeff_eq_0 coeff_in_coeffs intro: Gcd_fin_dvd)
       
  2881 qed
       
  2882 
       
  2883 lemma content_dvd_coeffs: "c \<in> set (coeffs p) \<Longrightarrow> content p dvd c"
       
  2884   by (simp add: content_def Gcd_fin_dvd)
       
  2885 
       
  2886 lemma normalize_content [simp]: "normalize (content p) = content p"
       
  2887   by (simp add: content_def)
       
  2888 
       
  2889 lemma is_unit_content_iff [simp]: "is_unit (content p) \<longleftrightarrow> content p = 1"
       
  2890 proof
       
  2891   assume "is_unit (content p)"
       
  2892   then have "normalize (content p) = 1" by (simp add: is_unit_normalize del: normalize_content)
       
  2893   then show "content p = 1" by simp
       
  2894 qed auto
       
  2895 
       
  2896 lemma content_smult [simp]: "content (smult c p) = normalize c * content p"
       
  2897   by (simp add: content_def coeffs_smult Gcd_fin_mult)
       
  2898 
       
  2899 lemma content_eq_zero_iff [simp]: "content p = 0 \<longleftrightarrow> p = 0"
       
  2900   by (auto simp: content_def simp: poly_eq_iff coeffs_def)
       
  2901 
       
  2902 definition primitive_part :: "'a :: semiring_gcd poly \<Rightarrow> 'a poly"
       
  2903   where "primitive_part p = map_poly (\<lambda>x. x div content p) p"
       
  2904 
       
  2905 lemma primitive_part_0 [simp]: "primitive_part 0 = 0"
       
  2906   by (simp add: primitive_part_def)
       
  2907 
       
  2908 lemma content_times_primitive_part [simp]: "smult (content p) (primitive_part p) = p"
       
  2909   for p :: "'a :: semiring_gcd poly"
       
  2910 proof (cases "p = 0")
       
  2911   case True
       
  2912   then show ?thesis by simp
       
  2913 next
       
  2914   case False
       
  2915   then show ?thesis
       
  2916   unfolding primitive_part_def
       
  2917   by (auto simp: smult_conv_map_poly map_poly_map_poly o_def content_dvd_coeffs
       
  2918       intro: map_poly_idI)
       
  2919 qed
       
  2920 
       
  2921 lemma primitive_part_eq_0_iff [simp]: "primitive_part p = 0 \<longleftrightarrow> p = 0"
       
  2922 proof (cases "p = 0")
       
  2923   case True
       
  2924   then show ?thesis by simp
       
  2925 next
       
  2926   case False
       
  2927   then have "primitive_part p = map_poly (\<lambda>x. x div content p) p"
       
  2928     by (simp add:  primitive_part_def)
       
  2929   also from False have "\<dots> = 0 \<longleftrightarrow> p = 0"
       
  2930     by (intro map_poly_eq_0_iff) (auto simp: dvd_div_eq_0_iff content_dvd_coeffs)
       
  2931   finally show ?thesis
       
  2932     using False by simp
       
  2933 qed
       
  2934 
       
  2935 lemma content_primitive_part [simp]:
       
  2936   assumes "p \<noteq> 0"
       
  2937   shows "content (primitive_part p) = 1"
       
  2938 proof -
       
  2939   have "p = smult (content p) (primitive_part p)"
       
  2940     by simp
       
  2941   also have "content \<dots> = content (primitive_part p) * content p"
       
  2942     by (simp del: content_times_primitive_part add: ac_simps)
       
  2943   finally have "1 * content p = content (primitive_part p) * content p"
       
  2944     by simp
       
  2945   then have "1 * content p div content p = content (primitive_part p) * content p div content p"
       
  2946     by simp
       
  2947   with assms show ?thesis
       
  2948     by simp
       
  2949 qed
       
  2950 
       
  2951 lemma content_decompose:
       
  2952   obtains p' :: "'a::semiring_gcd poly" where "p = smult (content p) p'" "content p' = 1"
       
  2953 proof (cases "p = 0")
       
  2954   case True
       
  2955   then show ?thesis by (intro that[of 1]) simp_all
       
  2956 next
       
  2957   case False
       
  2958   from content_dvd[of p] obtain r where r: "p = [:content p:] * r"
       
  2959     by (rule dvdE)
       
  2960   have "content p * 1 = content p * content r"
       
  2961     by (subst r) simp
       
  2962   with False have "content r = 1"
       
  2963     by (subst (asm) mult_left_cancel) simp_all
       
  2964   with r show ?thesis
       
  2965     by (intro that[of r]) simp_all
       
  2966 qed
       
  2967 
       
  2968 lemma content_dvd_contentI [intro]: "p dvd q \<Longrightarrow> content p dvd content q"
       
  2969   using const_poly_dvd_iff_dvd_content content_dvd dvd_trans by blast
       
  2970 
       
  2971 lemma primitive_part_const_poly [simp]: "primitive_part [:x:] = [:unit_factor x:]"
       
  2972   by (simp add: primitive_part_def map_poly_pCons)
       
  2973 
       
  2974 lemma primitive_part_prim: "content p = 1 \<Longrightarrow> primitive_part p = p"
       
  2975   by (auto simp: primitive_part_def)
       
  2976 
       
  2977 lemma degree_primitive_part [simp]: "degree (primitive_part p) = degree p"
       
  2978 proof (cases "p = 0")
       
  2979   case True
       
  2980   then show ?thesis by simp
       
  2981 next
       
  2982   case False
       
  2983   have "p = smult (content p) (primitive_part p)"
       
  2984     by simp
       
  2985   also from False have "degree \<dots> = degree (primitive_part p)"
       
  2986     by (subst degree_smult_eq) simp_all
       
  2987   finally show ?thesis ..
       
  2988 qed
       
  2989 
       
  2990 
       
  2991 subsection \<open>Division of polynomials\<close>
       
  2992 
       
  2993 subsubsection \<open>Division in general\<close>
       
  2994 
       
  2995 instantiation poly :: (idom_divide) idom_divide
       
  2996 begin
       
  2997 
       
  2998 fun divide_poly_main :: "'a \<Rightarrow> 'a poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> 'a poly"
       
  2999   where
       
  3000     "divide_poly_main lc q r d dr (Suc n) =
       
  3001       (let cr = coeff r dr; a = cr div lc; mon = monom a n in
       
  3002         if False \<or> a * lc = cr then (* False \<or> is only because of problem in function-package *)
       
  3003           divide_poly_main
       
  3004             lc
       
  3005             (q + mon)
       
  3006             (r - mon * d)
       
  3007             d (dr - 1) n else 0)"
       
  3008   | "divide_poly_main lc q r d dr 0 = q"
       
  3009 
       
  3010 definition divide_poly :: "'a poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly"
       
  3011   where "divide_poly f g =
       
  3012     (if g = 0 then 0
       
  3013      else
       
  3014       divide_poly_main (coeff g (degree g)) 0 f g (degree f)
       
  3015         (1 + length (coeffs f) - length (coeffs g)))"
       
  3016 
       
  3017 lemma divide_poly_main:
       
  3018   assumes d: "d \<noteq> 0" "lc = coeff d (degree d)"
       
  3019     and "degree (d * r) \<le> dr" "divide_poly_main lc q (d * r) d dr n = q'"
       
  3020     and "n = 1 + dr - degree d \<or> dr = 0 \<and> n = 0 \<and> d * r = 0"
       
  3021   shows "q' = q + r"
       
  3022   using assms(3-)
       
  3023 proof (induct n arbitrary: q r dr)
       
  3024   case (Suc n)
       
  3025   let ?rr = "d * r"
       
  3026   let ?a = "coeff ?rr dr"
       
  3027   let ?qq = "?a div lc"
       
  3028   define b where [simp]: "b = monom ?qq n"
       
  3029   let ?rrr =  "d * (r - b)"
       
  3030   let ?qqq = "q + b"
       
  3031   note res = Suc(3)
       
  3032   from Suc(4) have dr: "dr = n + degree d" by auto
       
  3033   from d have lc: "lc \<noteq> 0" by auto
       
  3034   have "coeff (b * d) dr = coeff b n * coeff d (degree d)"
       
  3035   proof (cases "?qq = 0")
       
  3036     case True
       
  3037     then show ?thesis by simp
       
  3038   next
       
  3039     case False
       
  3040     then have n: "n = degree b"
       
  3041       by (simp add: degree_monom_eq)
       
  3042     show ?thesis
       
  3043       unfolding n dr by (simp add: coeff_mult_degree_sum)
       
  3044   qed
       
  3045   also have "\<dots> = lc * coeff b n"
       
  3046     by (simp add: d)
       
  3047   finally have c2: "coeff (b * d) dr = lc * coeff b n" .
       
  3048   have rrr: "?rrr = ?rr - b * d"
       
  3049     by (simp add: field_simps)
       
  3050   have c1: "coeff (d * r) dr = lc * coeff r n"
       
  3051   proof (cases "degree r = n")
       
  3052     case True
       
  3053     with Suc(2) show ?thesis
       
  3054       unfolding dr using coeff_mult_degree_sum[of d r] d by (auto simp: ac_simps)
       
  3055   next
       
  3056     case False
       
  3057     from dr Suc(2) have "degree r \<le> n"
       
  3058       by auto
       
  3059         (metis add.commute add_le_cancel_left d(1) degree_0 degree_mult_eq
       
  3060           diff_is_0_eq diff_zero le_cases)
       
  3061     with False have r_n: "degree r < n"
       
  3062       by auto
       
  3063     then have right: "lc * coeff r n = 0"
       
  3064       by (simp add: coeff_eq_0)
       
  3065     have "coeff (d * r) dr = coeff (d * r) (degree d + n)"
       
  3066       by (simp add: dr ac_simps)
       
  3067     also from r_n have "\<dots> = 0"
       
  3068       by (metis False Suc.prems(1) add.commute add_left_imp_eq coeff_degree_mult coeff_eq_0
       
  3069         coeff_mult_degree_sum degree_mult_le dr le_eq_less_or_eq)
       
  3070     finally show ?thesis
       
  3071       by (simp only: right)
       
  3072   qed
       
  3073   have c0: "coeff ?rrr dr = 0"
       
  3074     and id: "lc * (coeff (d * r) dr div lc) = coeff (d * r) dr"
       
  3075     unfolding rrr coeff_diff c2
       
  3076     unfolding b_def coeff_monom coeff_smult c1 using lc by auto
       
  3077   from res[unfolded divide_poly_main.simps[of lc q] Let_def] id
       
  3078   have res: "divide_poly_main lc ?qqq ?rrr d (dr - 1) n = q'"
       
  3079     by (simp del: divide_poly_main.simps add: field_simps)
       
  3080   note IH = Suc(1)[OF _ res]
       
  3081   from Suc(4) have dr: "dr = n + degree d" by auto
       
  3082   from Suc(2) have deg_rr: "degree ?rr \<le> dr" by auto
       
  3083   have deg_bd: "degree (b * d) \<le> dr"
       
  3084     unfolding dr b_def by (rule order.trans[OF degree_mult_le]) (auto simp: degree_monom_le)
       
  3085   have "degree ?rrr \<le> dr"
       
  3086     unfolding rrr by (rule degree_diff_le[OF deg_rr deg_bd])
       
  3087   with c0 have deg_rrr: "degree ?rrr \<le> (dr - 1)"
       
  3088     by (rule coeff_0_degree_minus_1)
       
  3089   have "n = 1 + (dr - 1) - degree d \<or> dr - 1 = 0 \<and> n = 0 \<and> ?rrr = 0"
       
  3090   proof (cases dr)
       
  3091     case 0
       
  3092     with Suc(4) have 0: "dr = 0" "n = 0" "degree d = 0"
       
  3093       by auto
       
  3094     with deg_rrr have "degree ?rrr = 0"
       
  3095       by simp
       
  3096     from degree_eq_zeroE[OF this] obtain a where rrr: "?rrr = [:a:]"
       
  3097       by metis
       
  3098     show ?thesis
       
  3099       unfolding 0 using c0 unfolding rrr 0 by simp
       
  3100   next
       
  3101     case _: Suc
       
  3102     with Suc(4) show ?thesis by auto
       
  3103   qed
       
  3104   from IH[OF deg_rrr this] show ?case
       
  3105     by simp
       
  3106 next
       
  3107   case 0
       
  3108   show ?case
       
  3109   proof (cases "r = 0")
       
  3110     case True
       
  3111     with 0 show ?thesis by auto
       
  3112   next
       
  3113     case False
       
  3114     from d False have "degree (d * r) = degree d + degree r"
       
  3115       by (subst degree_mult_eq) auto
       
  3116     with 0 d show ?thesis by auto
       
  3117   qed
       
  3118 qed
       
  3119 
       
  3120 lemma divide_poly_main_0: "divide_poly_main 0 0 r d dr n = 0"
       
  3121 proof (induct n arbitrary: r d dr)
       
  3122   case 0
       
  3123   then show ?case by simp
       
  3124 next
       
  3125   case Suc
       
  3126   show ?case
       
  3127     unfolding divide_poly_main.simps[of _ _ r] Let_def
       
  3128     by (simp add: Suc del: divide_poly_main.simps)
       
  3129 qed
       
  3130 
       
  3131 lemma divide_poly:
       
  3132   assumes g: "g \<noteq> 0"
       
  3133   shows "(f * g) div g = (f :: 'a poly)"
       
  3134 proof -
       
  3135   have len: "length (coeffs f) = Suc (degree f)" if "f \<noteq> 0" for f :: "'a poly"
       
  3136     using that unfolding degree_eq_length_coeffs by auto
       
  3137   have "divide_poly_main (coeff g (degree g)) 0 (g * f) g (degree (g * f))
       
  3138     (1 + length (coeffs (g * f)) - length (coeffs  g)) = (f * g) div g"
       
  3139     by (simp add: divide_poly_def Let_def ac_simps)
       
  3140   note main = divide_poly_main[OF g refl le_refl this]
       
  3141   have "(f * g) div g = 0 + f"
       
  3142   proof (rule main, goal_cases)
       
  3143     case 1
       
  3144     show ?case
       
  3145     proof (cases "f = 0")
       
  3146       case True
       
  3147       with g show ?thesis
       
  3148         by (auto simp: degree_eq_length_coeffs)
       
  3149     next
       
  3150       case False
       
  3151       with g have fg: "g * f \<noteq> 0" by auto
       
  3152       show ?thesis
       
  3153         unfolding len[OF fg] len[OF g] by auto
       
  3154     qed
       
  3155   qed
       
  3156   then show ?thesis by simp
       
  3157 qed
       
  3158 
       
  3159 lemma divide_poly_0: "f div 0 = 0"
       
  3160   for f :: "'a poly"
       
  3161   by (simp add: divide_poly_def Let_def divide_poly_main_0)
       
  3162 
       
  3163 instance
       
  3164   by standard (auto simp: divide_poly divide_poly_0)
       
  3165 
       
  3166 end
       
  3167 
       
  3168 instance poly :: (idom_divide) algebraic_semidom ..
       
  3169 
       
  3170 lemma div_const_poly_conv_map_poly:
       
  3171   assumes "[:c:] dvd p"
       
  3172   shows "p div [:c:] = map_poly (\<lambda>x. x div c) p"
       
  3173 proof (cases "c = 0")
       
  3174   case True
       
  3175   then show ?thesis
       
  3176     by (auto intro!: poly_eqI simp: coeff_map_poly)
       
  3177 next
       
  3178   case False
       
  3179   from assms obtain q where p: "p = [:c:] * q" by (rule dvdE)
       
  3180   moreover {
       
  3181     have "smult c q = [:c:] * q"
       
  3182       by simp
       
  3183     also have "\<dots> div [:c:] = q"
       
  3184       by (rule nonzero_mult_div_cancel_left) (use False in auto)
       
  3185     finally have "smult c q div [:c:] = q" .
       
  3186   }
       
  3187   ultimately show ?thesis by (intro poly_eqI) (auto simp: coeff_map_poly False)
       
  3188 qed
       
  3189 
       
  3190 lemma is_unit_monom_0:
       
  3191   fixes a :: "'a::field"
       
  3192   assumes "a \<noteq> 0"
       
  3193   shows "is_unit (monom a 0)"
       
  3194 proof
       
  3195   from assms show "1 = monom a 0 * monom (inverse a) 0"
       
  3196     by (simp add: mult_monom)
       
  3197 qed
       
  3198 
       
  3199 lemma is_unit_triv: "a \<noteq> 0 \<Longrightarrow> is_unit [:a:]"
       
  3200   for a :: "'a::field"
       
  3201   by (simp add: is_unit_monom_0 monom_0 [symmetric])
       
  3202 
       
  3203 lemma is_unit_iff_degree:
       
  3204   fixes p :: "'a::field poly"
       
  3205   assumes "p \<noteq> 0"
       
  3206   shows "is_unit p \<longleftrightarrow> degree p = 0"
       
  3207     (is "?lhs \<longleftrightarrow> ?rhs")
       
  3208 proof
       
  3209   assume ?rhs
       
  3210   then obtain a where "p = [:a:]"
       
  3211     by (rule degree_eq_zeroE)
       
  3212   with assms show ?lhs
       
  3213     by (simp add: is_unit_triv)
       
  3214 next
       
  3215   assume ?lhs
       
  3216   then obtain q where "q \<noteq> 0" "p * q = 1" ..
       
  3217   then have "degree (p * q) = degree 1"
       
  3218     by simp
       
  3219   with \<open>p \<noteq> 0\<close> \<open>q \<noteq> 0\<close> have "degree p + degree q = 0"
       
  3220     by (simp add: degree_mult_eq)
       
  3221   then show ?rhs by simp
       
  3222 qed
       
  3223 
       
  3224 lemma is_unit_pCons_iff: "is_unit (pCons a p) \<longleftrightarrow> p = 0 \<and> a \<noteq> 0"
       
  3225   for p :: "'a::field poly"
       
  3226   by (cases "p = 0") (auto simp: is_unit_triv is_unit_iff_degree)
       
  3227 
       
  3228 lemma is_unit_monom_trival: "is_unit p \<Longrightarrow> monom (coeff p (degree p)) 0 = p"
       
  3229   for p :: "'a::field poly"
       
  3230   by (cases p) (simp_all add: monom_0 is_unit_pCons_iff)
       
  3231 
       
  3232 lemma is_unit_const_poly_iff: "[:c:] dvd 1 \<longleftrightarrow> c dvd 1"
       
  3233   for c :: "'a::{comm_semiring_1,semiring_no_zero_divisors}"
       
  3234   by (auto simp: one_poly_def)
       
  3235 
       
  3236 lemma is_unit_polyE:
       
  3237   fixes p :: "'a :: {comm_semiring_1,semiring_no_zero_divisors} poly"
       
  3238   assumes "p dvd 1"
       
  3239   obtains c where "p = [:c:]" "c dvd 1"
       
  3240 proof -
       
  3241   from assms obtain q where "1 = p * q"
       
  3242     by (rule dvdE)
       
  3243   then have "p \<noteq> 0" and "q \<noteq> 0"
       
  3244     by auto
       
  3245   from \<open>1 = p * q\<close> have "degree 1 = degree (p * q)"
       
  3246     by simp
       
  3247   also from \<open>p \<noteq> 0\<close> and \<open>q \<noteq> 0\<close> have "\<dots> = degree p + degree q"
       
  3248     by (simp add: degree_mult_eq)
       
  3249   finally have "degree p = 0" by simp
       
  3250   with degree_eq_zeroE obtain c where c: "p = [:c:]" .
       
  3251   with \<open>p dvd 1\<close> have "c dvd 1"
       
  3252     by (simp add: is_unit_const_poly_iff)
       
  3253   with c show thesis ..
       
  3254 qed
       
  3255 
       
  3256 lemma is_unit_polyE':
       
  3257   fixes p :: "'a::field poly"
       
  3258   assumes "is_unit p"
       
  3259   obtains a where "p = monom a 0" and "a \<noteq> 0"
       
  3260 proof -
       
  3261   obtain a q where "p = pCons a q"
       
  3262     by (cases p)
       
  3263   with assms have "p = [:a:]" and "a \<noteq> 0"
       
  3264     by (simp_all add: is_unit_pCons_iff)
       
  3265   with that show thesis by (simp add: monom_0)
       
  3266 qed
       
  3267 
       
  3268 lemma is_unit_poly_iff: "p dvd 1 \<longleftrightarrow> (\<exists>c. p = [:c:] \<and> c dvd 1)"
       
  3269   for p :: "'a::{comm_semiring_1,semiring_no_zero_divisors} poly"
       
  3270   by (auto elim: is_unit_polyE simp add: is_unit_const_poly_iff)
       
  3271 
       
  3272 
       
  3273 subsubsection \<open>Pseudo-Division\<close>
       
  3274 
       
  3275 text \<open>This part is by René Thiemann and Akihisa Yamada.\<close>
       
  3276 
       
  3277 fun pseudo_divmod_main ::
       
  3278   "'a :: comm_ring_1  \<Rightarrow> 'a poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly \<Rightarrow> nat \<Rightarrow> nat \<Rightarrow> 'a poly \<times> 'a poly"
       
  3279   where
       
  3280     "pseudo_divmod_main lc q r d dr (Suc n) =
       
  3281       (let
       
  3282         rr = smult lc r;
       
  3283         qq = coeff r dr;
       
  3284         rrr = rr - monom qq n * d;
       
  3285         qqq = smult lc q + monom qq n
       
  3286        in pseudo_divmod_main lc qqq rrr d (dr - 1) n)"
       
  3287   | "pseudo_divmod_main lc q r d dr 0 = (q,r)"
       
  3288 
       
  3289 definition pseudo_divmod :: "'a :: comm_ring_1 poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly \<times> 'a poly"
       
  3290   where "pseudo_divmod p q \<equiv>
       
  3291     if q = 0 then (0, p)
       
  3292     else
       
  3293       pseudo_divmod_main (coeff q (degree q)) 0 p q (degree p)
       
  3294         (1 + length (coeffs p) - length (coeffs q))"
       
  3295 
       
  3296 lemma pseudo_divmod_main:
       
  3297   assumes d: "d \<noteq> 0" "lc = coeff d (degree d)"
       
  3298     and "degree r \<le> dr" "pseudo_divmod_main lc q r d dr n = (q',r')"
       
  3299     and "n = 1 + dr - degree d \<or> dr = 0 \<and> n = 0 \<and> r = 0"
       
  3300   shows "(r' = 0 \<or> degree r' < degree d) \<and> smult (lc^n) (d * q + r) = d * q' + r'"
       
  3301   using assms(3-)
       
  3302 proof (induct n arbitrary: q r dr)
       
  3303   case 0
       
  3304   then show ?case by auto
       
  3305 next
       
  3306   case (Suc n)
       
  3307   let ?rr = "smult lc r"
       
  3308   let ?qq = "coeff r dr"
       
  3309   define b where [simp]: "b = monom ?qq n"
       
  3310   let ?rrr = "?rr - b * d"
       
  3311   let ?qqq = "smult lc q + b"
       
  3312   note res = Suc(3)
       
  3313   from res[unfolded pseudo_divmod_main.simps[of lc q] Let_def]
       
  3314   have res: "pseudo_divmod_main lc ?qqq ?rrr d (dr - 1) n = (q',r')"
       
  3315     by (simp del: pseudo_divmod_main.simps)
       
  3316   from Suc(4) have dr: "dr = n + degree d" by auto
       
  3317   have "coeff (b * d) dr = coeff b n * coeff d (degree d)"
       
  3318   proof (cases "?qq = 0")
       
  3319     case True
       
  3320     then show ?thesis by auto
       
  3321   next
       
  3322     case False
       
  3323     then have n: "n = degree b"
       
  3324       by (simp add: degree_monom_eq)
       
  3325     show ?thesis
       
  3326       unfolding n dr by (simp add: coeff_mult_degree_sum)
       
  3327   qed
       
  3328   also have "\<dots> = lc * coeff b n" by (simp add: d)
       
  3329   finally have "coeff (b * d) dr = lc * coeff b n" .
       
  3330   moreover have "coeff ?rr dr = lc * coeff r dr"
       
  3331     by simp
       
  3332   ultimately have c0: "coeff ?rrr dr = 0"
       
  3333     by auto
       
  3334   from Suc(4) have dr: "dr = n + degree d" by auto
       
  3335   have deg_rr: "degree ?rr \<le> dr"
       
  3336     using Suc(2) degree_smult_le dual_order.trans by blast
       
  3337   have deg_bd: "degree (b * d) \<le> dr"
       
  3338     unfolding dr by (rule order.trans[OF degree_mult_le]) (auto simp: degree_monom_le)
       
  3339   have "degree ?rrr \<le> dr"
       
  3340     using degree_diff_le[OF deg_rr deg_bd] by auto
       
  3341   with c0 have deg_rrr: "degree ?rrr \<le> (dr - 1)"
       
  3342     by (rule coeff_0_degree_minus_1)
       
  3343   have "n = 1 + (dr - 1) - degree d \<or> dr - 1 = 0 \<and> n = 0 \<and> ?rrr = 0"
       
  3344   proof (cases dr)
       
  3345     case 0
       
  3346     with Suc(4) have 0: "dr = 0" "n = 0" "degree d = 0" by auto
       
  3347     with deg_rrr have "degree ?rrr = 0" by simp
       
  3348     then have "\<exists>a. ?rrr = [:a:]"
       
  3349       by (metis degree_pCons_eq_if old.nat.distinct(2) pCons_cases)
       
  3350     from this obtain a where rrr: "?rrr = [:a:]"
       
  3351       by auto
       
  3352     show ?thesis
       
  3353       unfolding 0 using c0 unfolding rrr 0 by simp
       
  3354   next
       
  3355     case _: Suc
       
  3356     with Suc(4) show ?thesis by auto
       
  3357   qed
       
  3358   note IH = Suc(1)[OF deg_rrr res this]
       
  3359   show ?case
       
  3360   proof (intro conjI)
       
  3361     from IH show "r' = 0 \<or> degree r' < degree d"
       
  3362       by blast
       
  3363     show "smult (lc ^ Suc n) (d * q + r) = d * q' + r'"
       
  3364       unfolding IH[THEN conjunct2,symmetric]
       
  3365       by (simp add: field_simps smult_add_right)
       
  3366   qed
       
  3367 qed
       
  3368 
       
  3369 lemma pseudo_divmod:
       
  3370   assumes g: "g \<noteq> 0"
       
  3371     and *: "pseudo_divmod f g = (q,r)"
       
  3372   shows "smult (coeff g (degree g) ^ (Suc (degree f) - degree g)) f = g * q + r"  (is ?A)
       
  3373     and "r = 0 \<or> degree r < degree g"  (is ?B)
       
  3374 proof -
       
  3375   from *[unfolded pseudo_divmod_def Let_def]
       
  3376   have "pseudo_divmod_main (coeff g (degree g)) 0 f g (degree f)
       
  3377       (1 + length (coeffs f) - length (coeffs g)) = (q, r)"
       
  3378     by (auto simp: g)
       
  3379   note main = pseudo_divmod_main[OF _ _ _ this, OF g refl le_refl]
       
  3380   from g have "1 + length (coeffs f) - length (coeffs g) = 1 + degree f - degree g \<or>
       
  3381     degree f = 0 \<and> 1 + length (coeffs f) - length (coeffs g) = 0 \<and> f = 0"
       
  3382     by (cases "f = 0"; cases "coeffs g") (auto simp: degree_eq_length_coeffs)
       
  3383   note main' = main[OF this]
       
  3384   then show "r = 0 \<or> degree r < degree g" by auto
       
  3385   show "smult (coeff g (degree g) ^ (Suc (degree f) - degree g)) f = g * q + r"
       
  3386     by (subst main'[THEN conjunct2, symmetric], simp add: degree_eq_length_coeffs,
       
  3387         cases "f = 0"; cases "coeffs g", use g in auto)
       
  3388 qed
       
  3389 
       
  3390 definition "pseudo_mod_main lc r d dr n = snd (pseudo_divmod_main lc 0 r d dr n)"
       
  3391 
       
  3392 lemma snd_pseudo_divmod_main:
       
  3393   "snd (pseudo_divmod_main lc q r d dr n) = snd (pseudo_divmod_main lc q' r d dr n)"
       
  3394   by (induct n arbitrary: q q' lc r d dr) (simp_all add: Let_def)
       
  3395 
       
  3396 definition pseudo_mod :: "'a::{comm_ring_1,semiring_1_no_zero_divisors} poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly"
       
  3397   where "pseudo_mod f g = snd (pseudo_divmod f g)"
       
  3398 
       
  3399 lemma pseudo_mod:
       
  3400   fixes f g :: "'a::{comm_ring_1,semiring_1_no_zero_divisors} poly"
       
  3401   defines "r \<equiv> pseudo_mod f g"
       
  3402   assumes g: "g \<noteq> 0"
       
  3403   shows "\<exists>a q. a \<noteq> 0 \<and> smult a f = g * q + r" "r = 0 \<or> degree r < degree g"
       
  3404 proof -
       
  3405   let ?cg = "coeff g (degree g)"
       
  3406   let ?cge = "?cg ^ (Suc (degree f) - degree g)"
       
  3407   define a where "a = ?cge"
       
  3408   from r_def[unfolded pseudo_mod_def] obtain q where pdm: "pseudo_divmod f g = (q, r)"
       
  3409     by (cases "pseudo_divmod f g") auto
       
  3410   from pseudo_divmod[OF g pdm] have id: "smult a f = g * q + r" and "r = 0 \<or> degree r < degree g"
       
  3411     by (auto simp: a_def)
       
  3412   show "r = 0 \<or> degree r < degree g" by fact
       
  3413   from g have "a \<noteq> 0"
       
  3414     by (auto simp: a_def)
       
  3415   with id show "\<exists>a q. a \<noteq> 0 \<and> smult a f = g * q + r"
       
  3416     by auto
       
  3417 qed
       
  3418 
       
  3419 lemma fst_pseudo_divmod_main_as_divide_poly_main:
       
  3420   assumes d: "d \<noteq> 0"
       
  3421   defines lc: "lc \<equiv> coeff d (degree d)"
       
  3422   shows "fst (pseudo_divmod_main lc q r d dr n) =
       
  3423     divide_poly_main lc (smult (lc^n) q) (smult (lc^n) r) d dr n"
       
  3424 proof (induct n arbitrary: q r dr)
       
  3425   case 0
       
  3426   then show ?case by simp
       
  3427 next
       
  3428   case (Suc n)
       
  3429   note lc0 = leading_coeff_neq_0[OF d, folded lc]
       
  3430   then have "pseudo_divmod_main lc q r d dr (Suc n) =
       
  3431     pseudo_divmod_main lc (smult lc q + monom (coeff r dr) n)
       
  3432       (smult lc r - monom (coeff r dr) n * d) d (dr - 1) n"
       
  3433     by (simp add: Let_def ac_simps)
       
  3434   also have "fst \<dots> = divide_poly_main lc
       
  3435      (smult (lc^n) (smult lc q + monom (coeff r dr) n))
       
  3436      (smult (lc^n) (smult lc r - monom (coeff r dr) n * d))
       
  3437      d (dr - 1) n"
       
  3438     by (simp only: Suc[unfolded divide_poly_main.simps Let_def])
       
  3439   also have "\<dots> = divide_poly_main lc (smult (lc ^ Suc n) q) (smult (lc ^ Suc n) r) d dr (Suc n)"
       
  3440     unfolding smult_monom smult_distribs mult_smult_left[symmetric]
       
  3441     using lc0 by (simp add: Let_def ac_simps)
       
  3442   finally show ?case .
       
  3443 qed
       
  3444 
       
  3445 
       
  3446 subsubsection \<open>Division in polynomials over fields\<close>
       
  3447 
       
  3448 lemma pseudo_divmod_field:
       
  3449   fixes g :: "'a::field poly"
       
  3450   assumes g: "g \<noteq> 0"
       
  3451     and *: "pseudo_divmod f g = (q,r)"
       
  3452   defines "c \<equiv> coeff g (degree g) ^ (Suc (degree f) - degree g)"
       
  3453   shows "f = g * smult (1/c) q + smult (1/c) r"
       
  3454 proof -
       
  3455   from leading_coeff_neq_0[OF g] have c0: "c \<noteq> 0"
       
  3456     by (auto simp: c_def)
       
  3457   from pseudo_divmod(1)[OF g *, folded c_def] have "smult c f = g * q + r"
       
  3458     by auto
       
  3459   also have "smult (1 / c) \<dots> = g * smult (1 / c) q + smult (1 / c) r"
       
  3460     by (simp add: smult_add_right)
       
  3461   finally show ?thesis
       
  3462     using c0 by auto
       
  3463 qed
       
  3464 
       
  3465 lemma divide_poly_main_field:
       
  3466   fixes d :: "'a::field poly"
       
  3467   assumes d: "d \<noteq> 0"
       
  3468   defines lc: "lc \<equiv> coeff d (degree d)"
       
  3469   shows "divide_poly_main lc q r d dr n =
       
  3470     fst (pseudo_divmod_main lc (smult ((1 / lc)^n) q) (smult ((1 / lc)^n) r) d dr n)"
       
  3471   unfolding lc by (subst fst_pseudo_divmod_main_as_divide_poly_main) (auto simp: d power_one_over)
       
  3472 
       
  3473 lemma divide_poly_field:
       
  3474   fixes f g :: "'a::field poly"
       
  3475   defines "f' \<equiv> smult ((1 / coeff g (degree g)) ^ (Suc (degree f) - degree g)) f"
       
  3476   shows "f div g = fst (pseudo_divmod f' g)"
       
  3477 proof (cases "g = 0")
       
  3478   case True
       
  3479   show ?thesis
       
  3480     unfolding divide_poly_def pseudo_divmod_def Let_def f'_def True
       
  3481     by (simp add: divide_poly_main_0)
       
  3482 next
       
  3483   case False
       
  3484   from leading_coeff_neq_0[OF False] have "degree f' = degree f"
       
  3485     by (auto simp: f'_def)
       
  3486   then show ?thesis
       
  3487     using length_coeffs_degree[of f'] length_coeffs_degree[of f]
       
  3488     unfolding divide_poly_def pseudo_divmod_def Let_def
       
  3489       divide_poly_main_field[OF False]
       
  3490       length_coeffs_degree[OF False]
       
  3491       f'_def
       
  3492     by force
       
  3493 qed
       
  3494 
       
  3495 instantiation poly :: ("{semidom_divide_unit_factor,idom_divide}") normalization_semidom
       
  3496 begin
       
  3497 
       
  3498 definition unit_factor_poly :: "'a poly \<Rightarrow> 'a poly"
       
  3499   where "unit_factor_poly p = [:unit_factor (lead_coeff p):]"
       
  3500 
       
  3501 definition normalize_poly :: "'a poly \<Rightarrow> 'a poly"
       
  3502   where "normalize p = p div [:unit_factor (lead_coeff p):]"
       
  3503 
       
  3504 instance
       
  3505 proof
       
  3506   fix p :: "'a poly"
       
  3507   show "unit_factor p * normalize p = p"
       
  3508   proof (cases "p = 0")
       
  3509     case True
       
  3510     then show ?thesis
       
  3511       by (simp add: unit_factor_poly_def normalize_poly_def)
       
  3512   next
       
  3513     case False
       
  3514     then have "lead_coeff p \<noteq> 0"
       
  3515       by simp
       
  3516     then have *: "unit_factor (lead_coeff p) \<noteq> 0"
       
  3517       using unit_factor_is_unit [of "lead_coeff p"] by auto
       
  3518     then have "unit_factor (lead_coeff p) dvd 1"
       
  3519       by (auto intro: unit_factor_is_unit)
       
  3520     then have **: "unit_factor (lead_coeff p) dvd c" for c
       
  3521       by (rule dvd_trans) simp
       
  3522     have ***: "unit_factor (lead_coeff p) * (c div unit_factor (lead_coeff p)) = c" for c
       
  3523     proof -
       
  3524       from ** obtain b where "c = unit_factor (lead_coeff p) * b" ..
       
  3525       with False * show ?thesis by simp
       
  3526     qed
       
  3527     have "p div [:unit_factor (lead_coeff p):] =
       
  3528       map_poly (\<lambda>c. c div unit_factor (lead_coeff p)) p"
       
  3529       by (simp add: const_poly_dvd_iff div_const_poly_conv_map_poly **)
       
  3530     then show ?thesis
       
  3531       by (simp add: normalize_poly_def unit_factor_poly_def
       
  3532         smult_conv_map_poly map_poly_map_poly o_def ***)
       
  3533   qed
       
  3534 next
       
  3535   fix p :: "'a poly"
       
  3536   assume "is_unit p"
       
  3537   then obtain c where p: "p = [:c:]" "c dvd 1"
       
  3538     by (auto simp: is_unit_poly_iff)
       
  3539   then show "unit_factor p = p"
       
  3540     by (simp add: unit_factor_poly_def monom_0 is_unit_unit_factor)
       
  3541 next
       
  3542   fix p :: "'a poly"
       
  3543   assume "p \<noteq> 0"
       
  3544   then show "is_unit (unit_factor p)"
       
  3545     by (simp add: unit_factor_poly_def monom_0 is_unit_poly_iff unit_factor_is_unit)
       
  3546 qed (simp_all add: normalize_poly_def unit_factor_poly_def monom_0 lead_coeff_mult unit_factor_mult)
       
  3547 
       
  3548 end
       
  3549 
       
  3550 lemma normalize_poly_eq_map_poly: "normalize p = map_poly (\<lambda>x. x div unit_factor (lead_coeff p)) p"
       
  3551 proof -
       
  3552   have "[:unit_factor (lead_coeff p):] dvd p"
       
  3553     by (metis unit_factor_poly_def unit_factor_self)
       
  3554   then show ?thesis
       
  3555     by (simp add: normalize_poly_def div_const_poly_conv_map_poly)
       
  3556 qed
       
  3557 
       
  3558 lemma coeff_normalize [simp]:
       
  3559   "coeff (normalize p) n = coeff p n div unit_factor (lead_coeff p)"
       
  3560   by (simp add: normalize_poly_eq_map_poly coeff_map_poly)
       
  3561 
       
  3562 class field_unit_factor = field + unit_factor +
       
  3563   assumes unit_factor_field [simp]: "unit_factor = id"
       
  3564 begin
       
  3565 
       
  3566 subclass semidom_divide_unit_factor
       
  3567 proof
       
  3568   fix a
       
  3569   assume "a \<noteq> 0"
       
  3570   then have "1 = a * inverse a" by simp
       
  3571   then have "a dvd 1" ..
       
  3572   then show "unit_factor a dvd 1" by simp
       
  3573 qed simp_all
       
  3574 
       
  3575 end
       
  3576 
       
  3577 lemma unit_factor_pCons:
       
  3578   "unit_factor (pCons a p) = (if p = 0 then [:unit_factor a:] else unit_factor p)"
       
  3579   by (simp add: unit_factor_poly_def)
       
  3580 
       
  3581 lemma normalize_monom [simp]: "normalize (monom a n) = monom (normalize a) n"
       
  3582   by (cases "a = 0") (simp_all add: map_poly_monom normalize_poly_eq_map_poly degree_monom_eq)
       
  3583 
       
  3584 lemma unit_factor_monom [simp]: "unit_factor (monom a n) = [:unit_factor a:]"
       
  3585   by (cases "a = 0") (simp_all add: unit_factor_poly_def degree_monom_eq)
       
  3586 
       
  3587 lemma normalize_const_poly: "normalize [:c:] = [:normalize c:]"
       
  3588   by (simp add: normalize_poly_eq_map_poly map_poly_pCons)
       
  3589 
       
  3590 lemma normalize_smult: "normalize (smult c p) = smult (normalize c) (normalize p)"
       
  3591 proof -
       
  3592   have "smult c p = [:c:] * p" by simp
       
  3593   also have "normalize \<dots> = smult (normalize c) (normalize p)"
       
  3594     by (subst normalize_mult) (simp add: normalize_const_poly)
       
  3595   finally show ?thesis .
       
  3596 qed
       
  3597 
       
  3598 lemma smult_content_normalize_primitive_part [simp]:
       
  3599   "smult (content p) (normalize (primitive_part p)) = normalize p"
       
  3600 proof -
       
  3601   have "smult (content p) (normalize (primitive_part p)) =
       
  3602       normalize ([:content p:] * primitive_part p)"
       
  3603     by (subst normalize_mult) (simp_all add: normalize_const_poly)
       
  3604   also have "[:content p:] * primitive_part p = p" by simp
       
  3605   finally show ?thesis .
       
  3606 qed
       
  3607 
       
  3608 inductive eucl_rel_poly :: "'a::field poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly \<times> 'a poly \<Rightarrow> bool"
       
  3609   where
       
  3610     eucl_rel_poly_by0: "eucl_rel_poly x 0 (0, x)"
       
  3611   | eucl_rel_poly_dividesI: "y \<noteq> 0 \<Longrightarrow> x = q * y \<Longrightarrow> eucl_rel_poly x y (q, 0)"
       
  3612   | eucl_rel_poly_remainderI:
       
  3613       "y \<noteq> 0 \<Longrightarrow> degree r < degree y \<Longrightarrow> x = q * y + r \<Longrightarrow> eucl_rel_poly x y (q, r)"
       
  3614 
       
  3615 lemma eucl_rel_poly_iff:
       
  3616   "eucl_rel_poly x y (q, r) \<longleftrightarrow>
       
  3617     x = q * y + r \<and> (if y = 0 then q = 0 else r = 0 \<or> degree r < degree y)"
       
  3618   by (auto elim: eucl_rel_poly.cases
       
  3619       intro: eucl_rel_poly_by0 eucl_rel_poly_dividesI eucl_rel_poly_remainderI)
       
  3620 
       
  3621 lemma eucl_rel_poly_0: "eucl_rel_poly 0 y (0, 0)"
       
  3622   by (simp add: eucl_rel_poly_iff)
       
  3623 
       
  3624 lemma eucl_rel_poly_by_0: "eucl_rel_poly x 0 (0, x)"
       
  3625   by (simp add: eucl_rel_poly_iff)
       
  3626 
       
  3627 lemma eucl_rel_poly_pCons:
       
  3628   assumes rel: "eucl_rel_poly x y (q, r)"
       
  3629   assumes y: "y \<noteq> 0"
       
  3630   assumes b: "b = coeff (pCons a r) (degree y) / coeff y (degree y)"
       
  3631   shows "eucl_rel_poly (pCons a x) y (pCons b q, pCons a r - smult b y)"
       
  3632     (is "eucl_rel_poly ?x y (?q, ?r)")
       
  3633 proof -
       
  3634   from assms have x: "x = q * y + r" and r: "r = 0 \<or> degree r < degree y"
       
  3635     by (simp_all add: eucl_rel_poly_iff)
       
  3636   from b x have "?x = ?q * y + ?r" by simp
       
  3637   moreover
       
  3638   have "?r = 0 \<or> degree ?r < degree y"
       
  3639   proof (rule eq_zero_or_degree_less)
       
  3640     show "degree ?r \<le> degree y"
       
  3641     proof (rule degree_diff_le)
       
  3642       from r show "degree (pCons a r) \<le> degree y"
       
  3643         by auto
       
  3644       show "degree (smult b y) \<le> degree y"
       
  3645         by (rule degree_smult_le)
       
  3646     qed
       
  3647     from \<open>y \<noteq> 0\<close> show "coeff ?r (degree y) = 0"
       
  3648       by (simp add: b)
       
  3649   qed
       
  3650   ultimately show ?thesis
       
  3651     unfolding eucl_rel_poly_iff using \<open>y \<noteq> 0\<close> by simp
       
  3652 qed
       
  3653 
       
  3654 lemma eucl_rel_poly_exists: "\<exists>q r. eucl_rel_poly x y (q, r)"
       
  3655   apply (cases "y = 0")
       
  3656    apply (fast intro!: eucl_rel_poly_by_0)
       
  3657   apply (induct x)
       
  3658    apply (fast intro!: eucl_rel_poly_0)
       
  3659   apply (fast intro!: eucl_rel_poly_pCons)
       
  3660   done
       
  3661 
       
  3662 lemma eucl_rel_poly_unique:
       
  3663   assumes 1: "eucl_rel_poly x y (q1, r1)"
       
  3664   assumes 2: "eucl_rel_poly x y (q2, r2)"
       
  3665   shows "q1 = q2 \<and> r1 = r2"
       
  3666 proof (cases "y = 0")
       
  3667   assume "y = 0"
       
  3668   with assms show ?thesis
       
  3669     by (simp add: eucl_rel_poly_iff)
       
  3670 next
       
  3671   assume [simp]: "y \<noteq> 0"
       
  3672   from 1 have q1: "x = q1 * y + r1" and r1: "r1 = 0 \<or> degree r1 < degree y"
       
  3673     unfolding eucl_rel_poly_iff by simp_all
       
  3674   from 2 have q2: "x = q2 * y + r2" and r2: "r2 = 0 \<or> degree r2 < degree y"
       
  3675     unfolding eucl_rel_poly_iff by simp_all
       
  3676   from q1 q2 have q3: "(q1 - q2) * y = r2 - r1"
       
  3677     by (simp add: algebra_simps)
       
  3678   from r1 r2 have r3: "(r2 - r1) = 0 \<or> degree (r2 - r1) < degree y"
       
  3679     by (auto intro: degree_diff_less)
       
  3680   show "q1 = q2 \<and> r1 = r2"
       
  3681   proof (rule classical)
       
  3682     assume "\<not> ?thesis"
       
  3683     with q3 have "q1 \<noteq> q2" and "r1 \<noteq> r2" by auto
       
  3684     with r3 have "degree (r2 - r1) < degree y" by simp
       
  3685     also have "degree y \<le> degree (q1 - q2) + degree y" by simp
       
  3686     also from \<open>q1 \<noteq> q2\<close> have "\<dots> = degree ((q1 - q2) * y)"
       
  3687       by (simp add: degree_mult_eq)
       
  3688     also from q3 have "\<dots> = degree (r2 - r1)"
       
  3689       by simp
       
  3690     finally have "degree (r2 - r1) < degree (r2 - r1)" .
       
  3691     then show ?thesis by simp
       
  3692   qed
       
  3693 qed
       
  3694 
       
  3695 lemma eucl_rel_poly_0_iff: "eucl_rel_poly 0 y (q, r) \<longleftrightarrow> q = 0 \<and> r = 0"
       
  3696   by (auto dest: eucl_rel_poly_unique intro: eucl_rel_poly_0)
       
  3697 
       
  3698 lemma eucl_rel_poly_by_0_iff: "eucl_rel_poly x 0 (q, r) \<longleftrightarrow> q = 0 \<and> r = x"
       
  3699   by (auto dest: eucl_rel_poly_unique intro: eucl_rel_poly_by_0)
       
  3700 
       
  3701 lemmas eucl_rel_poly_unique_div = eucl_rel_poly_unique [THEN conjunct1]
       
  3702 
       
  3703 lemmas eucl_rel_poly_unique_mod = eucl_rel_poly_unique [THEN conjunct2]
       
  3704 
       
  3705 instantiation poly :: (field) semidom_modulo
       
  3706 begin
       
  3707 
       
  3708 definition modulo_poly :: "'a poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly"
       
  3709   where mod_poly_def: "f mod g =
       
  3710     (if g = 0 then f else pseudo_mod (smult ((1 / lead_coeff g) ^ (Suc (degree f) - degree g)) f) g)"
       
  3711 
       
  3712 instance
       
  3713 proof
       
  3714   fix x y :: "'a poly"
       
  3715   show "x div y * y + x mod y = x"
       
  3716   proof (cases "y = 0")
       
  3717     case True
       
  3718     then show ?thesis
       
  3719       by (simp add: divide_poly_0 mod_poly_def)
       
  3720   next
       
  3721     case False
       
  3722     then have "pseudo_divmod (smult ((1 / lead_coeff y) ^ (Suc (degree x) - degree y)) x) y =
       
  3723         (x div y, x mod y)"
       
  3724       by (simp add: divide_poly_field mod_poly_def pseudo_mod_def)
       
  3725     with False pseudo_divmod [OF False this] show ?thesis
       
  3726       by (simp add: power_mult_distrib [symmetric] ac_simps)
       
  3727   qed
       
  3728 qed
       
  3729 
       
  3730 end
       
  3731 
       
  3732 lemma eucl_rel_poly: "eucl_rel_poly x y (x div y, x mod y)"
       
  3733   unfolding eucl_rel_poly_iff
       
  3734 proof
       
  3735   show "x = x div y * y + x mod y"
       
  3736     by (simp add: div_mult_mod_eq)
       
  3737   show "if y = 0 then x div y = 0 else x mod y = 0 \<or> degree (x mod y) < degree y"
       
  3738   proof (cases "y = 0")
       
  3739     case True
       
  3740     then show ?thesis by auto
       
  3741   next
       
  3742     case False
       
  3743     with pseudo_mod[OF this] show ?thesis
       
  3744       by (simp add: mod_poly_def)
       
  3745   qed
       
  3746 qed
       
  3747 
       
  3748 lemma div_poly_eq: "eucl_rel_poly x y (q, r) \<Longrightarrow> x div y = q"
       
  3749   for x :: "'a::field poly"
       
  3750   by (rule eucl_rel_poly_unique_div [OF eucl_rel_poly])
       
  3751 
       
  3752 lemma mod_poly_eq: "eucl_rel_poly x y (q, r) \<Longrightarrow> x mod y = r"
       
  3753   for x :: "'a::field poly"
       
  3754   by (rule eucl_rel_poly_unique_mod [OF eucl_rel_poly])
       
  3755 
       
  3756 instance poly :: (field) ring_div
       
  3757 proof
       
  3758   fix x y z :: "'a poly"
       
  3759   assume "y \<noteq> 0"
       
  3760   with eucl_rel_poly [of x y] have "eucl_rel_poly (x + z * y) y (z + x div y, x mod y)"
       
  3761     by (simp add: eucl_rel_poly_iff distrib_right)
       
  3762   then show "(x + z * y) div y = z + x div y"
       
  3763     by (rule div_poly_eq)
       
  3764 next
       
  3765   fix x y z :: "'a poly"
       
  3766   assume "x \<noteq> 0"
       
  3767   show "(x * y) div (x * z) = y div z"
       
  3768   proof (cases "y \<noteq> 0 \<and> z \<noteq> 0")
       
  3769     have "\<And>x::'a poly. eucl_rel_poly x 0 (0, x)"
       
  3770       by (rule eucl_rel_poly_by_0)
       
  3771     then have [simp]: "\<And>x::'a poly. x div 0 = 0"
       
  3772       by (rule div_poly_eq)
       
  3773     have "\<And>x::'a poly. eucl_rel_poly 0 x (0, 0)"
       
  3774       by (rule eucl_rel_poly_0)
       
  3775     then have [simp]: "\<And>x::'a poly. 0 div x = 0"
       
  3776       by (rule div_poly_eq)
       
  3777     case False
       
  3778     then show ?thesis by auto
       
  3779   next
       
  3780     case True
       
  3781     then have "y \<noteq> 0" and "z \<noteq> 0" by auto
       
  3782     with \<open>x \<noteq> 0\<close> have "\<And>q r. eucl_rel_poly y z (q, r) \<Longrightarrow> eucl_rel_poly (x * y) (x * z) (q, x * r)"
       
  3783       by (auto simp: eucl_rel_poly_iff algebra_simps) (rule classical, simp add: degree_mult_eq)
       
  3784     moreover from eucl_rel_poly have "eucl_rel_poly y z (y div z, y mod z)" .
       
  3785     ultimately have "eucl_rel_poly (x * y) (x * z) (y div z, x * (y mod z))" .
       
  3786     then show ?thesis
       
  3787       by (simp add: div_poly_eq)
       
  3788   qed
       
  3789 qed
       
  3790 
       
  3791 lemma div_pCons_eq:
       
  3792   "pCons a p div q =
       
  3793     (if q = 0 then 0
       
  3794      else pCons (coeff (pCons a (p mod q)) (degree q) / lead_coeff q) (p div q))"
       
  3795   using eucl_rel_poly_pCons [OF eucl_rel_poly _ refl, of q a p]
       
  3796   by (auto intro: div_poly_eq)
       
  3797 
       
  3798 lemma mod_pCons_eq:
       
  3799   "pCons a p mod q =
       
  3800     (if q = 0 then pCons a p
       
  3801      else pCons a (p mod q) - smult (coeff (pCons a (p mod q)) (degree q) / lead_coeff q) q)"
       
  3802   using eucl_rel_poly_pCons [OF eucl_rel_poly _ refl, of q a p]
       
  3803   by (auto intro: mod_poly_eq)
       
  3804 
       
  3805 lemma div_mod_fold_coeffs:
       
  3806   "(p div q, p mod q) =
       
  3807     (if q = 0 then (0, p)
       
  3808      else
       
  3809       fold_coeffs
       
  3810         (\<lambda>a (s, r).
       
  3811           let b = coeff (pCons a r) (degree q) / coeff q (degree q)
       
  3812           in (pCons b s, pCons a r - smult b q)) p (0, 0))"
       
  3813   by (rule sym, induct p) (auto simp: div_pCons_eq mod_pCons_eq Let_def)
       
  3814 
       
  3815 lemma degree_mod_less: "y \<noteq> 0 \<Longrightarrow> x mod y = 0 \<or> degree (x mod y) < degree y"
       
  3816   using eucl_rel_poly [of x y] unfolding eucl_rel_poly_iff by simp
       
  3817 
       
  3818 lemma degree_mod_less': "b \<noteq> 0 \<Longrightarrow> a mod b \<noteq> 0 \<Longrightarrow> degree (a mod b) < degree b"
       
  3819   using degree_mod_less[of b a] by auto
       
  3820 
       
  3821 lemma div_poly_less:
       
  3822   fixes x :: "'a::field poly"
       
  3823   assumes "degree x < degree y"
       
  3824   shows "x div y = 0"
       
  3825 proof -
       
  3826   from assms have "eucl_rel_poly x y (0, x)"
       
  3827     by (simp add: eucl_rel_poly_iff)
       
  3828   then show "x div y = 0"
       
  3829     by (rule div_poly_eq)
       
  3830 qed
       
  3831 
       
  3832 lemma mod_poly_less:
       
  3833   assumes "degree x < degree y"
       
  3834   shows "x mod y = x"
       
  3835 proof -
       
  3836   from assms have "eucl_rel_poly x y (0, x)"
       
  3837     by (simp add: eucl_rel_poly_iff)
       
  3838   then show "x mod y = x"
       
  3839     by (rule mod_poly_eq)
       
  3840 qed
       
  3841 
       
  3842 lemma eucl_rel_poly_smult_left:
       
  3843   "eucl_rel_poly x y (q, r) \<Longrightarrow> eucl_rel_poly (smult a x) y (smult a q, smult a r)"
       
  3844   by (simp add: eucl_rel_poly_iff smult_add_right)
       
  3845 
       
  3846 lemma div_smult_left: "(smult a x) div y = smult a (x div y)"
       
  3847   for x y :: "'a::field poly"
       
  3848   by (rule div_poly_eq, rule eucl_rel_poly_smult_left, rule eucl_rel_poly)
       
  3849 
       
  3850 lemma mod_smult_left: "(smult a x) mod y = smult a (x mod y)"
       
  3851   by (rule mod_poly_eq, rule eucl_rel_poly_smult_left, rule eucl_rel_poly)
       
  3852 
       
  3853 lemma poly_div_minus_left [simp]: "(- x) div y = - (x div y)"
       
  3854   for x y :: "'a::field poly"
       
  3855   using div_smult_left [of "- 1::'a"] by simp
       
  3856 
       
  3857 lemma poly_mod_minus_left [simp]: "(- x) mod y = - (x mod y)"
       
  3858   for x y :: "'a::field poly"
       
  3859   using mod_smult_left [of "- 1::'a"] by simp
       
  3860 
       
  3861 lemma eucl_rel_poly_add_left:
       
  3862   assumes "eucl_rel_poly x y (q, r)"
       
  3863   assumes "eucl_rel_poly x' y (q', r')"
       
  3864   shows "eucl_rel_poly (x + x') y (q + q', r + r')"
       
  3865   using assms unfolding eucl_rel_poly_iff
       
  3866   by (auto simp: algebra_simps degree_add_less)
       
  3867 
       
  3868 lemma poly_div_add_left: "(x + y) div z = x div z + y div z"
       
  3869   for x y z :: "'a::field poly"
       
  3870   using eucl_rel_poly_add_left [OF eucl_rel_poly eucl_rel_poly]
       
  3871   by (rule div_poly_eq)
       
  3872 
       
  3873 lemma poly_mod_add_left: "(x + y) mod z = x mod z + y mod z"
       
  3874   for x y z :: "'a::field poly"
       
  3875   using eucl_rel_poly_add_left [OF eucl_rel_poly eucl_rel_poly]
       
  3876   by (rule mod_poly_eq)
       
  3877 
       
  3878 lemma poly_div_diff_left: "(x - y) div z = x div z - y div z"
       
  3879   for x y z :: "'a::field poly"
       
  3880   by (simp only: diff_conv_add_uminus poly_div_add_left poly_div_minus_left)
       
  3881 
       
  3882 lemma poly_mod_diff_left: "(x - y) mod z = x mod z - y mod z"
       
  3883   for x y z :: "'a::field poly"
       
  3884   by (simp only: diff_conv_add_uminus poly_mod_add_left poly_mod_minus_left)
       
  3885 
       
  3886 lemma eucl_rel_poly_smult_right:
       
  3887   "a \<noteq> 0 \<Longrightarrow> eucl_rel_poly x y (q, r) \<Longrightarrow> eucl_rel_poly x (smult a y) (smult (inverse a) q, r)"
       
  3888   by (simp add: eucl_rel_poly_iff)
       
  3889 
       
  3890 lemma div_smult_right: "a \<noteq> 0 \<Longrightarrow> x div (smult a y) = smult (inverse a) (x div y)"
       
  3891   for x y :: "'a::field poly"
       
  3892   by (rule div_poly_eq, erule eucl_rel_poly_smult_right, rule eucl_rel_poly)
       
  3893 
       
  3894 lemma mod_smult_right: "a \<noteq> 0 \<Longrightarrow> x mod (smult a y) = x mod y"
       
  3895   by (rule mod_poly_eq, erule eucl_rel_poly_smult_right, rule eucl_rel_poly)
       
  3896 
       
  3897 lemma poly_div_minus_right [simp]: "x div (- y) = - (x div y)"
       
  3898   for x y :: "'a::field poly"
       
  3899   using div_smult_right [of "- 1::'a"] by (simp add: nonzero_inverse_minus_eq)
       
  3900 
       
  3901 lemma poly_mod_minus_right [simp]: "x mod (- y) = x mod y"
       
  3902   for x y :: "'a::field poly"
       
  3903   using mod_smult_right [of "- 1::'a"] by simp
       
  3904 
       
  3905 lemma eucl_rel_poly_mult:
       
  3906   "eucl_rel_poly x y (q, r) \<Longrightarrow> eucl_rel_poly q z (q', r') \<Longrightarrow>
       
  3907     eucl_rel_poly x (y * z) (q', y * r' + r)"
       
  3908   apply (cases "z = 0", simp add: eucl_rel_poly_iff)
       
  3909   apply (cases "y = 0", simp add: eucl_rel_poly_by_0_iff eucl_rel_poly_0_iff)
       
  3910   apply (cases "r = 0")
       
  3911    apply (cases "r' = 0")
       
  3912     apply (simp add: eucl_rel_poly_iff)
       
  3913    apply (simp add: eucl_rel_poly_iff field_simps degree_mult_eq)
       
  3914   apply (cases "r' = 0")
       
  3915    apply (simp add: eucl_rel_poly_iff degree_mult_eq)
       
  3916   apply (simp add: eucl_rel_poly_iff field_simps)
       
  3917   apply (simp add: degree_mult_eq degree_add_less)
       
  3918   done
       
  3919 
       
  3920 lemma poly_div_mult_right: "x div (y * z) = (x div y) div z"
       
  3921   for x y z :: "'a::field poly"
       
  3922   by (rule div_poly_eq, rule eucl_rel_poly_mult, (rule eucl_rel_poly)+)
       
  3923 
       
  3924 lemma poly_mod_mult_right: "x mod (y * z) = y * (x div y mod z) + x mod y"
       
  3925   for x y z :: "'a::field poly"
       
  3926   by (rule mod_poly_eq, rule eucl_rel_poly_mult, (rule eucl_rel_poly)+)
       
  3927 
       
  3928 lemma mod_pCons:
       
  3929   fixes a :: "'a::field"
       
  3930     and x y :: "'a::field poly"
       
  3931   assumes y: "y \<noteq> 0"
       
  3932   defines "b \<equiv> coeff (pCons a (x mod y)) (degree y) / coeff y (degree y)"
       
  3933   shows "(pCons a x) mod y = pCons a (x mod y) - smult b y"
       
  3934   unfolding b_def
       
  3935   by (rule mod_poly_eq, rule eucl_rel_poly_pCons [OF eucl_rel_poly y refl])
       
  3936 
       
  3937 
       
  3938 subsubsection \<open>List-based versions for fast implementation\<close>
       
  3939 (* Subsection by:
       
  3940       Sebastiaan Joosten
       
  3941       René Thiemann
       
  3942       Akihisa Yamada
       
  3943     *)
       
  3944 fun minus_poly_rev_list :: "'a :: group_add list \<Rightarrow> 'a list \<Rightarrow> 'a list"
       
  3945   where
       
  3946     "minus_poly_rev_list (x # xs) (y # ys) = (x - y) # (minus_poly_rev_list xs ys)"
       
  3947   | "minus_poly_rev_list xs [] = xs"
       
  3948   | "minus_poly_rev_list [] (y # ys) = []"
       
  3949 
       
  3950 fun pseudo_divmod_main_list ::
       
  3951   "'a::comm_ring_1 \<Rightarrow> 'a list \<Rightarrow> 'a list \<Rightarrow> 'a list \<Rightarrow> nat \<Rightarrow> 'a list \<times> 'a list"
       
  3952   where
       
  3953     "pseudo_divmod_main_list lc q r d (Suc n) =
       
  3954       (let
       
  3955         rr = map (op * lc) r;
       
  3956         a = hd r;
       
  3957         qqq = cCons a (map (op * lc) q);
       
  3958         rrr = tl (if a = 0 then rr else minus_poly_rev_list rr (map (op * a) d))
       
  3959        in pseudo_divmod_main_list lc qqq rrr d n)"
       
  3960   | "pseudo_divmod_main_list lc q r d 0 = (q, r)"
       
  3961 
       
  3962 fun pseudo_mod_main_list :: "'a::comm_ring_1 \<Rightarrow> 'a list \<Rightarrow> 'a list \<Rightarrow> nat \<Rightarrow> 'a list"
       
  3963   where
       
  3964     "pseudo_mod_main_list lc r d (Suc n) =
       
  3965       (let
       
  3966         rr = map (op * lc) r;
       
  3967         a = hd r;
       
  3968         rrr = tl (if a = 0 then rr else minus_poly_rev_list rr (map (op * a) d))
       
  3969        in pseudo_mod_main_list lc rrr d n)"
       
  3970   | "pseudo_mod_main_list lc r d 0 = r"
       
  3971 
       
  3972 
       
  3973 fun divmod_poly_one_main_list ::
       
  3974     "'a::comm_ring_1 list \<Rightarrow> 'a list \<Rightarrow> 'a list \<Rightarrow> nat \<Rightarrow> 'a list \<times> 'a list"
       
  3975   where
       
  3976     "divmod_poly_one_main_list q r d (Suc n) =
       
  3977       (let
       
  3978         a = hd r;
       
  3979         qqq = cCons a q;
       
  3980         rr = tl (if a = 0 then r else minus_poly_rev_list r (map (op * a) d))
       
  3981        in divmod_poly_one_main_list qqq rr d n)"
       
  3982   | "divmod_poly_one_main_list q r d 0 = (q, r)"
       
  3983 
       
  3984 fun mod_poly_one_main_list :: "'a::comm_ring_1 list \<Rightarrow> 'a list \<Rightarrow> nat \<Rightarrow> 'a list"
       
  3985   where
       
  3986     "mod_poly_one_main_list r d (Suc n) =
       
  3987       (let
       
  3988         a = hd r;
       
  3989         rr = tl (if a = 0 then r else minus_poly_rev_list r (map (op * a) d))
       
  3990        in mod_poly_one_main_list rr d n)"
       
  3991   | "mod_poly_one_main_list r d 0 = r"
       
  3992 
       
  3993 definition pseudo_divmod_list :: "'a::comm_ring_1 list \<Rightarrow> 'a list \<Rightarrow> 'a list \<times> 'a list"
       
  3994   where "pseudo_divmod_list p q =
       
  3995     (if q = [] then ([], p)
       
  3996      else
       
  3997       (let rq = rev q;
       
  3998         (qu,re) = pseudo_divmod_main_list (hd rq) [] (rev p) rq (1 + length p - length q)
       
  3999        in (qu, rev re)))"
       
  4000 
       
  4001 definition pseudo_mod_list :: "'a::comm_ring_1 list \<Rightarrow> 'a list \<Rightarrow> 'a list"
       
  4002   where "pseudo_mod_list p q =
       
  4003     (if q = [] then p
       
  4004      else
       
  4005       (let
       
  4006         rq = rev q;
       
  4007         re = pseudo_mod_main_list (hd rq) (rev p) rq (1 + length p - length q)
       
  4008        in rev re))"
       
  4009 
       
  4010 lemma minus_zero_does_nothing: "minus_poly_rev_list x (map (op * 0) y) = x"
       
  4011   for x :: "'a::ring list"
       
  4012   by (induct x y rule: minus_poly_rev_list.induct) auto
       
  4013 
       
  4014 lemma length_minus_poly_rev_list [simp]: "length (minus_poly_rev_list xs ys) = length xs"
       
  4015   by (induct xs ys rule: minus_poly_rev_list.induct) auto
       
  4016 
       
  4017 lemma if_0_minus_poly_rev_list:
       
  4018   "(if a = 0 then x else minus_poly_rev_list x (map (op * a) y)) =
       
  4019     minus_poly_rev_list x (map (op * a) y)"
       
  4020   for a :: "'a::ring"
       
  4021   by(cases "a = 0") (simp_all add: minus_zero_does_nothing)
       
  4022 
       
  4023 lemma Poly_append: "Poly (a @ b) = Poly a + monom 1 (length a) * Poly b"
       
  4024   for a :: "'a::comm_semiring_1 list"
       
  4025   by (induct a) (auto simp: monom_0 monom_Suc)
       
  4026 
       
  4027 lemma minus_poly_rev_list: "length p \<ge> length q \<Longrightarrow>
       
  4028   Poly (rev (minus_poly_rev_list (rev p) (rev q))) =
       
  4029     Poly p - monom 1 (length p - length q) * Poly q"
       
  4030   for p q :: "'a :: comm_ring_1 list"
       
  4031 proof (induct "rev p" "rev q" arbitrary: p q rule: minus_poly_rev_list.induct)
       
  4032   case (1 x xs y ys)
       
  4033   then have "length (rev q) \<le> length (rev p)"
       
  4034     by simp
       
  4035   from this[folded 1(2,3)] have ys_xs: "length ys \<le> length xs"
       
  4036     by simp
       
  4037   then have *: "Poly (rev (minus_poly_rev_list xs ys)) =
       
  4038       Poly (rev xs) - monom 1 (length xs - length ys) * Poly (rev ys)"
       
  4039     by (subst "1.hyps"(1)[of "rev xs" "rev ys", unfolded rev_rev_ident length_rev]) auto
       
  4040   have "Poly p - monom 1 (length p - length q) * Poly q =
       
  4041     Poly (rev (rev p)) - monom 1 (length (rev (rev p)) - length (rev (rev q))) * Poly (rev (rev q))"
       
  4042     by simp
       
  4043   also have "\<dots> =
       
  4044       Poly (rev (x # xs)) - monom 1 (length (x # xs) - length (y # ys)) * Poly (rev (y # ys))"
       
  4045     unfolding 1(2,3) by simp
       
  4046   also from ys_xs have "\<dots> =
       
  4047     Poly (rev xs) + monom x (length xs) -
       
  4048       (monom 1 (length xs - length ys) * Poly (rev ys) + monom y (length xs))"
       
  4049     by (simp add: Poly_append distrib_left mult_monom smult_monom)
       
  4050   also have "\<dots> = Poly (rev (minus_poly_rev_list xs ys)) + monom (x - y) (length xs)"
       
  4051     unfolding * diff_monom[symmetric] by simp
       
  4052   finally show ?case
       
  4053     by (simp add: 1(2,3)[symmetric] smult_monom Poly_append)
       
  4054 qed auto
       
  4055 
       
  4056 lemma smult_monom_mult: "smult a (monom b n * f) = monom (a * b) n * f"
       
  4057   using smult_monom [of a _ n] by (metis mult_smult_left)
       
  4058 
       
  4059 lemma head_minus_poly_rev_list:
       
  4060   "length d \<le> length r \<Longrightarrow> d \<noteq> [] \<Longrightarrow>
       
  4061     hd (minus_poly_rev_list (map (op * (last d)) r) (map (op * (hd r)) (rev d))) = 0"
       
  4062   for d r :: "'a::comm_ring list"
       
  4063 proof (induct r)
       
  4064   case Nil
       
  4065   then show ?case by simp
       
  4066 next
       
  4067   case (Cons a rs)
       
  4068   then show ?case by (cases "rev d") (simp_all add: ac_simps)
       
  4069 qed
       
  4070 
       
  4071 lemma Poly_map: "Poly (map (op * a) p) = smult a (Poly p)"
       
  4072 proof (induct p)
       
  4073   case Nil
       
  4074   then show ?case by simp
       
  4075 next
       
  4076   case (Cons x xs)
       
  4077   then show ?case by (cases "Poly xs = 0") auto
       
  4078 qed
       
  4079 
       
  4080 lemma last_coeff_is_hd: "xs \<noteq> [] \<Longrightarrow> coeff (Poly xs) (length xs - 1) = hd (rev xs)"
       
  4081   by (simp_all add: hd_conv_nth rev_nth nth_default_nth nth_append)
       
  4082 
       
  4083 lemma pseudo_divmod_main_list_invar:
       
  4084   assumes leading_nonzero: "last d \<noteq> 0"
       
  4085     and lc: "last d = lc"
       
  4086     and "d \<noteq> []"
       
  4087     and "pseudo_divmod_main_list lc q (rev r) (rev d) n = (q', rev r')"
       
  4088     and "n = 1 + length r - length d"
       
  4089   shows "pseudo_divmod_main lc (monom 1 n * Poly q) (Poly r) (Poly d) (length r - 1) n =
       
  4090     (Poly q', Poly r')"
       
  4091   using assms(4-)
       
  4092 proof (induct n arbitrary: r q)
       
  4093   case (Suc n)
       
  4094   from Suc.prems have *: "\<not> Suc (length r) \<le> length d"
       
  4095     by simp
       
  4096   with \<open>d \<noteq> []\<close> have "r \<noteq> []"
       
  4097     using Suc_leI length_greater_0_conv list.size(3) by fastforce
       
  4098   let ?a = "(hd (rev r))"
       
  4099   let ?rr = "map (op * lc) (rev r)"
       
  4100   let ?rrr = "rev (tl (minus_poly_rev_list ?rr (map (op * ?a) (rev d))))"
       
  4101   let ?qq = "cCons ?a (map (op * lc) q)"
       
  4102   from * Suc(3) have n: "n = (1 + length r - length d - 1)"
       
  4103     by simp
       
  4104   from * have rr_val:"(length ?rrr) = (length r - 1)"
       
  4105     by auto
       
  4106   with \<open>r \<noteq> []\<close> * have rr_smaller: "(1 + length r - length d - 1) = (1 + length ?rrr - length d)"
       
  4107     by auto
       
  4108   from * have id: "Suc (length r) - length d = Suc (length r - length d)"
       
  4109     by auto
       
  4110   from Suc.prems *
       
  4111   have "pseudo_divmod_main_list lc ?qq (rev ?rrr) (rev d) (1 + length r - length d - 1) = (q', rev r')"
       
  4112     by (simp add: Let_def if_0_minus_poly_rev_list id)
       
  4113   with n have v: "pseudo_divmod_main_list lc ?qq (rev ?rrr) (rev d) n = (q', rev r')"
       
  4114     by auto
       
  4115   from * have sucrr:"Suc (length r) - length d = Suc (length r - length d)"
       
  4116     using Suc_diff_le not_less_eq_eq by blast
       
  4117   from Suc(3) \<open>r \<noteq> []\<close> have n_ok : "n = 1 + (length ?rrr) - length d"
       
  4118     by simp
       
  4119   have cong: "\<And>x1 x2 x3 x4 y1 y2 y3 y4. x1 = y1 \<Longrightarrow> x2 = y2 \<Longrightarrow> x3 = y3 \<Longrightarrow> x4 = y4 \<Longrightarrow>
       
  4120       pseudo_divmod_main lc x1 x2 x3 x4 n = pseudo_divmod_main lc y1 y2 y3 y4 n"
       
  4121     by simp
       
  4122   have hd_rev: "coeff (Poly r) (length r - Suc 0) = hd (rev r)"
       
  4123     using last_coeff_is_hd[OF \<open>r \<noteq> []\<close>] by simp
       
  4124   show ?case
       
  4125     unfolding Suc.hyps(1)[OF v n_ok, symmetric] pseudo_divmod_main.simps Let_def
       
  4126   proof (rule cong[OF _ _ refl], goal_cases)
       
  4127     case 1
       
  4128     show ?case
       
  4129       by (simp add: monom_Suc hd_rev[symmetric] smult_monom Poly_map)
       
  4130   next
       
  4131     case 2
       
  4132     show ?case
       
  4133     proof (subst Poly_on_rev_starting_with_0, goal_cases)
       
  4134       show "hd (minus_poly_rev_list (map (op * lc) (rev r)) (map (op * (hd (rev r))) (rev d))) = 0"
       
  4135         by (fold lc, subst head_minus_poly_rev_list, insert * \<open>d \<noteq> []\<close>, auto)
       
  4136       from * have "length d \<le> length r"
       
  4137         by simp
       
  4138       then show "smult lc (Poly r) - monom (coeff (Poly r) (length r - 1)) n * Poly d =
       
  4139           Poly (rev (minus_poly_rev_list (map (op * lc) (rev r)) (map (op * (hd (rev r))) (rev d))))"
       
  4140         by (fold rev_map) (auto simp add: n smult_monom_mult Poly_map hd_rev [symmetric]
       
  4141             minus_poly_rev_list)
       
  4142     qed
       
  4143   qed simp
       
  4144 qed simp
       
  4145 
       
  4146 lemma pseudo_divmod_impl [code]:
       
  4147   "pseudo_divmod f g = map_prod poly_of_list poly_of_list (pseudo_divmod_list (coeffs f) (coeffs g))"
       
  4148     for f g :: "'a::comm_ring_1 poly"
       
  4149 proof (cases "g = 0")
       
  4150   case False
       
  4151   then have "last (coeffs g) \<noteq> 0"
       
  4152     and "last (coeffs g) = lead_coeff g"
       
  4153     and "coeffs g \<noteq> []"
       
  4154     by (simp_all add: last_coeffs_eq_coeff_degree)
       
  4155   moreover obtain q r where qr: "pseudo_divmod_main_list
       
  4156     (last (coeffs g)) (rev [])
       
  4157     (rev (coeffs f)) (rev (coeffs g))
       
  4158     (1 + length (coeffs f) -
       
  4159     length (coeffs g)) = (q, rev (rev r))"
       
  4160     by force
       
  4161   ultimately have "(Poly q, Poly (rev r)) = pseudo_divmod_main (lead_coeff g) 0 f g
       
  4162     (length (coeffs f) - Suc 0) (Suc (length (coeffs f)) - length (coeffs g))"
       
  4163     by (subst pseudo_divmod_main_list_invar [symmetric]) auto
       
  4164   moreover have "pseudo_divmod_main_list
       
  4165     (hd (rev (coeffs g))) []
       
  4166     (rev (coeffs f)) (rev (coeffs g))
       
  4167     (1 + length (coeffs f) -
       
  4168     length (coeffs g)) = (q, r)"
       
  4169     using qr hd_rev [OF \<open>coeffs g \<noteq> []\<close>] by simp
       
  4170   ultimately show ?thesis
       
  4171     by (auto simp: degree_eq_length_coeffs pseudo_divmod_def pseudo_divmod_list_def Let_def)
       
  4172 next
       
  4173   case True
       
  4174   then show ?thesis
       
  4175     by (auto simp add: pseudo_divmod_def pseudo_divmod_list_def)
       
  4176 qed
       
  4177 
       
  4178 lemma pseudo_mod_main_list:
       
  4179   "snd (pseudo_divmod_main_list l q xs ys n) = pseudo_mod_main_list l xs ys n"
       
  4180   by (induct n arbitrary: l q xs ys) (auto simp: Let_def)
       
  4181 
       
  4182 lemma pseudo_mod_impl[code]: "pseudo_mod f g = poly_of_list (pseudo_mod_list (coeffs f) (coeffs g))"
       
  4183 proof -
       
  4184   have snd_case: "\<And>f g p. snd ((\<lambda>(x,y). (f x, g y)) p) = g (snd p)"
       
  4185     by auto
       
  4186   show ?thesis
       
  4187     unfolding pseudo_mod_def pseudo_divmod_impl pseudo_divmod_list_def
       
  4188       pseudo_mod_list_def Let_def
       
  4189     by (simp add: snd_case pseudo_mod_main_list)
       
  4190 qed
       
  4191 
       
  4192 
       
  4193 subsubsection \<open>Improved Code-Equations for Polynomial (Pseudo) Division\<close>
       
  4194 
       
  4195 lemma pdivmod_pdivmodrel: "eucl_rel_poly p q (r, s) \<longleftrightarrow> (p div q, p mod q) = (r, s)"
       
  4196   by (metis eucl_rel_poly eucl_rel_poly_unique)
       
  4197 
       
  4198 lemma pdivmod_via_pseudo_divmod:
       
  4199   "(f div g, f mod g) =
       
  4200     (if g = 0 then (0, f)
       
  4201      else
       
  4202       let
       
  4203         ilc = inverse (coeff g (degree g));
       
  4204         h = smult ilc g;
       
  4205         (q,r) = pseudo_divmod f h
       
  4206       in (smult ilc q, r))"
       
  4207   (is "?l = ?r")
       
  4208 proof (cases "g = 0")
       
  4209   case True
       
  4210   then show ?thesis by simp
       
  4211 next
       
  4212   case False
       
  4213   define lc where "lc = inverse (coeff g (degree g))"
       
  4214   define h where "h = smult lc g"
       
  4215   from False have h1: "coeff h (degree h) = 1" and lc: "lc \<noteq> 0"
       
  4216     by (auto simp: h_def lc_def)
       
  4217   then have h0: "h \<noteq> 0"
       
  4218     by auto
       
  4219   obtain q r where p: "pseudo_divmod f h = (q, r)"
       
  4220     by force
       
  4221   from False have id: "?r = (smult lc q, r)"
       
  4222     by (auto simp: Let_def h_def[symmetric] lc_def[symmetric] p)
       
  4223   from pseudo_divmod[OF h0 p, unfolded h1]
       
  4224   have f: "f = h * q + r" and r: "r = 0 \<or> degree r < degree h"
       
  4225     by auto
       
  4226   from f r h0 have "eucl_rel_poly f h (q, r)"
       
  4227     by (auto simp: eucl_rel_poly_iff)
       
  4228   then have "(f div h, f mod h) = (q, r)"
       
  4229     by (simp add: pdivmod_pdivmodrel)
       
  4230   with lc have "(f div g, f mod g) = (smult lc q, r)"
       
  4231     by (auto simp: h_def div_smult_right[OF lc] mod_smult_right[OF lc])
       
  4232   with id show ?thesis
       
  4233     by auto
       
  4234 qed
       
  4235 
       
  4236 lemma pdivmod_via_pseudo_divmod_list:
       
  4237   "(f div g, f mod g) =
       
  4238     (let cg = coeffs g in
       
  4239       if cg = [] then (0, f)
       
  4240       else
       
  4241         let
       
  4242           cf = coeffs f;
       
  4243           ilc = inverse (last cg);
       
  4244           ch = map (op * ilc) cg;
       
  4245           (q, r) = pseudo_divmod_main_list 1 [] (rev cf) (rev ch) (1 + length cf - length cg)
       
  4246         in (poly_of_list (map (op * ilc) q), poly_of_list (rev r)))"
       
  4247 proof -
       
  4248   note d = pdivmod_via_pseudo_divmod pseudo_divmod_impl pseudo_divmod_list_def
       
  4249   show ?thesis
       
  4250   proof (cases "g = 0")
       
  4251     case True
       
  4252     with d show ?thesis by auto
       
  4253   next
       
  4254     case False
       
  4255     define ilc where "ilc = inverse (coeff g (degree g))"
       
  4256     from False have ilc: "ilc \<noteq> 0"
       
  4257       by (auto simp: ilc_def)
       
  4258     with False have id: "g = 0 \<longleftrightarrow> False" "coeffs g = [] \<longleftrightarrow> False"
       
  4259       "last (coeffs g) = coeff g (degree g)"
       
  4260       "coeffs (smult ilc g) = [] \<longleftrightarrow> False"
       
  4261       by (auto simp: last_coeffs_eq_coeff_degree)
       
  4262     have id2: "hd (rev (coeffs (smult ilc g))) = 1"
       
  4263       by (subst hd_rev, insert id ilc, auto simp: coeffs_smult, subst last_map, auto simp: id ilc_def)
       
  4264     have id3: "length (coeffs (smult ilc g)) = length (coeffs g)"
       
  4265       "rev (coeffs (smult ilc g)) = rev (map (op * ilc) (coeffs g))"
       
  4266       unfolding coeffs_smult using ilc by auto
       
  4267     obtain q r where pair:
       
  4268       "pseudo_divmod_main_list 1 [] (rev (coeffs f)) (rev (map (op * ilc) (coeffs g)))
       
  4269         (1 + length (coeffs f) - length (coeffs g)) = (q, r)"
       
  4270       by force
       
  4271     show ?thesis
       
  4272       unfolding d Let_def id if_False ilc_def[symmetric] map_prod_def[symmetric] id2
       
  4273       unfolding id3 pair map_prod_def split
       
  4274       by (auto simp: Poly_map)
       
  4275   qed
       
  4276 qed
       
  4277 
       
  4278 lemma pseudo_divmod_main_list_1: "pseudo_divmod_main_list 1 = divmod_poly_one_main_list"
       
  4279 proof (intro ext, goal_cases)
       
  4280   case (1 q r d n)
       
  4281   have *: "map (op * 1) xs = xs" for xs :: "'a list"
       
  4282     by (induct xs) auto
       
  4283   show ?case
       
  4284     by (induct n arbitrary: q r d) (auto simp: * Let_def)
       
  4285 qed
       
  4286 
       
  4287 fun divide_poly_main_list :: "'a::idom_divide \<Rightarrow> 'a list \<Rightarrow> 'a list \<Rightarrow> 'a list \<Rightarrow> nat \<Rightarrow> 'a list"
       
  4288   where
       
  4289     "divide_poly_main_list lc q r d (Suc n) =
       
  4290       (let
       
  4291         cr = hd r
       
  4292         in if cr = 0 then divide_poly_main_list lc (cCons cr q) (tl r) d n else let
       
  4293         a = cr div lc;
       
  4294         qq = cCons a q;
       
  4295         rr = minus_poly_rev_list r (map (op * a) d)
       
  4296        in if hd rr = 0 then divide_poly_main_list lc qq (tl rr) d n else [])"
       
  4297   | "divide_poly_main_list lc q r d 0 = q"
       
  4298 
       
  4299 lemma divide_poly_main_list_simp [simp]:
       
  4300   "divide_poly_main_list lc q r d (Suc n) =
       
  4301     (let
       
  4302       cr = hd r;
       
  4303       a = cr div lc;
       
  4304       qq = cCons a q;
       
  4305       rr = minus_poly_rev_list r (map (op * a) d)
       
  4306      in if hd rr = 0 then divide_poly_main_list lc qq (tl rr) d n else [])"
       
  4307   by (simp add: Let_def minus_zero_does_nothing)
       
  4308 
       
  4309 declare divide_poly_main_list.simps(1)[simp del]
       
  4310 
       
  4311 definition divide_poly_list :: "'a::idom_divide poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly"
       
  4312   where "divide_poly_list f g =
       
  4313     (let cg = coeffs g in
       
  4314       if cg = [] then g
       
  4315       else
       
  4316         let
       
  4317           cf = coeffs f;
       
  4318           cgr = rev cg
       
  4319         in poly_of_list (divide_poly_main_list (hd cgr) [] (rev cf) cgr (1 + length cf - length cg)))"
       
  4320 
       
  4321 lemmas pdivmod_via_divmod_list = pdivmod_via_pseudo_divmod_list[unfolded pseudo_divmod_main_list_1]
       
  4322 
       
  4323 lemma mod_poly_one_main_list: "snd (divmod_poly_one_main_list q r d n) = mod_poly_one_main_list r d n"
       
  4324   by (induct n arbitrary: q r d) (auto simp: Let_def)
       
  4325 
       
  4326 lemma mod_poly_code [code]:
       
  4327   "f mod g =
       
  4328     (let cg = coeffs g in
       
  4329       if cg = [] then f
       
  4330       else
       
  4331         let
       
  4332           cf = coeffs f;
       
  4333           ilc = inverse (last cg);
       
  4334           ch = map (op * ilc) cg;
       
  4335           r = mod_poly_one_main_list (rev cf) (rev ch) (1 + length cf - length cg)
       
  4336         in poly_of_list (rev r))"
       
  4337   (is "_ = ?rhs")
       
  4338 proof -
       
  4339   have "snd (f div g, f mod g) = ?rhs"
       
  4340     unfolding pdivmod_via_divmod_list Let_def mod_poly_one_main_list [symmetric, of _ _ _ Nil]
       
  4341     by (auto split: prod.splits)
       
  4342   then show ?thesis by simp
       
  4343 qed
       
  4344 
       
  4345 definition div_field_poly_impl :: "'a :: field poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly"
       
  4346   where "div_field_poly_impl f g =
       
  4347     (let cg = coeffs g in
       
  4348       if cg = [] then 0
       
  4349       else
       
  4350         let
       
  4351           cf = coeffs f;
       
  4352           ilc = inverse (last cg);
       
  4353           ch = map (op * ilc) cg;
       
  4354           q = fst (divmod_poly_one_main_list [] (rev cf) (rev ch) (1 + length cf - length cg))
       
  4355         in poly_of_list ((map (op * ilc) q)))"
       
  4356 
       
  4357 text \<open>We do not declare the following lemma as code equation, since then polynomial division
       
  4358   on non-fields will no longer be executable. However, a code-unfold is possible, since
       
  4359   \<open>div_field_poly_impl\<close> is a bit more efficient than the generic polynomial division.\<close>
       
  4360 lemma div_field_poly_impl[code_unfold]: "op div = div_field_poly_impl"
       
  4361 proof (intro ext)
       
  4362   fix f g :: "'a poly"
       
  4363   have "fst (f div g, f mod g) = div_field_poly_impl f g"
       
  4364     unfolding div_field_poly_impl_def pdivmod_via_divmod_list Let_def
       
  4365     by (auto split: prod.splits)
       
  4366   then show "f div g =  div_field_poly_impl f g"
       
  4367     by simp
       
  4368 qed
       
  4369 
       
  4370 lemma divide_poly_main_list:
       
  4371   assumes lc0: "lc \<noteq> 0"
       
  4372     and lc: "last d = lc"
       
  4373     and d: "d \<noteq> []"
       
  4374     and "n = (1 + length r - length d)"
       
  4375   shows "Poly (divide_poly_main_list lc q (rev r) (rev d) n) =
       
  4376     divide_poly_main lc (monom 1 n * Poly q) (Poly r) (Poly d) (length r - 1) n"
       
  4377   using assms(4-)
       
  4378 proof (induct "n" arbitrary: r q)
       
  4379   case (Suc n)
       
  4380   from Suc.prems have ifCond: "\<not> Suc (length r) \<le> length d"
       
  4381     by simp
       
  4382   with d have r: "r \<noteq> []"
       
  4383     using Suc_leI length_greater_0_conv list.size(3) by fastforce
       
  4384   then obtain rr lcr where r: "r = rr @ [lcr]"
       
  4385     by (cases r rule: rev_cases) auto
       
  4386   from d lc obtain dd where d: "d = dd @ [lc]"
       
  4387     by (cases d rule: rev_cases) auto
       
  4388   from Suc(2) ifCond have n: "n = 1 + length rr - length d"
       
  4389     by (auto simp: r)
       
  4390   from ifCond have len: "length dd \<le> length rr"
       
  4391     by (simp add: r d)
       
  4392   show ?case
       
  4393   proof (cases "lcr div lc * lc = lcr")
       
  4394     case False
       
  4395     with r d show ?thesis
       
  4396       unfolding Suc(2)[symmetric]
       
  4397       by (auto simp add: Let_def nth_default_append)
       
  4398   next
       
  4399     case True
       
  4400     with r d have id:
       
  4401       "?thesis \<longleftrightarrow>
       
  4402         Poly (divide_poly_main_list lc (cCons (lcr div lc) q)
       
  4403           (rev (rev (minus_poly_rev_list (rev rr) (rev (map (op * (lcr div lc)) dd))))) (rev d) n) =
       
  4404           divide_poly_main lc
       
  4405             (monom 1 (Suc n) * Poly q + monom (lcr div lc) n)
       
  4406             (Poly r - monom (lcr div lc) n * Poly d)
       
  4407             (Poly d) (length rr - 1) n"
       
  4408       by (cases r rule: rev_cases; cases "d" rule: rev_cases)
       
  4409         (auto simp add: Let_def rev_map nth_default_append)
       
  4410     have cong: "\<And>x1 x2 x3 x4 y1 y2 y3 y4. x1 = y1 \<Longrightarrow> x2 = y2 \<Longrightarrow> x3 = y3 \<Longrightarrow> x4 = y4 \<Longrightarrow>
       
  4411         divide_poly_main lc x1 x2 x3 x4 n = divide_poly_main lc y1 y2 y3 y4 n"
       
  4412       by simp
       
  4413     show ?thesis
       
  4414       unfolding id
       
  4415     proof (subst Suc(1), simp add: n,
       
  4416         subst minus_poly_rev_list, force simp: len, rule cong[OF _ _ refl], goal_cases)
       
  4417       case 2
       
  4418       have "monom lcr (length rr) = monom (lcr div lc) (length rr - length dd) * monom lc (length dd)"
       
  4419         by (simp add: mult_monom len True)
       
  4420       then show ?case unfolding r d Poly_append n ring_distribs
       
  4421         by (auto simp: Poly_map smult_monom smult_monom_mult)
       
  4422     qed (auto simp: len monom_Suc smult_monom)
       
  4423   qed
       
  4424 qed simp
       
  4425 
       
  4426 lemma divide_poly_list[code]: "f div g = divide_poly_list f g"
       
  4427 proof -
       
  4428   note d = divide_poly_def divide_poly_list_def
       
  4429   show ?thesis
       
  4430   proof (cases "g = 0")
       
  4431     case True
       
  4432     show ?thesis by (auto simp: d True)
       
  4433   next
       
  4434     case False
       
  4435     then obtain cg lcg where cg: "coeffs g = cg @ [lcg]"
       
  4436       by (cases "coeffs g" rule: rev_cases) auto
       
  4437     with False have id: "(g = 0) = False" "(cg @ [lcg] = []) = False"
       
  4438       by auto
       
  4439     from cg False have lcg: "coeff g (degree g) = lcg"
       
  4440       using last_coeffs_eq_coeff_degree last_snoc by force
       
  4441     with False have "lcg \<noteq> 0" by auto
       
  4442     from cg Poly_coeffs [of g] have ltp: "Poly (cg @ [lcg]) = g"
       
  4443       by auto
       
  4444     show ?thesis
       
  4445       unfolding d cg Let_def id if_False poly_of_list_def
       
  4446       by (subst divide_poly_main_list, insert False cg \<open>lcg \<noteq> 0\<close>)
       
  4447         (auto simp: lcg ltp, simp add: degree_eq_length_coeffs)
       
  4448   qed
       
  4449 qed
       
  4450 
       
  4451 no_notation cCons (infixr "##" 65)
       
  4452 
       
  4453 end