# Theory Polynomials.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›

"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

instantiation mpoly :: (zero) zero
begin

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

instance ..

end

begin

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

instance

end

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

begin

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

instance

end

begin

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

instance

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",
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

"monom m (a + b) = monom m a + monom m b"

lemma monom_uminus:
"monom m (- a) = - monom m a"

lemma monom_diff:
"monom m (a - b) = monom m a - monom m b"

lemma monom_numeral [simp]:
"monom 0 (numeral n) = numeral n"

lemma monom_of_nat [simp]:
"monom 0 (of_nat n) = of_nat n"

lemma of_nat_monom:
"of_nat = monom 0 ∘ of_nat"

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)"

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"

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"

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)

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 (rule finite_mult_not_eq_zero_rightI)
apply simp
apply (rule finite_mult_not_eq_zero_rightI)
apply simp
done

"insertion f (p + q) = insertion f p + insertion f q"

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"
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"
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)"
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)
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))"
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))"
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"

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"

lemma mult_smult_left: "smult s p * q = smult s (p * q)"

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)"

lemma pmod_code:
"p pmod q = snd (pdivmod p q)"

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
```