Binary XGCD

[erik-3milabs]

Continuation of the work started in #755

[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 😅 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 😅

[erik-3milabs]

@kayabaNerve I think I found the culprit 🙈 As it turns out, in rare circumstances (like yours) this does not work 😅

I knew this was sketchy, and did it anyway. Next time, I'll think twice before committing something I know is sketchy.

[erik-3milabs]

FYI: the inputs you provided trigger a particular situation somewhere along the way where a > b but a.compact() < b.compact. As a result, this matrix multiplication should be doing something different than it is now.

[kayabaNerve]

Fun!

Is the fix as simple as reverting back to the top/bottom row negation, instead of the entire matrix negation, or is it unfortunately a bit more annoying to navigate?

[erik-3milabs]

Sadly, solving this problem is more complicated than that ☹️

Reverting to the old version is not really an option, since that version was not capable of computing the xgcd for Uints that had their msb set: I needed that top bit to represent the sign of the matrix elements. I solved this using an interesting trick I previously described here. It turns out that the input you presented is an exception to "the trick" I described there.

We can, however, use the fact that a > b and a.compact() < b.compact() to conclude that the top K-1 bits of both a and b are the same. This implies that just subtracting b from a already allows the latter to shrink by K-1 bits. We can, therefore, replace the outcome of the (..., update_matrix) = partial_xgcd(...) with a simple ((1, -1),(0,1)) matrix and still be guaranteed that the algorithm is done after all steps.

Sadly, there are some intricate details I skipped over that make it a tad more complicated 😅

1. The mirrored situation can also take place: a < b and a.compact() > b.compact()
2. Subtracting b from a can make the latter even, which means we have to swap them.
3. We want the solution to be fast to execute.

[erik-3milabs]

@kayabaNerve, the new commits I added should fix the bug. Please try it out and check whether it works on your end as well!

[erik-3milabs]

I'm happy to report the patch has a negligible impact on performance: ~1% 👌

[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 🙈 )

[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 🤔

[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