Binary XGCD

[erik-3milabs]

Continuation of the work started in #755

[erik-3milabs]

Hmmm. It looks like the miri target is on the slow side...

It is probably getting stuck on the randomized tests. Any proposals on how to approach this such that we don't overload miri?

[tarcieri]

I believe you can put #[cfg(not(miri))] on anything that's particularly slow

[kayabaNerve]

For some 200-bit inputs, after hundreds? of calls, I get invalid x, y coefficients. Specifically, both are negative, for the positive gcd = 2. This is with Uint::binxgcd(a, b) (so a, b are positive and at least either x or y should also be). I'm unsure what trap my inputs are following into (as I have tests for small numbers with gcd != 1) and can't determine their pattern, solely their infrequency.

[erik-3milabs]

For some 200-bit inputs, after hundreds? of calls, I get invalid x, y coefficients. Specifically, both are negative, for the positive gcd = 2. This is with Uint::binxgcd(a, b) (so a, b are positive and at least either x or y should also be). I'm unsure what trap my inputs are following into (as I have tests for small numbers with gcd != 1) and can't determine their pattern, solely their infrequency.

@kayabaNerve Hmmm, that is odd indeed. Thank you for letting me know! I don't have a clue about what is going on here, so I'll have to investigate this :)

If you happen to find an input pair (x,y) that gives you this issue, I would love to hear about them; such a pair will make my life 100x easier :heart:

Also, can you tell me what number of limbs you're using? Are you using U256s? Are you using an 32-bit or 64-bit platform?

[kayabaNerve]

I was experimenting with this (understanding it still is to-be-merged and not expecting it to be final) when my code calling stopped working. I was reviewing the exact issue for a good bit before I noticed this :sweat_smile: Since I didn't see it discussed, I wanted to ensure it was raised.

I'll post the pair in a moment.

[kayabaNerve]

    let (gcd, u, v) =   (crypto_bigint_xgcd::U512::from_be_hex("0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001B5DFB3BA1D549DFAF611B8D4C"))
        .extended_gcd(crypto_bigint_xgcd::U512::from_be_hex(
          "0000000000000000000000000000000000000000000000000000000000000000000000000000345EAEDFA8CA03C1F0F5B578A787FE2D23B82A807F178B37FD8E",
        ));

u is also not the multiplicative inverse of (a / g) % (b / g). I spent so long trying to determine why u was invalid when I noticed the more glaring issue of them both being negative, and figured that'd be more notable to flag (ideally leading to a singular root cause).

This is on a 64-bit platform with the latest commit on your branch (from 13 hours ago). My snippet is to my wrapper fn (as I prior had defined a wrapper around gcd and was experimenting with moving it to bingcd), yet calling Uint::binxgcd with those values should produce the results I observed.

[kayabaNerve]

Sorry, I think that may be the wrong pair. If so, please try:

self = Uint(0x00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001A0DEEF6F3AC2566149D925044)

    other = Uint(0x00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000072B69C9DD0AA15F135675EA9C5180CF8FF0A59298CFC92E87FA)

Those serialization are their 1024-bit form yet it has the same issue even with just a 512-bit Uint.

[erik-3milabs]

Sorry, I think that may be the wrong pair.

The first pair you sent is already producing an incorrect output on my end. Still, thanks for the second pair!

[kayabaNerve]

Happy to hear I was able to help! Really excited for this as someone who spends 50% of my runtime on safegcd :sweat_smile:

[kayabaNerve]

The new commit, "Fix bug", which I have not personally reviewed, passes my personal code! Thank you!

[kayabaNerve]

crypto-bigint/src/modular/bingcd/xgcd.rs:309:13:
b is never negative

This is from the "Fix bug" commit I prior commented worked for my prior noted issues. I'm actually using it in the same context, solely running a different test suite, so I'm unsure why this current test suite is managing to trip this low-level bug when my prior one didn't. I'll try to find a reproducible test case, but at least having noted this ideally lets some review be started on.

EDIT: binxgcd on the following Uints, which are of size of the size of their serializations. 64-bit host.

Uint(0x00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000007FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD03641424EB38E6AC0E34DE2F34BFAF22DE683E1F4B92847B6871C780488D797042229E1)

Uint(0x0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFD755DB9CD5E9140777FA4BD19A06C82839D671CD581C69BC5E697F5E45BCD07C52EC373A8BDC598B4493F50A1380E1281)

[erik-3milabs]

@kayabaNerve, thank you for being a persistent tester! I hope to find some time later this week to debug the issue you presented. I appreciate your patience!

(I'm delighted I added in those .excepts: I now know exactly where to look for the issue :see_no_evil: )

[kayabaNerve]

Also just hit a is never negative. Unfortunately, that seems to be much more infrequent for the numbers I'm generating, so I may not be able to provide a vector.

[erik-3milabs]

crypto-bigint/src/modular/bingcd/xgcd.rs:309:13:
b is never negative

This is from the "Fix bug" commit I prior commented worked for my prior noted issues. I'm actually using it in the same context, solely running a different test suite, so I'm unsure why this current test suite is managing to trip this low-level bug when my prior one didn't. I'll try to find a reproducible test case, but at least having noted this ideally lets some review be started on.

EDIT: binxgcd on the following Uints, which are of size of the size of their serializations. 64-bit host.

Uint(0x00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000007FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD03641424EB38E6AC0E34DE2F34BFAF22DE683E1F4B92847B6871C780488D797042229E1)

Uint(0x0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFD755DB9CD5E9140777FA4BD19A06C82839D671CD581C69BC5E697F5E45BCD07C52EC373A8BDC598B4493F50A1380E1281)

@kayabaNerve I had a quick look today. Your counterexample demonstrates there are more annoying ways in which the compact_a/compact_b variables misrepresent the actual a and b. I will have to think a bit about how to resolve this :thinking:

[kayabaNerve]

I did try poking at this myself to see if I could fix it, but it really wasn't immediately familiar/obvious to me. I ended up just making a fork which never uses the optimized algorithm as to not block my work.

Pros: Works. Cons: My runtime on GCDs went from 2s to ~20s when I have a deadline for an academic paper in just a few days 😬

I may try to work on this again if I myself have free time before the deadline, but any further help from you would truly be appreciated! This truly is a massive improvement for crypto-bigint on this topic and it's fundamentally what makes my current research topic go from 'initial impl for reference which isn't discussable' to 'initial impl which yields feasible results for practical deployment'.

EDIT: The Bezout coefficients are correct or their additive inverse modulo the other value divided by their GCD, for this specific test case. From these incorrect values, it's trivial to calculate the proper values with a post-processing run.

I tried to do so on every result, to restore functioning despite not having a proper fix, but that yields other errors so this isn't a universal fact.

[erik-3milabs]

I've worked on this here-and-there for the past couple days, without much luck.

The problem

The failing input presented by @kayabaNerve illustrates that the trick does not always apply.

Solutions

I've thought of a couple solutions:

1) Using Int's in BinXgcdMatrix.

Not an option, as integer overflows can happen when the input values have their most significant bit set.

2) Using ExtendedInt's in BinXgcdMatrix.

While possible on paper, this makes the multiplication between matrix and update_matrix and absolute nightmare.

3) Adding an second pattern attribute to BinXgcdMatrix.

In particular, the idea would be to have a top_pattern and bottom_pattern attribute. The first would indicate whether the signs in the first row of the matrix are either [ - + ] or [+ - ] and the bottom_pattern would do the same for the bottom row.

This also does not seem to work, as there are exceptional cases where the signs in a row are BOTH positive. Specifically, this seems to happen for inputs a, b such that for some small x and y it holds that x * a + y * b is divisible by 2^k >> x, y.

4) Solution 1, but now with an extra limb.

That more-or-less requires us to keep UNSAT_LIMBS around which is exactly what we're trying to avoid.

Now what?

I don't know yet. If you're working with Uints that do not have their top bit set, you could opt for the first solution. But that is a temporary solution at best.

I'll be away for the next couple weeks; I hope to return to this issue afterwards.

[tarcieri]

Solution 1, but now with an extra limb.

That more-or-less requires us to keep UNSAT_LIMBS around which is exactly what we're trying to avoid.

@erik-3milabs if it's just a single extra limb, with some effort you could always make a struct that adds that extra limb which avoids type-level size calculations

[tarcieri]

@erik-3milabs looks like this one needs a rebase

[erik-3milabs]

Alright, it took some time, but I think I managed to fix the bug.

@erik-3milabs if it's just a single extra limb, with some effort you could always make a struct that adds that extra limb which avoids type-level size calculations

@tarcieri this was indeed the solution. I thought this couldn't work, because you'd have to multiply two of these ExtraLimbInts, but it turns out that is not necessary!

@erik-3milabs looks like this one needs a rebase Because we merged bingcd, this now needs a lot more than just a rebase :sweat_smile:

To ease the review process, I propose to split this PR into (at least) two parts:

1. classic_binxgcd + all the boilerplate code to provide it to Uint, NonZeroUint, OddUint, Int, NonZeroInt and OddInt; this one won't be all that complicated, since all the parts are already in place for this.
2. optimized_binxgcd, which is a rather complicated algorithm that needs many eyes to go over it. Let's not bog the reviewer down with boiler plate code for that PR.

As such, this PR will remain a draft for a bit longer.

[erik-3milabs]

The content of this PR is moved to #854 and #856

[erik-3milabs]

Closing this PR; it is superseded by #854 and #856