Opened 11 years ago
Closed 10 years ago
#1482 closed defect (fixed)
[with patch, with 2 positive reviews] xgcd suboptimal output
Reported by: | nbruin | Owned by: | dmharvey |
---|---|---|---|
Priority: | major | Milestone: | sage-2.9.2 |
Component: | basic arithmetic | Keywords: | |
Cc: | Merged in: | ||
Authors: | Reviewers: | ||
Report Upstream: | Work issues: | ||
Branch: | Commit: | ||
Dependencies: | Stopgaps: |
Description
I was expecting xgcd(x,y) to produce u,v of minimal absolute value such that u*x+v*y=gcd(x,y). However, presently it doesn't in the edge case where x divides y or vice versa (i.e., where (u,v)=(1,0) or (u,v)=(0,1) would be valid)
sage: xgcd(2,4) (2, -1, 1) sage: xgcd(4,2) (2, 1, -1) sage: xgcd(12,2) (2, 1, -5)
This is not a bug in the strictest sense of the word, since the documentation does not promise u,v to be minimal, but for a lot of users minimal (u,v) would be the expected behaviour.
As far as I can see, this is directly the answer of mpz_gcdext. Possible solutions:
- fix mpz_gcdext
- test in sage.rings.integer.Integer._xgcd if the gcd equals self.abs() or right.abs() and return (u,v)=(+-1,0) or (u,v)=(0,+-1) as appropriate (this is cython code
- do this test in the python wrapper sage.rings.integer.Integer.xgcd instead.
Since xgcds are so crucial, I am hesitant to fix it myself. Can someone who knows the code well see what the appropriate solution is?
Attachments (1)
Change History (9)
comment:1 Changed 11 years ago by
- Milestone set to sage-2.9
comment:2 Changed 10 years ago by
- Milestone changed from sage-2.9.1 to sage-2.9.2
- Owner changed from somebody to dmharvey
- Status changed from new to assigned
comment:3 Changed 10 years ago by
From Nils
I don't think anybody should care about the signs. Given the close 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)
comment:4 Changed 10 years ago by
Earlier comments in the above thread:
Changed 10 years ago by
comment:5 Changed 10 years ago by
- Summary changed from xgcd suboptimal output to [with patch, needs review] xgcd suboptimal output
comment:6 Changed 10 years ago by
- Summary changed from [with patch, needs review] xgcd suboptimal output to [with patch, with positive review] xgcd suboptimal output
Interesting...
This looks good, and good documentation too.
comment:7 Changed 10 years ago by
- Summary changed from [with patch, with positive review] xgcd suboptimal output to [with patch, with 2 positive reviews] xgcd suboptimal output
comment:8 Changed 10 years ago by
- Resolution set to fixed
- Status changed from assigned to closed
Merged in 2.9.2.rc1.
I've started investigating. The issue is that since we included the fast gcd code of Niels Möller (see #424), the behaviour of
mpn_gcdext
has changed slightly, and the changes flows ontompz_gcdext
. Note that this is not a bug in the code, since neither the mpz nor mpn documentation claim any minimality condition on the cofactors.Here is some mpn-level code which produces different answers under GMP 4.2.1 and the patched version:
Vanilla GMP 4.2.1 prints -1, patched GMP prints 1.
I plan to fix this as follows. I think it would be good for Sage to guarantee minimality. I'm going to make the xgcd function call mpz_gcdext, then check for the "obvious" non-minimality issues and make corrections for those, and finally fall back on a (slower) but guaranteed algorithm if that doesn't work. Hopefully the latter should never actually be called, since it will need to divide by something somewhere.