Keeping only the primes with smaller valuation

Stream: Is there code for X?

[Riccardo Brasca]

Let n and m be natural numbers, written as n = ∏ pᵢ^nᵢand m = ∏ pᵢ^mᵢ (product over the primes). What is the best way of writing the number ∏ pᵢ^kᵢ, where kᵢ = nᵢ if nᵢ < mᵢ and 0 otherwise? I can do by hand, but I feel it's better to use docs#Nat.factorization in a clever way.

For contest: I am proving that if x and y are commuting elements in a group with orderOf x = n and orderOf y = m then there is an element of order lcm n m: it's possible to give the actual element, but it involves the number above.

[Kevin Buzzard]

Re the context: my guess is that you can use docs#Nat.xgcd and docs#Nat.gcdA etc to do this without having to factor anything.

[Kevin Buzzard]

(my instinct is that the element you're after is x^gcdA m n * y^gcdB m n or something like that: note that I didn't check anything here)

[Riccardo Brasca]

Let me think about that after lunch. The point is writing the lcm as the product of a divisor of m and a divisor of n, but the two divisors must be coprime.

[Riccardo Brasca]

I don't see how to make something like this work: the only case where we can compute the order of the product is when the orders are coprime, so I have the impression bezout is useless here.

[Johan Commelin]

I think we don't have a function

def minL (m n : Nat) : Nat :=
  if m < n then m else 0

[Johan Commelin]

Otherwise you could map that over the exponents

[Riccardo Brasca]

Yeah, I guess that doing it by hand is not that bad.

[Kevin Buzzard]

My instinct is still very much that one doesn't ever need to factorise anything and that your suggested approach is going to be more painful than a direct approach.

If $L$ and $G$ are the LCM and GCD of $m,n$ then $x':=x^{n/G}$ and $y':=y^{m/G}$ both have order $G$, so $y'=x'^u$ for some $u$ coprime to $G$. Then something like $x^uy$ should work? Even if I've got the details wrong (which I could well have, I'm in a bit of a rush), something like this seems to me to be a far less painful approach than factoring anything.

[Antoine Chambert-Loir]

Each prime gives you forbidden congruences for $u$: if the exponent of $p$ in $m$ is greater than that of $m$, $p$ should not divide $u$, if it's equal, one value is forbidden mod $p$ to the exponent, otherwise it's ok. How can you get $u$ that satisfies all of these without factoring?

[Michael Stoll]

Maybe one can avoid using factorizations explicitly by going via docs#induction_on_primes ?
(Just a vague idea; I didn't actually try it.)

[Riccardo Brasca]

Kevin Buzzard said:

My instinct is still very much that one doesn't ever need to factorise anything and that your suggested approach is going to be more painful than a direct approach.

If $L$ and $G$ are the LCM and GCD of $m,n$ then $x':=x^{n/G}$ and $y':=y^{m/G}$ both have order $G$, so $y'=x'^u$ for some $u$ coprime to $G$. Then something like $x^uy$ should work? Even if I've got the details wrong (which I could well have, I'm in a bit of a rush), something like this seems to me to be a far less painful approach than factoring anything.

...both have order $G$, so $y'=x'^u$ for some $u$ coprime to $G$. How can you find a relation between x and y? Note that I am not in a cyclic group (where the proof would be easier).

[Riccardo Brasca]

Michael Stoll said:

Maybe one can avoid using factorizations explicitly by going via docs#induction_on_primes ?
(Just a vague idea; I didn't actually try it.)

It can be useful, but only to prove the existence, and I would like to give the precise element.

I guess I will go with the factorization, let's see how painful it is.

[Riccardo Brasca]

Well, it's quite painful :sweat_smile:

[Riccardo Brasca]

I am not sure what induction principle to use to define the function, it seems we don't have any way of producing a function ℕ → ℕ → ℕ using some induction related to the factorization

[Kevin Buzzard]

Aah I was only thinking about the cyclic case! I was imagining that this was a question coming from my comment on your PR about roots of unity, and I had assumed we were in a cyclic group, but I see now that you don't need cyclicity for the result.

[Riccardo Brasca]

Yes, it comes from that PR, but since I am at this I would like to do the general case

[Kevin Buzzard]

Yeah then I retract my comments about it all being Bezout in disguise. Sorry for the noise!

[Junyan Xu]

This may be easier: for each prime pᵢ you can get an element of order pᵢ^max nᵢ mᵢ; it will be a power of x or a power of y depending on whether nᵢ or mᵢ is larger. If you take the product over all primes you get an element with desired order lcm n m (should be easy to prove by induction). (You'd also need to prove commutativity by induction, but that seems manageable.)

[Riccardo Brasca]

Yes, this can work. But now I am a little angry against Lean that I am not able to do the other approach, that is very direct mathematically, so I would really like to make it work...

[Junyan Xu]

Hmm it doesn't actually work because docs#Finsupp.prod requires CommMonoid ... we'd need an analogue of docs#Multiset.noncommProd or use docs#List.prod.

[Riccardo Brasca]

I really think this should not be hard:
given a b : ℕ, write Nat.lcm a b = m * n where:

• m ∣ a
• n ∣ b
• m and n are coprime.

[Riccardo Brasca]

We even have an explicit formula:

import Mathlib

open Nat

variable (a b : ℕ)

def fst_big :=
    (a.primeFactors.filter (fun p ↦ b.factorization p ≤ a.factorization p)).prod
      (fun p ↦ p ^ a.factorization p)

def snd_big :=
    (b.primeFactors.filter (fun p ↦ a.factorization p < b.factorization p)).prod
      (fun p ↦ p ^ b.factorization p)

variable {a b}

lemma foo (ha : a ≠ 0) (hb : b ≠ 0) : (fst_big a b) * (snd_big a b) = Nat.lcm a b := by
  revert hb
  refine recOnPosPrimePosCoprime (fun p n hp hn _ ↦ ?_) ?_ ?_ ?_ b
  · refine recOnPosPrimePosCoprime (fun q m hq hm ↦ ?_) ?_ ?_ ?_ a
    by_cases hpq : p = q
    · subst hpq
      by_cases hnm : n ≤ m
      · have : Nat.lcm (p ^ m) (p ^ n) = p ^ m := by
          obtain ⟨k, hk⟩ := pow_dvd_pow p hnm
          rw [← mul_one (p ^ n), hk, Nat.lcm_mul_left]
          simp
        simp_all [fst_big, primeFactors_prime_pow hm.ne', snd_big, primeFactors_prime_pow hn.ne']
      · have : Nat.lcm (p ^ m) (p ^ n) = p ^ n := by
          obtain ⟨k, hk⟩ := pow_dvd_pow p <| le_of_not_ge hnm
          rw [← mul_one (p ^ m), hk, Nat.lcm_mul_left]
          simp
        simp_all [fst_big, primeFactors_prime_pow hm.ne', snd_big, primeFactors_prime_pow hn.ne']
    · sorry
    · sorry
    · sorry
    · sorry
  · sorry
  · sorry
  · sorry

this is surely going to work, but it's very brutal.

[Riccardo Brasca]

For example

simp_all [fst_big, primeFactors_prime_pow hm.ne', snd_big, primeFactors_prime_pow hn.ne',
        Ne.symm hpq, Coprime.lcm_eq_mul <| coprime_pow_primes m n hq hp (Ne.symm hpq)]

proves the first sorry

[Arend Mellendijk]

This reminds me of a problem I ran into with arithmetic functions where the lcm didn't behave well. Would an argument like in #7787 work better than going through induction? This feels pretty promising:

import Mathlib

open Nat

variable (a b : ℕ)

def fst_big :=
    ((a*b).primeFactors.filter (fun p ↦ b.factorization p ≤ a.factorization p)).prod
      (fun p ↦ p ^ a.factorization p)

def snd_big :=
    ((a*b).primeFactors.filter (fun p ↦ a.factorization p < b.factorization p)).prod
      (fun p ↦ p ^ b.factorization p)

variable {a b}

lemma foo (ha : a ≠ 0) (hb : b ≠ 0) : (fst_big a b) * (snd_big a b) = Nat.lcm a b := by
  rw [fst_big, snd_big]
  have hab : a.lcm b ≠ 0 := by exact lcm_ne_zero ha hb
  erw [Nat.multiplicative_factorization (fun x ↦ x) _ _ hab,
    Finsupp.prod_of_support_subset _ _ _ (fun _ _ => sorry)
    (s:=(a*b).primeFactors), Finset.prod_filter, Finset.prod_filter,
    ← Finset.prod_mul_distrib]

  sorry
  · sorry
  · intros; rfl
  · rfl

[Arend Mellendijk]

(note I slightly changed the definitions to make the argument smoother)

[Riccardo Brasca]

Ah, using (a*b).primeFactors is surely a good idea, thanks!

[Riccardo Brasca]

But I really don't see how you can avoid using some kind of induction with the primes.

[Arend Mellendijk]

Like this? (this needs some cleaning up, of course)

import Mathlib

open Nat BigOperators

variable (a b : ℕ)

def fst_big :=
    ((a*b).primeFactors.filter (fun p ↦ b.factorization p ≤ a.factorization p)).prod
      (fun p ↦ p ^ a.factorization p)

def snd_big :=
    ((a*b).primeFactors.filter (fun p ↦ a.factorization p < b.factorization p)).prod
      (fun p ↦ p ^ b.factorization p)

variable {a b}

lemma fst_big_eq :
    fst_big a b = (∏ p in (a*b).primeFactors, p ^ (if b.factorization p ≤ a.factorization p then a.factorization p else 0)) := by
  rw [fst_big, Finset.prod_filter]
  congr with p
  split_ifs <;> simp

lemma snd_big_eq :
    snd_big a b = (∏ p in (a*b).primeFactors, p ^ (if a.factorization p < b.factorization p then b.factorization p else 0)) := by
  rw [snd_big, Finset.prod_filter]
  congr with p
  split_ifs <;> simp

lemma foo (ha : a ≠ 0) (hb : b ≠ 0) : (fst_big a b) * (snd_big a b) = Nat.lcm a b := by
  rw [fst_big_eq, snd_big_eq]
  have hab : a.lcm b ≠ 0 := by exact lcm_ne_zero ha hb
  erw [Nat.multiplicative_factorization (fun x ↦ x) _ _ hab,
    Finsupp.prod_of_support_subset _ _ _ (fun _ _ ↦ rfl)
    (s:=(a*b).primeFactors),
    ← Finset.prod_mul_distrib]
  simp_rw [←pow_add]
  apply Finset.prod_congr rfl
  intro p _
  congr
  rw [Nat.factorization_lcm ha hb]
  simp only [Finsupp.sup_apply]
  by_cases h : Nat.factorization b p ≤ Nat.factorization a p
  · rw [if_pos h, if_neg (Nat.not_lt.mpr h), sup_of_le_left h, add_zero]
  · rw [if_neg h, if_pos (Nat.not_le.mp h),
      sup_of_le_right (Nat.le_of_not_ge h), zero_add]
  · rw [Nat.factorization_lcm ha hb, Finsupp.support_sup, primeFactors_mul ha hb, Nat.support_factorization , Nat.support_factorization]
  · intros; rfl
  · rfl

[Riccardo Brasca]

Wow, thanks!! You can even change the definition to what fst_big_eq gives (I plan to use a notation)

[Michael Stoll]

Here is an attempt using prime-induction. (Two sorries left, but I have to cook dinner...)

import Mathlib.RingTheory.Int.Basic

lemma foo (a b : ℕ) : ∃ m n, Nat.lcm a b = m * n ∧ m ∣ a ∧ n ∣ b ∧ m.Coprime n := by
  induction' a using induction_on_primes with p a hp ih
  · exact ⟨0, 1, by simp, dvd_rfl, one_dvd _, rfl⟩
  · refine ⟨1, b, by simp, dvd_rfl, dvd_rfl, Nat.gcd_one_left b⟩
  · obtain ⟨m, n, h₁, h₂, h₃, h₄⟩ := ih
    by_cases H : p ∣ m
    · have : ¬ p ∣ n := by
        contrapose! h₄
        exact Nat.Prime.not_coprime_iff_dvd.mpr ⟨p, hp, H, h₄⟩
      refine ⟨p * m, n, dvd_antisymm ?_ ?_, Nat.mul_dvd_mul_left p h₂, h₃, ?_⟩
      · rw [mul_assoc, ← h₁]
        exact Nat.lcm_dvd (Nat.mul_dvd_mul_left p <| Nat.dvd_lcm_left a b) <| Dvd.dvd.mul_left h₂ p
      · sorry
      · exact Nat.Coprime.mul ((Nat.Prime.coprime_iff_not_dvd hp).mpr this) h₄
    · refine ⟨m, n, dvd_antisymm ?_ ?_, Dvd.dvd.mul_left h₂ p, h₃, h₄⟩
      · rw [← h₁]
        refine Nat.lcm_dvd ?_ <| Nat.dvd_lcm_right a b
        sorry
      · exact h₁ ▸ Nat.lcm_dvd ((Nat.dvd_mul_left a p).trans (Nat.dvd_lcm_left (p * a) b) <|
          Nat.dvd_lcm_right ..

[Junyan Xu]

Another way to define fst_big / snd_big:

import Mathlib

open Nat

variable (a b : ℕ)

def fst_big : ℕ := (Nat.lcm a b).factorization.prod fun p n ↦
  if b.factorization p ≤ a.factorization p then p ^ n else 1

def snd_big : ℕ := (Nat.lcm a b).factorization.prod fun p n ↦
  if b.factorization p ≤ a.factorization p then 1 else p ^ n

variable (ha : a ≠ 0) (hb : b ≠ 0)

lemma foo : (fst_big a b) * (snd_big a b) = Nat.lcm a b := by
  conv_rhs => rw [← factorization_prod_pow_eq_self (lcm_ne_zero ha hb)]
  rw [fst_big, snd_big, ← Finsupp.prod_mul]
  congr; ext p n; split_ifs <;> simp

lemma fst_big_dvd : fst_big a b ∣ a := by
  conv_rhs => rw [← factorization_prod_pow_eq_self ha]
  rw [Finsupp.prod_of_support_subset (s := (Nat.lcm a b).factorization.support)]
  · apply Finset.prod_dvd_prod_of_dvd; rintro p -; dsimp only; split_ifs with le
    · rw [factorization_lcm ha hb]; apply pow_dvd_pow; exact sup_le le_rfl le
    · apply one_dvd
  · intro p hp; rw [Finsupp.mem_support_iff] at hp ⊢
    rw [factorization_lcm ha hb]; exact (lt_sup_iff.mpr <| .inl <| Nat.pos_of_ne_zero hp).ne'
  · intros; rw [pow_zero]

[Riccardo Brasca]

This is really a wonderful community, now I only have to choose what is my preferred answer :)

[Riccardo Brasca]

The fact that the proof is two lines long

import Mathlib

open Nat

variable (a b : ℕ)

def fst_big : ℕ := (Nat.lcm a b).factorization.prod fun p n ↦
  if b.factorization p ≤ a.factorization p then p ^ n else 1

def snd_big : ℕ := (Nat.lcm a b).factorization.prod fun p n ↦
  if b.factorization p ≤ a.factorization p then 1 else p ^ n

variable (ha : a ≠ 0) (hb : b ≠ 0)

lemma foo : (fst_big a b) * (snd_big a b) = Nat.lcm a b := by
  rw [← factorization_prod_pow_eq_self (lcm_ne_zero ha hb), fst_big, snd_big, ← Finsupp.prod_mul]
  congr; ext p n; split_ifs <;> simp

seems a very good sign. Let me see how to prove that they're coprime.

[Junyan Xu]

Maybe I've shuffled all difficulties into proving coprimality ... I now think we should make use of docs#Nat.factorizationEquiv

[Riccardo Brasca]

I have to stop for today and I will not have time tomorrow, but if any of you want to PR this I will be very happy to review. (Note that having this result about Nat will make the result about the order completely trivial)

[Riccardo Brasca]

Anyway it seems that by definition, for any prime prime p, the valuation of one of the two numbers is 1, and I guess we know (in some form) that this implies coprimality

[Riccardo Brasca]

Ops, I confused being equal to 1 and having valuation 0

[Michael Stoll]

Upon further thought, I think my approach does not work the way I tried it (at some point it may be necessary to move the maximal p-power from n to m).

[Michael Stoll]

(BTW, the efficient way to implement this would be to use coprime factorization : given positive natural numbers $a$ and $b$, one can efficiently compute pairwise coprime numbers $p_1, \ldots, p_k$ greater than $1$ such that $a$ and $b$ can be written as power products of the $p_j$.)

[Riccardo Brasca]

I think @Junyan Xu 's approach is a good one, maybe going through docs#Nat.factorizationEquiv_inv_apply to compute the factorization of the two numbers

[Michael Stoll]

This should be generalizable to GCD monoids or at least factorization monoids, I assume?

[Junyan Xu]

Coprimality is actually not hard:

lemma coprime : (fst_big a b).Coprime (snd_big a b) := by
  rw [fst_big, snd_big]
  refine coprime_prod_left_iff.mpr fun p hp ↦ coprime_prod_right_iff.mpr fun q hq ↦ ?_
  dsimp only; split_ifs with h h'
  any_goals apply coprime_one_left
  · apply coprime_one_right
  refine coprime_pow_primes _ _ (prime_of_mem_primeFactors hp) (prime_of_mem_primeFactors hq) ?_
  contrapose! h'; rwa [← h']

I tried to use

def fst_big : ℕ+ := factorizationEquiv.symm
  ⟨(Nat.lcm a b).factorization.filter (fun p ↦ b.factorization p ≤ a.factorization p),
    fun _ hp ↦ prime_of_mem_primeFactors (Finset.filter_subset _ _ hp)⟩

but it doesn't simplify things (and it's annoying that Lean can't rewrite the definition using docs#Nat.factorizationEquiv_inv_apply because it uses Subtype.val rather than PNat.val / coercion ... we should fix the lemma)

[Riccardo Brasca]

Great!

[Junyan Xu]

Michael Stoll said:

This should be generalizable to GCD monoids or at least factorization monoids, I assume?

I can prove it for UFDs but I don't know if it generalizes to GCD monoids:

import Mathlib

variable {α} [CancelCommMonoidWithZero α] [UniqueFactorizationMonoid α]

lemma IsRelPrime.exists_isRelPrime_mul_eq {a b : α} (h : IsRelPrime a b) (c : α) :
    c ≠ 0 → ∃ ca cb, IsRelPrime ca a ∧ IsRelPrime cb b ∧ IsRelPrime ca cb ∧ ca * cb = c := by
  refine UniqueFactorizationMonoid.induction_on_coprime c (fun h ↦ (h rfl).elim) ?_ ?_ ?_
  · rintro x hu -
    exact ⟨1, x, isRelPrime_one_left, hu.isRelPrime_left, isRelPrime_one_left, one_mul _⟩
  · rintro p i hp -
    obtain dvd | hpa := hp.irreducible.dvd_or_isRelPrime (n := a)
    · exact ⟨1, p ^ i, isRelPrime_one_left, (h.of_dvd_left dvd).pow_left, isRelPrime_one_left, one_mul _⟩
    · exact ⟨p ^ i, 1, hpa.pow_left, isRelPrime_one_left, isRelPrime_one_right, mul_one _⟩
  intro x y hp hx hy h0
  obtain ⟨xa, xb, hxa, hxb, hxab, rfl⟩ := hx (left_ne_zero_of_mul h0)
  obtain ⟨ya, yb, hya, hyb, hyab, rfl⟩ := hy (right_ne_zero_of_mul h0)
  simp_rw [IsRelPrime.mul_left_iff, IsRelPrime.mul_right_iff] at hp
  refine ⟨xa * ya, xb * yb, hxa.mul_left hya, hxb.mul_left hyb,
    (hxab.mul_right hp.1.2).mul_left (hp.2.1.symm.mul_right hyab), by rw [mul_mul_mul_comm]⟩

theorem extract_gcd' {α} [CancelCommMonoidWithZero α] [GCDMonoid α] (x y : α) :
    ∃ x' y', x = gcd x y * x' ∧ y = gcd x y * y' ∧ IsRelPrime x' y' ∧
      Associated (lcm x y) (gcd x y * x' * y') := by
  obtain ⟨x', y', hx, hy, hu⟩ := extract_gcd x y
  rw [gcd_isUnit_iff_isRelPrime] at hu
  refine ⟨x', y', hx, hy, hu, ?_⟩
  obtain rfl | h0 := eq_or_ne x 0
  · rw [lcm_zero_left, ← hx, zero_mul]
  refine .of_mul_left ((gcd_mul_lcm x y).trans ?_) (.refl _) ?_
  · rw [← hx, mul_left_comm, ← hy]
  · exact fun h ↦ h0 ((gcd_eq_zero_iff x y).mp h).1

lemma UniqueFactorizationMonoid.exists_isRelPrime_mul_associated_lcm [GCDMonoid α] (a b : α) :
    ∃ ca cb, ca ∣ a ∧ cb ∣ b ∧ IsRelPrime ca cb ∧ ca * cb = lcm a b := by
  obtain ⟨x, y, hx, hy, hp, assoc⟩ := extract_gcd' a b
  obtain ⟨u, eq_lcm⟩ := assoc.symm
  obtain rfl | h0 := eq_or_ne a 0
  · exact ⟨0, 1, dvd_zero _, one_dvd _, isRelPrime_one_right,
      (zero_mul _).trans (lcm_zero_left _).symm⟩
  obtain ⟨ca, cb, ha, hb, hab, eq⟩ := hp.exists_isRelPrime_mul_eq _
    fun h ↦ h0 ((gcd_eq_zero_iff a b).mp h).1
  refine ⟨cb * x, ca * y * u, ⟨ca, ?_⟩, ⟨u⁻¹ * cb, ?_⟩, ?_, ?_⟩
  · rw [hx, ← eq, mul_assoc, mul_comm]
  · rw [hy, ← eq]; simp_rw [mul_assoc, u.mul_inv_cancel_left, mul_comm]
  · rw [isRelPrime_mul_unit_right_right u.isUnit]
    exact (hab.mul_right ha).symm.mul_right (hb.mul_left hp)
  · simp_rw [← eq_lcm, ← eq, ← mul_assoc]; congr 2; rw [mul_comm _ ca, mul_assoc]

[Junyan Xu]

asked as https://math.stackexchange.com/questions/4875914/factorization-into-coprimes-subordinate-to-two-given-coprimes-in-gcd-domain

[Riccardo Brasca]

We want the version for N anyway: everything is computable there.

[Junyan Xu]

BTW is there a tactic to close equality goals in commutative monoids such as those in the lines

  · rw [hx, ← eq, mul_assoc, mul_comm]
  · rw [hy, ← eq]; simp_rw [mul_assoc, u.mul_inv_cancel_left, mul_comm]

? I tried group and abel but they don't work, and ring requires addition.

[Riccardo Brasca]

I think abel is supposed to do it.

import Mathlib

example (M : Type*) [AddCommMonoid M] (ca cb x : M) : ca + cb + x = cb + x + ca := by
  abel

this works, but not if I remove the Add

[Riccardo Brasca]

Anyway, do you want to PR these results?

[Riccardo Brasca]

Indeed the docstring for abel says

/-!
# The `abel` tactic

Evaluate expressions in the language of additive, commutative monoids and groups.

-/

[Riccardo Brasca]

@Junyan Xu do you mind if I PR your work for N? I need it for an application to cyclotomic fields.

[Junyan Xu]

Please always feel free to!

[Riccardo Brasca]

Thanks!

[Riccardo Brasca]

#11235