Jacobi symbol

Stream: maths

[Michael Stoll]

First of all, thanks to everybody @Johan Commelin @Riccardo Brasca @Junyan Xu @Kevin Buzzard for your help with additive and multiplicative characters, Gauss sums, and the Legendre symbol!
I have now written code that defines the Jacobi symbol and proves its most important properties.

[Kevin Buzzard]

This is great. Do you think cubic reciprocity is in sight? Could it be an MSc project?

[Michael Stoll]

The PR is #16290.

[Michael Stoll]

I had to make some design decisions that I would like to mention here.
The Jacobi symbol is defined as a product over Legendre symbols (running over the prime factors of the number below). I am basically using nat.factor and a product over the resulting list. However, since the Legendre symbol (and the relevant properties we want to use) requires its second (actually the first, as implemented) argument to be a prime, I need to use a modified version of nat.factors, whose elements know that they are primes. For this purpose, I have introduced a variant nat.prime_factors with values in list primes. I had to duplicate part of the API for nat.factors to be able to use it easily.

[Michael Stoll]

Given this, the definition of the Jacobi symbol is simply

def jacobi_sym (a : ℤ) (b : ℕ) : ℤ :=
(b.prime_factors.map (λ (p : primes), @legendre_sym p.1 ⟨p.2⟩ a)).prod

Note that this definition returns something for all natural numbers b, even though the definition really only makes sense when b is odd. This requirement shows up later (in the form of an argument (hb : odd b)) in some of the lemmas.

[Riccardo Brasca]

Thanks!! A first trivial comment: the PR is pretty long, I think you can split if, PRing first of all the material about nat.prime_factors, maybe in another file :)

[Michael Stoll]

I introduce the notation [a | b]ⱼ for the symbol (localized to number_theory_symbols, which I anticipate will later also contain notation for the Hilbert symbol). Unfortunately, there is no upper-case subscipt J available (at least not via \_J in VSCode). (I noticed earlier when playing around with notation for the Hilbert symbol that round parentheses don't seem to work well, probably because they are used for a lot of other notation.)

[Michael Stoll]

For the proofs, I first show that the symbol is multiplicative in the second argument (which more or less follows directly from the definition) and then I use "multiplicative induction" via nat.rec_on_mul. For some of the proofs (equalities with multiplicative right hand side), one can abstract this into a custom proof principle (jacobi_sym_value), which results in one-line proofs of [-1 | b}ⱼ = χ₄ b and the similar versions for 2 and -2 and also in a fairly short proof of quadratic reciprocity.
In a number of the earlier proofs, it seems that I have to spell out the proof each time again, because the statements are too different in shape. (This contributes somewhat to the length of the file.) Ideas for simplifications are appreciated!

[Michael Stoll]

Riccardo Brasca said:

Thanks!! A first trivial comment: the PR is pretty long, I think you can split if, PRing first of all the material about nat.prime_factors, maybe in another file :)

I think the material on nat.prime_factors should be moved to data.nat.prime anyway, but I wanted to wait to see if more API is needed before making the move.
One could split the material on the Jacobi symbol into two parts: one containing the sections on "Definition" and "Properties" and one containing the sections on "Values" and "Quadratic Reciprocity". Would that be helpful?

[Michael Stoll]

Kevin Buzzard said:

This is great. Do you think cubic reciprocity is in sight? Could it be an MSc project?

Possibly. You'd have to do cubic characters, cubic residue symbols on Eisenstein integers (I don't know how much of the theory of Eisenstein integers is already available), Jacobi sums (which are a likely topic for my seminar, but I'll have to see how many students will actually be interested), at least this is what I see when I flip through Ireland&Rosen.

[Alex J. Best]

Did you try using docs#list.attach instead of defining a new function nat.prime_factors? That should allow you to use the fact that a member of prime factors is prime whenever needed

[Ruben Van de Velde]

There's already about seven ways to talk about the prime factors of a natural, adding more doesn't sound that great

[Michael Stoll]

Alex J. Best said:

Did you try using docs#list.attach instead of defining a new function nat.prime_factors? That should allow you to use the fact that a member of prime factors is prime whenever needed

nat.prime_factors is defined using list.attach. Using it directly seems to be a pain -- if l1 = l2, you cannot rewrite l1.attach into l2.attach, since l1.attach = l2.attach doesn't even type-check...

[Michael Stoll]

Ruben Van de Velde said:

There's already about seven ways to talk about the prime factors of a natural, adding more doesn't sound that great

Can you be more precise? I.e., which ways are there, and which one do you think could be used here most easily?

[Michael Stoll]

I have now shortened the lengthy proofs somewhat and also added text to the module docstring.

[Junyan Xu]

Michael Stoll said:

nat.prime_factors is defined using list.attach. Using it directly seems to be a pain -- if l1 = l2, you cannot rewrite l1.attach into l2.attach, since l1.attach = l2.attach doesn't even type-check...

Can you use docs#list.pmap? There is docs#list.pmap_congr.

[Michael Stoll]

I'll try to see if docs#list.pmap can be made to work in a reasonable way.
On the other hand, I would think that having a function prime_factors : ℕ → list primes is quite natural. After all, what is the type primes for if not for an application like this?

[Yakov Pechersky]

Why do you need a list? Why can't you do (n.factors.map _).prod _?

[Michael Stoll]

n.factors is a list, so it is not clear to me what you mean here.

[Yakov Pechersky]

Sorry, I was confused with the UFD factors

[Michael Stoll]

I cannot simply do (n.factors.map (λ p, legendre_sym p a)).prod, because legendre_sym needs to know that p is a prime.

[Michael Stoll]

With docs#list.pmap, the definition would be

def jacobi_sym' (a : ℤ) (b : ℕ) : ℤ :=
(b.factors.pmap (λ (p : ℕ) (pp : p.prime), @legendre_sym p ⟨pp⟩ a)
                (λ p (pf : p ∈ b.factors), prime_of_mem_factors pf)).prod

which is a bit more involved. I'll now try to use it...

[Alex J. Best]

I'm not sure what the type primes is for to be honest, it doesn't have too many meaningful algebraic structures on it, so it does seem easier to stick to the is_prime predicate where possible, and use subtypes when needed

[Michael Stoll]

(primes is a subtype...)

[Yakov Pechersky]

import number_theory.legendre_symbol.quadratic_reciprocity

open nat zmod

/-- The Jacobi symbol of `a` and `b` -/
def jacobi_sym (a : ℤ) (b : ℕ) : ℤ :=
(b.factors.pmap (λ n hn, @legendre_sym n ⟨hn⟩ a) (λ n hn, (prime_of_mem_factors hn : n.prime))).prod

[Alex J. Best]

Exactly, so giving it a name is just adding some extra layer of obfuscation to applying lemmas about subtype to it

[Yakov Pechersky]

Ah I see you got it

[Michael Stoll]

I'm missing docs#list.pmap_append (the analogue of docs#list.map_append). Is it (or something equivalent) available, perhaps under a different name?

[Michael Stoll]

Apart from this (which is needed for the fundamental multiplicativity in the second argument), the alternative definition seems to work well.

[Sebastien Gouezel]

Michael Stoll said:

I cannot simply do (n.factors.map (λ p, legendre_sym p a)).prod, because legendre_sym needs to know that p is a prime.

To me, this looks like an indication that the definition of legendre_sym is not good: it should be 1 (or 0, depending on what works best) if p is not prime, and the usual thing when p is prime.

[Michael Stoll]

That would introduce ites, and anyway, the problem would only be shifted, since all relevant properties of the Legendre symbol are only valid when p is prime.

[Sebastien Gouezel]

That's the concept of junk value: introduce total functions, even at places where they are not well defined mathematically (just like division by zero is defined as a total function), and prove the lemmas assuming the additional properties when needed. So many lemmas about legendre_sym would indeed have an assumption prime p, but the definition wouldn't. Just like the definition of the square root of a real number doesn't need the number to be nonnegative, but many lemmas on the square root have nonnegativity as an assumption.

This may look crazy at first sight, but experimentally this is a much better design than partial functions (as you are seeing just now).

[Michael Stoll]

The problem is that as soon as you want to use properties of the Legendre symbol to prove properties of the Jacobi symbol, you will have to know that the elements of b.factors are primes. So this only moves the problem to another place.
(And BTW, I do define the Jacobi symbol for all natural numbers b.)

[Yaël Dillies]

In fact, you will need to know many things, but you don't bake them into the definition!

[Michael Stoll]

Anyway, here is a version of list.pmap_append:

lemma list.pmap_append {α β : Type*} {p : α → Prop} (f : Π (a : α), p a → β) (l₁ l₂ : list α)
  (h₁ : ∀ (a : α), a ∈ l₁ → p a) (h₂ : ∀ (a : α), a ∈ l₂ → p a)
  (h : ∀ (a : α), a ∈ l₁ ++ l₂ → p a) :
  (l₁ ++ l₂).pmap f h = l₁.pmap f h₁ ++ l₂.pmap f h₂ :=
begin
  induction l₁ with hd tl ih,
  { simp only [list.pmap, list.nil_append], },
  { simp only [list.pmap, list.cons_append, eq_self_iff_true, true_and],
    apply ih, }
end

(Of course, h follows from h₁ and h₂, but it is not clear to me how to factor that in in a way that allows for easy use later.)

[Yaël Dillies]

You can give it a default value.

[Michael Stoll]

How do I spell this out?

[Johan Commelin]

Michael Stoll said:

you will have to know that the elements of b.factors are primes

That's known right? It's just that it is not available when doing a finset.prod because finset.prod takes a total function as argument.

[Michael Stoll]

Johan Commelin said:

Michael Stoll said:

you will have to know that the elements of b.factors are primes

That's known right? It's just that it is not available when doing a finset.prod because finset.prod takes a total function as argument.

Exactly.

[Johan Commelin]

So extending legendre_sym to return a junk value of 1 makes sense to me.

[Michael Stoll]

With list.pmap_appendavailable, I can do what I need using docs#list.pmap.

[Yaël Dillies]

But more likely you want to remove h_0 and h_1 and deduce them from h, so that you can rewrite without having to prove everything over and over again.

[Michael Stoll]

I.e., I can get rid of nat.prime_factors, and everything will work fine.

[Michael Stoll]

Yaël Dillies said:

But more likely you want to remove h_0 and h_1 and deduce them from h, so that you can rewrite without having to prove everything over and over again.

In my use case, you know h₁ and h₂ already, but you don't immediately have h, so I think the variant with building h out of the other two seems preferable.

[Michael Stoll]

(Because the lists are of the form n.factors, so you know they consist of primes.)

[Michael Stoll]

But maybe you are right. I'll try it out...

[Michael Stoll]

lemma list.pmap_append' {α β : Type*} {p : α → Prop} (f : Π (a : α), p a → β) (l₁ l₂ : list α)
  (h : ∀ (a : α), a ∈ l₁ ++ l₂ → p a) :
  (l₁ ++ l₂).pmap f h = l₁.pmap f (λ a ha, h a (list.mem_append_left l₂ ha))
                         ++ l₂.pmap f (λ a ha, h a (list.mem_append_right l₁ ha)) :=
begin
  induction l₁ with hd tl ih,
  { simp only [list.pmap, list.nil_append], },
  { simp only [list.pmap, list.cons_append, eq_self_iff_true, true_and, ih], }
end

[Moritz Doll]

Johan Commelin said:

So extending legendre_sym to return a junk value of 1 makes sense to me.

I would think that legendre_sym p should be the identity for p not prime, that would make even more stuff true in general.

[Michael Stoll]

This works in the same way in my application, so perhaps the second version is preferable.

[Michael Stoll]

The Legendre symbol is defined as the quadratic character of zmod p (composed with the natural map from the integers to zmod p), and the quadratic character is defined for (finite) fields.
My opinion is that once you have the Jacobi symbol, you don't have to consider Legendre symbols anyway, so in a sense the implementation of the latter is fairly irrelevant.

[Yaël Dillies]

Michael Stoll said:

In my use case, you know h₁ and h₂ already, but you don't immediately have h, so I think the variant with building h out of the other two seems preferable.

There's always the option of having both versions.

[Yaël Dillies]

I wish there was a way to tell Lean to infer a set of assumptions from another and back without having to write several lemmas.

[Michael Stoll]

The new version is now on github (#16290).

[Michael Stoll]

I think the following should be tagged @[simp]:
jacobi_sym_zero_right, jacobi_sym_one_right, jacobi_sym_one_left, jacobi_sym_zero_left (although the latter may rarely fire because of the condition on b). Should there be more simp lemmas?

[Alex J. Best]

I don't think making the last one a simp lemma is really worth it then, seeing as it will be so unlikely to fire, the others look good though

[Michael Stoll]

I think #16290 should be good now, except for moving the auxiliary lemmas at the beginning of the files to appropriate locations.
I have made some suggestions, but in some cases, it is not clear to me where to put them:

• nat.two_pow_mul_odd needs imports from data.nat.factorization.basic and data.nat.parity, but there does not seem to be an appropriate file that imports both
• int.sign_submonoid needs group_theory.submonoid.basic and the algebraic structure of the integers
• zmod.eq_zero_iff_gcd_ne_one and friends need data.zmod.basic and ring_theory.int.basic. data.zmod.quotient is a file that imports both, but has a different topic.

Suggestions? @Junyan Xu @Johan Commelin @Yaël Dillies

[Michael Stoll]

I have now moved the auxiliary lemmas to places I think are reasonable, For the three lemmas regarding the vanishing of a mod p, I have created a new file, since there was no suitable existing one with the right imports.
I hope this can now be merged. @Johan Commelin @Riccardo Brasca

[Johan Commelin]

I still had several comments. The PR is now well over 500 lines, and touches 8 files. I think the non-Jacobi files are pretty uncontroversial. If you move the Jacobi file to a new PR, then the rest can be merged quickly.

[Michael Stoll]

@Johan Commelin Thanks for the comments! I have fixed the documentation now. The Jacobi symbol file has been fairly thoroughly reviewed by @Junyan Xu, so my impression was that it should also be pretty uncontroversial right now. But if it helps, I can split it off, of course. (The PR was actually quite a bit longer initially, but thanks to Junyan's efforts at shortening proofs, it has shrunk over time.)

[Johan Commelin]

Well, it helps me review PRs in little gaps of free time if they are bit-sized and do only 1 thing.

[Michael Stoll]

OK; I'll remove the file jacobi_symbol.lean, wait until the rest gets merged (hopefully soon) and then PR it again.

[Michael Stoll]

@Johan Commelin #16290 has passed CI again.

[Riccardo Brasca]

Can you please update the title and the description of the PR? Thanks!

[Michael Stoll]

Done. (Plus the suggested simpa simplification.)

[Johan Commelin]

Michael Stoll said:

OK; I'll remove the file jacobi_symbol.lean, wait until the rest gets merged (hopefully soon) and then PR it again.

Thanks. Experience shows that review time is not linear in the total number of lines. Especially if many files change, then the changes should be trivial.

[Michael Stoll]

@Johan Commelin @Riccardo Brasca Can you reply to my question in the comments to #16290? Thanks!

[Riccardo Brasca]

I am busy this weekend, I will have a look on Monday, sorry.

[Michael Stoll]

As a teaser (in the hope that this will provide even more motivation to push along with #16290 and the follow-up Jacobi symbol PRs), I can now do the following in my set-up.

example : ¬ is_square (2345 : zmod 6789) :=
@nonsquare_of_jacobi_sym_eq_neg_one 2345 6789 (by norm_num)

The Jacobi symbol is evaluated using Quadratic Reciprocity, of course.

[Johan Commelin]

That's great!

[Johan Commelin]

It would be cool if we could have a norm_num extension that simplifies is_square x y to true or false.

[Johan Commelin]

Do you think you could do that?

[Johan Commelin]

It would involve a bit of meta-programming.

[Michael Stoll]

I have done a bit of meta-programming already for the norm_num extension (not that I would claim that I fully understand everything I'm doing).

[Michael Stoll]

But deciding whether a is a square mod b or not is at least as hard as factoring b, if I remember correctly.

[Michael Stoll]

You can quickly detect non-squares with the Jacobi symbol, but when b is not prime, there will be non-squares on which the symbol takes the value 1.

[Michael Stoll]

Of course, norm_num can evaluate b.factors, so one can use that, and it would work for whatever numbers the factorization can be done in this way.

[Johan Commelin]

Well, if it only works for numbers < 10^5 that's still nice.

[Johan Commelin]

We can work on fast factoring algorithms when we need them.

[Johan Commelin]

And we could even off-load it to an external program if we really need to.

[Michael Stoll]

I have to say that with larger numbers, I had issues with recursion depth, which I find surprising, since the proof terms produced by norm_num for p.prime can be much larger than what I get in the example above for the Jacobi symbol.

[Johan Commelin]

$ factor 987654321
987654321: 3 3 17 17 379721

have aux : 987654321 = 3 * 3 * 17 * 17 * 379721, rw aux,

for a very hands-on offloading of the factoring.

[Johan Commelin]

Hmmz, recursion depth shouldn't be a problem, unless your numbers get really large.

[Michael Stoll]

Johan Commelin said:

And we could even off-load it to an external program if we really need to.

We'd need certificates for the lagre primes that Lean can check efficiently.

[Michael Stoll]

Johan Commelin said:

Hmmz, recursion depth shouldn't be a problem, unless your numbers get really large.

I'll try again...

[Johan Commelin]

Sure, but what is the current limit for norm_num?

[Johan Commelin]

I guess it can handle primes of size ~ 10^6 easily? I have no idea.

[Michael Stoll]

example : [102334155 | 165580141]ⱼ= -1 := by norm_num

leads to deep recursion was detected at 'replace' (potential solution: increase stack space in your system)

[Michael Stoll]

The proof term I get (show_term {norm_num} does work) is roughly two-thirds the size (in lines or characters) of the one for

example : nat.prime 6791 := by norm_num

which does not run into any problems...

[Johan Commelin]

Hmm, I'm afraid we need someone with actual expertise in metaprogramming to debug this... not me...

[Johan Commelin]

Btw, note that we have the lucas primality test in mathlib. But I'm not sure if norm_num uses that.

[Michael Stoll]

It doesn't. norm_num does just trial factorization.

[Kevin Buzzard]

Michael Stoll said:

But deciding whether a is a square mod b or not is at least as hard as factoring b, if I remember correctly.

If you want to start thinking seriously about this sort of thing then I guess you should switch to Lean 4. Any interest in recruiting some people to work on writing the algorithms in Cohen's "A Course in Computational Algebraic Number Theory" in Lean 4, to see how competitive it is? The proofs they're correct can come later.

[Johan Commelin]

We could certainly get https://en.wikipedia.org/wiki/Primality_certificate#Pratt_certificates pretty easily.

[Michael Stoll]

Four-digit numbers seem to be about the limit for the Jacobi symbol.

[Michael Stoll]

example : nat.prime 1000003 := by norm_num works.

[Michael Stoll]

My proof term above is actually only one-third the size of the other one (I copied the infoview window, which showed the term twice).

[Johan Commelin]

Michael Stoll said:

Four-digit numbers seem to be about the limit for the Jacobi symbol.

Hmm, any idea what the computational bottle neck is?

[Michael Stoll]

No. It does not take long, and it uses less space than nat.prime 6791 (in terms of "allocated objects" and "allocated closures" as printed by the profiler). It is just the "deep recursion" that seems to be triggered at some point.

[Johan Commelin]

Hmm, weird. Because I would imagine there's at most 5 or 6 nested calls?

[Michael Stoll]

Here is the term for 2345 and 6789:

(id_tag ()
   (((λ (a a_1 : ℤ) (e_1 : a = a_1) (ᾰ ᾰ_1 : ℤ) (e_2 : ᾰ = ᾰ_1), congr (congr_arg eq e_1) e_2)
       [2345 | 6789]ⱼ
       (-1)
       (jacobi_sym_mod_left 2345 6789
          (ne_of_gt
             (bit1_pos'
                (bit0_pos
                   (bit1_pos'
                      (bit0_pos
                         (bit0_pos
                            (bit0_pos
                               (bit0_pos (bit1_pos' (bit0_pos (bit1_pos' (bit0_pos (bit1_pos' zero_lt_one')))))))))))))
          (-1)
          (jacobi_sym_qr₁ 2345 6789
             (norm_num.nat_mod 2345 4 586 1 2344
                (norm_num.mul_bit0_bit0 293 2 586 (norm_num.mul_bit0' 293 1 293 (mul_one 293)))
                (norm_num.one_add 2344 2345 (norm_num.bit0_succ 1172))
                (norm_num.lt_one_bit0 2 (norm_num.le_one_bit0 1 (le_refl 1))))
             (norm_num.nat_mod 6789 2 3394 1 6788 (norm_num.mul_bit0_bit0 1697 1 1697 (mul_one 1697))
                (norm_num.one_add 6788 6789 (norm_num.bit0_succ 3394))
                (norm_num.lt_one_bit0 1 (le_refl 1)))
             (-1)
             (jacobi_sym_nat_mod_left 6789 2345 (-1)
                (jacobi_sym_qr₁' 2099 2345
                   (norm_num.nat_mod 2099 2 1049 1 2098 (norm_num.mul_bit0' 1049 1 1049 (mul_one 1049))
                      (norm_num.one_add 2098 2099 (norm_num.bit0_succ 1049))
                      (norm_num.lt_one_bit0 1 (le_refl 1)))
                   (norm_num.nat_mod 2345 4 586 1 2344
                      (norm_num.mul_bit0_bit0 293 2 586 (norm_num.mul_bit0' 293 1 293 (mul_one 293)))
                      (norm_num.one_add 2344 2345 (norm_num.bit0_succ 1172))
                      (norm_num.lt_one_bit0 2 (norm_num.le_one_bit0 1 (le_refl 1))))
                   (-1)
                   (jacobi_sym_nat_mod_left 2345 2099 (-1)
                      (jacobi_sym_even_odd₃ 246 2099
                         (norm_num.nat_mod 246 2 123 0 246 (norm_num.mul_bit0' 123 1 123 (mul_one 123)) (zero_add 246)
                            (bit0_pos zero_lt_one'))
                         (norm_num.nat_mod 2099 8 262 3 2096
                            (norm_num.mul_bit0_bit0 131 4 524
                               (norm_num.mul_bit0' 131 2 262 (norm_num.mul_bit0' 131 1 131 (mul_one 131))))
                            (norm_num.add_bit1_bit0 1 1048 1049 (norm_num.one_add 1048 1049 (norm_num.bit0_succ 524)))
                            (norm_num.lt_bit1_bit0 1 4 (norm_num.sle_one_bit0 2 (norm_num.le_one_bit0 1 (le_refl 1)))))
                         1
                         (jacobi_sym_qr₃ 123 2099
                            (norm_num.nat_mod 123 4 30 3 120
                               (norm_num.mul_bit0_bit0 15 2 30 (norm_num.mul_bit0' 15 1 15 (mul_one 15)))
                               (norm_num.add_bit1_bit0 1 60 61 (norm_num.one_add 60 61 (norm_num.bit0_succ 30)))
                               (norm_num.lt_bit1_bit0 1 2 (norm_num.sle_one_bit0 1 (le_refl 1))))
                            (norm_num.nat_mod 2099 4 524 3 2096
                               (norm_num.mul_bit0_bit0 262 2 524 (norm_num.mul_bit0_bit0 131 1 131 (mul_one 131)))
                               (norm_num.add_bit1_bit0 1 1048 1049
                                  (norm_num.one_add 1048 1049 (norm_num.bit0_succ 524)))
                               (norm_num.lt_bit1_bit0 1 2 (norm_num.sle_one_bit0 1 (le_refl 1))))
                            (-1)
                            (jacobi_sym_nat_mod_left 2099 123 (-1)
                               (jacobi_sym_even_odd₃ 8 123
                                  (norm_num.nat_mod 8 2 4 0 8 (norm_num.mul_bit0_bit0 2 1 2 (mul_one 2)) (zero_add 8)
                                     (bit0_pos zero_lt_one'))
                                  (norm_num.nat_mod 123 8 15 3 120
                                     (norm_num.mul_bit0' 15 4 60
                                        (norm_num.mul_bit0' 15 2 30 (norm_num.mul_bit0' 15 1 15 (mul_one 15))))
                                     (norm_num.add_bit1_bit0 1 60 61 (norm_num.one_add 60 61 (norm_num.bit0_succ 30)))
                                     (norm_num.lt_bit1_bit0 1 4
                                        (norm_num.sle_one_bit0 2 (norm_num.le_one_bit0 1 (le_refl 1)))))
                                  1
                                  (jacobi_sym_even_odd₃ 4 123
                                     (norm_num.nat_mod 4 2 2 0 4 (norm_num.mul_bit0_bit0 1 1 1 (one_mul 1)) (zero_add 4)
                                        (bit0_pos zero_lt_one'))
                                     (norm_num.nat_mod 123 8 15 3 120
                                        (norm_num.mul_bit0' 15 4 60
                                           (norm_num.mul_bit0' 15 2 30 (norm_num.mul_bit0' 15 1 15 (mul_one 15))))
                                        (norm_num.add_bit1_bit0 1 60 61
                                           (norm_num.one_add 60 61 (norm_num.bit0_succ 30)))
                                        (norm_num.lt_bit1_bit0 1 4
                                           (norm_num.sle_one_bit0 2 (norm_num.le_one_bit0 1 (le_refl 1)))))
                                     (-1)
                                     (jacobi_sym_even_odd₃ 2 123
                                        (norm_num.nat_mod 2 2 1 0 2 (one_mul 2) (zero_add 2) (bit0_pos zero_lt_one'))
                                        (norm_num.nat_mod 123 8 15 3 120
                                           (norm_num.mul_bit0' 15 4 60
                                              (norm_num.mul_bit0' 15 2 30 (norm_num.mul_bit0' 15 1 15 (mul_one 15))))
                                           (norm_num.add_bit1_bit0 1 60 61
                                              (norm_num.one_add 60 61 (norm_num.bit0_succ 30)))
                                           (norm_num.lt_bit1_bit0 1 4
                                              (norm_num.sle_one_bit0 2 (norm_num.le_one_bit0 1 (le_refl 1)))))
                                        1
                                        (jacobi_sym_qr₁ 1 123
                                           (norm_num.nat_mod 1 4 0 1 0 (zero_mul 4) (add_zero 1)
                                              (norm_num.lt_one_bit0 2 (norm_num.le_one_bit0 1 (le_refl 1))))
                                           (norm_num.nat_mod 123 2 61 1 122 (norm_num.mul_bit0' 61 1 61 (mul_one 61))
                                              (norm_num.one_add 122 123 (norm_num.bit0_succ 61))
                                              (norm_num.lt_one_bit0 1 (le_refl 1)))
                                           1
                                           (jacobi_sym_nat_one_right 123)))))))))))))
       (-1)
       (-1)
       (eq.refl (-1))).trans
      (eq_true_intro (eq.refl (-1))))).mpr
  trivial

It's not that deeply nested (much less so than the primality proof term I was referring to).

[Johan Commelin]

Yeah, looks like ~10 nestings of jacobi_XYZ

[Michael Stoll]

The term for the large (8-digit) numbers has nesting depth 30-40.

[Michael Stoll]

I have now a new version (of the norm_num extension) that has no problems with

example : jacobi_sym_nat 102334155 165580141 = -1 := by norm_num

say, where jacobi_sym_nat is a version of the Jacobi symbol restricted to natural numbers.
However, the conversion from the Jacobi symbol to the restricted version still has a problem. This is done via

jacobi_sym a b = jacobi_sym_nat (a % b) b

and the issue seems to be some problem with matching types (along the lines of

example : 100000 % 3 = 1 := rfl

which runs into the deep recursion problem). I suspect some nat vs int problem.

Here is a question: How do I produce a proof of (a : int) = b within the framework of the norm_num machinery when I have a numeral expression ea for a : nat and a numeral expression eb for b : int (and of course a and b are equal)? @Mario Carneiro

[Mario Carneiro]

use prove_nat_uncast

[Michael Stoll]

The relevant part of my code looks like the following.

lemma jacobi_sym_mod_left (a : ℤ) (b : ℕ) (ab r : ℤ) (hab : a % b = ab)
  (hr : jacobi_sym_nat ab.to_nat b = r) :
  jacobi_sym a b = r := ...

/-- This evaluates `r := jacobi_sym a b` and produces a proof term for the equality.
This is done by reducing to `r := jacobi_sym_nat (a % b) b`. -/
meta def prove_jacobi_sym :
  instance_cache → instance_cache → expr → expr
    → tactic (instance_cache × instance_cache × expr × expr)
| zc nc ea eb := do
  a ← ea.to_int,
  b ← eb.to_nat,
  match b with -- deal with simple cases right away
  | 0 := pure (zc, nc, `(1 : ℤ), `(jacobi_sym_zero_right).mk_app [ea])
  | 1 := pure (zc, nc, `(1 : ℤ), `(jacobi_sym_one_right).mk_app [ea])
  | _ := do -- Now `1 < b`. Compute `jacobi_sym_nat (a % b) b` instead.
    (zc, eb') ← zc.of_int (b : ℤ),
    (zc, amb, phab) ← norm_num.prove_div_mod zc ea eb' tt, -- compute `a % b`
    let a' := amb.to_int,
    match a' with -- not sure why this is necessary, but `int.to_nat a'` gives an error.
    | none := failed -- we'll never get here
    | (some aa) := do
      (nc, amb') ← nc.of_nat (int.to_nat aa), -- `a % b` as a natural number
      (zc, nc, er, p) ← prove_jacobi_sym_nat zc nc amb' eb, -- compute `jacobi_sym_nat (a % b) b`
      pure (zc, nc, er, `(jacobi_sym_mod_left).mk_app [ea, eb, amb, er, phab, p])
    end
  end


/-- This is the `norm_num` plug-in that evaluates Jacobi symbols. -/
@[norm_num] meta def eval_jacobi_sym : expr → tactic (expr × expr)
| `(jacobi_sym %%ea %%eb) := do
    zc ← mk_instance_cache `(ℤ),
    nc ← mk_instance_cache `(ℕ),
    (prod.snd ∘ prod.snd) <$> prove_jacobi_sym zc nc ea eb
| _ := failed

[Michael Stoll]

Mario Carneiro said:

use prove_nat_uncast

I tried that, but it didn't help. Maybe I should try again...

[Michael Stoll]

(The code has changed a bit since then...)

[Mario Carneiro]

your lemma should look something like:

lemma jacobi_sym_mod_left (a : ℤ) (b c : ℕ) (ab r : ℤ) (hab : a % b = ab) (hc : (c : ℤ) = ab)
  (hr : jacobi_sym_nat c b = r) :
  jacobi_sym a b = r := ...

and then generate c and hc using prove_nat_uncast applied to ab

[Michael Stoll]

I had tried something like this, but that was before I re-wrote the jacobi_sym_nat part. I'll see if it works now...

[Michael Stoll]

With this:

meta def prove_jacobi_sym :
  instance_cache → instance_cache → expr → expr
    → tactic (instance_cache × instance_cache × expr × expr)
| zc nc ea eb := do
  a ← ea.to_int,
  b ← eb.to_nat,
  match b with -- deal with simple cases right away
  | 0 := pure (zc, nc, `(1 : ℤ), `(jacobi_sym_zero_right).mk_app [ea])
  | 1 := pure (zc, nc, `(1 : ℤ), `(jacobi_sym_one_right).mk_app [ea])
  | _ := do -- Now `1 < b`. Compute `jacobi_sym_nat (a % b) b` instead.
    (zc, eb') ← zc.of_int (b : ℤ),
    (zc, amb, phab) ← norm_num.prove_div_mod zc ea eb' tt, -- compute `a % b`
    (zc, nc, amb', phab') ← norm_num.prove_nat_uncast zc nc amb,
    (zc, nc, er, p) ← prove_jacobi_sym_nat zc nc amb' eb, -- compute `jacobi_sym_nat (a % b) b`
    pure (zc, nc, er, `(jacobi_sym_mod_left).mk_app [ea, eb, amb', amb, er, phab, phab', p])
  end

it still gives me "deep recursion" on larger numbers. Should I do something different to obtain eb'?

[Michael Stoll]

I really have no clue what the problem might be. The types of the sub-terms of the proof term that norm_num generates seem to match the signature of (the new version of) jacobi_sym_mod_left.

[Michael Stoll]

I have some idea how to work around the problem, but it is getting too late today...

[Mario Carneiro]

Your a % b = ab subgoal has coercions in it

[Mario Carneiro]

it should either be all-int or all-nat

[Mario Carneiro]

it's hard to write exact examples without the MWE

[Mario Carneiro]

what is the type of jacobi_sym and jacobi_sym_nat?

[Mario Carneiro]

You have to be really careful to not have type errors between ints and coerced nats, because they cause kernel timeouts

[Michael Stoll]

OK; here is an #mwe that works:

import tactic.norm_num
import tactic.zify

namespace tactic
namespace norm_num

def test (a : ℤ) (b : ℕ) : ℕ := if b = 0 then 0 else (a % b).to_nat

lemma test_def_zero (a : ℤ) : test a 0 = 0 := rfl

lemma test_def_pos (a : ℤ) (b : ℕ) (hb : b ≠ 0) : (test a b : ℤ) = a % b :=
begin
  rw [test],
  split_ifs with h, exact false.elim (hb h),
  refine (int.to_nat_of_nonneg (int.mod_nonneg a _)),
  exact_mod_cast hb,
end

lemma test_nn (a : ℤ) (b r : ℕ) (ab b' : ℤ) (hb₀ : b ≠ 0)
  (hb : (b : ℤ) = b' ) (hr : (r : ℤ) = ab) (hab : a % b' = ab) :
  test a b = r :=
begin
  zify,
  rwa [test_def_pos a b hb₀, hb, hr],
end

meta def prove_test : instance_cache → instance_cache → expr → expr →
  tactic (instance_cache × instance_cache × expr × expr)
| zc nc ea eb := do
  b ← eb.to_nat,
  match b with
  | 0 := pure (zc, nc, `(0 : ℕ), `(test_def_zero).mk_app [ea])
  | _ := do
    (nc, phb) ← norm_num.prove_ne nc eb `(0 : ℕ) b 0, -- proof of `b ≠ 0`
    (zc, eb') ← zc.of_int (b : ℤ),
    (zc, nc, eb₁, pb') ← norm_num.prove_nat_uncast zc nc eb',
    (zc, amb, phab) ← norm_num.prove_div_mod zc ea eb' tt, -- compute `a % b`
    (zc, nc, amb', phab') ← norm_num.prove_nat_uncast zc nc amb,
    pure (zc, nc, amb', `(test_nn).mk_app [ea, eb₁, amb', amb, eb', phb, pb', phab', phab])
  end

@[norm_num] meta def eval_test : expr → tactic (expr × expr)
| `(test %%ea %%eb) := do
    nc ← mk_instance_cache `(ℕ),
    zc ← mk_instance_cache `(ℤ),
    (prod.snd ∘ prod.snd) <$> prove_test zc nc ea eb
| _ := failed

example : test 3 0 = 0 := by norm_num
example : test 100000000 33333333 = 1 := by {norm_num,}

end norm_num
end tactic

Using the same modification (i.e., introducing b' and the proof that (b : int) = b') then also works in the Jacobi symbol setting. Thanks! @Mario Carneiro

[Mario Carneiro]

you shouldn't need to use of_int there

[Mario Carneiro]

prove_nat_uncast returns the new expression

[Mario Carneiro]

Oh nvm, I see you are using prove_nat_uncast "backwards" (i.e. you are actually computing a cast and not an uncast). That hasn't come up before but I guess you need it to evaluate int.mod on nats

[Michael Stoll]

Yes, int.mod on an int and a nat seems to be causing the problems here.

[Michael Stoll]

I have now PRed the Jacobi symbol code again: #16395
@Johan Commelin @Riccardo Brasca

[Michael Stoll]

Next up: a norm_num extension that computes Jacobi (and Legendre) symbols.
This can now also deal with large numbers; the problems have been resolved thanks to @Mario Carneiro's help; see the discussion above.

[Michael Stoll]

On my machine, according to the profiler,

example : [58378362899022564339483801989973056405585914719065 |
           53974350278769849773003214636618718468638750007307]ⱼ = -1 := by norm_num

takes less than 0.4 seconds.

[Michael Stoll]

Here is a question regarding notation (see the comments in the file added in #16395).

After defining notation [a | b]ⱼ for the Jacobi symbol, it seems that, at least in some contexts, list notation [x, y, z] does not work anymore. (It does in tactic mode, and square brackets work for instance arguments, though.) This surprises me, since Lean has no problems with differentiating notations like {1, 2, 3}, ¸{x : nat | odd x} and {x : nat // odd x}`, which looks like a similar situation to me.

• Can somebody explain what the problem may be?
• Is there a fix that allows using both Jacobi symbol and list notation at the same time?
• If not, my suggestion would be to switch to ⱼ( a | b ) instead (round parens are closer to the usual notation), even though I don't like the ⱼ on the left very much. But at least it does not seem to lead to conflilcts...
• Alternatives?

[Mario Carneiro]

Michael Stoll said:

This surprises me, since Lean has no problems with differentiating notations like {1, 2, 3}, ¸{x : nat | odd x} and {x : nat // odd x}, which looks like a similar situation to me.

Actually, lean can't handle that either. All notations starting with { are built in and not implemented via the usual notation framework.

[Michael Stoll]

Does this mean that the opening part of a bracket-like notation has to be unique?

[Mario Carneiro]

yes, I would recommend that

[Michael Stoll]

OK; then I'll change to ⱼ( a | b ) as per my proposal above. Unless somebody has a better suggestion.

[Mario Carneiro]

How about just J(a | b)?

[Michael Stoll]

OK, let's have a poll. What should the notation for the Jacobi symbol be?

[Mario Carneiro]

isn't the jacobi symbol usually written like the legendre symbol? (i.e. $\left(\dfrac{a}{b}\right)$)

[Michael Stoll]

Yes. But this is hard to replicate in unicode.

[Mario Carneiro]

maybe something more like a division e.g. a /_j b

[Michael Stoll]

/poll What should be the notation for the Jacobi symbol?
ⱼ( a | b )
J(a | b)
a /_j b

[Michael Stoll]

NOTE: Everything starting with just ( or [ will not work. See the discussion above and below.

I'd prefer to have some kind of brackets.

[Michael Stoll]

The symbol has not much to do with division, other than that in a way, we consider a modulo b.

[Michael Stoll]

( a | b )ⱼ will run into even worse problems than [a | b]ⱼ, since round parentheses are ubiquitous...

[Mario Carneiro]

it's not even legal

[Adam Topaz]

Oh, that's a shame.

[Michael Stoll]

(That was what I tried at the very beginning :smile: )

[Adam Topaz]

I suppose Matt's suggestion (a |_j b) is also illegal

[Mario Carneiro]

Michael Stoll said:

The symbol has not much to do with division, other than that in a way, we consider a modulo b.

Sure, but you are asking about the notation. Maybe you should tell this to Legendre

[Adam Topaz]

Do we have notation for Legendre?

[Mario Carneiro]

not that I can see

[Michael Stoll]

No (other than legendre_sym).

[Adam Topaz]

Maybe we can find some fancy unicode parens to use for both Legendre and Jacobi (with an additional subscript j).

[Mario Carneiro]

FWIW I used /_L for the legendre / jacobi / kronecker symbol in metamath

[Michael Stoll]

I think that ?(a | b) looks like a horizontal version of $\left(\frac{a}{b}\right)$, that's why I'm suggesting the vertical bar.

[Adam Topaz]

What's the unicode analogue of https://detexify.kirelabs.org/classify.html ?

[Matthew Ballard]

https://shapecatcher.com ?

[Mario Carneiro]

note that the vertical bar is also taken for a bunch of things in lean

[Mario Carneiro]

which is why dvd uses \| instead of |

[Michael Stoll]

Anyway, J(a | b) doesn't seem to give trouble, so I'd be willing to use that.

[Adam Topaz]

⟮ a | b ⟯
Those are unicode characters 0x27EE and 0x27EF

[Yuyang Zhao]

Does _j in choices means ⱼ?

[Yuyang Zhao]

Or they are different things?

[Michael Stoll]

Adam Topaz said:

⟮ a | b ⟯
Those are unicode characters 0x27EE and 0x27EF

Then we'd need an easy way of typing them in.

[Adam Topaz]

Michael Stoll said:

Adam Topaz said:

⟮ a | b ⟯
Those are unicode characters 0x27EE and 0x27EF

Then we'd need an easy way of typing them in.

Right... but it's fairly easy to add symbols to the vscode extension

[Adam Topaz]

For example we did a similar trick for the Lie bracket
https://github.com/leanprover-community/mathlib/blob/3c227fc0bf65ba9eb3048eb7845f8ce2a6084d45/src/data/bracket.lean#L37

[Adam Topaz]

We might be the only people in the world who use ⁅

[Mario Carneiro]

⟮ is already fairly overloaded in mathlib, I see it used for field extensions and manifold sets

[Michael Stoll]

VSCode says type ⟮ using \([ and `type ⟯ using \]).

[Adam Topaz]

There are other nice looking parens as well

[Michael Stoll]

Mario Carneiro said:

⟮ is already fairly overloaded in mathlib, I see it used for field extensions and manifold sets

That may not be a big problem, as the notation is localized and does not have to be used.

[Adam Topaz]

http://shapecatcher.com/unicode/info/10647

[Mario Carneiro]

not that this is a big problem since the notation is localized, but I guess it means you can't use jacobi symbols and field extensions at the same time

[Michael Stoll]

You can certainly use the symbols, but you'd have to write jacobi_sym a b instead.

[Mario Carneiro]

personally I prefer to just use no notation here. It's not like this is an algebraic operation that gets deeply nested - most equations only use the symbol once or twice

[Mario Carneiro]

and using a notation means that people have to look it up if it's not super obvious

[Mario Carneiro]

and [a | b], is certainly not obvious

[Michael Stoll]

So you vote for "nothing" in the poll above?

[Mario Carneiro]

I could, but I'm not sufficiently invested in this part of the library to care so I'll abstain

[Adam Topaz]

I suppose the notation situation would be (theoretically) better with Lean4? If so, we could leave it as "nothing" for now, with a "plan" to do something about it after the transition to Lean4?

[Mario Carneiro]

yes, lean 4's parser allows for backtracking so [a | b] is theoretically possible (although it is probably ambiguous with other stuff)

[Michael Stoll]

Maybe I'll wait to see if there are more votes, and if no clear majority (for a legal notation) emerges, I'll remove the notation for the time being.

[Mario Carneiro]

now we just need to get it to support

       [ a ]
#check [---]
       [ b ]

[Adam Topaz]

Isn't there some ITP that lets you do stuff like this?

[Mario Carneiro]

I have a vague recollection of such but I couldn't say which

[Yaël Dillies]

If you make a /_j b the notation, then you can also write it (a /_j b), right?

[Yaël Dillies]

And it makes the parentheses mean the same thing as elsewhere, which sounds desirable.

[Kevin Buzzard]

Ha ha you could also give it super high binding power which would force people to use brackets

[Michael Stoll]

I would think that super high binding power would usually make brackets unnecessary.

[Michael Stoll]

We should also think of the (quadratic) Hilbert symbol, which should be something like (a b | v) or perhaps (a, b | v) with some decorations to make the parser happy.

[Michael Stoll]

J(a | b) and H(a b | v) would be reasonable.

[Michael Stoll]

I'm less convinced by (a /ⱼ b) and (a b /ₕ v) or (a |ⱼ b) and (a b |ₕ v), but am willing to get more convinced...

[Johan Commelin]

I think In Lean 3 H(a b | v) is not possible. You will need a ,.

[Michael Stoll]

A comma is fine or perhaps even better, since the math notation is $\left(\frac{a,b}{v}\right)$.

[Michael Stoll]

I have now changed notation to J(a | b). It is easy to read (and type!) and unlikely to lead to confusion.

Perhaps somebody can look through #16395 carefully and check that it is OK for merging. @Johan Commelin @Riccardo Brasca

[Johan Commelin]

I am unfortunately (for all practical purposes) off-line till Monday.

[Riccardo Brasca]

I've left a couple of small comments, but LGTM!

[Michael Stoll]

I have now prepared a PR that moves zmod.{legendre|jacobi}_sym* to the root namespace: #16461.
To disambiguate with the version of QR for the Jacobi symbol, I have moved zmod.quadratic_reciprocity* to legendre_sym_quadratic_reciprocity*.
This also moves the results related to Gauss' and Eisenstein's lemmas from legendre_symbol. to zmod..
@Riccardo Brasca @Johan Commelin

[Michael Stoll]

@Riccardo Brasca suggested, and I had myself also thought about it, to change {legendre|jacobi}_sym_stuff to {legendre|jacobi}_sym.stuff, i.e., to move the various lemmas into a namespace legendre_sym or jacobi_sym. One possible problem is that replacing jacobi_sym_two with jacobi_sym.two (there are a few other similar names) looks a bit strange, in particular if locally referred to simply as two. I would propose changing this to jacobi_sym.value_at_two; then the move should be fine and I'm inclined to do it. But I'd like to put this up for discussion first, so please comment! (Here or at #16461.)

[Alex J. Best]

You can use protected lemma to ensure that a specific lemma with an overly general name (when not namespaced) can't be referred to without the full prefix

[Damiano Testa]

Alternatives could be jacobi_sym_two, avoiding namespacing or jacobi_sym.at_two, shorter than value_at_two and in the namespace.

[Michael Stoll]

I think I'll mostly go with Damiano's suggestion.

[Michael Stoll]

#16461 has now been merged (thanks @Riccardo Brasca), so the next step is to PR the norm_num extension that (provably) computes the value of the Jacobi (and Legendre) symbol. This is about 500 lines of code, so it should go in a separate file. Would number_theory.legendre_symbol.norm_num be a reasonable name for that file?

[Michael Stoll]

Since nobody complained, I've named the file as suggested. See #16519. @Mario Carneiro @Johan Commelin @Riccardo Brasca

[Michael Stoll]

@Mario Carneiro has approved the PR. Should somebody else look at it before it can be merged?

[Michael Stoll]

@Mario Carneiro thanks!

[Michael Stoll]

I noticed that after changing the notation to J(a | b), the parenthetical remark "(Unfortunately, there is no subscript "J" in unicode.)" in the docstring for the notation no longer makes sense. I've PRed a fix (#16541). Can somebody merge it after it has passed CI? @Riccardo Brasca @Johan Commelin

[Michael Stoll]

@Riccardo Brasca thanks!

[Riccardo Brasca]

@Michael Stoll Can you please check if this is still needed with latest master? It's possible that #16349 fixes it.

[Michael Stoll]

@Riccardo Brasca Did you mean #16463? #16349 was closed, but refers to it. Anyway, it looks like with the current state of affairs, this instance is no longer needed. So I'm removing it in #16659.

[Riccardo Brasca]

Yes sorry. Thanks

[Riccardo Brasca]

I am away from the computer, but feel free to request my review

[Michael Stoll]

Done.