Opened 10 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)

1482.hg (2.4 KB) - added by dmharvey 10 years ago.

Download all attachments as: .zip

Change History (9)

comment:1 Changed 10 years ago by nbruin

  • Milestone set to sage-2.9

comment:2 Changed 10 years ago by dmharvey

  • Milestone changed from sage-2.9.1 to sage-2.9.2
  • Owner changed from somebody to dmharvey
  • Status changed from new to assigned

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 onto mpz_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:

   mp_limb_t s1p = 4;
   mp_limb_t s2p = 2;
   mp_limb_t r2p[2];
   mp_size_t r2n;
   mp_limb_t r1p[2];
   mp_size_t r1n;
   
   r1n = mpn_gcdext(r1p, r2p, &r2n, &s1p, 1, &s2p, 1);
   
   printf("first cofactor = %ld\n", (long)(r2p[0]) * ((r2n > 0) ? 1 : -1));

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.

comment:3 Changed 10 years ago by was

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)

Changed 10 years ago by dmharvey

comment:5 Changed 10 years ago by dmharvey

  • Summary changed from xgcd suboptimal output to [with patch, needs review] xgcd suboptimal output

comment:6 Changed 10 years ago by robertwb

  • 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 was

  • 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 mabshoff

  • Resolution set to fixed
  • Status changed from assigned to closed

Merged in 2.9.2.rc1.

Note: See TracTickets for help on using tickets.