Annotation of OpenXM_contrib/gmp/mpz/gcdext.c, Revision 1.1.1.3
1.1 maekawa 1: /* mpz_gcdext(g, s, t, a, b) -- Set G to gcd(a, b), and S and T such that
2: g = as + bt.
3:
1.1.1.3 ! ohara 4: Copyright 1991, 1993, 1994, 1995, 1996, 1997, 2000, 2001 Free Software
1.1.1.2 maekawa 5: Foundation, Inc.
1.1 maekawa 6:
7: This file is part of the GNU MP Library.
8:
9: The GNU MP Library is free software; you can redistribute it and/or modify
1.1.1.2 maekawa 10: it under the terms of the GNU Lesser General Public License as published by
11: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1 maekawa 12: option) any later version.
13:
14: The GNU MP Library is distributed in the hope that it will be useful, but
15: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.1.1.2 maekawa 16: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 17: License for more details.
18:
1.1.1.2 maekawa 19: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 20: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
21: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
22: MA 02111-1307, USA. */
23:
1.1.1.2 maekawa 24: #include <stdio.h> /* for NULL */
1.1 maekawa 25: #include "gmp.h"
26: #include "gmp-impl.h"
27:
28: void
29: mpz_gcdext (mpz_ptr g, mpz_ptr s, mpz_ptr t, mpz_srcptr a, mpz_srcptr b)
30: {
1.1.1.2 maekawa 31: mp_size_t asize, bsize, usize, vsize;
32: mp_srcptr ap, bp;
33: mp_ptr up, vp;
34: mp_size_t gsize, ssize, tmp_ssize;
35: mp_ptr gp, sp, tmp_gp, tmp_sp;
36: mpz_srcptr u, v;
37: mpz_ptr ss, tt;
38: __mpz_struct stmp, gtmp;
39: TMP_DECL (marker);
40:
41: TMP_MARK (marker);
42:
43: /* mpn_gcdext requires that U >= V. Therefore, we often have to swap U and
44: V. This in turn leads to a lot of complications. The computed cofactor
45: will be the wrong one, so we have to fix that up at the end. */
46:
47: asize = ABS (SIZ (a));
48: bsize = ABS (SIZ (b));
49: ap = PTR (a);
50: bp = PTR (b);
51: if (asize > bsize || (asize == bsize && mpn_cmp (ap, bp, asize) > 0))
52: {
53: usize = asize;
54: vsize = bsize;
55: up = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB);
56: vp = (mp_ptr) TMP_ALLOC ((vsize + 1) * BYTES_PER_MP_LIMB);
57: MPN_COPY (up, ap, usize);
58: MPN_COPY (vp, bp, vsize);
59: u = a;
60: v = b;
61: ss = s;
62: tt = t;
63: }
64: else
65: {
66: usize = bsize;
67: vsize = asize;
68: up = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB);
69: vp = (mp_ptr) TMP_ALLOC ((vsize + 1) * BYTES_PER_MP_LIMB);
70: MPN_COPY (up, bp, usize);
71: MPN_COPY (vp, ap, vsize);
72: u = b;
73: v = a;
74: ss = t;
75: tt = s;
76: }
1.1 maekawa 77:
1.1.1.2 maekawa 78: tmp_gp = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB);
79: tmp_sp = (mp_ptr) TMP_ALLOC ((usize + 1) * BYTES_PER_MP_LIMB);
1.1 maekawa 80:
1.1.1.2 maekawa 81: if (vsize == 0)
1.1 maekawa 82: {
1.1.1.2 maekawa 83: tmp_sp[0] = 1;
84: tmp_ssize = 1;
85: MPN_COPY (tmp_gp, up, usize);
86: gsize = usize;
1.1 maekawa 87: }
1.1.1.2 maekawa 88: else
89: gsize = mpn_gcdext (tmp_gp, tmp_sp, &tmp_ssize, up, usize, vp, vsize);
90: ssize = ABS (tmp_ssize);
1.1 maekawa 91:
1.1.1.2 maekawa 92: PTR (>mp) = tmp_gp;
93: SIZ (>mp) = gsize;
94:
95: PTR (&stmp) = tmp_sp;
96: SIZ (&stmp) = (tmp_ssize ^ SIZ (u)) >= 0 ? ssize : -ssize;
97:
98: if (tt != NULL)
1.1 maekawa 99: {
1.1.1.2 maekawa 100: if (SIZ (v) == 0)
101: SIZ (tt) = 0;
1.1 maekawa 102: else
1.1.1.2 maekawa 103: {
104: mpz_t x;
105: MPZ_TMP_INIT (x, ssize + usize + 1);
106: mpz_mul (x, &stmp, u);
107: mpz_sub (x, >mp, x);
108: mpz_tdiv_q (tt, x, v);
109: }
1.1 maekawa 110: }
1.1.1.2 maekawa 111:
112: if (ss != NULL)
1.1 maekawa 113: {
1.1.1.2 maekawa 114: if (ALLOC (ss) < ssize)
115: _mpz_realloc (ss, ssize);
116: sp = PTR (ss);
117: MPN_COPY (sp, tmp_sp, ssize);
118: SIZ (ss) = SIZ (&stmp);
1.1 maekawa 119: }
120:
1.1.1.2 maekawa 121: if (ALLOC (g) < gsize)
122: _mpz_realloc (g, gsize);
123: gp = PTR (g);
124: MPN_COPY (gp, tmp_gp, gsize);
125: SIZ (g) = gsize;
126:
127: TMP_FREE (marker);
1.1 maekawa 128: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>