mpq_cmp
Niels Möller
nisse at lysator.liu.se
Thu May 26 21:47:28 UTC 2016
Hi,
I'm not on the gmp-discuss list, but I've read the thread at the
archive, and I've tried to understand the issues. Please cc me and/or
gmp-devel on any replies.
So we want to compare a/b to c/d, i.e., compute the sign of the cross
difference
a d - c b
We assume gcd(a,b) = gcd(c,d) = 1, i.e., fully reduced inputs. Now, if
it happens that g = gcd(a,c) is large, it can be more efficient to
compute
(a/g) d - (c/g) b
because the gcd is cheap (few euclid or lehmer steps) and divexact with
small quotient is cheap too.
Similarly, if g = gcd(b, d) is large,
a (d/g) - c (b/g)
can be efficient. And we can do both, divide out common factors of the
numerators and common factors of the denominators.
And the reason that mpq_sgn (mpq_sub(x, y)) is faster than mpq_cmp in
some cases is that mpq_sub always computes gcd(b, d) and divides it out,
because that's the natural way to get a fully reduced result.
mpq_cmp can't really do that, because gcd in general is more expensive
than the full product, and there's no simple way to detect a large gcd,
I'm not aware of any better way than to attempt to compute it, and give
up if it doesn't find the gcd in a small numer of steps.
Do I understand things correctly, so far?
Let's make things a little more abstract, and observe that what we're
computing is the sign of the determinant of the matrix (a, b ; c, d).
It's clear that we can do any row or column operation on the matrix,
without changing the result, or multiply it from the left or right by
any matrix M with det M > 0.
E.g., computing the gcd of b, d , with b > d, the first step corresponds
to setting q = floor (b/d) and updating the matrix as
(a, b ; c, d) <--- (a - q c, b - q d; c, d)
The quotient sequence of a/c and b/d start with the same quotient, iff
we also get 0 <= a - qc < c.
And we can determine the sign quickly if we can find a matrix M with
small elements, det > 0, so that multiplying the matrix (a, c; b, d) by
M results in a matrix with some zero element, since then the sign of the
determinant is determined by only the location of that zero and the
signs of the remaining non-zero elements.
For mpq_cmp, I think we can agree that a reasonable approach is to first
check for equality (cheap thanks to fully reduced form). Then rule out
the far from equal case by first comparing bit sizes, and then do a
single word approximation and compare using some error bounds. This will
include lots of cases with large gcd.
The work so far is O(1) work, and if sgn (ad - bc) isn't yet determined,
then the numbers are pretty close, relative difference is bounded by
some small number, 2^{-32} or so.
If we have two distinct but close rationals, how likely is it to have a
large gcd? Is it worth trying to optimize for?
One trick that might be worth trying is to chop off the top two
(normalized) limbs of a and c, and compute M = hgcd2(HI(a), HI(c)). If hgcd2
succeeds, compute M (c, d), and the result should give an indication if
the continued fraction expansion differs early.
Or if hgcd2 fails, that means that either a and c differ a lot in size,
or there's cancellation in a - c.
In the formar case, there's a large quotient, which we may want to take
out. Can we do anything clever in the latter case? Hmm, in the case of
a/b close to c/d, we're going to have similar cancellation also in the
pair b - d, right?
To compare a/b and (a+d1)/(b+d2), with d1 and d2 small, cancellation
means that we need to compute only
sgn (a d2 - b d1)
Hmm, so I'd suggest the following algorithm. As I said above, I have no
idea if the hairier parts makes any real differece, except for specially
chosen inputs.
mpq_cmp (a/b, c/d):
if (b == d)
return sgn (a - c);
else if (a == c)
return sgn (d - b)
bit_diff = bit_size(a) + bit_size(d) - bit_size(b) - bit_size(c)
if (abs(bit_diff) is large) // current test, >= 2 maybe?
return sgn(bit_diff)
high_diff = mul_high_limb(a, d) - mul_high_limb (b, c)
if (abs(high_diff) > max_error)
return sign(high_diff)
while max(a, b, c, d) several limbs larger than min(a, b, c, d)
Find the smallest of a, b, c, d, and the largest of the
(non-diagonal) neighbors. Say this pair is b > d.
q <-- floor(b/d)
b <-- b - qd
a <-- a - qc
if (a < 0 || a >= c)
// Note that we'll be pretty close, since we have already ruled out
// far-from-equal numbers.
quotient sequence differs, terminate.
Try M = hgcd2(HI(a), HI(c))
If hgcd2 fails, that means that either a and c are of very different
size, or that there's a lot of cancellation in a - c. Maybe do
something clever.
If hgcd2 succeeded, set
(a; c) <-- M (a; c) // Positive, and reduced by about one limb
(b; d) <-- M (b; d)
If b, d aren't similarly reduced by this operation, that means that
the quotient sequence differs, so figure out how and terminate.
If (a == c || b == d)
// A large gcd found
Compare the un-equal pair of elements and terminate.
// Compute full products. If practical, high limb first with early
// termination.
return sgn (a d - b c)
Regards,
/Niels
--
Niels Möller. PGP-encrypted email is preferred. Keyid C0B98E26.
Internet email is subject to wholesale government surveillance.
More information about the gmp-discuss
mailing list