Theory MPoly_Type

(* Author: Andreas Lochbihler, ETH Zurich
   Author: Florian Haftmann, TU Muenchen
*)

section ‹An abstract type for multivariate polynomials›

theory MPoly_Type
imports "HOL-Library.Poly_Mapping"
begin

subsection ‹Abstract type definition›

typedef (overloaded) 'a mpoly =
  "UNIV :: ((nat ⇒⇩0 nat) ⇒⇩0 'a::zero) set"
  morphisms mapping_of MPoly
 ..

setup_lifting type_definition_mpoly

(* these theorems are automatically generated by setup_lifting... *)
thm mapping_of_inverse   thm MPoly_inverse
thm mapping_of_inject    thm MPoly_inject
thm mapping_of_induct    thm MPoly_induct
thm mapping_of_cases     thm MPoly_cases


subsection ‹Additive structure›

instantiation mpoly :: (zero) zero
begin

lift_definition zero_mpoly :: "'a mpoly"
  is "0 :: (nat ⇒⇩0 nat) ⇒⇩0 'a" .

instance ..

end

instantiation mpoly :: (monoid_add) monoid_add
begin

lift_definition plus_mpoly :: "'a mpoly ⇒ 'a mpoly ⇒ 'a mpoly"
  is "Groups.plus :: ((nat ⇒⇩0 nat) ⇒⇩0 'a) ⇒ _" .

instance
  by intro_classes (transfer, simp add: fun_eq_iff add.assoc)+

end

instance mpoly :: (comm_monoid_add) comm_monoid_add
  by intro_classes (transfer, simp add: fun_eq_iff ac_simps)+

instantiation mpoly :: (cancel_comm_monoid_add) cancel_comm_monoid_add
begin

lift_definition minus_mpoly :: "'a mpoly ⇒ 'a mpoly ⇒ 'a mpoly"
  is "Groups.minus :: ((nat ⇒⇩0 nat) ⇒⇩0 'a) ⇒ _" .

instance
  by intro_classes (transfer, simp add: fun_eq_iff diff_diff_add)+

end

instantiation mpoly :: (ab_group_add) ab_group_add
begin

lift_definition uminus_mpoly :: "'a mpoly ⇒ 'a mpoly"
  is "Groups.uminus :: ((nat ⇒⇩0 nat) ⇒⇩0 'a) ⇒ _" .


instance
  by intro_classes (transfer, simp add: fun_eq_iff add_uminus_conv_diff)+

end


subsection ‹Multiplication by a coefficient›
(* ?do we need inc_power on abstract polynomials? *)

lift_definition smult :: "'a::{times,zero} ⇒ 'a mpoly ⇒ 'a mpoly"
  is "λa. Poly_Mapping.map (Groups.times a) :: ((nat ⇒⇩0 nat) ⇒⇩0 'a) ⇒ _" .

(* left lemmas in subsection ‹Pseudo-division of polynomials›,
   because I couldn't disentangle them and the notion of monomials. *)

subsection ‹Multiplicative structure›

instantiation mpoly :: (zero_neq_one) zero_neq_one
begin

lift_definition one_mpoly :: "'a mpoly"
  is "1 :: ((nat ⇒⇩0 nat) ⇒⇩0 'a)" .

instance
  by intro_classes (transfer, simp)

end

instantiation mpoly :: (semiring_0) semiring_0
begin

lift_definition times_mpoly :: "'a mpoly ⇒ 'a mpoly ⇒ 'a mpoly"
  is "Groups.times :: ((nat ⇒⇩0 nat) ⇒⇩0 'a) ⇒ _" .

instance
  by intro_classes (transfer, simp add: algebra_simps)+

end

instance mpoly :: (comm_semiring_0) comm_semiring_0
  by intro_classes (transfer, simp add: algebra_simps)+

instance mpoly :: (semiring_0_cancel) semiring_0_cancel
  ..

instance mpoly :: (comm_semiring_0_cancel) comm_semiring_0_cancel
  ..

instance mpoly :: (semiring_1) semiring_1
  by intro_classes (transfer, simp)+

instance mpoly :: (comm_semiring_1) comm_semiring_1
  by intro_classes (transfer, simp)+

instance mpoly :: (semiring_1_cancel) semiring_1_cancel
  ..

(*instance mpoly :: (comm_semiring_1_cancel) comm_semiring_1_cancel
  .. FIXME unclear whether this holds *)

instance mpoly :: (ring) ring
  ..

instance mpoly :: (comm_ring) comm_ring
  ..

instance mpoly :: (ring_1) ring_1
  ..

instance mpoly :: (comm_ring_1) comm_ring_1
  ..


subsection ‹Monomials›

text ‹
  Terminology is not unique here, so we use the notions as follows:
  A "monomial" and a "coefficient" together give a "term".
  These notions are significant in connection with "leading",
  "leading term", "leading coefficient" and "leading monomial",
  which all rely on a monomial order.
›

lift_definition monom :: "(nat ⇒⇩0 nat) ⇒ 'a::zero ⇒ 'a mpoly"
  is "Poly_Mapping.single :: (nat ⇒⇩0 nat) ⇒ _" .

lemma mapping_of_monom [simp]:
  "mapping_of (monom m a) = Poly_Mapping.single m a"
  by(fact monom.rep_eq)

lemma monom_zero [simp]:
  "monom 0 0 = 0"
  by transfer simp

lemma monom_one [simp]:
  "monom 0 1 = 1"
  by transfer simp

lemma monom_add:
  "monom m (a + b) = monom m a + monom m b"
  by transfer (simp add: single_add)

lemma monom_uminus:
  "monom m (- a) = - monom m a"
  by transfer (simp add: single_uminus)

lemma monom_diff:
  "monom m (a - b) = monom m a - monom m b"
  by transfer (simp add: single_diff)

lemma monom_numeral [simp]:
  "monom 0 (numeral n) = numeral n"
  by (induct n) (simp_all only: numeral.simps numeral_add monom_zero monom_one monom_add)

lemma monom_of_nat [simp]:
  "monom 0 (of_nat n) = of_nat n"
  by (induct n) (simp_all add: monom_add)

lemma of_nat_monom:
  "of_nat = monom 0 ∘ of_nat"
  by (simp add: fun_eq_iff)

lemma inj_monom [iff]:
  "inj (monom m)"
proof (rule injI, transfer)
  fix a b :: 'a and m :: "nat ⇒⇩0 nat"
  assume "Poly_Mapping.single m a = Poly_Mapping.single m b"
  with injD [of "Poly_Mapping.single m" a b]
  show "a = b" by simp
qed

lemma mult_monom: "monom x a * monom y b = monom (x + y) (a * b)"
  by transfer' (simp add: Poly_Mapping.mult_single)

instance mpoly :: (semiring_char_0) semiring_char_0
  by intro_classes (auto simp add: of_nat_monom inj_of_nat intro: inj_compose)

instance mpoly :: (ring_char_0) ring_char_0
  ..

lemma monom_of_int [simp]:
  "monom 0 (of_int k) = of_int k"
  apply (cases k)
  apply simp_all
  unfolding monom_diff monom_uminus
  apply simp
  done

subsection ‹Constants and Indeterminates›

text ‹Embedding of indeterminates and constants in type-class polynomials,
  can be used as constructors.›

definition Var⇩0 :: "'a ⇒ ('a ⇒⇩0 nat) ⇒⇩0 'b::{one,zero}" where
  "Var⇩0 n ≡ Poly_Mapping.single (Poly_Mapping.single n 1) 1"
definition Const⇩0 :: "'b ⇒ ('a ⇒⇩0 nat) ⇒⇩0 'b::zero" where "Const⇩0 c ≡ Poly_Mapping.single 0 c"

lemma Const⇩0_one: "Const⇩0 1 = 1"
  by (simp add: Const⇩0_def)

lemma Const⇩0_numeral: "Const⇩0 (numeral x) = numeral x"
  by (auto intro!: poly_mapping_eqI simp: Const⇩0_def lookup_numeral)

lemma Const⇩0_minus: "Const⇩0 (- x) = - Const⇩0 x"
  by (simp add: Const⇩0_def single_uminus)

lemma Const⇩0_zero: "Const⇩0 0 = 0"
  by (auto intro!: poly_mapping_eqI simp: Const⇩0_def)

lemma Var⇩0_power: "Var⇩0 v ^ n = Poly_Mapping.single (Poly_Mapping.single v n) 1"
  by (induction n) (auto simp: Var⇩0_def mult_single single_add[symmetric])

lift_definition Var::"nat ⇒ 'b::{one,zero} mpoly" is Var⇩0 .
lift_definition Const::"'b::zero ⇒ 'b mpoly" is Const⇩0 .


subsection ‹Integral domains›

instance mpoly :: (ring_no_zero_divisors) ring_no_zero_divisors
  by intro_classes (transfer, simp)

instance mpoly :: (ring_1_no_zero_divisors) ring_1_no_zero_divisors
  ..

instance mpoly :: (idom) idom
  ..


subsection ‹Monom coefficient lookup›

definition coeff :: "'a::zero mpoly ⇒ (nat ⇒⇩0 nat) ⇒ 'a"
where
  "coeff p = Poly_Mapping.lookup (mapping_of p)"


subsection ‹Insertion morphism›

definition insertion_fun_natural :: "(nat ⇒ 'a) ⇒ ((nat ⇒ nat) ⇒ 'a) ⇒ 'a::comm_semiring_1"
where
  "insertion_fun_natural f p = (∑m. p m * (∏v. f v ^ m v))"

definition insertion_fun :: "(nat ⇒ 'a) ⇒ ((nat ⇒⇩0 nat) ⇒ 'a) ⇒ 'a::comm_semiring_1"
where
  "insertion_fun f p = (∑m. p m * (∏v. f v ^ Poly_Mapping.lookup m v))"

text ‹N.b. have been unable to relate this to @{const insertion_fun_natural} using lifting!›

lift_definition insertion_aux :: "(nat ⇒ 'a) ⇒ ((nat ⇒⇩0 nat) ⇒⇩0 'a) ⇒ 'a::comm_semiring_1"
  is "insertion_fun" .

lift_definition insertion :: "(nat ⇒ 'a) ⇒ 'a mpoly ⇒ 'a::comm_semiring_1"
  is "insertion_aux" .

lemma aux:
  "Poly_Mapping.lookup f = (λ_. 0) ⟷ f = 0"
  apply transfer apply simp done

lemma insertion_trivial [simp]:
  "insertion (λ_. 0) p = coeff p 0"
proof -
  { fix f :: "(nat ⇒⇩0 nat) ⇒⇩0 'a"
    have "insertion_aux (λ_. 0) f = Poly_Mapping.lookup f 0"
      apply (simp add: insertion_aux_def insertion_fun_def power_Sum_any [symmetric])
      apply (simp add: zero_power_eq mult_when aux)
      done
  }
  then show ?thesis by (simp add: coeff_def insertion_def)
qed

lemma insertion_zero [simp]:
  "insertion f 0 = 0"
  by transfer (simp add: insertion_aux_def insertion_fun_def)

lemma insertion_fun_add:
  fixes f p q
  shows "insertion_fun f (Poly_Mapping.lookup (p + q)) =
    insertion_fun f (Poly_Mapping.lookup p) +
      insertion_fun f (Poly_Mapping.lookup q)"
  unfolding insertion_fun_def
  apply (subst Sum_any.distrib [symmetric])
  apply (simp_all add: plus_poly_mapping.rep_eq algebra_simps)
  apply (rule finite_mult_not_eq_zero_rightI)
  apply simp
  apply (rule finite_mult_not_eq_zero_rightI)
  apply simp
  done

lemma insertion_add:
  "insertion f (p + q) = insertion f p + insertion f q"
  by transfer (simp add: insertion_aux_def insertion_fun_add)

lemma insertion_one [simp]:
  "insertion f 1 = 1"
  by transfer (simp add: insertion_aux_def insertion_fun_def one_poly_mapping.rep_eq when_mult)

lemma insertion_fun_mult:
  fixes f p q
  shows "insertion_fun f (Poly_Mapping.lookup (p * q)) =
    insertion_fun f (Poly_Mapping.lookup p) *
      insertion_fun f (Poly_Mapping.lookup q)"
proof -
  { fix m :: "nat ⇒⇩0 nat"
    have "finite {v. Poly_Mapping.lookup m v ≠ 0}"
      by simp
    then have "finite {v. f v ^ Poly_Mapping.lookup m v ≠ 1}"
      by (rule rev_finite_subset) (auto intro: ccontr)
  }
  moreover define g where "g m = (∏v. f v ^ Poly_Mapping.lookup m v)" for m
  ultimately have *: "⋀a b. g (a + b) = g a * g b"
    by (simp add: plus_poly_mapping.rep_eq power_add Prod_any.distrib)
  have bij: "bij (λ(l, n, m). (m, l, n))"
    by (auto intro!: bijI injI simp add: image_def)
  let ?P = "{l. Poly_Mapping.lookup p l ≠ 0}"
  let ?Q = "{n. Poly_Mapping.lookup q n ≠ 0}"
  let ?PQ = "{l + n | l n. l ∈ Poly_Mapping.keys p ∧ n ∈ Poly_Mapping.keys q}"
  have "finite {l + n | l n. Poly_Mapping.lookup p l ≠ 0 ∧ Poly_Mapping.lookup q n ≠ 0}"
    by (rule finite_not_eq_zero_sumI) simp_all
  then have fin_PQ: "finite ?PQ"
    by (simp add: in_keys_iff)
  have "(∑m. Poly_Mapping.lookup (p * q) m * g m) =
    (∑m. (∑l. Poly_Mapping.lookup p l * (∑n. Poly_Mapping.lookup q n when m = l + n)) * g m)"
    by (simp add: times_poly_mapping.rep_eq prod_fun_def)
  also have "… = (∑m. (∑l. (∑n. g m * (Poly_Mapping.lookup p l * Poly_Mapping.lookup q n) when m = l + n)))"
    apply (subst Sum_any_left_distrib)
    apply (auto intro: finite_mult_not_eq_zero_rightI)
    apply (subst Sum_any_right_distrib)
    apply (auto intro: finite_mult_not_eq_zero_rightI)
    apply (subst Sum_any_left_distrib)
    apply (auto intro: finite_mult_not_eq_zero_leftI)
    apply (simp add: ac_simps mult_when)
    done
  also have "… = (∑m. (∑(l, n). g m * (Poly_Mapping.lookup p l * Poly_Mapping.lookup q n) when m = l + n))"
    apply (subst (2) Sum_any.cartesian_product [of "?P × ?Q"])
    apply (auto dest!: mult_not_zero)
    done
  also have "… = (∑(m, l, n). g m * (Poly_Mapping.lookup p l * Poly_Mapping.lookup q n) when m = l + n)"
    apply (subst Sum_any.cartesian_product [of "?PQ × (?P × ?Q)"])
      apply (auto dest!: mult_not_zero simp add: fin_PQ)
    apply (auto simp: in_keys_iff)
    done
  also have "… = (∑(l, n, m). g m * (Poly_Mapping.lookup p l * Poly_Mapping.lookup q n) when m = l + n)"
    using bij by (rule Sum_any.reindex_cong [of "λ(l, n, m). (m, l, n)"]) (simp add: fun_eq_iff)
  also have "… = (∑(l, n). ∑m. g m * (Poly_Mapping.lookup p l * Poly_Mapping.lookup q n) when m = l + n)"
    apply (subst Sum_any.cartesian_product2 [of "(?P × ?Q) × ?PQ"])
    apply (auto dest!: mult_not_zero simp add: fin_PQ )
    apply (auto simp: in_keys_iff)
    done
  also have "… = (∑(l, n). (g l * g n) * (Poly_Mapping.lookup p l * Poly_Mapping.lookup q n))"
    by (simp add: *)
  also have "… = (∑l. ∑n. (g l * g n) * (Poly_Mapping.lookup p l * Poly_Mapping.lookup q n))"
    apply (subst Sum_any.cartesian_product [of "?P × ?Q"])
    apply (auto dest!: mult_not_zero)
    done
  also have "… = (∑l. ∑n. (Poly_Mapping.lookup p l * g l) * (Poly_Mapping.lookup q n * g n))"
    by (simp add: ac_simps)
  also have "… =
    (∑m. Poly_Mapping.lookup p m * g m) *
    (∑m. Poly_Mapping.lookup q m * g m)"
    by (rule Sum_any_product [symmetric]) (auto intro: finite_mult_not_eq_zero_rightI)
  finally show ?thesis by (simp add: insertion_fun_def g_def)
qed

lemma insertion_mult:
  "insertion f (p * q) = insertion f p * insertion f q"
  by transfer (simp add: insertion_aux_def insertion_fun_mult)


subsection ‹Degree›

lift_definition degree :: "'a::zero mpoly ⇒ nat ⇒ nat"
is "λp v. Max (insert 0 ((λm. Poly_Mapping.lookup m v) ` Poly_Mapping.keys p))" .


lift_definition total_degree :: "'a::zero mpoly ⇒ nat"
is "λp. Max (insert 0 ((λm. sum (Poly_Mapping.lookup m) (Poly_Mapping.keys m)) ` Poly_Mapping.keys p))" .

lemma degree_zero [simp]:
  "degree 0 v = 0"
  by transfer simp

lemma total_degree_zero [simp]:
  "total_degree 0 = 0"
  by transfer simp
(*
value [code] "total_degree (0 :: int mpoly)" (***)
*)

lemma degree_one [simp]:
  "degree 1 v = 0"
  by transfer simp

lemma total_degree_one [simp]:
  "total_degree 1 = 0"
  by transfer simp

subsection ‹Pseudo-division of polynomials›

lemma smult_conv_mult: "smult s p = monom 0 s * p"
by transfer (simp add: mult_map_scale_conv_mult)

lemma smult_monom [simp]:
  fixes c :: "_ :: mult_zero"
  shows "smult c (monom x c') = monom x (c * c')"
by transfer simp

lemma smult_0 [simp]:
  fixes p :: "_ :: mult_zero mpoly"
  shows "smult 0 p = 0"
by transfer(simp add: map_eq_zero_iff)

lemma mult_smult_left: "smult s p * q = smult s (p * q)"
by(simp add: smult_conv_mult mult.assoc)

lift_definition sdiv :: "'a::euclidean_ring ⇒ 'a mpoly ⇒ 'a mpoly"
  is "λa. Poly_Mapping.map (λb. b div a) :: ((nat ⇒⇩0 nat) ⇒⇩0 'a) ⇒ _"
.
text ‹
  \qt{Polynomial division} is only possible on univariate polynomials ‹K[x]›
  over a field ‹K›, all other kinds of polynomials only allow pseudo-division
  [1]p.40/41":

  ‹∀x y :: 'a mpoly. y ≠ 0 ⇒ ∃a q r. smult a x = q * y + r›

  The introduction of pseudo-division below generalises @{file ‹~~/src/HOL/Computational_Algebra/Polynomial.thy›}.
  [1] Winkler, Polynomial Algorithms, 1996.
  The generalisation raises issues addressed by Wenda Li and commented below.
  Florian replied to the issues conjecturing, that the abstract mpoly needs not
  be aware of the issues, in case these are only concerned with executability.
›

definition pseudo_divmod_rel
  :: "'a::euclidean_ring => 'a mpoly => 'a mpoly => 'a mpoly => 'a mpoly => bool"
where
  "pseudo_divmod_rel a x y q r ⟷
    smult a x = q * y + r ∧ (if y = 0 then q = 0 else r = 0 ∨ degree r < degree y)"
(* the notion of degree resigns a main variable in multivariate polynomials;
   also, there are infinitely many tuples (a,q,r) such that a x = q y + r *)

definition pdiv :: "'a::euclidean_ring mpoly ⇒ 'a mpoly ⇒ ('a × 'a mpoly)" (infixl "pdiv" 70)
where
  "x pdiv y = (THE (a, q). ∃r. pseudo_divmod_rel a x y q r)"

definition pmod :: "'a::euclidean_ring mpoly ⇒ 'a mpoly ⇒ 'a mpoly" (infixl "pmod" 70)
where
  "x pmod y = (THE r. ∃a q. pseudo_divmod_rel a x y q r)"

definition pdivmod :: "'a::euclidean_ring mpoly ⇒ 'a mpoly ⇒ ('a × 'a mpoly) × 'a mpoly"
where
  "pdivmod p q = (p pdiv q, p pmod q)"

(* "_code" seems inappropriate, since "THE" in definitions pdiv and pmod is not unique,
   see comment to definition pseudo_divmod_rel; so this cannot be executable ... *)
lemma pdiv_code:
  "p pdiv q = fst (pdivmod p q)"
  by (simp add: pdivmod_def)

lemma pmod_code:
  "p pmod q = snd (pdivmod p q)"
  by (simp add: pdivmod_def)

subsection ‹Primitive poly, etc›

lift_definition coeffs :: "'a :: zero mpoly ⇒ 'a set"
is "Poly_Mapping.range :: ((nat ⇒⇩0 nat) ⇒⇩0 'a) ⇒ _" .

lemma finite_coeffs [simp]: "finite (coeffs p)"
by transfer simp

text ‹[1]p.82
  A "primitive'" polynomial has coefficients with GCD equal to 1.
  A polynomial is factored into "content" and "primitive part"
  for many different purposes.›

definition primitive :: "'a::{euclidean_ring,semiring_Gcd} mpoly ⇒ bool"
where
  "primitive p ⟷ Gcd (coeffs p) = 1"

definition content_primitive :: "'a::{euclidean_ring,GCD.Gcd} mpoly ⇒ 'a × 'a mpoly"
where
  "content_primitive p = (
    let d = Gcd (coeffs p)
    in (d, sdiv d p))"

value "let p = M [1,2,3] (4::int) + M [2,0,4] 6 + M [2,0,5] 8
  in content_primitive p"


end