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. _______________________________________________ gmp-devel mailing list gmp-devel@gmplib.org https://gmplib.org/mailman/listinfo/gmp-devel