mpq_cmp

[Niels Möller]

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.