Fast modexp

Stream: general

[Bolton Bailey]

For the Miller-Rabin project with Sean, it would be nice if there were an efficient way to compute whether a number is an efficient Miller-Rabin witness to a given base. Currently, it seems the bottleneck is that modular exponentiation is naively evaluated:

import data.zmod.basic
#eval (1234565432124 : zmod 10000000000009) * (12327261237214456) -- fast
#eval (1234324 : zmod 1000000009) ^ (1239382) -- about 1 sec
#eval (1234324 : zmod 1000000009) ^ (12393825) -- about 10 sec
#eval (1234324 : zmod 1000000009) ^ (123938256) -- excessive memory consumption detected at 'vm'

What do you all think will be the best approach for getting #eval to do a fast modexp?

[Kevin Buzzard]

There's a Miller-Rabin project?? You know #eval can't be used in proofs, right?

[Bolton Bailey]

Yes, I know eval can't be used in proofs. It's a research apprenticeship program, and Sean is my mentee. For the sake of the poster session at the end of term, it would be nice to make a graph showing the relative performance of the Miller-Rabin test we're implementing and norm_num

[Bolton Bailey]

The PR is here #12254

[Mario Carneiro]

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)

#eval (1234565432124 : zmod 10000000000009) * (12327261237214456) -- fast
#eval binpow (1234324 : zmod 1000000009) (1239382) -- fast
#eval binpow (1234324 : zmod 1000000009) (12393825) -- fast
#eval binpow (1234324 : zmod 1000000009) (123938256) -- fast

[Bolton Bailey]

Thanks!

[Simon Hudon]

The simplest solution is to write a modexp function with the right performance characteristics. If you're working in Lean 4 (or plan to), I believe (I haven't used it yet) writing csimp lemmas will let you specify rewrite rules for the compiler to transform your ^ expression into a call to your optimized function

[Mario Carneiro]

the issue is that the default implementation of monoid.pow is npow_rec which works by unary recursion. The main downside of making binpow the default is that it has different defeqs

[Mario Carneiro]

By the way, norm_num always uses binary recursion for implementing ^, so it might even be able to beat the #eval implementation for large enough exponent

[Mario Carneiro]

@Bolton Bailey, I updated the code snippet to fix a bug

[Gabriel Ebner]

We should probably give fin and zmod (and nat and int) the fast npow implementation then.

[Mario Carneiro]

should we make it the default for all monoids?

[Gabriel Ebner]

That seems like a good idea.

[Eric Wieser]

Is docs#fin.comm_ring using the default power? (Edit: yes)

[Eric Wieser]

The implementation I expected was the one that makes docs#fin.coe_pow defeq, but that lemma doesn't even exist!