Binary GCD

[erik-3milabs]

This PR introduces an implementation of Optimized Binary GCD. Ref: Pornin, Algorithm 2.

Upsides to this technique:

• it is up to 27x faster than the gcd algorithm currently implemented in crypto_bigint (see below) (really, it just has a different complexity bound).
• does not need UNSAT_LIMBS
• it is actually constant time (unlike the current implementation, which is sneakily vartime in the maximum of the bitsizes of the two operands).

Benchmark results

Word = u64

limbs	gcd (vt)	gcd (ct)	new_gcd (ct)
2	10.687 µs	20.619 µs	3.6090 µs
4	29.121 µs	56.433 µs	7.1124 µs
8	99.819 µs	195.02 µs	16.184 µs
16	359.39 µs	710.26 µs	44.294 µs
32	1.6804 ms	3.3097 ms	136.49 µs
64	6.9717 ms	13.028 ms	494.16 µs
128	29.099 ms	57.325 ms	2.3335 ms
256	143.22 ms	244.89 ms	8.7722 ms

[dolevmu]

Looks great overall. I find the implementation constant time. Offered a few annotations and documentations to clarify the code after going over it, and suggested some changes to boost performance. I think it is worth reading and considering, but I approve this.

Great work thanks @erik-3milabs !

[dolevmu]

Good question, but looks like you answered it, your problem is that the answer may depend on the machine. Suppose I gave you an answer, do you know what to do in the code? If I say "this is optimal for 64-bit architecture and that for 32-bit".. If so, let's check on another machine.

[dolevmu]

According to Algorithm 2, you multiply by 2 always, but I suppose you know what you are doing and he is wrong on this edge case.

[dolevmu]

Also, it multiplies the second row and not the first one (f1,g1) is multiplied by two.

[erik-3milabs]

1. This is a modification I made to the algorithm. It prevents integer overflows in an upcoming PR. I've updated the docstring of this function to indicate this.
2. yes, but this code has a and b swapped, so I should be doubling the first row ;-)

Note: that doubling problem was fixed because I reverted the (a, b) swap.

[erik-3milabs]

Turns out, this is wrong: this leads to buggy behavior for specific pairs of numbers. I'm refactoring some code to bypass this issue.

[erik-3milabs]

Fixed in #761

[dolevmu]

In light of my previous comment, it is worth to consider if we can efficiently select between two vectors or matrices, and similarly swap columns/rows efficiently (which you did implement).

[erik-3milabs]

@tarcieri do you have any ideas on improving the performance of this and other swapping operations?

[erik-3milabs]

@tarcieri what do you think of this line? Previously, it was like this:

a = Uint::select(&a, &a.wrapping_sub(&b), a_odd);

The current code is 25-10% faster for Uints with few limbs (1, 2, 3, etc.)

[tarcieri]

I'm surprised there's that much of a difference. Are you sure it's always faster or is it faster depending on a?

[erik-3milabs]

I'm pretty confident it is always faster.

The reason I think this is faster, is that we are now selecting between a constant and a variable, instead of two variables. Given that select is const and loves to be inlined, the compiler can now optimize the select operation.

Recall, Uint::select calls Limb::select, which in turn calls

impl ConstChoice {
    /// Return `b` if `self` is truthy, otherwise return `a`.
    #[inline]
    pub(crate) const fn select_word(&self, a: Word, b: Word) -> Word {
        a ^ (self.0 & (a ^ b))
    }
}

When a is the constant ZERO, this can be optimized as:

        self.0 & b

saving two XOR operations, or 2/3's of this operation.

Returning to the gcd subroutine, this select is in the hot loop of this algorithm. In total, the loop executes:

• Uint::is_odd (1 op)
• Uint::lt (2 ops/word),
• ConstChoice::and (1 op),
• Uint::wrapping_sub (4 ops/word),
• Uint::select (3 ops/word -> 1 ops/word)
• Uint::shr (3 ops/word)

So, there is a reduction from 12 to 10 ops/word, or a 17% improvement.

[tarcieri]

When a is the constant ZERO

Aah, ok

[tarcieri]

@erik-3milabs I just set master to v0.7.0-pre in #765

This means you can now make breaking changes, such as removing the existing safegcd implementation and changing trait impls like Gcd and InvMod to use bingcd instead

[erik-3milabs]

@erik-3milabs I just set master to v0.7.0-pre in #765

This means you can now make breaking changes, such as removing the existing safegcd implementation and changing trait impls like Gcd and InvMod to use bingcd instead

@tarcieri While I could still modify the Gcd trait to use this algorithm, this PR does not yet introduce the tools necessary to replace InvMod.

Aside from that, what else would be required to see this PR merged?

[tarcieri]

Aah, lack of invmod support would definitely be a problem. Is it something you plan on addressing eventually? My understanding is, like safegcd, that invmod is a big part of binary GCD's intended usage.

It seems a little weird to have multiple implementations of GCD algorithms which effectively do the same thing, though to completely replace safegcd in addition to invmod you'd also need support for boxed types (though with const_mut_refs stable it should be a lot easier to share an implementation).

Also seems it needs a rebase due to upstream changes.

[erik-3milabs]

Aah, lack of invmod support would definitely be a problem. Is it something you plan on addressing eventually? My understanding is, like safegcd, that invmod is a big part of binary GCD's intended usage.

You're right. This PR only introduces the gcd algorithm; PR #761 extends the algorithm into xgcd. Stripping some things from the xgcd algorithm gives invmod. Given that I don't need invmod myself, I am not too keen on implementing it 🙈

It seems a little weird to have multiple implementations of GCD algorithms which effectively do the same thing, though to completely replace safegcd in addition to invmod you'd also need support for boxed types (though with const_mut_refs stable it should be a lot easier to share an implementation).

I agree that having two algorithms is overkill. Let me see about implementing this for Boxed<X> as well.

Also seems it needs a rebase due to upstream changes.

Yeah, you're right. Let me address that right away.