diff options
author | Pedro Alvarez <pedro.alvarez@codethink.co.uk> | 2016-05-27 17:39:31 +0100 |
---|---|---|
committer | Pedro Alvarez <pedro.alvarez@codethink.co.uk> | 2016-05-27 17:53:32 +0100 |
commit | 26c75cf8267919f81a1759c9c965a52c660233f9 (patch) | |
tree | cf2a39cf56c2c8ac45760854413ab233e6263974 /gmp/mpn/generic/gcdext_1.c | |
parent | 56892c1d217baea02092b51a09bbc924130ca84c (diff) | |
download | gcc-tarball-26c75cf8267919f81a1759c9c965a52c660233f9.tar.gz |
go to gmp 4.3.2baserock/pedroalvarez/gcc-5.3.0-gmp432
Diffstat (limited to 'gmp/mpn/generic/gcdext_1.c')
-rw-r--r-- | gmp/mpn/generic/gcdext_1.c | 265 |
1 files changed, 9 insertions, 256 deletions
diff --git a/gmp/mpn/generic/gcdext_1.c b/gmp/mpn/generic/gcdext_1.c index ea46cceb72..f1dd9ee963 100644 --- a/gmp/mpn/generic/gcdext_1.c +++ b/gmp/mpn/generic/gcdext_1.c @@ -1,273 +1,27 @@ /* mpn_gcdext -- Extended Greatest Common Divisor. -Copyright 1996, 1998, 2000-2005, 2008, 2009 Free Software Foundation, Inc. +Copyright 1996, 1998, 2000, 2001, 2002, 2003, 2004, 2005, 2008, 2009 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 either: - - * the GNU Lesser General Public License as published by the Free - Software Foundation; either version 3 of the License, or (at your - option) any later version. - -or - - * the GNU General Public License as published by the Free Software - Foundation; either version 2 of the License, or (at your option) any - later version. - -or both in parallel, as here. +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation; either version 3 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 General Public License -for more details. +or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public +License for more details. -You should have received copies of the GNU General Public License and the -GNU Lesser General Public License along with the GNU MP Library. If not, -see https://www.gnu.org/licenses/. */ +You should have received a copy of the GNU Lesser General Public License +along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" -#ifndef GCDEXT_1_USE_BINARY -#define GCDEXT_1_USE_BINARY 0 -#endif - -#ifndef GCDEXT_1_BINARY_METHOD -#define GCDEXT_1_BINARY_METHOD 2 -#endif - -#ifndef USE_ZEROTAB -#define USE_ZEROTAB 1 -#endif - -#if GCDEXT_1_USE_BINARY - -#if USE_ZEROTAB -static unsigned char zerotab[0x40] = { - 6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, - 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, - 5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, - 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0 -}; -#endif - -mp_limb_t -mpn_gcdext_1 (mp_limb_signed_t *sp, mp_limb_signed_t *tp, - mp_limb_t u, mp_limb_t v) -{ - /* Maintain - - U = t1 u + t0 v - V = s1 u + s0 v - - where U, V are the inputs (without any shared power of two), - and the matrix has determinant ± 2^{shift}. - */ - mp_limb_t s0 = 1; - mp_limb_t t0 = 0; - mp_limb_t s1 = 0; - mp_limb_t t1 = 1; - mp_limb_t ug; - mp_limb_t vg; - mp_limb_t ugh; - mp_limb_t vgh; - unsigned zero_bits; - unsigned shift; - unsigned i; -#if GCDEXT_1_BINARY_METHOD == 2 - mp_limb_t det_sign; -#endif - - ASSERT (u > 0); - ASSERT (v > 0); - - count_trailing_zeros (zero_bits, u | v); - u >>= zero_bits; - v >>= zero_bits; - - if ((u & 1) == 0) - { - count_trailing_zeros (shift, u); - u >>= shift; - t1 <<= shift; - } - else if ((v & 1) == 0) - { - count_trailing_zeros (shift, v); - v >>= shift; - s0 <<= shift; - } - else - shift = 0; - -#if GCDEXT_1_BINARY_METHOD == 1 - while (u != v) - { - unsigned count; - if (u > v) - { - u -= v; -#if USE_ZEROTAB - count = zerotab [u & 0x3f]; - u >>= count; - if (UNLIKELY (count == 6)) - { - unsigned c; - do - { - c = zerotab[u & 0x3f]; - u >>= c; - count += c; - } - while (c == 6); - } -#else - count_trailing_zeros (count, u); - u >>= count; -#endif - t0 += t1; t1 <<= count; - s0 += s1; s1 <<= count; - } - else - { - v -= u; -#if USE_ZEROTAB - count = zerotab [v & 0x3f]; - v >>= count; - if (UNLIKELY (count == 6)) - { - unsigned c; - do - { - c = zerotab[v & 0x3f]; - v >>= c; - count += c; - } - while (c == 6); - } -#else - count_trailing_zeros (count, v); - v >>= count; -#endif - t1 += t0; t0 <<= count; - s1 += s0; s0 <<= count; - } - shift += count; - } -#else -# if GCDEXT_1_BINARY_METHOD == 2 - u >>= 1; - v >>= 1; - - det_sign = 0; - - while (u != v) - { - unsigned count; - mp_limb_t d = u - v; - mp_limb_t vgtu = LIMB_HIGHBIT_TO_MASK (d); - mp_limb_t sx; - mp_limb_t tx; - - /* When v <= u (vgtu == 0), the updates are: - - (u; v) <-- ( (u - v) >> count; v) (det = +(1<<count) for corr. M factor) - (t1, t0) <-- (t1 << count, t0 + t1) - - and when v > 0, the updates are - - (u; v) <-- ( (v - u) >> count; u) (det = -(1<<count)) - (t1, t0) <-- (t0 << count, t0 + t1) - - and similarly for s1, s0 - */ - - /* v <-- min (u, v) */ - v += (vgtu & d); - - /* u <-- |u - v| */ - u = (d ^ vgtu) - vgtu; - - /* Number of trailing zeros is the same no matter if we look at - * d or u, but using d gives more parallelism. */ -#if USE_ZEROTAB - count = zerotab[d & 0x3f]; - if (UNLIKELY (count == 6)) - { - unsigned c = 6; - do - { - d >>= c; - c = zerotab[d & 0x3f]; - count += c; - } - while (c == 6); - } -#else - count_trailing_zeros (count, d); -#endif - det_sign ^= vgtu; - - tx = vgtu & (t0 - t1); - sx = vgtu & (s0 - s1); - t0 += t1; - s0 += s1; - t1 += tx; - s1 += sx; - - count++; - u >>= count; - t1 <<= count; - s1 <<= count; - shift += count; - } - u = (u << 1) + 1; -# else /* GCDEXT_1_BINARY_METHOD == 2 */ -# error Unknown GCDEXT_1_BINARY_METHOD -# endif -#endif - - /* Now u = v = g = gcd (u,v). Compute U/g and V/g */ - ug = t0 + t1; - vg = s0 + s1; - - ugh = ug/2 + (ug & 1); - vgh = vg/2 + (vg & 1); - - /* Now ±2^{shift} g = s0 U - t0 V. Get rid of the power of two, using - s0 U - t0 V = (s0 + V/g) U - (t0 + U/g) V. */ - for (i = 0; i < shift; i++) - { - mp_limb_t mask = - ( (s0 | t0) & 1); - - s0 /= 2; - t0 /= 2; - s0 += mask & vgh; - t0 += mask & ugh; - } - /* FIXME: Try simplifying this condition. */ - if ( (s0 > 1 && 2*s0 >= vg) || (t0 > 1 && 2*t0 >= ug) ) - { - s0 -= vg; - t0 -= ug; - } -#if GCDEXT_1_BINARY_METHOD == 2 - /* Conditional negation. */ - s0 = (s0 ^ det_sign) - det_sign; - t0 = (t0 ^ det_sign) - det_sign; -#endif - *sp = s0; - *tp = -t0; - - return u << zero_bits; -} - -#else /* !GCDEXT_1_USE_BINARY */ - /* FIXME: Takes two single-word limbs. It could be extended to a * function that accepts a bignum for the first input, and only @@ -325,4 +79,3 @@ mpn_gcdext_1 (mp_limb_signed_t *up, mp_limb_signed_t *vp, v1 -= q * v0; } } -#endif /* !GCDEXT_1_USE_BINARY */ |