version 1.1.1.1, 2000/01/10 15:35:27 |
version 1.1.1.2, 2000/09/09 14:12:52 |
|
|
/* mpz_gcdext(g, s, t, a, b) -- Set G to gcd(a, b), and S and T such that |
/* mpz_gcdext(g, s, t, a, b) -- Set G to gcd(a, b), and S and T such that |
g = as + bt. |
g = as + bt. |
|
|
Copyright (C) 1991, 1993, 1994, 1995 Free Software Foundation, Inc. |
Copyright (C) 1991, 1993, 1994, 1995, 1996, 1997, 2000 Free Software |
|
Foundation, Inc. |
|
|
This file is part of the GNU MP Library. |
This file is part of the GNU MP Library. |
|
|
The GNU MP Library is free software; you can redistribute it and/or modify |
The GNU MP Library is free software; you can redistribute it and/or modify |
it under the terms of the GNU Library General Public License as published by |
it under the terms of the GNU Lesser General Public License as published by |
the Free Software Foundation; either version 2 of the License, or (at your |
the Free Software Foundation; either version 2.1 of the License, or (at your |
option) any later version. |
option) any later version. |
|
|
The GNU MP Library is distributed in the hope that it will be useful, but |
The GNU MP Library is distributed in the hope that it will be useful, but |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |
License for more details. |
License for more details. |
|
|
You should have received a copy of the GNU Library General Public License |
You should have received a copy of the GNU Lesser General Public License |
along with the GNU MP Library; see the file COPYING.LIB. If not, write to |
along with the GNU MP Library; see the file COPYING.LIB. If not, write to |
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
MA 02111-1307, USA. */ |
MA 02111-1307, USA. */ |
|
|
|
#include <stdio.h> /* for NULL */ |
#include "gmp.h" |
#include "gmp.h" |
#include "gmp-impl.h" |
#include "gmp-impl.h" |
|
|
/* Botch: SLOW! */ |
|
|
|
void |
void |
#if __STDC__ |
#if __STDC__ |
mpz_gcdext (mpz_ptr g, mpz_ptr s, mpz_ptr t, mpz_srcptr a, mpz_srcptr b) |
mpz_gcdext (mpz_ptr g, mpz_ptr s, mpz_ptr t, mpz_srcptr a, mpz_srcptr b) |
Line 37 mpz_gcdext (g, s, t, a, b) |
|
Line 37 mpz_gcdext (g, s, t, a, b) |
|
mpz_srcptr b; |
mpz_srcptr b; |
#endif |
#endif |
{ |
{ |
mpz_t s0, s1, q, r, x, d0, d1; |
mp_size_t asize, bsize, usize, vsize; |
|
mp_srcptr ap, bp; |
|
mp_ptr up, vp; |
|
mp_size_t gsize, ssize, tmp_ssize; |
|
mp_ptr gp, sp, tmp_gp, tmp_sp; |
|
mpz_srcptr u, v; |
|
mpz_ptr ss, tt; |
|
__mpz_struct stmp, gtmp; |
|
TMP_DECL (marker); |
|
|
mpz_init_set_ui (s0, 1L); |
TMP_MARK (marker); |
mpz_init_set_ui (s1, 0L); |
|
mpz_init (q); |
|
mpz_init (r); |
|
mpz_init (x); |
|
mpz_init_set (d0, a); |
|
mpz_init_set (d1, b); |
|
|
|
while (d1->_mp_size != 0) |
/* mpn_gcdext requires that U >= V. Therefore, we often have to swap U and |
|
V. This in turn leads to a lot of complications. The computed cofactor |
|
will be the wrong one, so we have to fix that up at the end. */ |
|
|
|
asize = ABS (SIZ (a)); |
|
bsize = ABS (SIZ (b)); |
|
ap = PTR (a); |
|
bp = PTR (b); |
|
if (asize > bsize || (asize == bsize && mpn_cmp (ap, bp, asize) > 0)) |
{ |
{ |
mpz_tdiv_qr (q, r, d0, d1); |
usize = asize; |
mpz_set (d0, d1); |
vsize = bsize; |
mpz_set (d1, r); |
up = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB); |
|
vp = (mp_ptr) TMP_ALLOC ((vsize + 1) * BYTES_PER_MP_LIMB); |
|
MPN_COPY (up, ap, usize); |
|
MPN_COPY (vp, bp, vsize); |
|
u = a; |
|
v = b; |
|
ss = s; |
|
tt = t; |
|
} |
|
else |
|
{ |
|
usize = bsize; |
|
vsize = asize; |
|
up = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB); |
|
vp = (mp_ptr) TMP_ALLOC ((vsize + 1) * BYTES_PER_MP_LIMB); |
|
MPN_COPY (up, bp, usize); |
|
MPN_COPY (vp, ap, vsize); |
|
u = b; |
|
v = a; |
|
ss = t; |
|
tt = s; |
|
} |
|
|
mpz_mul (x, s1, q); |
tmp_gp = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB); |
mpz_sub (x, s0, x); |
tmp_sp = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB); |
mpz_set (s0, s1); |
|
mpz_set (s1, x); |
if (vsize == 0) |
|
{ |
|
tmp_sp[0] = 1; |
|
tmp_ssize = 1; |
|
MPN_COPY (tmp_gp, up, usize); |
|
gsize = usize; |
} |
} |
|
else |
|
gsize = mpn_gcdext (tmp_gp, tmp_sp, &tmp_ssize, up, usize, vp, vsize); |
|
ssize = ABS (tmp_ssize); |
|
|
if (t != NULL) |
PTR (>mp) = tmp_gp; |
|
SIZ (>mp) = gsize; |
|
|
|
PTR (&stmp) = tmp_sp; |
|
SIZ (&stmp) = (tmp_ssize ^ SIZ (u)) >= 0 ? ssize : -ssize; |
|
|
|
if (tt != NULL) |
{ |
{ |
mpz_mul (x, s0, a); |
if (SIZ (v) == 0) |
mpz_sub (x, d0, x); |
SIZ (tt) = 0; |
if (b->_mp_size == 0) |
|
t->_mp_size = 0; |
|
else |
else |
mpz_tdiv_q (t, x, b); |
{ |
|
mpz_t x; |
|
MPZ_TMP_INIT (x, ssize + usize + 1); |
|
mpz_mul (x, &stmp, u); |
|
mpz_sub (x, >mp, x); |
|
mpz_tdiv_q (tt, x, v); |
|
} |
} |
} |
mpz_set (s, s0); |
|
mpz_set (g, d0); |
if (ss != NULL) |
if (g->_mp_size < 0) |
|
{ |
{ |
g->_mp_size = -g->_mp_size; |
if (ALLOC (ss) < ssize) |
s->_mp_size = -s->_mp_size; |
_mpz_realloc (ss, ssize); |
if (t != NULL) |
sp = PTR (ss); |
t->_mp_size = -t->_mp_size; |
MPN_COPY (sp, tmp_sp, ssize); |
|
SIZ (ss) = SIZ (&stmp); |
} |
} |
|
|
mpz_clear (s0); |
if (ALLOC (g) < gsize) |
mpz_clear (s1); |
_mpz_realloc (g, gsize); |
mpz_clear (q); |
gp = PTR (g); |
mpz_clear (r); |
MPN_COPY (gp, tmp_gp, gsize); |
mpz_clear (x); |
SIZ (g) = gsize; |
mpz_clear (d0); |
|
mpz_clear (d1); |
TMP_FREE (marker); |
} |
} |