Moduli of towers of powers

Stream: general

[Bolton Bailey]

Kenny Lau said:

it seems like we don't have the simp lemma Coprime a n -> a^b%n = a^(b%Phi(n))%n

Incidentally, I actually proved this while working on olympiad problems for Numina:

import Mathlib

open Nat

-- Nat.ModEq.pow_totient for raw modulus
lemma Nat.pow_totient_mod_eq_one {x n : ℕ} (hn : 1 < n) (h : x.Coprime n) : (x ^ n.totient) % n = 1 := by
  exact mod_eq_of_modEq (Nat.ModEq.pow_totient h) hn

lemma Nat.pow_add_totient_mod_eq_one {x k n : ℕ} (hn : 1 < n) (h : x.Coprime n) : (x ^ (k + n.totient)) % n = (x ^ k) % n := by
  rw [Nat.pow_add]
  rw [Nat.mul_mod]
  rw [Nat.pow_totient_mod_eq_one hn h]
  simp

lemma Nat.pow_add_mul_totient_mod_eq_one {x k l n : ℕ} (hn : 1 < n) (h : x.Coprime n) : (x ^ (k + l * n.totient)) % n = (x ^ k) % n := by
  induction l with
  | zero => simp
  | succ l ih =>
    rw [add_mul, one_mul]
    rw [<-add_assoc]
    rw [Nat.pow_add_totient_mod_eq_one hn h]
    exact ih

lemma Nat.pow_totient_mod {x k n : ℕ} (hn : 1 < n) (h : x.Coprime n) :
    x ^ k % n = x ^ (k % n.totient) % n := by
  let d := k / n.totient
  let r := k % n.totient
  have keq : k = d * n.totient + r := by exact Eq.symm (div_add_mod' k (φ n))
  have rle : r < n.totient := by
    refine mod_lt k ?_
    refine totient_pos.mpr ?_
    linarith
  rw [keq] at *
  rw [add_comm]
  simp [Nat.pow_add_mul_totient_mod_eq_one hn h]
  have : r % φ n = r := by
    exact mod_mod k (φ n)
  rw [this]

[Kenny Lau]

Please PR!

[Bolton Bailey]

(I was thinking that I should also write something that simplifies this in the case where a and n are not Coprime, and perhaps learn enough about simprocs to write one that could reduce a more general class of expressions of this type)

[Kenny Lau]

I would state it as a ^ (φ n + b) % n = a^(φ n + b % φ n) % n; and I would be more careful about tagging it with simp because by itself it's a loop and there has to be some way around this that I don't know yet

[Kenny Lau]

perhaps just something as simple as adding (h : φ n < b := by norm_num) as hypothesis

[Kenny Lau]

update: nope, that doesn't work, don't really know how to make it work

[Yaël Dillies]

Make it a simproc that only triggers when all values are numerals :wink:

[Kenny Lau]

@[simp] lemma Nat.pow_totient_mod' {x k n : ℕ} (hn : 1 < n := by simp)
    (h : x.Coprime n := by decide) (h₂ : n.totient < k := by simp [totient]) :
    x ^ k % n = x ^ (k % n.totient) % n :=
  Nat.pow_totient_mod hn h

lemma fooa : φ 10 = 4 := by decide -- yes!
lemma foob : φ 10 = 4 := by simp -- nope
lemma fooc : φ 10 = 4 := by norm_num -- nope

lemma foo1 : 3^3^3^3^3^3%10 = 4 := by
  simp
  rw [Nat.pow_totient_mod']

(attach to the code above)

[Kenny Lau]

what I'm thinking is something like this, but it doesn't work

[Kenny Lau]

we need to teach norm_num or simp how to evaluate φ, for example

[Kenny Lau]

Yaël Dillies said:

Make it a simproc that only triggers when all values are numerals :wink:

that wouldn't work, i want to trigger it for 3 ^ 3 ^ 3 ^ 7625597484987 % 10 = 4

[Kenny Lau]

do we need a custom bigpow tactic lol

[Aaron Liu]

We already have reduce_mod_char which does modular fastexp

[Kenny Lau]

does it work for my example?

[Aaron Liu]

Probably not, but let's see

[Bolton Bailey]

(We seem to have gotten somewhat off-topic, I'm going to try migrating this thread to a new topic.)

[Notification Bot]

16 messages were moved here from #general > How to formalize "compute this expression" exercises? by Bolton Bailey.

[Aaron Liu]

import Mathlib

example : 3 ^ 3 ^ 3 ^ 7625597484987 % 10 = 4 := by
  rw [show 4 = 4 % 10 from rfl, ← ZMod.natCast_eq_natCast_iff',
    Nat.cast_pow, Nat.cast_ofNat]
  -- takes forever
  reduce_mod_char

[Aaron Liu]

We could extend reduce_mod_char to recursively reduce the exponents too

[Kenny Lau]

maybe we need to define something like fpow : ZMod n -> ZMod (φ n) -> ZMod n (but with special support to deal with cases when the base is not coprime to n)

[Aaron Liu]

What would that do?

[Kenny Lau]

but the auxiliary problem we have is this: can someone write a tactic to prove 12345 < 3 ^ 3 ^ 3 ^ 7625597484987?

[Aaron Liu]

I feel like that's a separate problem

[Kenny Lau]

Aaron Liu said:

What would that do?

further reduce that to a more finite problem, so for example (using a+nZ to mean the element of Z/nZ), we would have ↑(3^100%10) = (3+10ℤ)^100 = (3+10ℤ)^(100+4ℤ)

[Kenny Lau]

Aaron Liu said:

I feel like that's a separate problem

it's really not, if the base is not coprime to n then it's only "semiperiodic"

[Kenny Lau]

semiperiodic means a->b->c->b->c->b->c...

[Kenny Lau]

like the decimal expansion of a rational number (is that a theorem? XD)

[Aaron Liu]

or the continued fraction of a root of a quadratic equation with integer coefficients (is that a theorem?)

[Aaron Liu]

A tactic for proving _ < _ ^ _ ^ ... would probably work by repeatedly logging both sides

[Aaron Liu]

using docs#Nat.lt_pow_iff_log_lt or similar

[Aaron Liu]

Are we assuming that the modulus is relatively small here?

[Bolton Bailey]

At the very least, we have to assume that we can compute n, totient n, totient (totient n), ... for about as many totients as there are exponents.

[Kenny Lau]

btw here is the 'manual' proof for reference:

lemma foo : 3^3^3^3^3^3%10 = 7 := by
  rw [pow_totient_mod, show φ 10 = 4 by decide,
    pow_totient_mod (n := 4), show φ 4 = 2 by decide,
    pow_totient_mod (n := 2), show φ 2 = 1 by decide,
    mod_one, pow_zero, show 1 % 2 = 1 by decide,
    show 3 ^ 1 = 3 by decide, show 3 % 4 = 3 by decide,
    show 3 ^ 3 % 10 = 7 by decide]

(attach to the code at the top)

[Aaron Liu]

What method are we using to compute totient n? Can we get in a @[csimp] lemma for it?

[Bolton Bailey]

Well, totient computation requires factorizing the number at hand. So if the modulus turns out to have large prime factors, we have an issue.

[Kenny Lau]

it's about time we had a tactic to factorise numbers :melt:

[Bolton Bailey]

I'm not too familiar with csimp, is there some mechanism where we can tell it to give up if it takes too long to factorize?

[Aaron Liu]

I did some experiments before and found that out of all the candidate implementations this one was the fastest. Unfortunately it was also the hardest to prove stuff about and I was never able to show it equals totient.

import Mathlib

def Nat.fastTotient (n : Nat) : Nat :=
  if n = 0 then 0 else
  (n.primeFactorsList.foldl
    (fun (prev, acc) curr =>
      if curr = prev then (curr, acc * curr) else (curr, acc * (curr - 1)))
    (0, 1)).2

[Aaron Liu]

Actually, maybe it was the other one that was faster? I forgot. I'll have to experiment again.

[Eric Wieser]

Kenny Lau said:

it's about time we had a tactic to factorise numbers :melt:

I think simp can do it? docs#Nat.primeFactorsList_ofNat

[Aaron Liu]

Yes, this one seems to have been the fastest:

import Mathlib

def Nat.fastTotient (n : Nat) : Nat :=
  (n.primeFactorsList.foldl
    (fun (prev, acc) curr =>
      if curr = prev then (curr, acc) else (curr, acc * (curr - 1) / curr))
    (0, n)).2

[Bolton Bailey]

Kenny Lau said:

Please PR!

I haven't had time to do the non-coprime case, but I PRed these lemmas just now #26694

[Bolton Bailey]

Kenny Lau said:

I would state it as a ^ (φ n + b) % n = a^(φ n + b % φ n) % n; and I would be more careful about tagging it with simp because by itself it's a loop and there has to be some way around this that I don't know yet

The thing that makes me nervous about using this approach to simplify (b ^ e) % n is: How do we determine, if e is itself a potentially complicated expression involving exponents, whether n < e. This seems necessary so that we know whether to apply the lemma or just evaluate b ^ e normally, but I'm not sure how to do it without introducing a more complicated decision procedure.

Might another approach be: Determine the factorization of n, use CRT to reduce the problem of computing (b ^ e) % (p ^ a) for all p^a in the factorization, reexpress b = p^k * b' and then evaluate (p^k ^ e) % (p^a) and b' % (p ^ a) separately?

[Aaron Liu]

I could write a tactic

[Bolton Bailey]

What I really want is for @Yaël Dillies and @Paul Lezeau to publish more of their simprocs for the working mathematician series so we can make this a simproc :grinning:

[Aaron Liu]

A simproc might be challenging

[Aaron Liu]

since you have to do it somewhat recursively

[Bolton Bailey]

Nevertheless, I think this feels more appropriate as a simproc, because I would like it to automatically happen when someone calls simp without anyone having to learn the name and scope of a new tactic.

[Aaron Liu]

It feels somewhat outside the scope of simp

[Aaron Liu]

simp isn't for computing numeric expression, that's norm_num

[Aaron Liu]

And norm_num doesn't really work on this either, since you have to take into account a modulus

[Bolton Bailey]

Can simprocs be added to norm_num?

[Aaron Liu]

norm_num calls simp

[Aaron Liu]

It's also extensible

[Aaron Liu]

So you can register a "norm_num extension" to whatever operation you want

[Aaron Liu]

For example in #26317 I added one for Nat.clog

[Bolton Bailey]

Right, ok, this norm_num extension thing sounds like the most-in-scope home for something like this.

[Aaron Liu]

But you don't get a modulus to work with

[Aaron Liu]

Something like reduce_mod_char seems like a much better fit for this kind of simplification, but reduce_mod_char is kind of broken and just about the only thing it's good for is fast modexp

[Bolton Bailey]

Not sure what you mean, is it not possible to write a norm_num extension that does something like "check if the term looks like (b ^ <expression of literals> % n) and use this approach if so"?

[Aaron Liu]

It is possible, but I don't think it's a good idea

[Aaron Liu]

simp goes bottom up anyways, so if you run norm_num you'll probably end up simplifying the inside first

[Aaron Liu]

But I wouldn't know since I haven't tried it

[Bolton Bailey]

If I have an expression like this and I think "I want to normalize this into an explicit number", I think that the next thing I am likely to think is "maybe this is something that norm_num does" and try it.

Frankly, looking at reduce_mod_char I am not sure why this exists as a separate tactic from norm num, couldn't the functionality of reduce_mod_char be folded into norm_num?

bottom up

I though there were ways to counteract this, like the little down-arrow version you see on reduceIte

[Aaron Liu]

Bolton Bailey said:

Frankly, looking at reduce_mod_char I am not sure why this exists as a separate tactic from norm num, couldn't the functionality of reduce_mod_char be folded into norm_num?

According to some comments around the implementation for reduce_mod_char:

This is not a norm_num plugin because it does not match on the syntax of e,
rather it matches on the type of e.

[Paul Lezeau]

Bolton Bailey said:

What I really want is for Yaël Dillies and Paul Lezeau to publish more of their simprocs for the working mathematician series so we can make this a simproc :grinning:

Yaël and I have both been distracted with other stuff lately but the plan is to try and finish another one up over the current month or so!

[Yaël Dillies]

Aaron Liu said:

since you have to do it somewhat recursively

That's no issue :cowboy:

[Aaron Liu]

Yaël Dillies said:

Aaron Liu said:

since you have to do it somewhat recursively

That's no issue :cowboy:

Can you explain

[Yaël Dillies]

Aaron Liu said:

simp goes bottom up anyways, so if you run norm_num you'll probably end up simplifying the inside first

That is factually incorrect, or at least requires a few disclaimers

[Aaron Liu]

I'll take a look at the code and see what it does

[Aaron Liu]

oh it does it twice?

[Aaron Liu]

oh there's pre and post norm_num extensions

[Aaron Liu]

that makes sense

[Yaël Dillies]

Aaron Liu said:

Can you explain

You can read the code of docs#existsAndEq or wait for Paul and I to publish the blogpost about it