Annotation of OpenXM_contrib/gmp/mpz/divegcd.c, Revision 1.1
1.1 ! ohara 1: /* mpz_divexact_gcd -- exact division optimized for GCDs.
! 2:
! 3: THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE AND ARE ALMOST CERTAIN TO
! 4: BE SUBJECT TO INCOMPATIBLE CHANGES IN FUTURE GNU MP RELEASES.
! 5:
! 6: Copyright 2000 Free Software Foundation, Inc.
! 7:
! 8: This file is part of the GNU MP Library.
! 9:
! 10: The GNU MP Library is free software; you can redistribute it and/or modify
! 11: it under the terms of the GNU Lesser General Public License as published by
! 12: the Free Software Foundation; either version 2.1 of the License, or (at your
! 13: option) any later version.
! 14:
! 15: The GNU MP Library is distributed in the hope that it will be useful, but
! 16: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 17: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 18: License for more details.
! 19:
! 20: You should have received a copy of the GNU Lesser General Public License
! 21: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 22: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 23: MA 02111-1307, USA. */
! 24:
! 25: #include "gmp.h"
! 26: #include "gmp-impl.h"
! 27: #include "longlong.h"
! 28:
! 29:
! 30: /* Set q to a/d, expecting d to be from a GCD and therefore usually small.
! 31:
! 32: The distribution of GCDs of random numbers can be found in Knuth volume 2
! 33: section 4.5.2 theorem D.
! 34:
! 35: GCD chance
! 36: 1 60.7%
! 37: 2^k 20.2% (1<=k<32)
! 38: 3*2^k 9.0% (1<=k<32)
! 39: other 10.1%
! 40:
! 41: Only the low limb is examined for optimizations, since GCDs bigger than
! 42: 2^32 (or 2^64) will occur very infrequently.
! 43:
! 44: Future: This could change to an mpn_divexact_gcd, possibly partly
! 45: inlined, if/when the relevant mpq functions change to an mpn based
! 46: implementation. */
! 47:
! 48:
! 49: static void
! 50: mpz_divexact_by3 (mpz_ptr q, mpz_srcptr a)
! 51: {
! 52: mp_size_t size = SIZ(a);
! 53: if (size == 0)
! 54: {
! 55: SIZ(q) = 0;
! 56: return;
! 57: }
! 58: else
! 59: {
! 60: mp_size_t abs_size = ABS(size);
! 61: mp_ptr qp;
! 62:
! 63: MPZ_REALLOC (q, abs_size);
! 64:
! 65: qp = PTR(q);
! 66: mpn_divexact_by3 (qp, PTR(a), abs_size);
! 67:
! 68: abs_size -= (qp[abs_size-1] == 0);
! 69: SIZ(q) = (size>0 ? abs_size : -abs_size);
! 70: }
! 71: }
! 72:
! 73: void
! 74: mpz_divexact_gcd (mpz_ptr q, mpz_srcptr a, mpz_srcptr d)
! 75: {
! 76: ASSERT (mpz_sgn (d) > 0);
! 77:
! 78: if (SIZ(d) == 1)
! 79: {
! 80: mp_limb_t dl = PTR(d)[0];
! 81: int twos;
! 82:
! 83: if (dl == 1)
! 84: {
! 85: if (q != a)
! 86: mpz_set (q, a);
! 87: return;
! 88: }
! 89: if (dl == 3)
! 90: {
! 91: mpz_divexact_by3 (q, a);
! 92: return;
! 93: }
! 94:
! 95: count_trailing_zeros (twos, dl);
! 96: dl >>= twos;
! 97:
! 98: if (dl == 1)
! 99: {
! 100: mpz_tdiv_q_2exp (q, a, twos);
! 101: return;
! 102: }
! 103: if (dl == 3)
! 104: {
! 105: mpz_tdiv_q_2exp (q, a, twos);
! 106: mpz_divexact_by3 (q, q);
! 107: return;
! 108: }
! 109: }
! 110:
! 111: mpz_divexact (q, a, d);
! 112: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>