X-Git-Url: https://oss.titaniummirror.com/gitweb?a=blobdiff_plain;f=gmp%2Fmpn%2Fgeneric%2Fgcdext_1.c;fp=gmp%2Fmpn%2Fgeneric%2Fgcdext_1.c;h=efade2b4c9e00542ad751b948d2c4a2e4f89f1cb;hb=6fed43773c9b0ce596dca5686f37ac3fc0fa11c0;hp=0000000000000000000000000000000000000000;hpb=27b11d56b743098deb193d510b337ba22dc52e5c;p=msp430-gcc.git diff --git a/gmp/mpn/generic/gcdext_1.c b/gmp/mpn/generic/gcdext_1.c new file mode 100644 index 00000000..efade2b4 --- /dev/null +++ b/gmp/mpn/generic/gcdext_1.c @@ -0,0 +1,319 @@ +/* mpn_gcdext -- Extended Greatest Common Divisor. + +Copyright 1996, 1998, 2000, 2001, 2002, 2003, 2004, 2005, 2008 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 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 Lesser General Public +License for more details. + +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/. */ + +/* Default to binary gcdext_1, since it is best on most current machines. + We should teach tuneup to choose the right gcdext_1. */ +#define GCDEXT_1_USE_BINARY 1 + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" + +#ifndef NULL +# define NULL ((void *) 0) +#endif + +/* FIXME: Takes two single-word limbs. It could be extended to a + * function that accepts a bignum for the first input, and only + * returns the first co-factor. */ + +/* Returns g, u and v such that g = u A - v B. There are three + different cases for the result: + + g = u A - v B, 0 < u < b, 0 < v < a + g = A u = 1, v = 0 + g = B u = B, v = A - 1 + + We always return with 0 < u <= b, 0 <= v < a. +*/ +#if GCDEXT_1_USE_BINARY + +static mp_limb_t +gcdext_1_odd (mp_limb_t *up, mp_limb_t *vp, mp_limb_t a, mp_limb_t b) +{ + mp_limb_t u0; + mp_limb_t v0; + mp_limb_t v1; + mp_limb_t u1; + + mp_limb_t B = b; + mp_limb_t A = a; + + /* Through out this function maintain + + a = u0 A - v0 B + b = u1 A - v1 B + + where A and B are odd. */ + + u0 = 1; v0 = 0; + u1 = b; v1 = a-1; + + if (A == 1) + { + *up = u0; *vp = v0; + return 1; + } + else if (B == 1) + { + *up = u1; *vp = v1; + return 1; + } + + while (a != b) + { + mp_limb_t mask; + + ASSERT (a % 2 == 1); + ASSERT (b % 2 == 1); + + ASSERT (0 < u0); ASSERT (u0 <= B); + ASSERT (0 < u1); ASSERT (u1 <= B); + + ASSERT (0 <= v0); ASSERT (v0 < A); + ASSERT (0 <= v1); ASSERT (v1 < A); + + if (a > b) + { + MP_LIMB_T_SWAP (a, b); + MP_LIMB_T_SWAP (u0, u1); + MP_LIMB_T_SWAP (v0, v1); + } + + ASSERT (a < b); + + /* Makes b even */ + b -= a; + + mask = - (mp_limb_t) (u1 < u0); + u1 += B & mask; + v1 += A & mask; + u1 -= u0; + v1 -= v0; + + ASSERT (b % 2 == 0); + + do + { + /* As b = u1 A + v1 B is even, while A and B are odd, + either both or none of u1, v1 is even */ + + ASSERT (u1 % 2 == v1 % 2); + + mask = -(u1 & 1); + u1 = u1 / 2 + ((B / 2) & mask) - mask; + v1 = v1 / 2 + ((A / 2) & mask) - mask; + + b /= 2; + } + while (b % 2 == 0); + } + + /* Now g = a = b */ + ASSERT (a == b); + ASSERT (u1 <= B); + ASSERT (v1 < A); + + ASSERT (A % a == 0); + ASSERT (B % a == 0); + ASSERT (u0 % (B/a) == u1 % (B/a)); + ASSERT (v0 % (A/a) == v1 % (A/a)); + + *up = u0; *vp = v0; + + return a; +} + +mp_limb_t +mpn_gcdext_1 (mp_limb_t *up, mp_limb_t *vp, mp_limb_t a, mp_limb_t b) +{ + unsigned shift = 0; + mp_limb_t g; + mp_limb_t u; + mp_limb_t v; + + /* We use unsigned values in the range 0, ... B - 1. As the values + are uniquely determined only modulo B, we can add B at will, to + get numbers in range or flip the least significant bit. */ + /* Deal with powers of two */ + while ((a | b) % 2 == 0) + { + a /= 2; b /= 2; shift++; + } + + if (b % 2 == 0) + { + unsigned k = 0; + + do { + b /= 2; k++; + } while (b % 2 == 0); + + g = gcdext_1_odd (&u, &v, a, b); + + while (k--) + { + /* We have g = u a + v b, and need to construct + g = u'a + v'(2b). + + If v is even, we can just set u' = u, v' = v/2 + If v is odd, we can set v' = (v + a)/2, u' = u + b + */ + + if (v % 2 == 0) + v /= 2; + else + { + u = u + b; + v = v/2 + a/2 + 1; + } + b *= 2; + } + } + else if (a % 2 == 0) + { + unsigned k = 0; + + do { + a /= 2; k++; + } while (a % 2 == 0); + + g = gcdext_1_odd (&u, &v, a, b); + + while (k--) + { + /* We have g = u a + v b, and need to construct + g = u'(2a) + v'b. + + If u is even, we can just set u' = u/2, v' = v. + If u is odd, we can set u' = (u + b)/2 + */ + + if (u % 2 == 0) + u /= 2; + else + { + u = u/2 + b/2 + 1; + v = v + a; + } + a *= 2; + } + } + else + /* Ok, both are odd */ + g = gcdext_1_odd (&u, &v, a, b); + + *up = u; + *vp = v; + + return g << shift; +} + +#else /* ! GCDEXT_1_USE_BINARY */ +static mp_limb_t +gcdext_1_u (mp_limb_t *up, mp_limb_t a, mp_limb_t b) +{ + /* Maintain + + a = u0 A mod B + b = - u1 A mod B + */ + mp_limb_t u0 = 1; + mp_limb_t u1 = 0; + mp_limb_t B = b; + + ASSERT (a >= b); + ASSERT (b > 0); + + for (;;) + { + mp_limb_t q; + + q = a / b; + a -= q * b; + + if (a == 0) + { + *up = B - u1; + return b; + } + u0 += q * u1; + + q = b / a; + b -= q * a; + + if (b == 0) + { + *up = u0; + return a; + } + u1 += q * u0; + } +} + +mp_limb_t +mpn_gcdext_1 (mp_limb_t *up, mp_limb_t *vp, mp_limb_t a, mp_limb_t b) +{ + /* Maintain + + a = u0 A - v0 B + b = - u1 A + v1 B = (B - u1) A - (A - v1) B + */ + mp_limb_t u0 = 1; + mp_limb_t v0 = 0; + mp_limb_t u1 = 0; + mp_limb_t v1 = 1; + + mp_limb_t A = a; + mp_limb_t B = b; + + ASSERT (a >= b); + ASSERT (b > 0); + + for (;;) + { + mp_limb_t q; + + q = a / b; + a -= q * b; + + if (a == 0) + { + *up = B - u1; + *vp = A - v1; + return b; + } + u0 += q * u1; + v0 += q * v1; + + q = b / a; + b -= q * a; + + if (b == 0) + { + *up = u0; + *vp = v0; + return a; + } + u1 += q * u0; + v1 += q * v0; + } +} +#endif /* ! GCDEXT_1_USE_BINARY */