large numbers computations

Stream: maths

[Jakob von Raumer]

Are there tactics to deal with goals like 10 ^ 18 < 3618502788666131213697322783095070105623107215331596699973092056135872020483?

[Riccardo Brasca]

norm_num is quite fast for me on this precise example

[Jakob von Raumer]

Thanks!

[Jakob von Raumer]

Do we have a fast prime checker as well?

[Riccardo Brasca]

This is really non trivial, right? I mean, even Pari/GP, or similar software, is very fast because it checks psudoprimality

[Jireh Loreaux]

We'll, no one ever really implements AKS, but PRIMES is in P.

[Riccardo Brasca]

Sure, but AFAIK AKS is totally useless in practice. I think that in practice if one as primality certificate (given by Pari or whatever) it shouldn't be too hard to formalize primality in Lean. But of course asking Lean to find the certificate is different.

[Jakob von Raumer]

Sure would be good to have a decision procedure that's at least faster than just brute force

[Riccardo Brasca]

norm_num also knows about primes, so you can use it, but it is not very fast.

[Riccardo Brasca]

On the other hand we have a good implementation of the Jacobi symbol, so checking for example if 2, 3 and 5 (say) are witness of nonprimality should be fast. And this is already a good pseudoprimality test.

[Riccardo Brasca]

Especially if one implements Miller-Rabin rather than Solovay–Strassen. According to Wikipedia, 2, 3 and 5 are enough up to 25326001.

[Eric Rodriguez]

Riccardo Brasca said:

Especially if one implements Miller-Rabin rather than Solovay–Strassen. According to Wikipedia, 2, 3 and 5 are enough up to 25326001.

sure, but then we have to prove that it's enough :)

[Eric Rodriguez]

we do have lucas lehmer if you want a mersenne prime though!

[Riccardo Brasca]

Yes, I was talking about psuedoprimality. If we want primality I am more or less sure the best we can do is prove that something like a Pratt certificate implies primality (this is surely almost already in mathlib) and then use something else to find the witness.

[Riccardo Brasca]

Like what polyrith does

[Eric Rodriguez]

ah, that's a good idea. indeed we do have docs#lucas_primality; in Lean4 I imagine this can be done quite fast

[Eric Rodriguez]

maybe with some cached tactics too

[Jakob von Raumer]

ECPP prime certificates would be an alternative. Might be a fun project to see what mathematics about elliptic curves is missing to implement a checker

[Jakob von Raumer]

Hasse's theorem seems to be the main result it relies on

[Kevin Buzzard]

We can't even prove the points form a group. There is a low-level proof of Hasse's theorem due to Manin I think, but only valid in characteristic >= 5 IIRC

[Jakob von Raumer]

Well characteristic >= 5 would be enough for the primality check, but it's probably better to go a high-level way and build a good library of elliptic curves first

[Bhavik Mehta]

I tried writing a verifier using norm_num and Pratt certificates, the summary was essentially that it works but doesn't seem to be much better than what norm_num already does - in particular it gets a bit further before timing out but not significantly so

[Riccardo Brasca]

Even using Pratt certificates recursively? I mean, to prove that all the factors are prime.

[Bolton Bailey]

@Bhavik Mehta does it at least work for primes like 15 * 2^27 + 1? It would be convenient for a project I'm working on to verify this prime in particular.

[Bhavik Mehta]

Riccardo Brasca said:

Even using Pratt certificates recursively? I mean, to prove that all the factors are prime.

Yes!

[Bhavik Mehta]

Bolton Bailey said:

Bhavik Mehta does it at least work for primes like 15 * 2^27 + 1? It would be convenient for a project I'm working on to verify this prime in particular.

I will check

[Bolton Bailey]

Another comment on primality testing, obviously there are several different approaches - brute force calculation, Pratt certificates, AKS, and all of these are going to have different performance on different number sizes (this is yet another thing that was first pointed out to me by @Mario Carneiro ). It would be cool if we could eventually have a primality tester that selectively applies different tests based on the size of the number being proven, in the same way that big-integer multiplication libraries use different algorithms on different sizes of integers.

[Jakob von Raumer]

@Bhavik Mehta I don't really get why it should be that slow? Because the modulo calculations are?

[Bhavik Mehta]

One factor to bear in mind is that in computation, Pratt certificates are an improvement because you replace many divisions in the trial division algorithm with multiplication and powmod operations. But when verifying, the cost of a natural division is practically the same as a multiplication, since a proof of division amounts to an inequality proof, addition proof and multiplication proof. So the usual benefit from Pratt certificates is lessened greatly

[Riccardo Brasca]

Bolton Bailey said:

Bhavik Mehta does it at least work for primes like 15 * 2^27 + 1? It would be convenient for a project I'm working on to verify this prime in particular.

This number is not that big. It would be sad if we're not able to prove it's prime...

[Riccardo Brasca]

If you really want this specific number I don't see any problem. I just tested and norm_num is able to compute 31^15, then one only have to take squares (I hope not by hand...). The prime factors of p-1 are super small so exploiting the condition nat.prime q → q ∣ p - 1 is probably easy (one may need a general result here, unrelated to this number).

[Riccardo Brasca]

I think it's more helpful to prove docs#lucas_primality using docs#nat.factorization

[Alex J. Best]

Is https://en.wikipedia.org/wiki/Proth%27s_theorem the general result you are thinking of?

[Riccardo Brasca]

Ah, this is even better, but I was just thinking that if 3*5*2^27 = n then we know all prime factors of n.

[Bolton Bailey]

Unfortunately, norm_num times out on the goal 2 ^ (2 ^ 27 * 15) = 1 in zmod (2 ^ 27 * 15 + 1). I think the issue is that exponentiation was not performant. Mario gave us a definition to get around this when we were working on Miller-Rabin (unfortunately that branch has been left hanging for a while, I don't have time to finish it up)

def binpow {M} [has_one M] [has_mul M] (m : M) : ℕ → M :=
nat.binary_rec 1 (λ b _ ih, let ih2 := ih * ih in cond b (m * ih2) ih2)

[Riccardo Brasca]

I think you can do it using pow_mul, helping Lean a little bit.

[Riccardo Brasca]

This is of course not the right way of doing it, but a least it should work in this precise case.

[Bhavik Mehta]

Bolton Bailey said:

Unfortunately, norm_num times out on the goal 2 ^ (2 ^ 27 * 15) = 1 in zmod (2 ^ 27 * 15 + 1). I think the issue is that exponentiation was not performant. Mario gave us a definition to get around this when we were working on Miller-Rabin (unfortunately that branch has been left hanging for a while, I don't have time to finish it up)

def binpow {M} [has_one M] [has_mul M] (m : M) : ℕ → M :=
nat.binary_rec 1 (λ b _ ih, let ih2 := ih * ih in cond b (m * ih2) ih2)

Yes, the first thing I did when trying to verify using Pratt certificates was convince norm_num to do efficient modular exponentiation

[Bhavik Mehta]

@Bolton Bailey update: I can prove that number is prime, and it's relatively fast - around 5 seconds

[Bolton Bailey]

Hehe, sweet!

[Bolton Bailey]

Is the proof fully automatic, or do you need to hand-solve some goals?

[Bhavik Mehta]

bad news is that the original setup I had to do primes in general was pretty dodgy - I adapted it from https://www.isa-afp.org/browser_info/current/AFP/Pratt_Certificate/outline.pdf, namely rules 1 and 2 from page 2, and it seems to be really inefficient in the case of a prime factor appearing multiple times (unsure if this is from the isabelle version or my translation of it)

[Bhavik Mehta]

The proof is:

lemma second_bit {x} : x ∈ nat.factors (15 * 2^27 + 1 - 1) ↔ x = 2 ∨ x = 3 ∨ x = 5 :=
begin
  norm_num,
end

lemma boltons_prime : prime (15 * 2^27 + 1) :=
begin
  refine pratt_rule_2' 31 (by norm_num) (by norm_num) _,
  simp only [pratt_predicate, second_bit, forall_eq_or_imp, forall_eq],
  change my_pow_mod _ _ _ ≠ 1 ∧ my_pow_mod _ _ _ ≠ 1 ∧ my_pow_mod _ _ _ ≠ 1,
  norm_num,
end

with about 100 lines of other stuff before it which I'd written a few months ago

[Bhavik Mehta]

what I have right now definitely isn't correctly designed to be usable, which is why I haven't tried to PR it anywhere!

[Bolton Bailey]

Still, very cool to have primality proofs for cryptographically important primes.

[Bhavik Mehta]

Thought about it a bit more, and the problem with what I had originally is precisely in the case where a prime factor q of p-1 has high multiplicity - the isabelle pratt certificate (and hence mine) do k powmod calculations if q^k | p-1, but it should be enough to only do 1 for each prime

[Bhavik Mehta]

(So my original attempt for your prime timed out, but the modified version above is shorter and doesn't time out!)

[Bhavik Mehta]

Here are the (average of 1) benchmarks for the four primes I've tested: 199999 takes 2s, 2013265921 (Bolton's prime) takes 3s, 1000000007 takes 5s (despite being smaller than Bolton's prime - maybe because this has the multiplicity problem) and 23509285402367 times out very early

[Bolton Bailey]

The "Goldilocks prime" 2^64 -2^32 + 1 is another prime people often use for the convenient hardware properties, which is one more than 18446744069414584320 = 2^32×3×5×17×257×65537. I imagine your approach be able to prove this relatively quickly as well?

[Mario Carneiro]

Do those numbers include finding the certificate?

[Bhavik Mehta]

Mario Carneiro said:

Do those numbers include finding the certificate?

No, I found the primitive root 31 here by using wolfram alpha, and the pratt_rule_2' 31 is where I input it

[Mario Carneiro]

so the calculation is just a powmod?

[Bhavik Mehta]

It's a number of powmods, but yeah - really there's not much more than the powmod implementation I showed you in messages a while back

[Mario Carneiro]

are the auxiliary computations all inexpensive?

[Mario Carneiro]

e.g. 31 is prime, those _ arguments (which are presumably divisions of the target number by powers of the factors or something)

[Bhavik Mehta]

I think so, but I haven't benchmarked - all the norm_num calls which aren't pow_mod are natural division (of a very large number by a very small prime factor), natural subtraction, and a natural non zero proof (or some combination)

[Mario Carneiro]

try splitting the powmods into their own lemmas

[Bhavik Mehta]

Sorry one second - I just realised I was using a worse powmod implementation one than what I actually did, let me check properly with the better version. It seems that 18446744069414584321 is too big

[Bhavik Mehta]

image.png Here's the breakdown for the goldilocks prime

[Bhavik Mehta]

In each case the 'tactic execution took' is very close to 80% of the elaboration time showed, so I only showed that

[Bhavik Mehta]

thing0 is where any benefit of 'recursive' certificates might give an improvement, but all the factors here are small so I didn't bother optimising that

[Bhavik Mehta]

Oh, changing thing2 to use 18446744069414584320 rather than 18446744069414584321 - 1 seems to speed it up by 2x... that's not what I expected!

[Bhavik Mehta]

Here's my scratch file: https://github.com/leanprover-community/mathlib/blob/pow_mod/src/pow_mod.lean

[Jakob von Raumer]

That's great, I'll try that on 0x800000000000011000000000000000000000000000000000000000000000001

[Bhavik Mehta]

Instead of proving the thing about nat.factors by norm_num like I did, for that prime you're probably better off showing that the particular list of factors does multiply to p-1, and that each element is a prime, because I guess norm_num might be too slow to verify the larger prime factors of p-1

[Bolton Bailey]

Riccardo Brasca said:

I think it's more helpful to prove docs#lucas_primality using docs#nat.factorization

I guess you mean it would be better if the condition docs#lucas_primality were stated as ∀ q ∈ p.factors, ..., I agree this would look a bit better, and I wouldn't object to someone changing it (although I think docs#factors should probably be named prime_factors).

[Bhavik Mehta]

Why not both? (also ∀ q ∈ (p - 1).factors)

[Riccardo Brasca]

I just mean a formulation such that it's easy to list, in Lean, all the cases.

[Bhavik Mehta]

With quite a bit of effort and broken up lemmas, I proved 47867742232066880047611079 is prime. This prime is interesting because it can't be written as the sum or difference of prime powers (from a 1999 result of Sun), and I verified it's not the sum or difference of a power of two and a prime (a 1975 result)

[Bhavik Mehta]

Perhaps also worth noting that p-1 here also has a large prime factor q = 90101681149415123, and q-1 also has a large prime factor 1123145497, and bare norm_num can't do this one either, so the certificate starts there

[Riccardo Brasca]

Nice work ! Is it somewhere?

[Kevin Buzzard]

That's a very philosophical question!

[Riccardo Brasca]

I mean, on some repository :sweat_smile:

[Bhavik Mehta]

Riccardo Brasca said:

Nice work ! Is it somewhere?

Not yet - if my primality verifying code from earlier is bad, this is atrocious :D