_basecase or _sec? [
Torbjorn Granlund
tg at gmplib.org
Fri May 3 14:40:49 CEST 2013
nisse at lysator.liu.se (Niels Möller) writes:
I think Newton analogues exist only when b is a power, not in general.
And the most important case is prime b.
I think it exists also for b can be factorised into prime powers...
> I am not familiar with the Jebelean (or Möller!) criteria for gcd, but
> my intuition suggests that there exists a O(k)-bit matrix which reduces
> k bits from the operands. And the same intuition suggests that it can
> be computed from O(k) high bits. Am I wrong?
You need to have roughly *2k* bits of each operand to make k bits
progress. And you get into trouble at least in the case that those 2k
bits are equal (you need to subtract one operand from the other, but
which depends on the other bits).
We don't need to insist on keeping operands positive.
Let me think aloud.
As a start, one should try to construct an algorithm which makes (at
least) two bits of progress per iteration, with all problematic special
cases handled in a side-channel silent manner. I looked into it briefly
before writing my current code, but I couldn't find anything easy.
Two bits per bignum operation will give a 100% speedup, and still leave
a lot of headroom for future GMP releases. :-)
And I'm afraid the table lookup even for just two bits of progress needs
5 bits from each operand, i.e., 1024 entries, which might leak
information through the cache instead.
My thought was not to have a table, but plain loop which executed a
limb-size dependent number of iterations, and has not internal branches.
Right-to-left algorithms will likely have easier bookkeeping than
left-to-right algorithms.
I have not looked at right-to-left algorithms for gcdext, at least not
in a very long time. Do they resolve gcd(a,b) = ax + by or something
like gcd(a,b)*2^t = ax + by?
Instead of a lehmer-like algorithms which eliminates k bits of each
operand using a vector/matrix multiply, it might be better (smaller
table) to eliminate k bits from the larger operand only. But then we
need to keep track of which operand is largest. I think this simplest
appproach was the one I looked a bit into, and it didn't seem very
promising. I think one can always eliminate two bits, by setting
m <-- min(a,b)
M <-- max(a,b)
a <-- (M + k m) / 4
b <-- m
where k = ± 1, which ever is needed to make M + k m a multiple of 4. And
then some additional logic for the case of even numbers (if M = 0 (mod
4), just set k = 0, and if M = 2 (mod 4), set k = 2). But to go beyond
two bits, one would need something like k-ary gcd, which might introduce
false factors.
Right, answering my question above. Cannot false factors of 2 be
handled?
The k \in {0, 1, 2, -1} could be a conditional subtraction and addition
1. Subtract: d <-- a - [a odd] b (if a odd)
2. Swap: b <-- a, d <-- -d (if underflow above, d < 0)
3. Add: a <-- a + 2 b (under certain conditions)
4. Shift: a <-- a / 4 Always
Current code does (1) and (2) only (and a shift by one bit).
If one wants a single addmul_1, one could choose k \in {0,1,2,3}
instead, but one still *needs* to take sign (a-b) into account; it's no
use to add the larger operand to the smaller. Analysis of the guaranteed
progress after n iterations is going to get more complicated with k > 0.
For any variant, I think the iteration needs to take into account:
Low bits of a and b (one can try to maintain b odd at all times, but
one needs to handle the case of even a). Or high bits for a
left-to-right algorithm, but I think that's going to be less
convenient.
Sign of (a - b).
The straight-forward way to determine the sign of (a - b) is to actually
compute the difference. The result of this subtraction is directly
useful for the "plain" binary algorithm, but a bit wasted if we try
something more sophisticated.
I suspect non-leaky code has to compute the sign as a full subtraction.
--
Torbjörn
More information about the gmp-devel
mailing list