aa9bc17601
BMakefiles and other bits will follow. Requested by: Andrey Chernov Made world by: Chuck Robey
442 lines
9.0 KiB
C
442 lines
9.0 KiB
C
/* mpn_gcdext -- Extended Greatest Common Divisor.
|
|
|
|
Copyright (C) 1996 Free Software Foundation, Inc.
|
|
|
|
This file is part of the GNU MP Library.
|
|
|
|
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
|
|
the Free Software Foundation; either version 2 of the License, or (at your
|
|
option) any later version.
|
|
|
|
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
|
|
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
|
|
License for more details.
|
|
|
|
You should have received a copy of the GNU Library General Public License
|
|
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,
|
|
MA 02111-1307, USA. */
|
|
|
|
#include "gmp.h"
|
|
#include "gmp-impl.h"
|
|
#include "longlong.h"
|
|
|
|
#ifndef EXTEND
|
|
#define EXTEND 1
|
|
#endif
|
|
|
|
#if STAT
|
|
int arr[BITS_PER_MP_LIMB];
|
|
#endif
|
|
|
|
#define SGN(A) (((A) < 0) ? -1 : ((A) > 0))
|
|
|
|
/* Idea 1: After we have performed a full division, don't shift operands back,
|
|
but instead account for the extra factors-of-2 thus introduced.
|
|
Idea 2: Simple generalization to use divide-and-conquer would give us an
|
|
algorithm that runs faster than O(n^2).
|
|
Idea 3: The input numbers need less space as the computation progresses,
|
|
while the s0 and s1 variables need more space. To save space, we
|
|
could make them share space, and have the latter variables grow
|
|
into the former. */
|
|
|
|
/* Precondition: U >= V. */
|
|
|
|
mp_size_t
|
|
#if EXTEND
|
|
#if __STDC__
|
|
mpn_gcdext (mp_ptr gp, mp_ptr s0p,
|
|
mp_ptr up, mp_size_t size, mp_ptr vp, mp_size_t vsize)
|
|
#else
|
|
mpn_gcdext (gp, s0p, up, size, vp, vsize)
|
|
mp_ptr gp;
|
|
mp_ptr s0p;
|
|
mp_ptr up;
|
|
mp_size_t size;
|
|
mp_ptr vp;
|
|
mp_size_t vsize;
|
|
#endif
|
|
#else
|
|
#if __STDC__
|
|
mpn_gcd (mp_ptr gp,
|
|
mp_ptr up, mp_size_t size, mp_ptr vp, mp_size_t vsize)
|
|
#else
|
|
mpn_gcd (gp, up, size, vp, vsize)
|
|
mp_ptr gp;
|
|
mp_ptr up;
|
|
mp_size_t size;
|
|
mp_ptr vp;
|
|
mp_size_t vsize;
|
|
#endif
|
|
#endif
|
|
{
|
|
mp_limb_t uh, vh;
|
|
mp_limb_signed_t A, B, C, D;
|
|
int cnt;
|
|
mp_ptr tp, wp;
|
|
#if RECORD
|
|
mp_limb_signed_t min = 0, max = 0;
|
|
#endif
|
|
#if EXTEND
|
|
mp_ptr s1p;
|
|
mp_ptr orig_s0p = s0p;
|
|
mp_size_t ssize, orig_size = size;
|
|
TMP_DECL (mark);
|
|
|
|
TMP_MARK (mark);
|
|
|
|
tp = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
|
|
wp = (mp_ptr) TMP_ALLOC ((size + 1) * BYTES_PER_MP_LIMB);
|
|
s1p = (mp_ptr) TMP_ALLOC (size * BYTES_PER_MP_LIMB);
|
|
|
|
MPN_ZERO (s0p, size);
|
|
MPN_ZERO (s1p, size);
|
|
|
|
s0p[0] = 1;
|
|
s1p[0] = 0;
|
|
ssize = 1;
|
|
#endif
|
|
|
|
if (size > vsize)
|
|
{
|
|
/* Normalize V (and shift up U the same amount). */
|
|
count_leading_zeros (cnt, vp[vsize - 1]);
|
|
if (cnt != 0)
|
|
{
|
|
mp_limb_t cy;
|
|
mpn_lshift (vp, vp, vsize, cnt);
|
|
cy = mpn_lshift (up, up, size, cnt);
|
|
up[size] = cy;
|
|
size += cy != 0;
|
|
}
|
|
|
|
mpn_divmod (up + vsize, up, size, vp, vsize);
|
|
#if EXTEND
|
|
/* This is really what it boils down to in this case... */
|
|
s0p[0] = 0;
|
|
s1p[0] = 1;
|
|
#endif
|
|
size = vsize;
|
|
if (cnt != 0)
|
|
{
|
|
mpn_rshift (up, up, size, cnt);
|
|
mpn_rshift (vp, vp, size, cnt);
|
|
}
|
|
{
|
|
mp_ptr xp;
|
|
xp = up; up = vp; vp = xp;
|
|
}
|
|
}
|
|
|
|
for (;;)
|
|
{
|
|
/* Figure out exact size of V. */
|
|
vsize = size;
|
|
MPN_NORMALIZE (vp, vsize);
|
|
if (vsize <= 1)
|
|
break;
|
|
|
|
/* Make UH be the most significant limb of U, and make VH be
|
|
corresponding bits from V. */
|
|
uh = up[size - 1];
|
|
vh = vp[size - 1];
|
|
count_leading_zeros (cnt, uh);
|
|
if (cnt != 0)
|
|
{
|
|
uh = (uh << cnt) | (up[size - 2] >> (BITS_PER_MP_LIMB - cnt));
|
|
vh = (vh << cnt) | (vp[size - 2] >> (BITS_PER_MP_LIMB - cnt));
|
|
}
|
|
|
|
#if 0
|
|
/* For now, only handle BITS_PER_MP_LIMB-1 bits. This makes
|
|
room for sign bit. */
|
|
uh >>= 1;
|
|
vh >>= 1;
|
|
#endif
|
|
A = 1;
|
|
B = 0;
|
|
C = 0;
|
|
D = 1;
|
|
|
|
for (;;)
|
|
{
|
|
mp_limb_signed_t q, T;
|
|
if (vh + C == 0 || vh + D == 0)
|
|
break;
|
|
|
|
q = (uh + A) / (vh + C);
|
|
if (q != (uh + B) / (vh + D))
|
|
break;
|
|
|
|
T = A - q * C;
|
|
A = C;
|
|
C = T;
|
|
T = B - q * D;
|
|
B = D;
|
|
D = T;
|
|
T = uh - q * vh;
|
|
uh = vh;
|
|
vh = T;
|
|
}
|
|
|
|
#if RECORD
|
|
min = MIN (A, min); min = MIN (B, min);
|
|
min = MIN (C, min); min = MIN (D, min);
|
|
max = MAX (A, max); max = MAX (B, max);
|
|
max = MAX (C, max); max = MAX (D, max);
|
|
#endif
|
|
|
|
if (B == 0)
|
|
{
|
|
mp_limb_t qh;
|
|
mp_size_t i;
|
|
|
|
/* This is quite rare. I.e., optimize something else! */
|
|
|
|
/* Normalize V (and shift up U the same amount). */
|
|
count_leading_zeros (cnt, vp[vsize - 1]);
|
|
if (cnt != 0)
|
|
{
|
|
mp_limb_t cy;
|
|
mpn_lshift (vp, vp, vsize, cnt);
|
|
cy = mpn_lshift (up, up, size, cnt);
|
|
up[size] = cy;
|
|
size += cy != 0;
|
|
}
|
|
|
|
qh = mpn_divmod (up + vsize, up, size, vp, vsize);
|
|
#if EXTEND
|
|
MPN_COPY (tp, s0p, ssize);
|
|
for (i = 0; i < size - vsize; i++)
|
|
{
|
|
mp_limb_t cy;
|
|
cy = mpn_addmul_1 (tp + i, s1p, ssize, up[vsize + i]);
|
|
if (cy != 0)
|
|
tp[ssize++] = cy;
|
|
}
|
|
if (qh != 0)
|
|
{
|
|
mp_limb_t cy;
|
|
abort ();
|
|
/* XXX since qh == 1, mpn_addmul_1 is overkill */
|
|
cy = mpn_addmul_1 (tp + size - vsize, s1p, ssize, qh);
|
|
if (cy != 0)
|
|
tp[ssize++] = cy;
|
|
}
|
|
#if 0
|
|
MPN_COPY (s0p, s1p, ssize); /* should be old ssize, kind of */
|
|
MPN_COPY (s1p, tp, ssize);
|
|
#else
|
|
{
|
|
mp_ptr xp;
|
|
xp = s0p; s0p = s1p; s1p = xp;
|
|
xp = s1p; s1p = tp; tp = xp;
|
|
}
|
|
#endif
|
|
#endif
|
|
size = vsize;
|
|
if (cnt != 0)
|
|
{
|
|
mpn_rshift (up, up, size, cnt);
|
|
mpn_rshift (vp, vp, size, cnt);
|
|
}
|
|
|
|
{
|
|
mp_ptr xp;
|
|
xp = up; up = vp; vp = xp;
|
|
}
|
|
MPN_NORMALIZE (up, size);
|
|
}
|
|
else
|
|
{
|
|
/* T = U*A + V*B
|
|
W = U*C + V*D
|
|
U = T
|
|
V = W */
|
|
|
|
if (SGN(A) == SGN(B)) /* should be different sign */
|
|
abort ();
|
|
if (SGN(C) == SGN(D)) /* should be different sign */
|
|
abort ();
|
|
#if STAT
|
|
{ mp_limb_t x;
|
|
x = ABS (A) | ABS (B) | ABS (C) | ABS (D);
|
|
count_leading_zeros (cnt, x);
|
|
arr[BITS_PER_MP_LIMB - cnt]++; }
|
|
#endif
|
|
if (A == 0)
|
|
{
|
|
if (B != 1) abort ();
|
|
MPN_COPY (tp, vp, size);
|
|
}
|
|
else
|
|
{
|
|
if (A < 0)
|
|
{
|
|
mpn_mul_1 (tp, vp, size, B);
|
|
mpn_submul_1 (tp, up, size, -A);
|
|
}
|
|
else
|
|
{
|
|
mpn_mul_1 (tp, up, size, A);
|
|
mpn_submul_1 (tp, vp, size, -B);
|
|
}
|
|
}
|
|
if (C < 0)
|
|
{
|
|
mpn_mul_1 (wp, vp, size, D);
|
|
mpn_submul_1 (wp, up, size, -C);
|
|
}
|
|
else
|
|
{
|
|
mpn_mul_1 (wp, up, size, C);
|
|
mpn_submul_1 (wp, vp, size, -D);
|
|
}
|
|
|
|
{
|
|
mp_ptr xp;
|
|
xp = tp; tp = up; up = xp;
|
|
xp = wp; wp = vp; vp = xp;
|
|
}
|
|
|
|
#if EXTEND
|
|
{ mp_limb_t cy;
|
|
MPN_ZERO (tp, orig_size);
|
|
if (A == 0)
|
|
{
|
|
if (B != 1) abort ();
|
|
MPN_COPY (tp, s1p, ssize);
|
|
}
|
|
else
|
|
{
|
|
if (A < 0)
|
|
{
|
|
cy = mpn_mul_1 (tp, s1p, ssize, B);
|
|
cy += mpn_addmul_1 (tp, s0p, ssize, -A);
|
|
}
|
|
else
|
|
{
|
|
cy = mpn_mul_1 (tp, s0p, ssize, A);
|
|
cy += mpn_addmul_1 (tp, s1p, ssize, -B);
|
|
}
|
|
if (cy != 0)
|
|
tp[ssize++] = cy;
|
|
}
|
|
MPN_ZERO (wp, orig_size);
|
|
if (C < 0)
|
|
{
|
|
cy = mpn_mul_1 (wp, s1p, ssize, D);
|
|
cy += mpn_addmul_1 (wp, s0p, ssize, -C);
|
|
}
|
|
else
|
|
{
|
|
cy = mpn_mul_1 (wp, s0p, ssize, C);
|
|
cy += mpn_addmul_1 (wp, s1p, ssize, -D);
|
|
}
|
|
if (cy != 0)
|
|
wp[ssize++] = cy;
|
|
}
|
|
{
|
|
mp_ptr xp;
|
|
xp = tp; tp = s0p; s0p = xp;
|
|
xp = wp; wp = s1p; s1p = xp;
|
|
}
|
|
#endif
|
|
#if 0 /* Is it a win to remove multiple zeros here? */
|
|
MPN_NORMALIZE (up, size);
|
|
#else
|
|
if (up[size - 1] == 0)
|
|
size--;
|
|
#endif
|
|
}
|
|
}
|
|
|
|
#if RECORD
|
|
printf ("min: %ld\n", min);
|
|
printf ("max: %ld\n", max);
|
|
#endif
|
|
|
|
if (vsize == 0)
|
|
{
|
|
if (gp != up)
|
|
MPN_COPY (gp, up, size);
|
|
#if EXTEND
|
|
if (orig_s0p != s0p)
|
|
MPN_COPY (orig_s0p, s0p, ssize);
|
|
#endif
|
|
TMP_FREE (mark);
|
|
return size;
|
|
}
|
|
else
|
|
{
|
|
mp_limb_t vl, ul, t;
|
|
#if EXTEND
|
|
mp_limb_t cy;
|
|
mp_size_t i;
|
|
#endif
|
|
vl = vp[0];
|
|
#if EXTEND
|
|
t = mpn_divmod_1 (wp, up, size, vl);
|
|
MPN_COPY (tp, s0p, ssize);
|
|
for (i = 0; i < size; i++)
|
|
{
|
|
cy = mpn_addmul_1 (tp + i, s1p, ssize, wp[i]);
|
|
if (cy != 0)
|
|
tp[ssize++] = cy;
|
|
}
|
|
#if 0
|
|
MPN_COPY (s0p, s1p, ssize);
|
|
MPN_COPY (s1p, tp, ssize);
|
|
#else
|
|
{
|
|
mp_ptr xp;
|
|
xp = s0p; s0p = s1p; s1p = xp;
|
|
xp = s1p; s1p = tp; tp = xp;
|
|
}
|
|
#endif
|
|
#else
|
|
t = mpn_mod_1 (up, size, vl);
|
|
#endif
|
|
ul = vl;
|
|
vl = t;
|
|
while (vl != 0)
|
|
{
|
|
mp_limb_t t;
|
|
#if EXTEND
|
|
mp_limb_t q, cy;
|
|
q = ul / vl;
|
|
t = ul - q*vl;
|
|
|
|
MPN_COPY (tp, s0p, ssize);
|
|
cy = mpn_addmul_1 (tp, s1p, ssize, q);
|
|
if (cy != 0)
|
|
tp[ssize++] = cy;
|
|
#if 0
|
|
MPN_COPY (s0p, s1p, ssize);
|
|
MPN_COPY (s1p, tp, ssize);
|
|
#else
|
|
{
|
|
mp_ptr xp;
|
|
xp = s0p; s0p = s1p; s1p = xp;
|
|
xp = s1p; s1p = tp; tp = xp;
|
|
}
|
|
#endif
|
|
|
|
#else
|
|
t = ul % vl;
|
|
#endif
|
|
ul = vl;
|
|
vl = t;
|
|
}
|
|
gp[0] = ul;
|
|
#if EXTEND
|
|
if (orig_s0p != s0p)
|
|
MPN_COPY (orig_s0p, s0p, ssize);
|
|
#endif
|
|
TMP_FREE (mark);
|
|
return 1;
|
|
}
|
|
}
|