src/HOL/Library/Fundamental_Theorem_Algebra.thy
author wenzelm
Mon Apr 28 23:43:13 2014 +0200 (2014-04-28)
changeset 56778 cb0929421ca6
parent 56776 309e1a61ee7c
child 56795 e8cce2bd23e5
permissions -rw-r--r--
tuned proofs;
     1 (* Author: Amine Chaieb, TU Muenchen *)
     2 
     3 header{*Fundamental Theorem of Algebra*}
     4 
     5 theory Fundamental_Theorem_Algebra
     6 imports Polynomial Complex_Main
     7 begin
     8 
     9 subsection {* Square root of complex numbers *}
    10 
    11 definition csqrt :: "complex \<Rightarrow> complex"
    12 where
    13   "csqrt z =
    14     (if Im z = 0 then
    15        if 0 \<le> Re z then Complex (sqrt(Re z)) 0
    16        else Complex 0 (sqrt(- Re z))
    17      else Complex (sqrt((cmod z + Re z) /2)) ((Im z / abs(Im z)) * sqrt((cmod z - Re z) /2)))"
    18 
    19 lemma csqrt[algebra]: "(csqrt z)\<^sup>2 = z"
    20 proof -
    21   obtain x y where xy: "z = Complex x y" by (cases z)
    22   {
    23     assume y0: "y = 0"
    24     {
    25       assume x0: "x \<ge> 0"
    26       then have ?thesis
    27         using y0 xy real_sqrt_pow2[OF x0]
    28         by (simp add: csqrt_def power2_eq_square)
    29     }
    30     moreover
    31     {
    32       assume "\<not> x \<ge> 0"
    33       then have x0: "- x \<ge> 0" by arith
    34       then have ?thesis
    35         using y0 xy real_sqrt_pow2[OF x0]
    36         by (simp add: csqrt_def power2_eq_square)
    37     }
    38     ultimately have ?thesis by blast
    39   }
    40   moreover
    41   {
    42     assume y0: "y \<noteq> 0"
    43     {
    44       fix x y
    45       let ?z = "Complex x y"
    46       from abs_Re_le_cmod[of ?z] have tha: "abs x \<le> cmod ?z"
    47         by auto
    48       then have "cmod ?z - x \<ge> 0" "cmod ?z + x \<ge> 0"
    49         by arith+
    50       then have "(sqrt (x * x + y * y) + x) / 2 \<ge> 0" "(sqrt (x * x + y * y) - x) / 2 \<ge> 0"
    51         by (simp_all add: power2_eq_square)
    52     }
    53     note th = this
    54     have sq4: "\<And>x::real. x\<^sup>2 / 4 = (x / 2)\<^sup>2"
    55       by (simp add: power2_eq_square)
    56     from th[of x y]
    57     have sq4': "sqrt (((sqrt (x * x + y * y) + x)\<^sup>2 / 4)) = (sqrt (x * x + y * y) + x) / 2"
    58       "sqrt (((sqrt (x * x + y * y) - x)\<^sup>2 / 4)) = (sqrt (x * x + y * y) - x) / 2"
    59       unfolding sq4 by simp_all
    60     then have th1: "sqrt ((sqrt (x * x + y * y) + x) * (sqrt (x * x + y * y) + x) / 4) -
    61         sqrt ((sqrt (x * x + y * y) - x) * (sqrt (x * x + y * y) - x) / 4) = x"
    62       unfolding power2_eq_square by simp
    63     have "sqrt 4 = sqrt (2\<^sup>2)"
    64       by simp
    65     then have sqrt4: "sqrt 4 = 2"
    66       by (simp only: real_sqrt_abs)
    67     have th2: "2 * (y * sqrt ((sqrt (x * x + y * y) - x) * (sqrt (x * x + y * y) + x) / 4)) / \<bar>y\<bar> = y"
    68       using iffD2[OF real_sqrt_pow2_iff sum_power2_ge_zero[of x y]] y0
    69       unfolding power2_eq_square
    70       by (simp add: algebra_simps real_sqrt_divide sqrt4)
    71     from y0 xy have ?thesis
    72       apply (simp add: csqrt_def power2_eq_square)
    73       apply (simp add: real_sqrt_sum_squares_mult_ge_zero[of x y]
    74         real_sqrt_pow2[OF th(1)[of x y], unfolded power2_eq_square]
    75         real_sqrt_pow2[OF th(2)[of x y], unfolded power2_eq_square]
    76         real_sqrt_mult[symmetric])
    77       using th1 th2  ..
    78   }
    79   ultimately show ?thesis by blast
    80 qed
    81 
    82 lemma csqrt_Complex: "x \<ge> 0 \<Longrightarrow> csqrt (Complex x 0) = Complex (sqrt x) 0"
    83   by (simp add: csqrt_def)
    84 
    85 lemma csqrt_0 [simp]: "csqrt 0 = 0"
    86   by (simp add: csqrt_def)
    87 
    88 lemma csqrt_1 [simp]: "csqrt 1 = 1"
    89   by (simp add: csqrt_def)
    90 
    91 lemma csqrt_principal: "0 < Re(csqrt(z)) | Re(csqrt(z)) = 0 & 0 \<le> Im(csqrt(z))"
    92 proof (cases z)
    93   case (Complex x y)
    94   then show ?thesis
    95     using real_sqrt_sum_squares_ge1 [of "x" y]
    96           real_sqrt_sum_squares_ge1 [of "-x" y]
    97           real_sqrt_sum_squares_eq_cancel [of x y]
    98     apply (auto simp: csqrt_def intro!: Rings.ordered_ring_class.split_mult_pos_le)
    99     apply (metis add_commute diff_add_cancel le_add_same_cancel1 real_sqrt_sum_squares_ge1)
   100     apply (metis add_commute less_eq_real_def power_minus_Bit0
   101             real_0_less_add_iff real_sqrt_sum_squares_eq_cancel)
   102     done
   103 qed
   104 
   105 lemma Re_csqrt: "0 \<le> Re(csqrt z)"
   106   by (metis csqrt_principal le_less)
   107 
   108 lemma csqrt_square: "0 < Re z \<or> Re z = 0 \<and> 0 \<le> Im z \<Longrightarrow> csqrt (z\<^sup>2) = z"
   109   using csqrt [of "z\<^sup>2"] csqrt_principal [of "z\<^sup>2"]
   110   by (cases z) (auto simp: power2_eq_iff)
   111 
   112 lemma csqrt_eq_0 [simp]: "csqrt z = 0 \<longleftrightarrow> z = 0"
   113   by auto (metis csqrt power_eq_0_iff)
   114 
   115 lemma csqrt_eq_1 [simp]: "csqrt z = 1 \<longleftrightarrow> z = 1"
   116   by auto (metis csqrt power2_eq_1_iff)
   117 
   118 
   119 subsection {* More lemmas about module of complex numbers *}
   120 
   121 lemma complex_of_real_power: "complex_of_real x ^ n = complex_of_real (x^n)"
   122   by (rule of_real_power [symmetric])
   123 
   124 text{* The triangle inequality for cmod *}
   125 lemma complex_mod_triangle_sub: "cmod w \<le> cmod (w + z) + norm z"
   126   using complex_mod_triangle_ineq2[of "w + z" "-z"] by auto
   127 
   128 
   129 subsection {* Basic lemmas about polynomials *}
   130 
   131 lemma poly_bound_exists:
   132   fixes p :: "'a::{comm_semiring_0,real_normed_div_algebra} poly"
   133   shows "\<exists>m. m > 0 \<and> (\<forall>z. norm z \<le> r \<longrightarrow> norm (poly p z) \<le> m)"
   134 proof (induct p)
   135   case 0
   136   then show ?case by (rule exI[where x=1]) simp
   137 next
   138   case (pCons c cs)
   139   from pCons.hyps obtain m where m: "\<forall>z. norm z \<le> r \<longrightarrow> norm (poly cs z) \<le> m"
   140     by blast
   141   let ?k = " 1 + norm c + \<bar>r * m\<bar>"
   142   have kp: "?k > 0" using abs_ge_zero[of "r*m"] norm_ge_zero[of c] by arith
   143   {
   144     fix z :: 'a
   145     assume H: "norm z \<le> r"
   146     from m H have th: "norm (poly cs z) \<le> m"
   147       by blast
   148     from H have rp: "r \<ge> 0" using norm_ge_zero[of z]
   149       by arith
   150     have "norm (poly (pCons c cs) z) \<le> norm c + norm (z* poly cs z)"
   151       using norm_triangle_ineq[of c "z* poly cs z"] by simp
   152     also have "\<dots> \<le> norm c + r * m"
   153       using mult_mono[OF H th rp norm_ge_zero[of "poly cs z"]]
   154       by (simp add: norm_mult)
   155     also have "\<dots> \<le> ?k"
   156       by simp
   157     finally have "norm (poly (pCons c cs) z) \<le> ?k" .
   158   }
   159   with kp show ?case by blast
   160 qed
   161 
   162 
   163 text{* Offsetting the variable in a polynomial gives another of same degree *}
   164 
   165 definition offset_poly :: "'a::comm_semiring_0 poly \<Rightarrow> 'a \<Rightarrow> 'a poly"
   166   where "offset_poly p h = fold_coeffs (\<lambda>a q. smult h q + pCons a q) p 0"
   167 
   168 lemma offset_poly_0: "offset_poly 0 h = 0"
   169   by (simp add: offset_poly_def)
   170 
   171 lemma offset_poly_pCons:
   172   "offset_poly (pCons a p) h =
   173     smult h (offset_poly p h) + pCons a (offset_poly p h)"
   174   by (cases "p = 0 \<and> a = 0") (auto simp add: offset_poly_def)
   175 
   176 lemma offset_poly_single: "offset_poly [:a:] h = [:a:]"
   177   by (simp add: offset_poly_pCons offset_poly_0)
   178 
   179 lemma poly_offset_poly: "poly (offset_poly p h) x = poly p (h + x)"
   180   apply (induct p)
   181   apply (simp add: offset_poly_0)
   182   apply (simp add: offset_poly_pCons algebra_simps)
   183   done
   184 
   185 lemma offset_poly_eq_0_lemma: "smult c p + pCons a p = 0 \<Longrightarrow> p = 0"
   186   by (induct p arbitrary: a) (simp, force)
   187 
   188 lemma offset_poly_eq_0_iff: "offset_poly p h = 0 \<longleftrightarrow> p = 0"
   189   apply (safe intro!: offset_poly_0)
   190   apply (induct p, simp)
   191   apply (simp add: offset_poly_pCons)
   192   apply (frule offset_poly_eq_0_lemma, simp)
   193   done
   194 
   195 lemma degree_offset_poly: "degree (offset_poly p h) = degree p"
   196   apply (induct p)
   197   apply (simp add: offset_poly_0)
   198   apply (case_tac "p = 0")
   199   apply (simp add: offset_poly_0 offset_poly_pCons)
   200   apply (simp add: offset_poly_pCons)
   201   apply (subst degree_add_eq_right)
   202   apply (rule le_less_trans [OF degree_smult_le])
   203   apply (simp add: offset_poly_eq_0_iff)
   204   apply (simp add: offset_poly_eq_0_iff)
   205   done
   206 
   207 definition "psize p = (if p = 0 then 0 else Suc (degree p))"
   208 
   209 lemma psize_eq_0_iff [simp]: "psize p = 0 \<longleftrightarrow> p = 0"
   210   unfolding psize_def by simp
   211 
   212 lemma poly_offset:
   213   fixes p :: "'a::comm_ring_1 poly"
   214   shows "\<exists>q. psize q = psize p \<and> (\<forall>x. poly q x = poly p (a + x))"
   215 proof (intro exI conjI)
   216   show "psize (offset_poly p a) = psize p"
   217     unfolding psize_def
   218     by (simp add: offset_poly_eq_0_iff degree_offset_poly)
   219   show "\<forall>x. poly (offset_poly p a) x = poly p (a + x)"
   220     by (simp add: poly_offset_poly)
   221 qed
   222 
   223 text{* An alternative useful formulation of completeness of the reals *}
   224 lemma real_sup_exists:
   225   assumes ex: "\<exists>x. P x"
   226     and bz: "\<exists>z. \<forall>x. P x \<longrightarrow> x < z"
   227   shows "\<exists>s::real. \<forall>y. (\<exists>x. P x \<and> y < x) \<longleftrightarrow> y < s"
   228 proof
   229   from bz have "bdd_above (Collect P)"
   230     by (force intro: less_imp_le)
   231   then show "\<forall>y. (\<exists>x. P x \<and> y < x) \<longleftrightarrow> y < Sup (Collect P)"
   232     using ex bz by (subst less_cSup_iff) auto
   233 qed
   234 
   235 subsection {* Fundamental theorem of algebra *}
   236 lemma  unimodular_reduce_norm:
   237   assumes md: "cmod z = 1"
   238   shows "cmod (z + 1) < 1 \<or> cmod (z - 1) < 1 \<or> cmod (z + ii) < 1 \<or> cmod (z - ii) < 1"
   239 proof -
   240   obtain x y where z: "z = Complex x y "
   241     by (cases z) auto
   242   from md z have xy: "x\<^sup>2 + y\<^sup>2 = 1"
   243     by (simp add: cmod_def)
   244   {
   245     assume C: "cmod (z + 1) \<ge> 1" "cmod (z - 1) \<ge> 1" "cmod (z + ii) \<ge> 1" "cmod (z - ii) \<ge> 1"
   246     from C z xy have "2 * x \<le> 1" "2 * x \<ge> -1" "2 * y \<le> 1" "2 * y \<ge> -1"
   247       by (simp_all add: cmod_def power2_eq_square algebra_simps)
   248     then have "abs (2 * x) \<le> 1" "abs (2 * y) \<le> 1"
   249       by simp_all
   250     then have "(abs (2 * x))\<^sup>2 \<le> 1\<^sup>2" "(abs (2 * y))\<^sup>2 \<le> 1\<^sup>2"
   251       by - (rule power_mono, simp, simp)+
   252     then have th0: "4 * x\<^sup>2 \<le> 1" "4 * y\<^sup>2 \<le> 1"
   253       by (simp_all add: power_mult_distrib)
   254     from add_mono[OF th0] xy have False by simp
   255   }
   256   then show ?thesis
   257     unfolding linorder_not_le[symmetric] by blast
   258 qed
   259 
   260 text{* Hence we can always reduce modulus of @{text "1 + b z^n"} if nonzero *}
   261 lemma reduce_poly_simple:
   262   assumes b: "b \<noteq> 0"
   263     and n: "n \<noteq> 0"
   264   shows "\<exists>z. cmod (1 + b * z^n) < 1"
   265   using n
   266 proof (induct n rule: nat_less_induct)
   267   fix n
   268   assume IH: "\<forall>m<n. m \<noteq> 0 \<longrightarrow> (\<exists>z. cmod (1 + b * z ^ m) < 1)"
   269   assume n: "n \<noteq> 0"
   270   let ?P = "\<lambda>z n. cmod (1 + b * z ^ n) < 1"
   271   {
   272     assume e: "even n"
   273     then have "\<exists>m. n = 2 * m"
   274       by presburger
   275     then obtain m where m: "n = 2 * m"
   276       by blast
   277     from n m have "m \<noteq> 0" "m < n"
   278       by presburger+
   279     with IH[rule_format, of m] obtain z where z: "?P z m"
   280       by blast
   281     from z have "?P (csqrt z) n" by (simp add: m power_mult csqrt)
   282     then have "\<exists>z. ?P z n" ..
   283   }
   284   moreover
   285   {
   286     assume o: "odd n"
   287     have th0: "cmod (complex_of_real (cmod b) / b) = 1"
   288       using b by (simp add: norm_divide)
   289     from o have "\<exists>m. n = Suc (2 * m)"
   290       by presburger+
   291     then obtain m where m: "n = Suc (2*m)"
   292       by blast
   293     from unimodular_reduce_norm[OF th0] o
   294     have "\<exists>v. cmod (complex_of_real (cmod b) / b + v^n) < 1"
   295       apply (cases "cmod (complex_of_real (cmod b) / b + 1) < 1", rule_tac x="1" in exI, simp)
   296       apply (cases "cmod (complex_of_real (cmod b) / b - 1) < 1", rule_tac x="-1" in exI, simp)
   297       apply (cases "cmod (complex_of_real (cmod b) / b + ii) < 1")
   298       apply (cases "even m", rule_tac x="ii" in exI, simp add: m power_mult)
   299       apply (rule_tac x="- ii" in exI, simp add: m power_mult)
   300       apply (cases "even m", rule_tac x="- ii" in exI, simp add: m power_mult)
   301       apply (auto simp add: m power_mult)
   302       apply (rule_tac x="ii" in exI)
   303       apply (auto simp add: m power_mult)
   304       done
   305     then obtain v where v: "cmod (complex_of_real (cmod b) / b + v^n) < 1"
   306       by blast
   307     let ?w = "v / complex_of_real (root n (cmod b))"
   308     from odd_real_root_pow[OF o, of "cmod b"]
   309     have th1: "?w ^ n = v^n / complex_of_real (cmod b)"
   310       by (simp add: power_divide complex_of_real_power)
   311     have th2:"cmod (complex_of_real (cmod b) / b) = 1"
   312       using b by (simp add: norm_divide)
   313     then have th3: "cmod (complex_of_real (cmod b) / b) \<ge> 0"
   314       by simp
   315     have th4: "cmod (complex_of_real (cmod b) / b) *
   316         cmod (1 + b * (v ^ n / complex_of_real (cmod b))) <
   317         cmod (complex_of_real (cmod b) / b) * 1"
   318       apply (simp only: norm_mult[symmetric] distrib_left)
   319       using b v
   320       apply (simp add: th2)
   321       done
   322     from mult_less_imp_less_left[OF th4 th3]
   323     have "?P ?w n" unfolding th1 .
   324     then have "\<exists>z. ?P z n" ..
   325   }
   326   ultimately show "\<exists>z. ?P z n" by blast
   327 qed
   328 
   329 text{* Bolzano-Weierstrass type property for closed disc in complex plane. *}
   330 
   331 lemma metric_bound_lemma: "cmod (x - y) \<le> \<bar>Re x - Re y\<bar> + \<bar>Im x - Im y\<bar>"
   332   using real_sqrt_sum_squares_triangle_ineq[of "Re x - Re y" 0 0 "Im x - Im y" ]
   333   unfolding cmod_def by simp
   334 
   335 lemma bolzano_weierstrass_complex_disc:
   336   assumes r: "\<forall>n. cmod (s n) \<le> r"
   337   shows "\<exists>f z. subseq f \<and> (\<forall>e >0. \<exists>N. \<forall>n \<ge> N. cmod (s (f n) - z) < e)"
   338 proof-
   339   from seq_monosub[of "Re \<circ> s"]
   340   obtain f where f: "subseq f" "monoseq (\<lambda>n. Re (s (f n)))"
   341     unfolding o_def by blast
   342   from seq_monosub[of "Im \<circ> s \<circ> f"]
   343   obtain g where g: "subseq g" "monoseq (\<lambda>n. Im (s (f (g n))))"
   344     unfolding o_def by blast
   345   let ?h = "f \<circ> g"
   346   from r[rule_format, of 0] have rp: "r \<ge> 0"
   347     using norm_ge_zero[of "s 0"] by arith
   348   have th: "\<forall>n. r + 1 \<ge> \<bar>Re (s n)\<bar>"
   349   proof
   350     fix n
   351     from abs_Re_le_cmod[of "s n"] r[rule_format, of n]
   352     show "\<bar>Re (s n)\<bar> \<le> r + 1" by arith
   353   qed
   354   have conv1: "convergent (\<lambda>n. Re (s (f n)))"
   355     apply (rule Bseq_monoseq_convergent)
   356     apply (simp add: Bseq_def)
   357     apply (metis gt_ex le_less_linear less_trans order.trans th)
   358     apply (rule f(2))
   359     done
   360   have th: "\<forall>n. r + 1 \<ge> \<bar>Im (s n)\<bar>"
   361   proof
   362     fix n
   363     from abs_Im_le_cmod[of "s n"] r[rule_format, of n]
   364     show "\<bar>Im (s n)\<bar> \<le> r + 1"
   365       by arith
   366   qed
   367 
   368   have conv2: "convergent (\<lambda>n. Im (s (f (g n))))"
   369     apply (rule Bseq_monoseq_convergent)
   370     apply (simp add: Bseq_def)
   371     apply (metis gt_ex le_less_linear less_trans order.trans th)
   372     apply (rule g(2))
   373     done
   374 
   375   from conv1[unfolded convergent_def] obtain x where "LIMSEQ (\<lambda>n. Re (s (f n))) x"
   376     by blast
   377   then have x: "\<forall>r>0. \<exists>n0. \<forall>n\<ge>n0. \<bar> Re (s (f n)) - x \<bar> < r"
   378     unfolding LIMSEQ_iff real_norm_def .
   379 
   380   from conv2[unfolded convergent_def] obtain y where "LIMSEQ (\<lambda>n. Im (s (f (g n)))) y"
   381     by blast
   382   then have y: "\<forall>r>0. \<exists>n0. \<forall>n\<ge>n0. \<bar> Im (s (f (g n))) - y \<bar> < r"
   383     unfolding LIMSEQ_iff real_norm_def .
   384   let ?w = "Complex x y"
   385   from f(1) g(1) have hs: "subseq ?h"
   386     unfolding subseq_def by auto
   387   {
   388     fix e :: real
   389     assume ep: "e > 0"
   390     then have e2: "e/2 > 0" by simp
   391     from x[rule_format, OF e2] y[rule_format, OF e2]
   392     obtain N1 N2 where N1: "\<forall>n\<ge>N1. \<bar>Re (s (f n)) - x\<bar> < e / 2"
   393       and N2: "\<forall>n\<ge>N2. \<bar>Im (s (f (g n))) - y\<bar> < e / 2" by blast
   394     {
   395       fix n
   396       assume nN12: "n \<ge> N1 + N2"
   397       then have nN1: "g n \<ge> N1" and nN2: "n \<ge> N2"
   398         using seq_suble[OF g(1), of n] by arith+
   399       from add_strict_mono[OF N1[rule_format, OF nN1] N2[rule_format, OF nN2]]
   400       have "cmod (s (?h n) - ?w) < e"
   401         using metric_bound_lemma[of "s (f (g n))" ?w] by simp
   402     }
   403     then have "\<exists>N. \<forall>n\<ge>N. cmod (s (?h n) - ?w) < e" by blast
   404   }
   405   with hs show ?thesis by blast
   406 qed
   407 
   408 text{* Polynomial is continuous. *}
   409 
   410 lemma poly_cont:
   411   fixes p :: "'a::{comm_semiring_0,real_normed_div_algebra} poly"
   412   assumes ep: "e > 0"
   413   shows "\<exists>d >0. \<forall>w. 0 < norm (w - z) \<and> norm (w - z) < d \<longrightarrow> norm (poly p w - poly p z) < e"
   414 proof -
   415   obtain q where q: "degree q = degree p" "\<And>x. poly q x = poly p (z + x)"
   416   proof
   417     show "degree (offset_poly p z) = degree p"
   418       by (rule degree_offset_poly)
   419     show "\<And>x. poly (offset_poly p z) x = poly p (z + x)"
   420       by (rule poly_offset_poly)
   421   qed
   422   have th: "\<And>w. poly q (w - z) = poly p w"
   423     using q(2)[of "w - z" for w] by simp
   424   show ?thesis unfolding th[symmetric]
   425   proof (induct q)
   426     case 0
   427     then show ?case
   428       using ep by auto
   429   next
   430     case (pCons c cs)
   431     from poly_bound_exists[of 1 "cs"]
   432     obtain m where m: "m > 0" "\<And>z. norm z \<le> 1 \<Longrightarrow> norm (poly cs z) \<le> m"
   433       by blast
   434     from ep m(1) have em0: "e/m > 0"
   435       by (simp add: field_simps)
   436     have one0: "1 > (0::real)"
   437       by arith
   438     from real_lbound_gt_zero[OF one0 em0]
   439     obtain d where d: "d > 0" "d < 1" "d < e / m"
   440       by blast
   441     from d(1,3) m(1) have dm: "d * m > 0" "d * m < e"
   442       by (simp_all add: field_simps)
   443     show ?case
   444     proof (rule ex_forward[OF real_lbound_gt_zero[OF one0 em0]], clarsimp simp add: norm_mult)
   445       fix d w
   446       assume H: "d > 0" "d < 1" "d < e/m" "w \<noteq> z" "norm (w - z) < d"
   447       then have d1: "norm (w-z) \<le> 1" "d \<ge> 0"
   448         by simp_all
   449       from H(3) m(1) have dme: "d*m < e"
   450         by (simp add: field_simps)
   451       from H have th: "norm (w - z) \<le> d"
   452         by simp
   453       from mult_mono[OF th m(2)[OF d1(1)] d1(2) norm_ge_zero] dme
   454       show "norm (w - z) * norm (poly cs (w - z)) < e"
   455         by simp
   456     qed
   457   qed
   458 qed
   459 
   460 text{* Hence a polynomial attains minimum on a closed disc
   461   in the complex plane. *}
   462 lemma poly_minimum_modulus_disc: "\<exists>z. \<forall>w. cmod w \<le> r \<longrightarrow> cmod (poly p z) \<le> cmod (poly p w)"
   463 proof -
   464   {
   465     assume "\<not> r \<ge> 0"
   466     then have ?thesis
   467       by (metis norm_ge_zero order.trans)
   468   }
   469   moreover
   470   {
   471     assume rp: "r \<ge> 0"
   472     from rp have "cmod 0 \<le> r \<and> cmod (poly p 0) = - (- cmod (poly p 0))"
   473       by simp
   474     then have mth1: "\<exists>x z. cmod z \<le> r \<and> cmod (poly p z) = - x"
   475       by blast
   476     {
   477       fix x z
   478       assume H: "cmod z \<le> r" "cmod (poly p z) = - x" "\<not> x < 1"
   479       then have "- x < 0 "
   480         by arith
   481       with H(2) norm_ge_zero[of "poly p z"] have False
   482         by simp
   483     }
   484     then have mth2: "\<exists>z. \<forall>x. (\<exists>z. cmod z \<le> r \<and> cmod (poly p z) = - x) \<longrightarrow> x < z"
   485       by blast
   486     from real_sup_exists[OF mth1 mth2] obtain s where
   487       s: "\<forall>y. (\<exists>x. (\<exists>z. cmod z \<le> r \<and> cmod (poly p z) = - x) \<and> y < x) \<longleftrightarrow> y < s" by blast
   488     let ?m = "- s"
   489     {
   490       fix y
   491       from s[rule_format, of "-y"]
   492       have "(\<exists>z x. cmod z \<le> r \<and> - (- cmod (poly p z)) < y) \<longleftrightarrow> ?m < y"
   493         unfolding minus_less_iff[of y ] equation_minus_iff by blast
   494     }
   495     note s1 = this[unfolded minus_minus]
   496     from s1[of ?m] have s1m: "\<And>z x. cmod z \<le> r \<Longrightarrow> cmod (poly p z) \<ge> ?m"
   497       by auto
   498     {
   499       fix n :: nat
   500       from s1[rule_format, of "?m + 1/real (Suc n)"]
   501       have "\<exists>z. cmod z \<le> r \<and> cmod (poly p z) < - s + 1 / real (Suc n)"
   502         by simp
   503     }
   504     then have th: "\<forall>n. \<exists>z. cmod z \<le> r \<and> cmod (poly p z) < - s + 1 / real (Suc n)" ..
   505     from choice[OF th] obtain g where
   506         g: "\<forall>n. cmod (g n) \<le> r" "\<forall>n. cmod (poly p (g n)) <?m + 1 /real(Suc n)"
   507       by blast
   508     from bolzano_weierstrass_complex_disc[OF g(1)]
   509     obtain f z where fz: "subseq f" "\<forall>e>0. \<exists>N. \<forall>n\<ge>N. cmod (g (f n) - z) < e"
   510       by blast
   511     {
   512       fix w
   513       assume wr: "cmod w \<le> r"
   514       let ?e = "\<bar>cmod (poly p z) - ?m\<bar>"
   515       {
   516         assume e: "?e > 0"
   517         then have e2: "?e/2 > 0" by simp
   518         from poly_cont[OF e2, of z p] obtain d where
   519             d: "d > 0" "\<forall>w. 0<cmod (w - z)\<and> cmod(w - z) < d \<longrightarrow> cmod(poly p w - poly p z) < ?e/2"
   520           by blast
   521         {
   522           fix w
   523           assume w: "cmod (w - z) < d"
   524           have "cmod(poly p w - poly p z) < ?e / 2"
   525             using d(2)[rule_format, of w] w e by (cases "w = z") simp_all
   526         }
   527         note th1 = this
   528 
   529         from fz(2) d(1) obtain N1 where N1: "\<forall>n\<ge>N1. cmod (g (f n) - z) < d"
   530           by blast
   531         from reals_Archimedean2[of "2/?e"] obtain N2 :: nat where N2: "2/?e < real N2"
   532           by blast
   533         have th2: "cmod (poly p (g (f (N1 + N2))) - poly p z) < ?e/2"
   534           using N1[rule_format, of "N1 + N2"] th1 by simp
   535         {
   536           fix a b e2 m :: real
   537           have "a < e2 \<Longrightarrow> \<bar>b - m\<bar> < e2 \<Longrightarrow> 2 * e2 \<le> \<bar>b - m\<bar> + a \<Longrightarrow> False"
   538             by arith
   539         }
   540         note th0 = this
   541         have ath: "\<And>m x e::real. m \<le> x \<Longrightarrow> x < m + e \<Longrightarrow> \<bar>x - m\<bar> < e"
   542           by arith
   543         from s1m[OF g(1)[rule_format]] have th31: "?m \<le> cmod(poly p (g (f (N1 + N2))))" .
   544         from seq_suble[OF fz(1), of "N1+N2"]
   545         have th00: "real (Suc (N1 + N2)) \<le> real (Suc (f (N1 + N2)))"
   546           by simp
   547         have th000: "0 \<le> (1::real)" "(1::real) \<le> 1" "real (Suc (N1 + N2)) > 0"
   548           using N2 by auto
   549         from frac_le[OF th000 th00]
   550         have th00: "?m +1 / real (Suc (f (N1 + N2))) \<le> ?m + 1 / real (Suc (N1 + N2))"
   551           by simp
   552         from g(2)[rule_format, of "f (N1 + N2)"]
   553         have th01:"cmod (poly p (g (f (N1 + N2)))) < - s + 1 / real (Suc (f (N1 + N2)))" .
   554         from order_less_le_trans[OF th01 th00]
   555         have th32: "cmod(poly p (g (f (N1 + N2)))) < ?m + (1/ real(Suc (N1 + N2)))" .
   556         from N2 have "2/?e < real (Suc (N1 + N2))"
   557           by arith
   558         with e2 less_imp_inverse_less[of "2/?e" "real (Suc (N1 + N2))"]
   559         have "?e/2 > 1/ real (Suc (N1 + N2))"
   560           by (simp add: inverse_eq_divide)
   561         with ath[OF th31 th32]
   562         have thc1: "\<bar>cmod(poly p (g (f (N1 + N2)))) - ?m\<bar> < ?e/2"
   563           by arith
   564         have ath2: "\<And>a b c m::real. \<bar>a - b\<bar> \<le> c \<Longrightarrow> \<bar>b - m\<bar> \<le> \<bar>a - m\<bar> + c"
   565           by arith
   566         have th22: "\<bar>cmod (poly p (g (f (N1 + N2)))) - cmod (poly p z)\<bar> \<le>
   567             cmod (poly p (g (f (N1 + N2))) - poly p z)"
   568           by (simp add: norm_triangle_ineq3)
   569         from ath2[OF th22, of ?m]
   570         have thc2: "2 * (?e/2) \<le>
   571             \<bar>cmod(poly p (g (f (N1 + N2)))) - ?m\<bar> + cmod (poly p (g (f (N1 + N2))) - poly p z)"
   572           by simp
   573         from th0[OF th2 thc1 thc2] have False .
   574       }
   575       then have "?e = 0"
   576         by auto
   577       then have "cmod (poly p z) = ?m"
   578         by simp
   579       with s1m[OF wr] have "cmod (poly p z) \<le> cmod (poly p w)"
   580         by simp
   581     }
   582     then have ?thesis by blast
   583   }
   584   ultimately show ?thesis by blast
   585 qed
   586 
   587 lemma "(rcis (sqrt (abs r)) (a/2))\<^sup>2 = rcis (abs r) a"
   588   unfolding power2_eq_square
   589   apply (simp add: rcis_mult)
   590   apply (simp add: power2_eq_square[symmetric])
   591   done
   592 
   593 lemma cispi: "cis pi = -1"
   594   by (simp add: cis_def)
   595 
   596 lemma "(rcis (sqrt (abs r)) ((pi + a)/2))\<^sup>2 = rcis (- abs r) a"
   597   unfolding power2_eq_square
   598   apply (simp add: rcis_mult add_divide_distrib)
   599   apply (simp add: power2_eq_square[symmetric] rcis_def cispi cis_mult[symmetric])
   600   done
   601 
   602 text {* Nonzero polynomial in z goes to infinity as z does. *}
   603 
   604 lemma poly_infinity:
   605   fixes p:: "'a::{comm_semiring_0,real_normed_div_algebra} poly"
   606   assumes ex: "p \<noteq> 0"
   607   shows "\<exists>r. \<forall>z. r \<le> norm z \<longrightarrow> d \<le> norm (poly (pCons a p) z)"
   608   using ex
   609 proof (induct p arbitrary: a d)
   610   case (pCons c cs a d)
   611   {
   612     assume H: "cs \<noteq> 0"
   613     with pCons.hyps obtain r where r: "\<forall>z. r \<le> norm z \<longrightarrow> d + norm a \<le> norm (poly (pCons c cs) z)"
   614       by blast
   615     let ?r = "1 + \<bar>r\<bar>"
   616     {
   617       fix z::'a
   618       assume h: "1 + \<bar>r\<bar> \<le> norm z"
   619       have r0: "r \<le> norm z" using h by arith
   620       from r[rule_format, OF r0] have th0: "d + norm a \<le> 1 * norm(poly (pCons c cs) z)"
   621         by arith
   622       from h have z1: "norm z \<ge> 1"
   623         by arith
   624       from order_trans[OF th0 mult_right_mono[OF z1 norm_ge_zero[of "poly (pCons c cs) z"]]]
   625       have th1: "d \<le> norm(z * poly (pCons c cs) z) - norm a"
   626         unfolding norm_mult by (simp add: algebra_simps)
   627       from norm_diff_ineq[of "z * poly (pCons c cs) z" a]
   628       have th2: "norm(z * poly (pCons c cs) z) - norm a \<le> norm (poly (pCons a (pCons c cs)) z)"
   629         by (simp add: algebra_simps)
   630       from th1 th2 have "d \<le> norm (poly (pCons a (pCons c cs)) z)" by arith
   631     }
   632     then have ?case by blast
   633   }
   634   moreover
   635   {
   636     assume cs0: "\<not> (cs \<noteq> 0)"
   637     with pCons.prems have c0: "c \<noteq> 0"
   638       by simp
   639     from cs0 have cs0': "cs = 0"
   640       by simp
   641     {
   642       fix z::'a
   643       assume h: "(\<bar>d\<bar> + norm a) / norm c \<le> norm z"
   644       from c0 have "norm c > 0"
   645         by simp
   646       from h c0 have th0: "\<bar>d\<bar> + norm a \<le> norm (z * c)"
   647         by (simp add: field_simps norm_mult)
   648       have ath: "\<And>mzh mazh ma. mzh \<le> mazh + ma \<Longrightarrow> \<bar>d\<bar> + ma \<le> mzh \<Longrightarrow> d \<le> mazh"
   649         by arith
   650       from norm_diff_ineq[of "z * c" a] have th1: "norm (z * c) \<le> norm (a + z * c) + norm a"
   651         by (simp add: algebra_simps)
   652       from ath[OF th1 th0] have "d \<le> norm (poly (pCons a (pCons c cs)) z)"
   653         using cs0' by simp
   654     }
   655     then have ?case  by blast
   656   }
   657   ultimately show ?case by blast
   658 qed simp
   659 
   660 text {* Hence polynomial's modulus attains its minimum somewhere. *}
   661 lemma poly_minimum_modulus: "\<exists>z.\<forall>w. cmod (poly p z) \<le> cmod (poly p w)"
   662 proof (induct p)
   663   case 0
   664   then show ?case by simp
   665 next
   666   case (pCons c cs)
   667   show ?case
   668   proof (cases "cs = 0")
   669     case False
   670     from poly_infinity[OF False, of "cmod (poly (pCons c cs) 0)" c]
   671     obtain r where r: "\<And>z. r \<le> cmod z \<Longrightarrow> cmod (poly (pCons c cs) 0) \<le> cmod (poly (pCons c cs) z)"
   672       by blast
   673     have ath: "\<And>z r. r \<le> cmod z \<or> cmod z \<le> \<bar>r\<bar>"
   674       by arith
   675     from poly_minimum_modulus_disc[of "\<bar>r\<bar>" "pCons c cs"]
   676     obtain v where v: "\<And>w. cmod w \<le> \<bar>r\<bar> \<Longrightarrow> cmod (poly (pCons c cs) v) \<le> cmod (poly (pCons c cs) w)"
   677       by blast
   678     {
   679       fix z
   680       assume z: "r \<le> cmod z"
   681       from v[of 0] r[OF z] have "cmod (poly (pCons c cs) v) \<le> cmod (poly (pCons c cs) z)"
   682         by simp
   683     }
   684     note v0 = this
   685     from v0 v ath[of r] show ?thesis
   686       by blast
   687   next
   688     case True
   689     with pCons.hyps show ?thesis by simp
   690   qed
   691 qed
   692 
   693 text{* Constant function (non-syntactic characterization). *}
   694 definition "constant f = (\<forall>x y. f x = f y)"
   695 
   696 lemma nonconstant_length: "\<not> constant (poly p) \<Longrightarrow> psize p \<ge> 2"
   697   by (induct p) (auto simp: constant_def psize_def)
   698 
   699 lemma poly_replicate_append:
   700   "poly (monom 1 n * p) (x::'a::{comm_ring_1}) = x^n * poly p x"
   701   by (simp add: poly_monom)
   702 
   703 text {* Decomposition of polynomial, skipping zero coefficients
   704   after the first.  *}
   705 
   706 lemma poly_decompose_lemma:
   707   assumes nz: "\<not> (\<forall>z. z \<noteq> 0 \<longrightarrow> poly p z = (0::'a::idom))"
   708   shows "\<exists>k a q. a \<noteq> 0 \<and> Suc (psize q + k) = psize p \<and>
   709     (\<forall>z. poly p z = z^k * poly (pCons a q) z)"
   710   unfolding psize_def
   711   using nz
   712 proof (induct p)
   713   case 0
   714   then show ?case by simp
   715 next
   716   case (pCons c cs)
   717   show ?case
   718   proof (cases "c = 0")
   719     case True
   720     from pCons.hyps pCons.prems True show ?thesis
   721       apply (auto)
   722       apply (rule_tac x="k+1" in exI)
   723       apply (rule_tac x="a" in exI, clarsimp)
   724       apply (rule_tac x="q" in exI)
   725       apply auto
   726       done
   727   next
   728     case False
   729     show ?thesis
   730       apply (rule exI[where x=0])
   731       apply (rule exI[where x=c], auto simp add: False)
   732       done
   733   qed
   734 qed
   735 
   736 lemma poly_decompose:
   737   assumes nc: "\<not> constant (poly p)"
   738   shows "\<exists>k a q. a \<noteq> (0::'a::idom) \<and> k \<noteq> 0 \<and>
   739                psize q + k + 1 = psize p \<and>
   740               (\<forall>z. poly p z = poly p 0 + z^k * poly (pCons a q) z)"
   741   using nc
   742 proof (induct p)
   743   case 0
   744   then show ?case
   745     by (simp add: constant_def)
   746 next
   747   case (pCons c cs)
   748   {
   749     assume C:"\<forall>z. z \<noteq> 0 \<longrightarrow> poly cs z = 0"
   750     {
   751       fix x y
   752       from C have "poly (pCons c cs) x = poly (pCons c cs) y"
   753         by (cases "x = 0") auto
   754     }
   755     with pCons.prems have False
   756       by (auto simp add: constant_def)
   757   }
   758   then have th: "\<not> (\<forall>z. z \<noteq> 0 \<longrightarrow> poly cs z = 0)" ..
   759   from poly_decompose_lemma[OF th]
   760   show ?case
   761     apply clarsimp
   762     apply (rule_tac x="k+1" in exI)
   763     apply (rule_tac x="a" in exI)
   764     apply simp
   765     apply (rule_tac x="q" in exI)
   766     apply (auto simp add: psize_def split: if_splits)
   767     done
   768 qed
   769 
   770 text{* Fundamental theorem of algebra *}
   771 
   772 lemma fundamental_theorem_of_algebra:
   773   assumes nc: "\<not> constant (poly p)"
   774   shows "\<exists>z::complex. poly p z = 0"
   775   using nc
   776 proof (induct "psize p" arbitrary: p rule: less_induct)
   777   case less
   778   let ?p = "poly p"
   779   let ?ths = "\<exists>z. ?p z = 0"
   780 
   781   from nonconstant_length[OF less(2)] have n2: "psize p \<ge> 2" .
   782   from poly_minimum_modulus obtain c where c: "\<forall>w. cmod (?p c) \<le> cmod (?p w)"
   783     by blast
   784 
   785   show ?ths
   786   proof (cases "?p c = 0")
   787     case True
   788     then show ?thesis by blast
   789   next
   790     case False
   791     note pc0 = this
   792     from poly_offset[of p c] obtain q where q: "psize q = psize p" "\<forall>x. poly q x = ?p (c + x)"
   793       by blast
   794     {
   795       assume h: "constant (poly q)"
   796       from q(2) have th: "\<forall>x. poly q (x - c) = ?p x" by auto
   797       {
   798         fix x y
   799         from th have "?p x = poly q (x - c)" by auto
   800         also have "\<dots> = poly q (y - c)"
   801           using h unfolding constant_def by blast
   802         also have "\<dots> = ?p y" using th by auto
   803         finally have "?p x = ?p y" .
   804       }
   805       with less(2) have False
   806         unfolding constant_def by blast
   807     }
   808     then have qnc: "\<not> constant (poly q)"
   809       by blast
   810     from q(2) have pqc0: "?p c = poly q 0"
   811       by simp
   812     from c pqc0 have cq0: "\<forall>w. cmod (poly q 0) \<le> cmod (?p w)"
   813       by simp
   814     let ?a0 = "poly q 0"
   815     from pc0 pqc0 have a00: "?a0 \<noteq> 0"
   816       by simp
   817     from a00 have qr: "\<forall>z. poly q z = poly (smult (inverse ?a0) q) z * ?a0"
   818       by simp
   819     let ?r = "smult (inverse ?a0) q"
   820     have lgqr: "psize q = psize ?r"
   821       using a00
   822       unfolding psize_def degree_def
   823       by (simp add: poly_eq_iff)
   824     {
   825       assume h: "\<And>x y. poly ?r x = poly ?r y"
   826       {
   827         fix x y
   828         from qr[rule_format, of x] have "poly q x = poly ?r x * ?a0"
   829           by auto
   830         also have "\<dots> = poly ?r y * ?a0"
   831           using h by simp
   832         also have "\<dots> = poly q y"
   833           using qr[rule_format, of y] by simp
   834         finally have "poly q x = poly q y" .
   835       }
   836       with qnc have False unfolding constant_def by blast
   837     }
   838     then have rnc: "\<not> constant (poly ?r)"
   839       unfolding constant_def by blast
   840     from qr[rule_format, of 0] a00 have r01: "poly ?r 0 = 1"
   841       by auto
   842     {
   843       fix w
   844       have "cmod (poly ?r w) < 1 \<longleftrightarrow> cmod (poly q w / ?a0) < 1"
   845         using qr[rule_format, of w] a00 by (simp add: divide_inverse mult_ac)
   846       also have "\<dots> \<longleftrightarrow> cmod (poly q w) < cmod ?a0"
   847         using a00 unfolding norm_divide by (simp add: field_simps)
   848       finally have "cmod (poly ?r w) < 1 \<longleftrightarrow> cmod (poly q w) < cmod ?a0" .
   849     }
   850     note mrmq_eq = this
   851     from poly_decompose[OF rnc] obtain k a s where
   852       kas: "a \<noteq> 0" "k \<noteq> 0" "psize s + k + 1 = psize ?r"
   853         "\<forall>z. poly ?r z = poly ?r 0 + z^k* poly (pCons a s) z" by blast
   854     {
   855       assume "psize p = k + 1"
   856       with kas(3) lgqr[symmetric] q(1) have s0: "s = 0"
   857         by auto
   858       {
   859         fix w
   860         have "cmod (poly ?r w) = cmod (1 + a * w ^ k)"
   861           using kas(4)[rule_format, of w] s0 r01 by (simp add: algebra_simps)
   862       }
   863       note hth = this [symmetric]
   864       from reduce_poly_simple[OF kas(1,2)] have "\<exists>w. cmod (poly ?r w) < 1"
   865         unfolding hth by blast
   866     }
   867     moreover
   868     {
   869       assume kn: "psize p \<noteq> k + 1"
   870       from kn kas(3) q(1) lgqr have k1n: "k + 1 < psize p"
   871         by simp
   872       have th01: "\<not> constant (poly (pCons 1 (monom a (k - 1))))"
   873         unfolding constant_def poly_pCons poly_monom
   874         using kas(1) apply simp
   875         apply (rule exI[where x=0])
   876         apply (rule exI[where x=1])
   877         apply simp
   878         done
   879       from kas(1) kas(2) have th02: "k + 1 = psize (pCons 1 (monom a (k - 1)))"
   880         by (simp add: psize_def degree_monom_eq)
   881       from less(1) [OF k1n [simplified th02] th01]
   882       obtain w where w: "1 + w^k * a = 0"
   883         unfolding poly_pCons poly_monom
   884         using kas(2) by (cases k) (auto simp add: algebra_simps)
   885       from poly_bound_exists[of "cmod w" s] obtain m where
   886         m: "m > 0" "\<forall>z. cmod z \<le> cmod w \<longrightarrow> cmod (poly s z) \<le> m" by blast
   887       have w0: "w \<noteq> 0" using kas(2) w
   888         by (auto simp add: power_0_left)
   889       from w have "(1 + w ^ k * a) - 1 = 0 - 1"
   890         by simp
   891       then have wm1: "w^k * a = - 1"
   892         by simp
   893       have inv0: "0 < inverse (cmod w ^ (k + 1) * m)"
   894         using norm_ge_zero[of w] w0 m(1)
   895         by (simp add: inverse_eq_divide zero_less_mult_iff)
   896       with real_lbound_gt_zero[OF zero_less_one] obtain t where
   897         t: "t > 0" "t < 1" "t < inverse (cmod w ^ (k + 1) * m)" by blast
   898       let ?ct = "complex_of_real t"
   899       let ?w = "?ct * w"
   900       have "1 + ?w^k * (a + ?w * poly s ?w) = 1 + ?ct^k * (w^k * a) + ?w^k * ?w * poly s ?w"
   901         using kas(1) by (simp add: algebra_simps power_mult_distrib)
   902       also have "\<dots> = complex_of_real (1 - t^k) + ?w^k * ?w * poly s ?w"
   903         unfolding wm1 by simp
   904       finally have "cmod (1 + ?w^k * (a + ?w * poly s ?w)) =
   905         cmod (complex_of_real (1 - t^k) + ?w^k * ?w * poly s ?w)"
   906         by metis
   907       with norm_triangle_ineq[of "complex_of_real (1 - t^k)" "?w^k * ?w * poly s ?w"]
   908       have th11: "cmod (1 + ?w^k * (a + ?w * poly s ?w)) \<le> \<bar>1 - t^k\<bar> + cmod (?w^k * ?w * poly s ?w)"
   909         unfolding norm_of_real by simp
   910       have ath: "\<And>x t::real. 0 \<le> x \<Longrightarrow> x < t \<Longrightarrow> t \<le> 1 \<Longrightarrow> \<bar>1 - t\<bar> + x < 1"
   911         by arith
   912       have "t * cmod w \<le> 1 * cmod w"
   913         apply (rule mult_mono)
   914         using t(1,2)
   915         apply auto
   916         done
   917       then have tw: "cmod ?w \<le> cmod w"
   918         using t(1) by (simp add: norm_mult)
   919       from t inv0 have "t * (cmod w ^ (k + 1) * m) < 1"
   920         by (simp add: inverse_eq_divide field_simps)
   921       with zero_less_power[OF t(1), of k] have th30: "t^k * (t* (cmod w ^ (k + 1) * m)) < t^k * 1"
   922         by (metis comm_mult_strict_left_mono)
   923       have "cmod (?w^k * ?w * poly s ?w) = t^k * (t* (cmod w ^ (k + 1) * cmod (poly s ?w)))"
   924         using w0 t(1)
   925         by (simp add: algebra_simps power_mult_distrib norm_power norm_mult)
   926       then have "cmod (?w^k * ?w * poly s ?w) \<le> t^k * (t* (cmod w ^ (k + 1) * m))"
   927         using t(1,2) m(2)[rule_format, OF tw] w0
   928         by auto
   929       with th30 have th120: "cmod (?w^k * ?w * poly s ?w) < t^k"
   930         by simp
   931       from power_strict_mono[OF t(2), of k] t(1) kas(2) have th121: "t^k \<le> 1"
   932         by auto
   933       from ath[OF norm_ge_zero[of "?w^k * ?w * poly s ?w"] th120 th121]
   934       have th12: "\<bar>1 - t^k\<bar> + cmod (?w^k * ?w * poly s ?w) < 1" .
   935       from th11 th12 have "cmod (1 + ?w^k * (a + ?w * poly s ?w)) < 1"
   936         by arith
   937       then have "cmod (poly ?r ?w) < 1"
   938         unfolding kas(4)[rule_format, of ?w] r01 by simp
   939       then have "\<exists>w. cmod (poly ?r w) < 1"
   940         by blast
   941     }
   942     ultimately have cr0_contr: "\<exists>w. cmod (poly ?r w) < 1"
   943       by blast
   944     from cr0_contr cq0 q(2) show ?thesis
   945       unfolding mrmq_eq not_less[symmetric] by auto
   946   qed
   947 qed
   948 
   949 text {* Alternative version with a syntactic notion of constant polynomial. *}
   950 
   951 lemma fundamental_theorem_of_algebra_alt:
   952   assumes nc: "\<not> (\<exists>a l. a \<noteq> 0 \<and> l = 0 \<and> p = pCons a l)"
   953   shows "\<exists>z. poly p z = (0::complex)"
   954   using nc
   955 proof (induct p)
   956   case 0
   957   then show ?case by simp
   958 next
   959   case (pCons c cs)
   960   show ?case
   961   proof (cases "c = 0")
   962     case True
   963     then show ?thesis by auto
   964   next
   965     case False
   966     {
   967       assume nc: "constant (poly (pCons c cs))"
   968       from nc[unfolded constant_def, rule_format, of 0]
   969       have "\<forall>w. w \<noteq> 0 \<longrightarrow> poly cs w = 0" by auto
   970       then have "cs = 0"
   971       proof (induct cs)
   972         case 0
   973         then show ?case by simp
   974       next
   975         case (pCons d ds)
   976         show ?case
   977         proof (cases "d = 0")
   978           case True
   979           then show ?thesis using pCons.prems pCons.hyps by simp
   980         next
   981           case False
   982           from poly_bound_exists[of 1 ds] obtain m where
   983             m: "m > 0" "\<forall>z. \<forall>z. cmod z \<le> 1 \<longrightarrow> cmod (poly ds z) \<le> m" by blast
   984           have dm: "cmod d / m > 0" using False m(1) by (simp add: field_simps)
   985           from real_lbound_gt_zero[OF dm zero_less_one] obtain x where
   986             x: "x > 0" "x < cmod d / m" "x < 1" by blast
   987           let ?x = "complex_of_real x"
   988           from x have cx: "?x \<noteq> 0"  "cmod ?x \<le> 1" by simp_all
   989           from pCons.prems[rule_format, OF cx(1)]
   990           have cth: "cmod (?x*poly ds ?x) = cmod d" by (simp add: eq_diff_eq[symmetric])
   991           from m(2)[rule_format, OF cx(2)] x(1)
   992           have th0: "cmod (?x*poly ds ?x) \<le> x*m"
   993             by (simp add: norm_mult)
   994           from x(2) m(1) have "x*m < cmod d" by (simp add: field_simps)
   995           with th0 have "cmod (?x*poly ds ?x) \<noteq> cmod d" by auto
   996           with cth show ?thesis by blast
   997         qed
   998       qed
   999     }
  1000     then have nc: "\<not> constant (poly (pCons c cs))" using pCons.prems False
  1001       by blast
  1002     from fundamental_theorem_of_algebra[OF nc] show ?thesis .
  1003   qed
  1004 qed
  1005 
  1006 
  1007 subsection{* Nullstellensatz, degrees and divisibility of polynomials *}
  1008 
  1009 lemma nullstellensatz_lemma:
  1010   fixes p :: "complex poly"
  1011   assumes "\<forall>x. poly p x = 0 \<longrightarrow> poly q x = 0"
  1012     and "degree p = n"
  1013     and "n \<noteq> 0"
  1014   shows "p dvd (q ^ n)"
  1015   using assms
  1016 proof (induct n arbitrary: p q rule: nat_less_induct)
  1017   fix n :: nat
  1018   fix p q :: "complex poly"
  1019   assume IH: "\<forall>m<n. \<forall>p q.
  1020                  (\<forall>x. poly p x = (0::complex) \<longrightarrow> poly q x = 0) \<longrightarrow>
  1021                  degree p = m \<longrightarrow> m \<noteq> 0 \<longrightarrow> p dvd (q ^ m)"
  1022     and pq0: "\<forall>x. poly p x = 0 \<longrightarrow> poly q x = 0"
  1023     and dpn: "degree p = n"
  1024     and n0: "n \<noteq> 0"
  1025   from dpn n0 have pne: "p \<noteq> 0" by auto
  1026   let ?ths = "p dvd (q ^ n)"
  1027   {
  1028     fix a
  1029     assume a: "poly p a = 0"
  1030     {
  1031       assume oa: "order a p \<noteq> 0"
  1032       let ?op = "order a p"
  1033       from pne have ap: "([:- a, 1:] ^ ?op) dvd p" "\<not> [:- a, 1:] ^ (Suc ?op) dvd p"
  1034         using order by blast+
  1035       note oop = order_degree[OF pne, unfolded dpn]
  1036       {
  1037         assume q0: "q = 0"
  1038         then have ?ths using n0
  1039           by (simp add: power_0_left)
  1040       }
  1041       moreover
  1042       {
  1043         assume q0: "q \<noteq> 0"
  1044         from pq0[rule_format, OF a, unfolded poly_eq_0_iff_dvd]
  1045         obtain r where r: "q = [:- a, 1:] * r" by (rule dvdE)
  1046         from ap(1) obtain s where s: "p = [:- a, 1:] ^ ?op * s"
  1047           by (rule dvdE)
  1048         have sne: "s \<noteq> 0" using s pne by auto
  1049         {
  1050           assume ds0: "degree s = 0"
  1051           from ds0 obtain k where kpn: "s = [:k:]"
  1052             by (cases s) (auto split: if_splits)
  1053           from sne kpn have k: "k \<noteq> 0" by simp
  1054           let ?w = "([:1/k:] * ([:-a,1:] ^ (n - ?op))) * (r ^ n)"
  1055           have "q ^ n = p * ?w"
  1056             apply (subst r, subst s, subst kpn)
  1057             using k oop [of a]
  1058             apply (subst power_mult_distrib, simp)
  1059             apply (subst power_add [symmetric], simp)
  1060             done
  1061           then have ?ths unfolding dvd_def by blast
  1062         }
  1063         moreover
  1064         {
  1065           assume ds0: "degree s \<noteq> 0"
  1066           from ds0 sne dpn s oa
  1067             have dsn: "degree s < n"
  1068               apply auto
  1069               apply (erule ssubst)
  1070               apply (simp add: degree_mult_eq degree_linear_power)
  1071               done
  1072             {
  1073               fix x assume h: "poly s x = 0"
  1074               {
  1075                 assume xa: "x = a"
  1076                 from h[unfolded xa poly_eq_0_iff_dvd] obtain u where u: "s = [:- a, 1:] * u"
  1077                   by (rule dvdE)
  1078                 have "p = [:- a, 1:] ^ (Suc ?op) * u"
  1079                   by (subst s, subst u, simp only: power_Suc mult_ac)
  1080                 with ap(2)[unfolded dvd_def] have False by blast
  1081               }
  1082               note xa = this
  1083               from h have "poly p x = 0" by (subst s) simp
  1084               with pq0 have "poly q x = 0" by blast
  1085               with r xa have "poly r x = 0"
  1086                 by auto
  1087             }
  1088             note impth = this
  1089             from IH[rule_format, OF dsn, of s r] impth ds0
  1090             have "s dvd (r ^ (degree s))" by blast
  1091             then obtain u where u: "r ^ (degree s) = s * u" ..
  1092             then have u': "\<And>x. poly s x * poly u x = poly r x ^ degree s"
  1093               by (simp only: poly_mult[symmetric] poly_power[symmetric])
  1094             let ?w = "(u * ([:-a,1:] ^ (n - ?op))) * (r ^ (n - degree s))"
  1095             from oop[of a] dsn have "q ^ n = p * ?w"
  1096               apply -
  1097               apply (subst s, subst r)
  1098               apply (simp only: power_mult_distrib)
  1099               apply (subst mult_assoc [where b=s])
  1100               apply (subst mult_assoc [where a=u])
  1101               apply (subst mult_assoc [where b=u, symmetric])
  1102               apply (subst u [symmetric])
  1103               apply (simp add: mult_ac power_add [symmetric])
  1104               done
  1105             then have ?ths unfolding dvd_def by blast
  1106         }
  1107         ultimately have ?ths by blast
  1108       }
  1109       ultimately have ?ths by blast
  1110     }
  1111     then have ?ths using a order_root pne by blast
  1112   }
  1113   moreover
  1114   {
  1115     assume exa: "\<not> (\<exists>a. poly p a = 0)"
  1116     from fundamental_theorem_of_algebra_alt[of p] exa
  1117     obtain c where ccs: "c \<noteq> 0" "p = pCons c 0"
  1118       by blast
  1119     then have pp: "\<And>x. poly p x = c"
  1120       by simp
  1121     let ?w = "[:1/c:] * (q ^ n)"
  1122     from ccs have "(q ^ n) = (p * ?w)"
  1123       by simp
  1124     then have ?ths
  1125       unfolding dvd_def by blast
  1126   }
  1127   ultimately show ?ths by blast
  1128 qed
  1129 
  1130 lemma nullstellensatz_univariate:
  1131   "(\<forall>x. poly p x = (0::complex) \<longrightarrow> poly q x = 0) \<longleftrightarrow>
  1132     p dvd (q ^ (degree p)) \<or> (p = 0 \<and> q = 0)"
  1133 proof -
  1134   {
  1135     assume pe: "p = 0"
  1136     then have eq: "(\<forall>x. poly p x = (0::complex) \<longrightarrow> poly q x = 0) \<longleftrightarrow> q = 0"
  1137       by (auto simp add: poly_all_0_iff_0)
  1138     {
  1139       assume "p dvd (q ^ (degree p))"
  1140       then obtain r where r: "q ^ (degree p) = p * r" ..
  1141       from r pe have False by simp
  1142     }
  1143     with eq pe have ?thesis by blast
  1144   }
  1145   moreover
  1146   {
  1147     assume pe: "p \<noteq> 0"
  1148     {
  1149       assume dp: "degree p = 0"
  1150       then obtain k where k: "p = [:k:]" "k \<noteq> 0" using pe
  1151         by (cases p) (simp split: if_splits)
  1152       then have th1: "\<forall>x. poly p x \<noteq> 0"
  1153         by simp
  1154       from k dp have "q ^ (degree p) = p * [:1/k:]"
  1155         by (simp add: one_poly_def)
  1156       then have th2: "p dvd (q ^ (degree p))" ..
  1157       from th1 th2 pe have ?thesis by blast
  1158     }
  1159     moreover
  1160     {
  1161       assume dp: "degree p \<noteq> 0"
  1162       then obtain n where n: "degree p = Suc n "
  1163         by (cases "degree p") auto
  1164       {
  1165         assume "p dvd (q ^ (Suc n))"
  1166         then obtain u where u: "q ^ (Suc n) = p * u" ..
  1167         {
  1168           fix x
  1169           assume h: "poly p x = 0" "poly q x \<noteq> 0"
  1170           then have "poly (q ^ (Suc n)) x \<noteq> 0"
  1171             by simp
  1172           then have False using u h(1)
  1173             by (simp only: poly_mult) simp
  1174         }
  1175       }
  1176       with n nullstellensatz_lemma[of p q "degree p"] dp
  1177       have ?thesis by auto
  1178     }
  1179     ultimately have ?thesis by blast
  1180   }
  1181   ultimately show ?thesis by blast
  1182 qed
  1183 
  1184 text{* Useful lemma *}
  1185 
  1186 lemma constant_degree:
  1187   fixes p :: "'a::{idom,ring_char_0} poly"
  1188   shows "constant (poly p) \<longleftrightarrow> degree p = 0" (is "?lhs = ?rhs")
  1189 proof
  1190   assume l: ?lhs
  1191   from l[unfolded constant_def, rule_format, of _ "0"]
  1192   have th: "poly p = poly [:poly p 0:]"
  1193     by auto
  1194   then have "p = [:poly p 0:]"
  1195     by (simp add: poly_eq_poly_eq_iff)
  1196   then have "degree p = degree [:poly p 0:]"
  1197     by simp
  1198   then show ?rhs
  1199     by simp
  1200 next
  1201   assume r: ?rhs
  1202   then obtain k where "p = [:k:]"
  1203     by (cases p) (simp split: if_splits)
  1204   then show ?lhs
  1205     unfolding constant_def by auto
  1206 qed
  1207 
  1208 lemma divides_degree:
  1209   assumes pq: "p dvd (q:: complex poly)"
  1210   shows "degree p \<le> degree q \<or> q = 0"
  1211   by (metis dvd_imp_degree_le pq)
  1212 
  1213 (* Arithmetic operations on multivariate polynomials.                        *)
  1214 
  1215 lemma mpoly_base_conv:
  1216   fixes x :: "'a::comm_ring_1"
  1217   shows "0 = poly 0 x" "c = poly [:c:] x" "x = poly [:0,1:] x"
  1218   by simp_all
  1219 
  1220 lemma mpoly_norm_conv:
  1221   fixes x :: "'a::comm_ring_1"
  1222   shows "poly [:0:] x = poly 0 x" "poly [:poly 0 y:] x = poly 0 x"
  1223   by simp_all
  1224 
  1225 lemma mpoly_sub_conv:
  1226   fixes x :: "'a::comm_ring_1"
  1227   shows "poly p x - poly q x = poly p x + -1 * poly q x"
  1228   by simp
  1229 
  1230 lemma poly_pad_rule: "poly p x = 0 \<Longrightarrow> poly (pCons 0 p) x = 0"
  1231   by simp
  1232 
  1233 lemma poly_cancel_eq_conv:
  1234   fixes x :: "'a::field"
  1235   shows "x = 0 \<Longrightarrow> a \<noteq> 0 \<Longrightarrow> (y = 0) = (a * y - b * x = 0)"
  1236   by auto
  1237 
  1238 lemma poly_divides_pad_rule:
  1239   fixes p:: "('a::comm_ring_1) poly"
  1240   assumes pq: "p dvd q"
  1241   shows "p dvd (pCons 0 q)"
  1242 proof -
  1243   have "pCons 0 q = q * [:0,1:]" by simp
  1244   then have "q dvd (pCons 0 q)" ..
  1245   with pq show ?thesis by (rule dvd_trans)
  1246 qed
  1247 
  1248 lemma poly_divides_conv0:
  1249   fixes p:: "'a::field poly"
  1250   assumes lgpq: "degree q < degree p"
  1251     and lq: "p \<noteq> 0"
  1252   shows "p dvd q \<longleftrightarrow> q = 0" (is "?lhs \<longleftrightarrow> ?rhs")
  1253 proof
  1254   assume r: ?rhs
  1255   then have "q = p * 0" by simp
  1256   then show ?lhs ..
  1257 next
  1258   assume l: ?lhs
  1259   show ?rhs
  1260   proof (cases "q = 0")
  1261     case True
  1262     then show ?thesis by simp
  1263   next
  1264     assume q0: "q \<noteq> 0"
  1265     from l q0 have "degree p \<le> degree q"
  1266       by (rule dvd_imp_degree_le)
  1267     with lgpq show ?thesis by simp
  1268   qed
  1269 qed
  1270 
  1271 lemma poly_divides_conv1:
  1272   fixes p :: "'a::field poly"
  1273   assumes a0: "a \<noteq> 0"
  1274     and pp': "p dvd p'"
  1275     and qrp': "smult a q - p' = r"
  1276   shows "p dvd q \<longleftrightarrow> p dvd r" (is "?lhs \<longleftrightarrow> ?rhs")
  1277 proof
  1278   from pp' obtain t where t: "p' = p * t" ..
  1279   {
  1280     assume l: ?lhs
  1281     then obtain u where u: "q = p * u" ..
  1282     have "r = p * (smult a u - t)"
  1283       using u qrp' [symmetric] t by (simp add: algebra_simps)
  1284     then show ?rhs ..
  1285   next
  1286     assume r: ?rhs
  1287     then obtain u where u: "r = p * u" ..
  1288     from u [symmetric] t qrp' [symmetric] a0
  1289     have "q = p * smult (1/a) (u + t)" by (simp add: algebra_simps)
  1290     then show ?lhs ..
  1291   }
  1292 qed
  1293 
  1294 lemma basic_cqe_conv1:
  1295   "(\<exists>x. poly p x = 0 \<and> poly 0 x \<noteq> 0) \<longleftrightarrow> False"
  1296   "(\<exists>x. poly 0 x \<noteq> 0) \<longleftrightarrow> False"
  1297   "(\<exists>x. poly [:c:] x \<noteq> 0) \<longleftrightarrow> c \<noteq> 0"
  1298   "(\<exists>x. poly 0 x = 0) \<longleftrightarrow> True"
  1299   "(\<exists>x. poly [:c:] x = 0) \<longleftrightarrow> c = 0"
  1300   by simp_all
  1301 
  1302 lemma basic_cqe_conv2:
  1303   assumes l:"p \<noteq> 0"
  1304   shows "(\<exists>x. poly (pCons a (pCons b p)) x = (0::complex))"
  1305 proof -
  1306   {
  1307     fix h t
  1308     assume h: "h \<noteq> 0" "t = 0" and "pCons a (pCons b p) = pCons h t"
  1309     with l have False by simp
  1310   }
  1311   then have th: "\<not> (\<exists> h t. h \<noteq> 0 \<and> t = 0 \<and> pCons a (pCons b p) = pCons h t)"
  1312     by blast
  1313   from fundamental_theorem_of_algebra_alt[OF th] show ?thesis
  1314     by auto
  1315 qed
  1316 
  1317 lemma  basic_cqe_conv_2b: "(\<exists>x. poly p x \<noteq> (0::complex)) \<longleftrightarrow> p \<noteq> 0"
  1318   by (metis poly_all_0_iff_0)
  1319 
  1320 lemma basic_cqe_conv3:
  1321   fixes p q :: "complex poly"
  1322   assumes l: "p \<noteq> 0"
  1323   shows "(\<exists>x. poly (pCons a p) x = 0 \<and> poly q x \<noteq> 0) \<longleftrightarrow> \<not> ((pCons a p) dvd (q ^ psize p))"
  1324 proof -
  1325   from l have dp: "degree (pCons a p) = psize p"
  1326     by (simp add: psize_def)
  1327   from nullstellensatz_univariate[of "pCons a p" q] l
  1328   show ?thesis
  1329     by (metis dp pCons_eq_0_iff)
  1330 qed
  1331 
  1332 lemma basic_cqe_conv4:
  1333   fixes p q :: "complex poly"
  1334   assumes h: "\<And>x. poly (q ^ n) x = poly r x"
  1335   shows "p dvd (q ^ n) \<longleftrightarrow> p dvd r"
  1336 proof -
  1337   from h have "poly (q ^ n) = poly r"
  1338     by auto
  1339   then have "(q ^ n) = r"
  1340     by (simp add: poly_eq_poly_eq_iff)
  1341   then show "p dvd (q ^ n) \<longleftrightarrow> p dvd r"
  1342     by simp
  1343 qed
  1344 
  1345 lemma poly_const_conv:
  1346   fixes x :: "'a::comm_ring_1"
  1347   shows "poly [:c:] x = y \<longleftrightarrow> c = y"
  1348   by simp
  1349 
  1350 end