src/HOL/Computational_Algebra/Polynomial_Factorial.thy
author haftmann
Sun Oct 08 22:28:21 2017 +0200 (20 months ago)
changeset 66808 1907167b6038
parent 66806 a4e82b58d833
child 66817 0b12755ccbb2
permissions -rw-r--r--
elementary definition of division on natural numbers
     1 (*  Title:      HOL/Computational_Algebra/Polynomial_Factorial.thy
     2     Author:     Manuel Eberl
     3 *)
     4 
     5 section \<open>Polynomials, fractions and rings\<close>
     6 
     7 theory Polynomial_Factorial
     8 imports 
     9   Complex_Main
    10   Polynomial
    11   Normalized_Fraction
    12 begin
    13 
    14 subsection \<open>Lifting elements into the field of fractions\<close>
    15 
    16 definition to_fract :: "'a :: idom \<Rightarrow> 'a fract"
    17   where "to_fract x = Fract x 1"
    18   \<comment> \<open>FIXME: more idiomatic name, abbreviation\<close>
    19 
    20 lemma to_fract_0 [simp]: "to_fract 0 = 0"
    21   by (simp add: to_fract_def eq_fract Zero_fract_def)
    22 
    23 lemma to_fract_1 [simp]: "to_fract 1 = 1"
    24   by (simp add: to_fract_def eq_fract One_fract_def)
    25 
    26 lemma to_fract_add [simp]: "to_fract (x + y) = to_fract x + to_fract y"
    27   by (simp add: to_fract_def)
    28 
    29 lemma to_fract_diff [simp]: "to_fract (x - y) = to_fract x - to_fract y"
    30   by (simp add: to_fract_def)
    31   
    32 lemma to_fract_uminus [simp]: "to_fract (-x) = -to_fract x"
    33   by (simp add: to_fract_def)
    34   
    35 lemma to_fract_mult [simp]: "to_fract (x * y) = to_fract x * to_fract y"
    36   by (simp add: to_fract_def)
    37 
    38 lemma to_fract_eq_iff [simp]: "to_fract x = to_fract y \<longleftrightarrow> x = y"
    39   by (simp add: to_fract_def eq_fract)
    40   
    41 lemma to_fract_eq_0_iff [simp]: "to_fract x = 0 \<longleftrightarrow> x = 0"
    42   by (simp add: to_fract_def Zero_fract_def eq_fract)
    43 
    44 lemma snd_quot_of_fract_nonzero [simp]: "snd (quot_of_fract x) \<noteq> 0"
    45   by transfer simp
    46 
    47 lemma Fract_quot_of_fract [simp]: "Fract (fst (quot_of_fract x)) (snd (quot_of_fract x)) = x"
    48   by transfer (simp del: fractrel_iff, subst fractrel_normalize_quot_left, simp)
    49 
    50 lemma to_fract_quot_of_fract:
    51   assumes "snd (quot_of_fract x) = 1"
    52   shows   "to_fract (fst (quot_of_fract x)) = x"
    53 proof -
    54   have "x = Fract (fst (quot_of_fract x)) (snd (quot_of_fract x))" by simp
    55   also note assms
    56   finally show ?thesis by (simp add: to_fract_def)
    57 qed
    58 
    59 lemma snd_quot_of_fract_Fract_whole:
    60   assumes "y dvd x"
    61   shows   "snd (quot_of_fract (Fract x y)) = 1"
    62   using assms by transfer (auto simp: normalize_quot_def Let_def gcd_proj2_if_dvd)
    63   
    64 lemma Fract_conv_to_fract: "Fract a b = to_fract a / to_fract b"
    65   by (simp add: to_fract_def)
    66 
    67 lemma quot_of_fract_to_fract [simp]: "quot_of_fract (to_fract x) = (x, 1)"
    68   unfolding to_fract_def by transfer (simp add: normalize_quot_def)
    69 
    70 lemma fst_quot_of_fract_eq_0_iff [simp]: "fst (quot_of_fract x) = 0 \<longleftrightarrow> x = 0"
    71   by transfer simp
    72  
    73 lemma snd_quot_of_fract_to_fract [simp]: "snd (quot_of_fract (to_fract x)) = 1"
    74   unfolding to_fract_def by (rule snd_quot_of_fract_Fract_whole) simp_all
    75 
    76 lemma coprime_quot_of_fract:
    77   "coprime (fst (quot_of_fract x)) (snd (quot_of_fract x))"
    78   by transfer (simp add: coprime_normalize_quot)
    79 
    80 lemma unit_factor_snd_quot_of_fract: "unit_factor (snd (quot_of_fract x)) = 1"
    81   using quot_of_fract_in_normalized_fracts[of x] 
    82   by (simp add: normalized_fracts_def case_prod_unfold)  
    83 
    84 lemma unit_factor_1_imp_normalized: "unit_factor x = 1 \<Longrightarrow> normalize x = x"
    85   by (subst (2) normalize_mult_unit_factor [symmetric, of x])
    86      (simp del: normalize_mult_unit_factor)
    87   
    88 lemma normalize_snd_quot_of_fract: "normalize (snd (quot_of_fract x)) = snd (quot_of_fract x)"
    89   by (intro unit_factor_1_imp_normalized unit_factor_snd_quot_of_fract)
    90 
    91 
    92 subsection \<open>Lifting polynomial coefficients to the field of fractions\<close>
    93 
    94 abbreviation (input) fract_poly 
    95   where "fract_poly \<equiv> map_poly to_fract"
    96 
    97 abbreviation (input) unfract_poly 
    98   where "unfract_poly \<equiv> map_poly (fst \<circ> quot_of_fract)"
    99   
   100 lemma fract_poly_smult [simp]: "fract_poly (smult c p) = smult (to_fract c) (fract_poly p)"
   101   by (simp add: smult_conv_map_poly map_poly_map_poly o_def)
   102 
   103 lemma fract_poly_0 [simp]: "fract_poly 0 = 0"
   104   by (simp add: poly_eqI coeff_map_poly)
   105 
   106 lemma fract_poly_1 [simp]: "fract_poly 1 = 1"
   107   by (simp add: map_poly_pCons)
   108 
   109 lemma fract_poly_add [simp]:
   110   "fract_poly (p + q) = fract_poly p + fract_poly q"
   111   by (intro poly_eqI) (simp_all add: coeff_map_poly)
   112 
   113 lemma fract_poly_diff [simp]:
   114   "fract_poly (p - q) = fract_poly p - fract_poly q"
   115   by (intro poly_eqI) (simp_all add: coeff_map_poly)
   116 
   117 lemma to_fract_sum [simp]: "to_fract (sum f A) = sum (\<lambda>x. to_fract (f x)) A"
   118   by (cases "finite A", induction A rule: finite_induct) simp_all 
   119 
   120 lemma fract_poly_mult [simp]:
   121   "fract_poly (p * q) = fract_poly p * fract_poly q"
   122   by (intro poly_eqI) (simp_all add: coeff_map_poly coeff_mult)
   123 
   124 lemma fract_poly_eq_iff [simp]: "fract_poly p = fract_poly q \<longleftrightarrow> p = q"
   125   by (auto simp: poly_eq_iff coeff_map_poly)
   126 
   127 lemma fract_poly_eq_0_iff [simp]: "fract_poly p = 0 \<longleftrightarrow> p = 0"
   128   using fract_poly_eq_iff[of p 0] by (simp del: fract_poly_eq_iff)
   129 
   130 lemma fract_poly_dvd: "p dvd q \<Longrightarrow> fract_poly p dvd fract_poly q"
   131   by (auto elim!: dvdE)
   132 
   133 lemma prod_mset_fract_poly: 
   134   "(\<Prod>x\<in>#A. map_poly to_fract (f x)) = fract_poly (prod_mset (image_mset f A))"
   135   by (induct A) (simp_all add: ac_simps)
   136   
   137 lemma is_unit_fract_poly_iff:
   138   "p dvd 1 \<longleftrightarrow> fract_poly p dvd 1 \<and> content p = 1"
   139 proof safe
   140   assume A: "p dvd 1"
   141   with fract_poly_dvd [of p 1] show "is_unit (fract_poly p)"
   142     by simp
   143   from A show "content p = 1"
   144     by (auto simp: is_unit_poly_iff normalize_1_iff)
   145 next
   146   assume A: "fract_poly p dvd 1" and B: "content p = 1"
   147   from A obtain c where c: "fract_poly p = [:c:]" by (auto simp: is_unit_poly_iff)
   148   {
   149     fix n :: nat assume "n > 0"
   150     have "to_fract (coeff p n) = coeff (fract_poly p) n" by (simp add: coeff_map_poly)
   151     also note c
   152     also from \<open>n > 0\<close> have "coeff [:c:] n = 0" by (simp add: coeff_pCons split: nat.splits)
   153     finally have "coeff p n = 0" by simp
   154   }
   155   hence "degree p \<le> 0" by (intro degree_le) simp_all
   156   with B show "p dvd 1" by (auto simp: is_unit_poly_iff normalize_1_iff elim!: degree_eq_zeroE)
   157 qed
   158   
   159 lemma fract_poly_is_unit: "p dvd 1 \<Longrightarrow> fract_poly p dvd 1"
   160   using fract_poly_dvd[of p 1] by simp
   161 
   162 lemma fract_poly_smult_eqE:
   163   fixes c :: "'a :: {idom_divide,ring_gcd} fract"
   164   assumes "fract_poly p = smult c (fract_poly q)"
   165   obtains a b 
   166     where "c = to_fract b / to_fract a" "smult a p = smult b q" "coprime a b" "normalize a = a"
   167 proof -
   168   define a b where "a = fst (quot_of_fract c)" and "b = snd (quot_of_fract c)"
   169   have "smult (to_fract a) (fract_poly q) = smult (to_fract b) (fract_poly p)"
   170     by (subst smult_eq_iff) (simp_all add: a_def b_def Fract_conv_to_fract [symmetric] assms)
   171   hence "fract_poly (smult a q) = fract_poly (smult b p)" by (simp del: fract_poly_eq_iff)
   172   hence "smult b p = smult a q" by (simp only: fract_poly_eq_iff)
   173   moreover have "c = to_fract a / to_fract b" "coprime b a" "normalize b = b"
   174     by (simp_all add: a_def b_def coprime_quot_of_fract gcd.commute
   175           normalize_snd_quot_of_fract Fract_conv_to_fract [symmetric])
   176   ultimately show ?thesis by (intro that[of a b])
   177 qed
   178 
   179 
   180 subsection \<open>Fractional content\<close>
   181 
   182 abbreviation (input) Lcm_coeff_denoms 
   183     :: "'a :: {semiring_Gcd,idom_divide,ring_gcd} fract poly \<Rightarrow> 'a"
   184   where "Lcm_coeff_denoms p \<equiv> Lcm (snd ` quot_of_fract ` set (coeffs p))"
   185   
   186 definition fract_content :: 
   187       "'a :: {factorial_semiring,semiring_Gcd,ring_gcd,idom_divide} fract poly \<Rightarrow> 'a fract" where
   188   "fract_content p = 
   189      (let d = Lcm_coeff_denoms p in Fract (content (unfract_poly (smult (to_fract d) p))) d)" 
   190 
   191 definition primitive_part_fract :: 
   192       "'a :: {factorial_semiring,semiring_Gcd,ring_gcd,idom_divide} fract poly \<Rightarrow> 'a poly" where
   193   "primitive_part_fract p = 
   194      primitive_part (unfract_poly (smult (to_fract (Lcm_coeff_denoms p)) p))"
   195 
   196 lemma primitive_part_fract_0 [simp]: "primitive_part_fract 0 = 0"
   197   by (simp add: primitive_part_fract_def)
   198 
   199 lemma fract_content_eq_0_iff [simp]:
   200   "fract_content p = 0 \<longleftrightarrow> p = 0"
   201   unfolding fract_content_def Let_def Zero_fract_def
   202   by (subst eq_fract) (auto simp: Lcm_0_iff map_poly_eq_0_iff)
   203 
   204 lemma content_primitive_part_fract [simp]: "p \<noteq> 0 \<Longrightarrow> content (primitive_part_fract p) = 1"
   205   unfolding primitive_part_fract_def
   206   by (rule content_primitive_part)
   207      (auto simp: primitive_part_fract_def map_poly_eq_0_iff Lcm_0_iff)  
   208 
   209 lemma content_times_primitive_part_fract:
   210   "smult (fract_content p) (fract_poly (primitive_part_fract p)) = p"
   211 proof -
   212   define p' where "p' = unfract_poly (smult (to_fract (Lcm_coeff_denoms p)) p)"
   213   have "fract_poly p' = 
   214           map_poly (to_fract \<circ> fst \<circ> quot_of_fract) (smult (to_fract (Lcm_coeff_denoms p)) p)"
   215     unfolding primitive_part_fract_def p'_def 
   216     by (subst map_poly_map_poly) (simp_all add: o_assoc)
   217   also have "\<dots> = smult (to_fract (Lcm_coeff_denoms p)) p"
   218   proof (intro map_poly_idI, unfold o_apply)
   219     fix c assume "c \<in> set (coeffs (smult (to_fract (Lcm_coeff_denoms p)) p))"
   220     then obtain c' where c: "c' \<in> set (coeffs p)" "c = to_fract (Lcm_coeff_denoms p) * c'"
   221       by (auto simp add: Lcm_0_iff coeffs_smult split: if_splits)
   222     note c(2)
   223     also have "c' = Fract (fst (quot_of_fract c')) (snd (quot_of_fract c'))"
   224       by simp
   225     also have "to_fract (Lcm_coeff_denoms p) * \<dots> = 
   226                  Fract (Lcm_coeff_denoms p * fst (quot_of_fract c')) (snd (quot_of_fract c'))"
   227       unfolding to_fract_def by (subst mult_fract) simp_all
   228     also have "snd (quot_of_fract \<dots>) = 1"
   229       by (intro snd_quot_of_fract_Fract_whole dvd_mult2 dvd_Lcm) (insert c(1), auto)
   230     finally show "to_fract (fst (quot_of_fract c)) = c"
   231       by (rule to_fract_quot_of_fract)
   232   qed
   233   also have "p' = smult (content p') (primitive_part p')" 
   234     by (rule content_times_primitive_part [symmetric])
   235   also have "primitive_part p' = primitive_part_fract p"
   236     by (simp add: primitive_part_fract_def p'_def)
   237   also have "fract_poly (smult (content p') (primitive_part_fract p)) = 
   238                smult (to_fract (content p')) (fract_poly (primitive_part_fract p))" by simp
   239   finally have "smult (to_fract (content p')) (fract_poly (primitive_part_fract p)) =
   240                       smult (to_fract (Lcm_coeff_denoms p)) p" .
   241   thus ?thesis
   242     by (subst (asm) smult_eq_iff)
   243        (auto simp add: Let_def p'_def Fract_conv_to_fract field_simps Lcm_0_iff fract_content_def)
   244 qed
   245 
   246 lemma fract_content_fract_poly [simp]: "fract_content (fract_poly p) = to_fract (content p)"
   247 proof -
   248   have "Lcm_coeff_denoms (fract_poly p) = 1"
   249     by (auto simp: set_coeffs_map_poly)
   250   hence "fract_content (fract_poly p) = 
   251            to_fract (content (map_poly (fst \<circ> quot_of_fract \<circ> to_fract) p))"
   252     by (simp add: fract_content_def to_fract_def fract_collapse map_poly_map_poly del: Lcm_1_iff)
   253   also have "map_poly (fst \<circ> quot_of_fract \<circ> to_fract) p = p"
   254     by (intro map_poly_idI) simp_all
   255   finally show ?thesis .
   256 qed
   257 
   258 lemma content_decompose_fract:
   259   fixes p :: "'a :: {factorial_semiring,semiring_Gcd,ring_gcd,idom_divide} fract poly"
   260   obtains c p' where "p = smult c (map_poly to_fract p')" "content p' = 1"
   261 proof (cases "p = 0")
   262   case True
   263   hence "p = smult 0 (map_poly to_fract 1)" "content 1 = 1" by simp_all
   264   thus ?thesis ..
   265 next
   266   case False
   267   thus ?thesis
   268     by (rule that[OF content_times_primitive_part_fract [symmetric] content_primitive_part_fract])
   269 qed
   270 
   271 
   272 subsection \<open>More properties of content and primitive part\<close>
   273 
   274 lemma lift_prime_elem_poly:
   275   assumes "prime_elem (c :: 'a :: semidom)"
   276   shows   "prime_elem [:c:]"
   277 proof (rule prime_elemI)
   278   fix a b assume *: "[:c:] dvd a * b"
   279   from * have dvd: "c dvd coeff (a * b) n" for n
   280     by (subst (asm) const_poly_dvd_iff) blast
   281   {
   282     define m where "m = (GREATEST m. \<not>c dvd coeff b m)"
   283     assume "\<not>[:c:] dvd b"
   284     hence A: "\<exists>i. \<not>c dvd coeff b i" by (subst (asm) const_poly_dvd_iff) blast
   285     have B: "\<forall>i. \<not>c dvd coeff b i \<longrightarrow> i \<le> degree b"
   286       by (auto intro: le_degree)
   287     have coeff_m: "\<not>c dvd coeff b m" unfolding m_def by (rule GreatestI_ex_nat[OF A B])
   288     have "i \<le> m" if "\<not>c dvd coeff b i" for i
   289       unfolding m_def by (rule Greatest_le_nat[OF that B])
   290     hence dvd_b: "c dvd coeff b i" if "i > m" for i using that by force
   291 
   292     have "c dvd coeff a i" for i
   293     proof (induction i rule: nat_descend_induct[of "degree a"])
   294       case (base i)
   295       thus ?case by (simp add: coeff_eq_0)
   296     next
   297       case (descend i)
   298       let ?A = "{..i+m} - {i}"
   299       have "c dvd coeff (a * b) (i + m)" by (rule dvd)
   300       also have "coeff (a * b) (i + m) = (\<Sum>k\<le>i + m. coeff a k * coeff b (i + m - k))"
   301         by (simp add: coeff_mult)
   302       also have "{..i+m} = insert i ?A" by auto
   303       also have "(\<Sum>k\<in>\<dots>. coeff a k * coeff b (i + m - k)) =
   304                    coeff a i * coeff b m + (\<Sum>k\<in>?A. coeff a k * coeff b (i + m - k))"
   305         (is "_ = _ + ?S")
   306         by (subst sum.insert) simp_all
   307       finally have eq: "c dvd coeff a i * coeff b m + ?S" .
   308       moreover have "c dvd ?S"
   309       proof (rule dvd_sum)
   310         fix k assume k: "k \<in> {..i+m} - {i}"
   311         show "c dvd coeff a k * coeff b (i + m - k)"
   312         proof (cases "k < i")
   313           case False
   314           with k have "c dvd coeff a k" by (intro descend.IH) simp
   315           thus ?thesis by simp
   316         next
   317           case True
   318           hence "c dvd coeff b (i + m - k)" by (intro dvd_b) simp
   319           thus ?thesis by simp
   320         qed
   321       qed
   322       ultimately have "c dvd coeff a i * coeff b m"
   323         by (simp add: dvd_add_left_iff)
   324       with assms coeff_m show "c dvd coeff a i"
   325         by (simp add: prime_elem_dvd_mult_iff)
   326     qed
   327     hence "[:c:] dvd a" by (subst const_poly_dvd_iff) blast
   328   }
   329   then show "[:c:] dvd a \<or> [:c:] dvd b" by blast
   330 next
   331   from assms show "[:c:] \<noteq> 0" and "\<not> [:c:] dvd 1"
   332     by (simp_all add: prime_elem_def is_unit_const_poly_iff)
   333 qed
   334 
   335 lemma prime_elem_const_poly_iff:
   336   fixes c :: "'a :: semidom"
   337   shows   "prime_elem [:c:] \<longleftrightarrow> prime_elem c"
   338 proof
   339   assume A: "prime_elem [:c:]"
   340   show "prime_elem c"
   341   proof (rule prime_elemI)
   342     fix a b assume "c dvd a * b"
   343     hence "[:c:] dvd [:a:] * [:b:]" by (simp add: mult_ac)
   344     from A and this have "[:c:] dvd [:a:] \<or> [:c:] dvd [:b:]" by (rule prime_elem_dvd_multD)
   345     thus "c dvd a \<or> c dvd b" by simp
   346   qed (insert A, auto simp: prime_elem_def is_unit_poly_iff)
   347 qed (auto intro: lift_prime_elem_poly)
   348 
   349 context
   350 begin
   351 
   352 private lemma content_1_mult:
   353   fixes f g :: "'a :: {semiring_Gcd,factorial_semiring} poly"
   354   assumes "content f = 1" "content g = 1"
   355   shows   "content (f * g) = 1"
   356 proof (cases "f * g = 0")
   357   case False
   358   from assms have "f \<noteq> 0" "g \<noteq> 0" by auto
   359 
   360   hence "f * g \<noteq> 0" by auto
   361   {
   362     assume "\<not>is_unit (content (f * g))"
   363     with False have "\<exists>p. p dvd content (f * g) \<and> prime p"
   364       by (intro prime_divisor_exists) simp_all
   365     then obtain p where "p dvd content (f * g)" "prime p" by blast
   366     from \<open>p dvd content (f * g)\<close> have "[:p:] dvd f * g"
   367       by (simp add: const_poly_dvd_iff_dvd_content)
   368     moreover from \<open>prime p\<close> have "prime_elem [:p:]" by (simp add: lift_prime_elem_poly)
   369     ultimately have "[:p:] dvd f \<or> [:p:] dvd g"
   370       by (simp add: prime_elem_dvd_mult_iff)
   371     with assms have "is_unit p" by (simp add: const_poly_dvd_iff_dvd_content)
   372     with \<open>prime p\<close> have False by simp
   373   }
   374   hence "is_unit (content (f * g))" by blast
   375   hence "normalize (content (f * g)) = 1" by (simp add: is_unit_normalize del: normalize_content)
   376   thus ?thesis by simp
   377 qed (insert assms, auto)
   378 
   379 lemma content_mult:
   380   fixes p q :: "'a :: {factorial_semiring, semiring_Gcd} poly"
   381   shows "content (p * q) = content p * content q"
   382 proof -
   383   from content_decompose[of p] guess p' . note p = this
   384   from content_decompose[of q] guess q' . note q = this
   385   have "content (p * q) = content p * content q * content (p' * q')"
   386     by (subst p, subst q) (simp add: mult_ac normalize_mult)
   387   also from p q have "content (p' * q') = 1" by (intro content_1_mult)
   388   finally show ?thesis by simp
   389 qed
   390 
   391 lemma fract_poly_dvdD:
   392   fixes p :: "'a :: {factorial_semiring,semiring_Gcd,ring_gcd,idom_divide} poly"
   393   assumes "fract_poly p dvd fract_poly q" "content p = 1"
   394   shows   "p dvd q"
   395 proof -
   396   from assms(1) obtain r where r: "fract_poly q = fract_poly p * r" by (erule dvdE)
   397   from content_decompose_fract[of r] guess c r' . note r' = this
   398   from r r' have eq: "fract_poly q = smult c (fract_poly (p * r'))" by simp  
   399   from fract_poly_smult_eqE[OF this] guess a b . note ab = this
   400   have "content (smult a q) = content (smult b (p * r'))" by (simp only: ab(2))
   401   hence eq': "normalize b = a * content q" by (simp add: assms content_mult r' ab(4))
   402   have "1 = gcd a (normalize b)" by (simp add: ab)
   403   also note eq'
   404   also have "gcd a (a * content q) = a" by (simp add: gcd_proj1_if_dvd ab(4))
   405   finally have [simp]: "a = 1" by simp
   406   from eq ab have "q = p * ([:b:] * r')" by simp
   407   thus ?thesis by (rule dvdI)
   408 qed
   409 
   410 end
   411 
   412 
   413 subsection \<open>Polynomials over a field are a Euclidean ring\<close>
   414 
   415 context
   416 begin
   417 
   418 interpretation field_poly: 
   419   unique_euclidean_ring where zero = "0 :: 'a :: field poly"
   420     and one = 1 and plus = plus
   421     and uminus = uminus and minus = minus
   422     and times = times
   423     and normalize = "\<lambda>p. smult (inverse (lead_coeff p)) p"
   424     and unit_factor = "\<lambda>p. [:lead_coeff p:]"
   425     and euclidean_size = "\<lambda>p. if p = 0 then 0 else 2 ^ degree p"
   426     and uniqueness_constraint = top
   427     and divide = divide and modulo = modulo
   428   rewrites "dvd.dvd (times :: 'a poly \<Rightarrow> _) = Rings.dvd"
   429     and "comm_monoid_mult.prod_mset times 1 = prod_mset"
   430     and "comm_semiring_1.irreducible times 1 0 = irreducible"
   431     and "comm_semiring_1.prime_elem times 1 0 = prime_elem"
   432 proof -
   433   show "dvd.dvd (times :: 'a poly \<Rightarrow> _) = Rings.dvd"
   434     by (simp add: dvd_dict)
   435   show "comm_monoid_mult.prod_mset times 1 = prod_mset"
   436     by (simp add: prod_mset_dict)
   437   show "comm_semiring_1.irreducible times 1 0 = irreducible"
   438     by (simp add: irreducible_dict)
   439   show "comm_semiring_1.prime_elem times 1 0 = prime_elem"
   440     by (simp add: prime_elem_dict)
   441   show "class.unique_euclidean_ring divide plus minus (0 :: 'a poly) times 1
   442     (\<lambda>p. [:lead_coeff p:]) (\<lambda>p. smult (inverse (lead_coeff p)) p) modulo
   443     (\<lambda>p. if p = 0 then 0 else 2 ^ degree p) uminus top"
   444   proof (standard, fold dvd_dict)
   445     fix p :: "'a poly"
   446     show "[:lead_coeff p:] * smult (inverse (lead_coeff p)) p = p"
   447       by (cases "p = 0") simp_all
   448   next
   449     fix p :: "'a poly" assume "is_unit p"
   450     then show "[:lead_coeff p:] = p"
   451       by (elim is_unit_polyE) (auto simp: monom_0 one_poly_def field_simps)
   452   next
   453     fix p :: "'a poly" assume "p \<noteq> 0"
   454     then show "is_unit [:lead_coeff p:]"
   455       by (simp add: is_unit_pCons_iff)
   456   next
   457     fix p q s :: "'a poly" assume "s \<noteq> 0"
   458     moreover assume "(if p = 0 then (0::nat) else 2 ^ degree p) < (if q = 0 then 0 else 2 ^ degree q)"
   459     ultimately show "(if p * s = 0 then (0::nat) else 2 ^ degree (p * s)) < (if q * s = 0 then 0 else 2 ^ degree (q * s))"
   460       by (auto simp add: degree_mult_eq)
   461   next
   462     fix p q r :: "'a poly"
   463     assume "p \<noteq> 0"
   464     with eucl_rel_poly [of r p] have "eucl_rel_poly (r + q * p) p (q + r div p, r mod p)"
   465       by (simp add: eucl_rel_poly_iff distrib_right)
   466     then have "(r + q * p) div p = q + r div p"
   467       by (rule div_poly_eq)
   468     moreover assume "(if r = 0 then (0::nat) else 2 ^ degree r)
   469       < (if p = 0 then 0 else 2 ^ degree p)"
   470     ultimately show "(q * p + r) div p = q"
   471       using \<open>p \<noteq> 0\<close>
   472       by (cases "r = 0") (simp_all add: div_poly_less ac_simps)
   473   qed (auto simp: lead_coeff_mult Rings.div_mult_mod_eq intro!: degree_mod_less' degree_mult_right_le)
   474 qed
   475 
   476 lemma field_poly_irreducible_imp_prime:
   477   "prime_elem p" if "irreducible p" for p :: "'a :: field poly"
   478   using that by (fact field_poly.irreducible_imp_prime_elem)
   479 
   480 lemma field_poly_prod_mset_prime_factorization:
   481   "prod_mset (field_poly.prime_factorization p) = smult (inverse (lead_coeff p)) p"
   482   if "p \<noteq> 0" for p :: "'a :: field poly"
   483   using that by (fact field_poly.prod_mset_prime_factorization)
   484 
   485 lemma field_poly_in_prime_factorization_imp_prime:
   486   "prime_elem p" if "p \<in># field_poly.prime_factorization x"
   487   for p :: "'a :: field poly"
   488   by (rule field_poly.prime_imp_prime_elem, rule field_poly.in_prime_factors_imp_prime)
   489     (fact that)
   490 
   491 
   492 subsection \<open>Primality and irreducibility in polynomial rings\<close>
   493 
   494 lemma nonconst_poly_irreducible_iff:
   495   fixes p :: "'a :: {factorial_semiring,semiring_Gcd,ring_gcd,idom_divide} poly"
   496   assumes "degree p \<noteq> 0"
   497   shows   "irreducible p \<longleftrightarrow> irreducible (fract_poly p) \<and> content p = 1"
   498 proof safe
   499   assume p: "irreducible p"
   500 
   501   from content_decompose[of p] guess p' . note p' = this
   502   hence "p = [:content p:] * p'" by simp
   503   from p this have "[:content p:] dvd 1 \<or> p' dvd 1" by (rule irreducibleD)
   504   moreover have "\<not>p' dvd 1"
   505   proof
   506     assume "p' dvd 1"
   507     hence "degree p = 0" by (subst p') (auto simp: is_unit_poly_iff)
   508     with assms show False by contradiction
   509   qed
   510   ultimately show [simp]: "content p = 1" by (simp add: is_unit_const_poly_iff)
   511   
   512   show "irreducible (map_poly to_fract p)"
   513   proof (rule irreducibleI)
   514     have "fract_poly p = 0 \<longleftrightarrow> p = 0" by (intro map_poly_eq_0_iff) auto
   515     with assms show "map_poly to_fract p \<noteq> 0" by auto
   516   next
   517     show "\<not>is_unit (fract_poly p)"
   518     proof
   519       assume "is_unit (map_poly to_fract p)"
   520       hence "degree (map_poly to_fract p) = 0"
   521         by (auto simp: is_unit_poly_iff)
   522       hence "degree p = 0" by (simp add: degree_map_poly)
   523       with assms show False by contradiction
   524    qed
   525  next
   526    fix q r assume qr: "fract_poly p = q * r"
   527    from content_decompose_fract[of q] guess cg q' . note q = this
   528    from content_decompose_fract[of r] guess cr r' . note r = this
   529    from qr q r p have nz: "cg \<noteq> 0" "cr \<noteq> 0" by auto
   530    from qr have eq: "fract_poly p = smult (cr * cg) (fract_poly (q' * r'))"
   531      by (simp add: q r)
   532    from fract_poly_smult_eqE[OF this] guess a b . note ab = this
   533    hence "content (smult a p) = content (smult b (q' * r'))" by (simp only:)
   534    with ab(4) have a: "a = normalize b" by (simp add: content_mult q r)
   535    hence "normalize b = gcd a b" by simp
   536    also from ab(3) have "\<dots> = 1" .
   537    finally have "a = 1" "is_unit b" by (simp_all add: a normalize_1_iff)
   538    
   539    note eq
   540    also from ab(1) \<open>a = 1\<close> have "cr * cg = to_fract b" by simp
   541    also have "smult \<dots> (fract_poly (q' * r')) = fract_poly (smult b (q' * r'))" by simp
   542    finally have "p = ([:b:] * q') * r'" by (simp del: fract_poly_smult)
   543    from p and this have "([:b:] * q') dvd 1 \<or> r' dvd 1" by (rule irreducibleD)
   544    hence "q' dvd 1 \<or> r' dvd 1" by (auto dest: dvd_mult_right simp del: mult_pCons_left)
   545    hence "fract_poly q' dvd 1 \<or> fract_poly r' dvd 1" by (auto simp: fract_poly_is_unit)
   546    with q r show "is_unit q \<or> is_unit r"
   547      by (auto simp add: is_unit_smult_iff dvd_field_iff nz)
   548  qed
   549 
   550 next
   551 
   552   assume irred: "irreducible (fract_poly p)" and primitive: "content p = 1"
   553   show "irreducible p"
   554   proof (rule irreducibleI)
   555     from irred show "p \<noteq> 0" by auto
   556   next
   557     from irred show "\<not>p dvd 1"
   558       by (auto simp: irreducible_def dest: fract_poly_is_unit)
   559   next
   560     fix q r assume qr: "p = q * r"
   561     hence "fract_poly p = fract_poly q * fract_poly r" by simp
   562     from irred and this have "fract_poly q dvd 1 \<or> fract_poly r dvd 1" 
   563       by (rule irreducibleD)
   564     with primitive qr show "q dvd 1 \<or> r dvd 1"
   565       by (auto simp:  content_prod_eq_1_iff is_unit_fract_poly_iff)
   566   qed
   567 qed
   568 
   569 private lemma irreducible_imp_prime_poly:
   570   fixes p :: "'a :: {factorial_semiring,semiring_Gcd,ring_gcd,idom_divide} poly"
   571   assumes "irreducible p"
   572   shows   "prime_elem p"
   573 proof (cases "degree p = 0")
   574   case True
   575   with assms show ?thesis
   576     by (auto simp: prime_elem_const_poly_iff irreducible_const_poly_iff
   577              intro!: irreducible_imp_prime_elem elim!: degree_eq_zeroE)
   578 next
   579   case False
   580   from assms False have irred: "irreducible (fract_poly p)" and primitive: "content p = 1"
   581     by (simp_all add: nonconst_poly_irreducible_iff)
   582   from irred have prime: "prime_elem (fract_poly p)" by (rule field_poly_irreducible_imp_prime)
   583   show ?thesis
   584   proof (rule prime_elemI)
   585     fix q r assume "p dvd q * r"
   586     hence "fract_poly p dvd fract_poly (q * r)" by (rule fract_poly_dvd)
   587     hence "fract_poly p dvd fract_poly q * fract_poly r" by simp
   588     from prime and this have "fract_poly p dvd fract_poly q \<or> fract_poly p dvd fract_poly r"
   589       by (rule prime_elem_dvd_multD)
   590     with primitive show "p dvd q \<or> p dvd r" by (auto dest: fract_poly_dvdD)
   591   qed (insert assms, auto simp: irreducible_def)
   592 qed
   593 
   594 lemma degree_primitive_part_fract [simp]:
   595   "degree (primitive_part_fract p) = degree p"
   596 proof -
   597   have "p = smult (fract_content p) (fract_poly (primitive_part_fract p))"
   598     by (simp add: content_times_primitive_part_fract)
   599   also have "degree \<dots> = degree (primitive_part_fract p)"
   600     by (auto simp: degree_map_poly)
   601   finally show ?thesis ..
   602 qed
   603 
   604 lemma irreducible_primitive_part_fract:
   605   fixes p :: "'a :: {idom_divide, ring_gcd, factorial_semiring, semiring_Gcd} fract poly"
   606   assumes "irreducible p"
   607   shows   "irreducible (primitive_part_fract p)"
   608 proof -
   609   from assms have deg: "degree (primitive_part_fract p) \<noteq> 0"
   610     by (intro notI) 
   611        (auto elim!: degree_eq_zeroE simp: irreducible_def is_unit_poly_iff dvd_field_iff)
   612   hence [simp]: "p \<noteq> 0" by auto
   613 
   614   note \<open>irreducible p\<close>
   615   also have "p = [:fract_content p:] * fract_poly (primitive_part_fract p)" 
   616     by (simp add: content_times_primitive_part_fract)
   617   also have "irreducible \<dots> \<longleftrightarrow> irreducible (fract_poly (primitive_part_fract p))"
   618     by (intro irreducible_mult_unit_left) (simp_all add: is_unit_poly_iff dvd_field_iff)
   619   finally show ?thesis using deg
   620     by (simp add: nonconst_poly_irreducible_iff)
   621 qed
   622 
   623 lemma prime_elem_primitive_part_fract:
   624   fixes p :: "'a :: {idom_divide, ring_gcd, factorial_semiring, semiring_Gcd} fract poly"
   625   shows "irreducible p \<Longrightarrow> prime_elem (primitive_part_fract p)"
   626   by (intro irreducible_imp_prime_poly irreducible_primitive_part_fract)
   627 
   628 lemma irreducible_linear_field_poly:
   629   fixes a b :: "'a::field"
   630   assumes "b \<noteq> 0"
   631   shows "irreducible [:a,b:]"
   632 proof (rule irreducibleI)
   633   fix p q assume pq: "[:a,b:] = p * q"
   634   also from pq assms have "degree \<dots> = degree p + degree q" 
   635     by (intro degree_mult_eq) auto
   636   finally have "degree p = 0 \<or> degree q = 0" using assms by auto
   637   with assms pq show "is_unit p \<or> is_unit q"
   638     by (auto simp: is_unit_const_poly_iff dvd_field_iff elim!: degree_eq_zeroE)
   639 qed (insert assms, auto simp: is_unit_poly_iff)
   640 
   641 lemma prime_elem_linear_field_poly:
   642   "(b :: 'a :: field) \<noteq> 0 \<Longrightarrow> prime_elem [:a,b:]"
   643   by (rule field_poly_irreducible_imp_prime, rule irreducible_linear_field_poly)
   644 
   645 lemma irreducible_linear_poly:
   646   fixes a b :: "'a::{idom_divide,ring_gcd,factorial_semiring,semiring_Gcd}"
   647   shows "b \<noteq> 0 \<Longrightarrow> coprime a b \<Longrightarrow> irreducible [:a,b:]"
   648   by (auto intro!: irreducible_linear_field_poly 
   649            simp:   nonconst_poly_irreducible_iff content_def map_poly_pCons)
   650 
   651 lemma prime_elem_linear_poly:
   652   fixes a b :: "'a::{idom_divide,ring_gcd,factorial_semiring,semiring_Gcd}"
   653   shows "b \<noteq> 0 \<Longrightarrow> coprime a b \<Longrightarrow> prime_elem [:a,b:]"
   654   by (rule irreducible_imp_prime_poly, rule irreducible_linear_poly)
   655 
   656  
   657 subsection \<open>Prime factorisation of polynomials\<close>   
   658 
   659 private lemma poly_prime_factorization_exists_content_1:
   660   fixes p :: "'a :: {factorial_semiring,semiring_Gcd,ring_gcd,idom_divide} poly"
   661   assumes "p \<noteq> 0" "content p = 1"
   662   shows   "\<exists>A. (\<forall>p. p \<in># A \<longrightarrow> prime_elem p) \<and> prod_mset A = normalize p"
   663 proof -
   664   let ?P = "field_poly.prime_factorization (fract_poly p)"
   665   define c where "c = prod_mset (image_mset fract_content ?P)"
   666   define c' where "c' = c * to_fract (lead_coeff p)"
   667   define e where "e = prod_mset (image_mset primitive_part_fract ?P)"
   668   define A where "A = image_mset (normalize \<circ> primitive_part_fract) ?P"
   669   have "content e = (\<Prod>x\<in>#field_poly.prime_factorization (map_poly to_fract p). 
   670                       content (primitive_part_fract x))"
   671     by (simp add: e_def content_prod_mset multiset.map_comp o_def)
   672   also have "image_mset (\<lambda>x. content (primitive_part_fract x)) ?P = image_mset (\<lambda>_. 1) ?P"
   673     by (intro image_mset_cong content_primitive_part_fract) auto
   674   finally have content_e: "content e = 1"
   675     by simp    
   676   
   677   from \<open>p \<noteq> 0\<close> have "fract_poly p = [:lead_coeff (fract_poly p):] * 
   678     smult (inverse (lead_coeff (fract_poly p))) (fract_poly p)"
   679     by simp 
   680   also have "[:lead_coeff (fract_poly p):] = [:to_fract (lead_coeff p):]" 
   681     by (simp add: monom_0 degree_map_poly coeff_map_poly)
   682   also from assms have "smult (inverse (lead_coeff (fract_poly p))) (fract_poly p) = prod_mset ?P" 
   683     by (subst field_poly_prod_mset_prime_factorization) simp_all
   684   also have "\<dots> = prod_mset (image_mset id ?P)" by simp
   685   also have "image_mset id ?P = 
   686                image_mset (\<lambda>x. [:fract_content x:] * fract_poly (primitive_part_fract x)) ?P"
   687     by (intro image_mset_cong) (auto simp: content_times_primitive_part_fract)
   688   also have "prod_mset \<dots> = smult c (fract_poly e)"
   689     by (subst prod_mset.distrib) (simp_all add: prod_mset_fract_poly prod_mset_const_poly c_def e_def)
   690   also have "[:to_fract (lead_coeff p):] * \<dots> = smult c' (fract_poly e)"
   691     by (simp add: c'_def)
   692   finally have eq: "fract_poly p = smult c' (fract_poly e)" .
   693   also obtain b where b: "c' = to_fract b" "is_unit b"
   694   proof -
   695     from fract_poly_smult_eqE[OF eq] guess a b . note ab = this
   696     from ab(2) have "content (smult a p) = content (smult b e)" by (simp only: )
   697     with assms content_e have "a = normalize b" by (simp add: ab(4))
   698     with ab have ab': "a = 1" "is_unit b" by (simp_all add: normalize_1_iff)
   699     with ab ab' have "c' = to_fract b" by auto
   700     from this and \<open>is_unit b\<close> show ?thesis by (rule that)
   701   qed
   702   hence "smult c' (fract_poly e) = fract_poly (smult b e)" by simp
   703   finally have "p = smult b e" by (simp only: fract_poly_eq_iff)
   704   hence "p = [:b:] * e" by simp
   705   with b have "normalize p = normalize e" 
   706     by (simp only: normalize_mult) (simp add: is_unit_normalize is_unit_poly_iff)
   707   also have "normalize e = prod_mset A"
   708     by (simp add: multiset.map_comp e_def A_def normalize_prod_mset)
   709   finally have "prod_mset A = normalize p" ..
   710   
   711   have "prime_elem p" if "p \<in># A" for p
   712     using that by (auto simp: A_def prime_elem_primitive_part_fract prime_elem_imp_irreducible 
   713                         dest!: field_poly_in_prime_factorization_imp_prime )
   714   from this and \<open>prod_mset A = normalize p\<close> show ?thesis
   715     by (intro exI[of _ A]) blast
   716 qed
   717 
   718 lemma poly_prime_factorization_exists:
   719   fixes p :: "'a :: {factorial_semiring,semiring_Gcd,ring_gcd,idom_divide} poly"
   720   assumes "p \<noteq> 0"
   721   shows   "\<exists>A. (\<forall>p. p \<in># A \<longrightarrow> prime_elem p) \<and> prod_mset A = normalize p"
   722 proof -
   723   define B where "B = image_mset (\<lambda>x. [:x:]) (prime_factorization (content p))"
   724   have "\<exists>A. (\<forall>p. p \<in># A \<longrightarrow> prime_elem p) \<and> prod_mset A = normalize (primitive_part p)"
   725     by (rule poly_prime_factorization_exists_content_1) (insert assms, simp_all)
   726   then guess A by (elim exE conjE) note A = this
   727   moreover from assms have "prod_mset B = [:content p:]"
   728     by (simp add: B_def prod_mset_const_poly prod_mset_prime_factorization)
   729   moreover have "\<forall>p. p \<in># B \<longrightarrow> prime_elem p"
   730     by (auto simp: B_def intro!: lift_prime_elem_poly dest: in_prime_factors_imp_prime)
   731   ultimately show ?thesis by (intro exI[of _ "B + A"]) auto
   732 qed
   733 
   734 end
   735 
   736 
   737 subsection \<open>Typeclass instances\<close>
   738 
   739 instance poly :: (factorial_ring_gcd) factorial_semiring
   740   by standard (rule poly_prime_factorization_exists)  
   741 
   742 instantiation poly :: (factorial_ring_gcd) factorial_ring_gcd
   743 begin
   744 
   745 definition gcd_poly :: "'a poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly" where
   746   [code del]: "gcd_poly = gcd_factorial"
   747 
   748 definition lcm_poly :: "'a poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly" where
   749   [code del]: "lcm_poly = lcm_factorial"
   750   
   751 definition Gcd_poly :: "'a poly set \<Rightarrow> 'a poly" where
   752  [code del]: "Gcd_poly = Gcd_factorial"
   753 
   754 definition Lcm_poly :: "'a poly set \<Rightarrow> 'a poly" where
   755  [code del]: "Lcm_poly = Lcm_factorial"
   756  
   757 instance by standard (simp_all add: gcd_poly_def lcm_poly_def Gcd_poly_def Lcm_poly_def)
   758 
   759 end
   760 
   761 instantiation poly :: ("{field,factorial_ring_gcd}") unique_euclidean_ring
   762 begin
   763 
   764 definition euclidean_size_poly :: "'a poly \<Rightarrow> nat"
   765   where "euclidean_size_poly p = (if p = 0 then 0 else 2 ^ degree p)"
   766 
   767 definition uniqueness_constraint_poly :: "'a poly \<Rightarrow> 'a poly \<Rightarrow> bool"
   768   where [simp]: "uniqueness_constraint_poly = top"
   769 
   770 instance proof
   771   show "(q * p + r) div p = q" if "p \<noteq> 0"
   772     and "euclidean_size r < euclidean_size p" for q p r :: "'a poly"
   773   proof -
   774     from \<open>p \<noteq> 0\<close> eucl_rel_poly [of r p] have "eucl_rel_poly (r + q * p) p (q + r div p, r mod p)"
   775       by (simp add: eucl_rel_poly_iff distrib_right)
   776     then have "(r + q * p) div p = q + r div p"
   777       by (rule div_poly_eq)
   778     with that show ?thesis
   779       by (cases "r = 0") (simp_all add: euclidean_size_poly_def div_poly_less ac_simps)
   780   qed
   781 qed (auto simp: euclidean_size_poly_def Rings.div_mult_mod_eq div_poly_less degree_mult_eq intro!: degree_mod_less' degree_mult_right_le
   782     split: if_splits)
   783 
   784 end
   785 
   786 instance poly :: ("{field,factorial_ring_gcd}") euclidean_ring_gcd
   787   by (rule euclidean_ring_gcd_class.intro, rule factorial_euclidean_semiring_gcdI)
   788     standard
   789 
   790   
   791 subsection \<open>Polynomial GCD\<close>
   792 
   793 lemma gcd_poly_decompose:
   794   fixes p q :: "'a :: factorial_ring_gcd poly"
   795   shows "gcd p q = 
   796            smult (gcd (content p) (content q)) (gcd (primitive_part p) (primitive_part q))"
   797 proof (rule sym, rule gcdI)
   798   have "[:gcd (content p) (content q):] * gcd (primitive_part p) (primitive_part q) dvd
   799           [:content p:] * primitive_part p" by (intro mult_dvd_mono) simp_all
   800   thus "smult (gcd (content p) (content q)) (gcd (primitive_part p) (primitive_part q)) dvd p"
   801     by simp
   802 next
   803   have "[:gcd (content p) (content q):] * gcd (primitive_part p) (primitive_part q) dvd
   804           [:content q:] * primitive_part q" by (intro mult_dvd_mono) simp_all
   805   thus "smult (gcd (content p) (content q)) (gcd (primitive_part p) (primitive_part q)) dvd q"
   806     by simp
   807 next
   808   fix d assume "d dvd p" "d dvd q"
   809   hence "[:content d:] * primitive_part d dvd 
   810            [:gcd (content p) (content q):] * gcd (primitive_part p) (primitive_part q)"
   811     by (intro mult_dvd_mono) auto
   812   thus "d dvd smult (gcd (content p) (content q)) (gcd (primitive_part p) (primitive_part q))"
   813     by simp
   814 qed (auto simp: normalize_smult)
   815   
   816 
   817 lemma gcd_poly_pseudo_mod:
   818   fixes p q :: "'a :: factorial_ring_gcd poly"
   819   assumes nz: "q \<noteq> 0" and prim: "content p = 1" "content q = 1"
   820   shows   "gcd p q = gcd q (primitive_part (pseudo_mod p q))"
   821 proof -
   822   define r s where "r = fst (pseudo_divmod p q)" and "s = snd (pseudo_divmod p q)"
   823   define a where "a = [:coeff q (degree q) ^ (Suc (degree p) - degree q):]"
   824   have [simp]: "primitive_part a = unit_factor a"
   825     by (simp add: a_def unit_factor_poly_def unit_factor_power monom_0)
   826   from nz have [simp]: "a \<noteq> 0" by (auto simp: a_def)
   827   
   828   have rs: "pseudo_divmod p q = (r, s)" by (simp add: r_def s_def)
   829   have "gcd (q * r + s) q = gcd q s"
   830     using gcd_add_mult[of q r s] by (simp add: gcd.commute add_ac mult_ac)
   831   with pseudo_divmod(1)[OF nz rs]
   832     have "gcd (p * a) q = gcd q s" by (simp add: a_def)
   833   also from prim have "gcd (p * a) q = gcd p q"
   834     by (subst gcd_poly_decompose)
   835        (auto simp: primitive_part_mult gcd_mult_unit1 primitive_part_prim 
   836              simp del: mult_pCons_right )
   837   also from prim have "gcd q s = gcd q (primitive_part s)"
   838     by (subst gcd_poly_decompose) (simp_all add: primitive_part_prim)
   839   also have "s = pseudo_mod p q" by (simp add: s_def pseudo_mod_def)
   840   finally show ?thesis .
   841 qed
   842 
   843 lemma degree_pseudo_mod_less:
   844   assumes "q \<noteq> 0" "pseudo_mod p q \<noteq> 0"
   845   shows   "degree (pseudo_mod p q) < degree q"
   846   using pseudo_mod(2)[of q p] assms by auto
   847 
   848 function gcd_poly_code_aux :: "'a :: factorial_ring_gcd poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly" where
   849   "gcd_poly_code_aux p q = 
   850      (if q = 0 then normalize p else gcd_poly_code_aux q (primitive_part (pseudo_mod p q)))" 
   851 by auto
   852 termination
   853   by (relation "measure ((\<lambda>p. if p = 0 then 0 else Suc (degree p)) \<circ> snd)")
   854      (auto simp: degree_pseudo_mod_less)
   855 
   856 declare gcd_poly_code_aux.simps [simp del]
   857 
   858 lemma gcd_poly_code_aux_correct:
   859   assumes "content p = 1" "q = 0 \<or> content q = 1"
   860   shows   "gcd_poly_code_aux p q = gcd p q"
   861   using assms
   862 proof (induction p q rule: gcd_poly_code_aux.induct)
   863   case (1 p q)
   864   show ?case
   865   proof (cases "q = 0")
   866     case True
   867     thus ?thesis by (subst gcd_poly_code_aux.simps) auto
   868   next
   869     case False
   870     hence "gcd_poly_code_aux p q = gcd_poly_code_aux q (primitive_part (pseudo_mod p q))"
   871       by (subst gcd_poly_code_aux.simps) simp_all
   872     also from "1.prems" False 
   873       have "primitive_part (pseudo_mod p q) = 0 \<or> 
   874               content (primitive_part (pseudo_mod p q)) = 1"
   875       by (cases "pseudo_mod p q = 0") auto
   876     with "1.prems" False 
   877       have "gcd_poly_code_aux q (primitive_part (pseudo_mod p q)) = 
   878               gcd q (primitive_part (pseudo_mod p q))"
   879       by (intro 1) simp_all
   880     also from "1.prems" False 
   881       have "\<dots> = gcd p q" by (intro gcd_poly_pseudo_mod [symmetric]) auto
   882     finally show ?thesis .
   883   qed
   884 qed
   885 
   886 definition gcd_poly_code 
   887     :: "'a :: factorial_ring_gcd poly \<Rightarrow> 'a poly \<Rightarrow> 'a poly" 
   888   where "gcd_poly_code p q = 
   889            (if p = 0 then normalize q else if q = 0 then normalize p else
   890               smult (gcd (content p) (content q)) 
   891                 (gcd_poly_code_aux (primitive_part p) (primitive_part q)))"
   892 
   893 lemma gcd_poly_code [code]: "gcd p q = gcd_poly_code p q"
   894   by (simp add: gcd_poly_code_def gcd_poly_code_aux_correct gcd_poly_decompose [symmetric])
   895 
   896 lemma lcm_poly_code [code]: 
   897   fixes p q :: "'a :: factorial_ring_gcd poly"
   898   shows "lcm p q = normalize (p * q) div gcd p q"
   899   by (fact lcm_gcd)
   900 
   901 lemmas Gcd_poly_set_eq_fold [code] = Gcd_set_eq_fold [where ?'a = "'a :: factorial_ring_gcd poly"]
   902 lemmas Lcm_poly_set_eq_fold [code] = Lcm_set_eq_fold [where ?'a = "'a :: factorial_ring_gcd poly"]
   903 
   904 text \<open>Example:
   905   @{lemma "Lcm {[:1, 2, 3:], [:2, 3, 4:]} = [:[:2:], [:7:], [:16:], [:17:], [:12 :: int:]:]" by eval}
   906 \<close>
   907   
   908 end