@robot-dreams Always nice to see new contributors in the repo. :) Just curious, are you simply helping or are you interested in getting this merged? If the latter, I'm further curious about your use case. Feel free just ignore my intrusive questions.

@real-or-random Thanks for the reply! I'm simply helping and try to learn about this repo, I don't have a personal use case.

@robot-dreams Thanks for your comments!

The failure with asm=arm is expected, as Apple M1 is 64-bit ARM, while the assembly code used here is 32-bit.

Made another change: in VERIFY mode, the upper bound on the number of posdivsteps iterations is reduced to ~750 (which is close to the median of how many iterations are needed for random input mod the field size), so that the uncomputable path is actually exercised in the tests.

@sipa Thanks for responding to the comments!

• Yes, 32-bit makes sense (in retrospect the fact that the file was called field_10x26 instead of field_5x52 should have given that away)
• Great idea to reduce the number of iterations for validation!

One potential concern is that with this change, there are now 4 slightly different implementations of a similar algorithm (32/64-bit, inv/jacobi). Do you think there's a nice way to consolidate the logic, at least among a given word size? (If there's no nice way or if it's too much trouble, I'm fine with leaving it as is and relying on having good tests.)

Arguably there are 6 implementations (32/64-bit, and const/var/varpos divsteps), and they do share some of the code (in particular the matrix application), but not everything. I don't think there is an obvious way of merging more code without reducing performance, but I also haven't tried very hard.

What was the motivation for using different approaches to clear the low-order bits of g between the 32 and 64 bit implementations (8-bit lookup table vs 4 and 6-bit tricks from Hacker's Delight)?

Just dumb benchmarking. The formula was faster on x86_64 systems we tried, and the table was faster on ARMv7 systems (it's explained in one of the commit messages that introduced it).

I agree that replacing f, g = g, -f with f, g = g, f ensures that all matrix entries are always positive. But aside from actually observing that it converges to f = 1 in practice, do you have any other intuitive justification for this? (The best I can come up with is that since you're constantly dividing g by 2 and swapping the roles of f and g, you're quickly moving towards a solution even if you replace g - f with g + f in the original GCD algorithm.)

I don't, no. I do have a proof that any f>0,g>0 converges to f=g somewhere in a finite number of steps, but the upper bound is something like f*g*log2(f*g), which is obviously useless for practical purposes. The technique explained in https://github.com/sipa/safegcd-bounds completely fails for posdivsteps (doesn't converge at all). With a lot of tricks I got it to converge, but in what appears quadratic (O(log2(f)^2)), and it's only tractable to evaluate for f up to 2^20 or so.

The difficulty is that there are both steps that take f and g further apart (g /= 2) and steps that bring them closer (g = (f+g)/2), and it's hard to infer tight bounds on how many times one can happen in between two executions of the other.

So in short, no; the actual argument why this works is purely empirical; all out of a 2.5 billion random inputs (with the field size modulus) need between 630 and 921 iterations.

Is there already a writeup describing this change? If not, would it make sense to add a brief section to https://github.com/bitcoin-core/secp256k1/blob/master/doc/safegcd_implementation.md? (If you think so, I'd be happy to give that a shot.)

I think that would be valuable, yes. I haven't done so yet, but was considering it. Feel free to take a shot at it.

Included @robot-dreams's update to the writeup.

ACK 53b7ce2

Rebased, and updated to only support coprime inputs in the secp256k1_secp256k1_jacobiXX_maybe_var functions. The return value of those functions for "unknown/timeout" has been changed from -2 to 0.

reACK mod nit

ACK ce3cfc7