#1482: xgcd suboptimal output
David Harvey
I've been looking at
http://www.sagetrac.org/sage_trac/ticket/1482
and approximately diagnosed the problem (see comments on the ticket),
but it's not clear to me exactly how to proceed.
The new underlying gcd code produces quite inscrutable output, for
example:
sage: xgcd(3, 3)
(3, 0, 1)
sage: xgcd(3, 6)
(3, -3, 2)
sage: xgcd(3, 9)
(3, -8, 3)
sage: xgcd(3, 12)
(3, -3, 1)
sage: xgcd(3, 15)
(3, -14, 3)
sage: xgcd(3, 18)
(3, -17, 3)
sage: xgcd(3, 21)
(3, -20, 3)
sage: xgcd(3, 24)
(3, -15, 2)
sage: xgcd(3, 27)
(3, -26, 3)
sage: xgcd(3, 30)
(3, -29, 3)
sage: xgcd(3, 33)
(3, -32, 3)
sage: xgcd(3, 36)
(3, -35, 3)
sage: xgcd(3, 39)
(3, -38, 3)
sage: xgcd(3, 42)
(3, -41, 3)
sage: xgcd(3, 45)
(3, -44, 3)
sage: xgcd(3, 48)
(3, -15, 1)
sage: xgcd(3, 51)
(3, -50, 3)
It's very hard to see any rhyme or reason in that. Certainly the GMP
documentation doesn't guarantee any properties of the cofactors,
apart from bezout, either at the mpz or mpn levels, so it's not a bug.
I'm inclined to do the following in Sage. Provide a "fast" flag which
just calls GMP and doesn't provide any guarantees about the output
(apart from satisfying bezout). If "fast" is False (default), then we
convert to a "standard" form, which should be unique, mathematically
sensible, and carefully documented.
But having thought about it for a while, it's not even clear to me
what exactly the right definition of "standard" should be. I don't
quite understand the minimality condition you proposed.
Apart from the kind of special case listed above, GMP does seem to
follow some "rules", but I'm not sure they are always satisfied, and
we can't rely on them since the GMP documentation doesn't promise
anything.
For example, if 0 < a < b, then xgcd(a, b)[1] always seems to be
negative and xgcd(a, b)[2] always positive, but it swaps around if 0
< b < a:
sage: xgcd(11, 17)
(1, -3, 2)
sage: xgcd(17, 11)
(1, 2, -3)
i.e. it looks like it swaps them to make a < b, and then swaps the
cofactors at the end.
It generally seems to handle input signs in the same way, i.e. it
first does xgcd(|a|, |b|), then throws the signs back in at the end:
sage: xgcd(11, 17)
(1, -3, 2)
sage: xgcd(-11, 17)
(1, 3, 2)
sage: xgcd(11, -17)
(1, -3, -2)
sage: xgcd(-11, -17)
(1, 3, -2)
But none of this is particularly canonical.
What's the "right" way to do this?
david
Nils Bruin
connection between continued fractions and Euclid's algorithm (which
does guarantee minimality), I guess you could try and see what signs
would be given back by a continued fractions approach.
It actually looks like they had a very good reason for departing from
returning minimal solutions in GMP's xgcd. It is nice to have an
option for minimality, but it should perhaps not be the default. The
way I ran into this problem was that I tested for one of the
parameters to be zero to detect if the gcd is equal to one of the
inputs. This example shows you're better off doing that by testing for
equality on the GCD and the inputs.
The "bug" should perhaps be the failure of the documentation to point
out that minimality is not guaranteed (people will naively expect
Euclid's algorithm)