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 6Copyright 2000, 2005, 2011, 2012 Free Software Foundation, Inc. 7 8This file is part of the GNU MP Library. 9 10The GNU MP Library is free software; you can redistribute it and/or modify 11it under the terms of either: 12 13 * the GNU Lesser General Public License as published by the Free 14 Software Foundation; either version 3 of the License, or (at your 15 option) any later version. 16 17or 18 19 * the GNU General Public License as published by the Free Software 20 Foundation; either version 2 of the License, or (at your option) any 21 later version. 22 23or both in parallel, as here. 24 25The GNU MP Library is distributed in the hope that it will be useful, but 26WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 27or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 28for more details. 29 30You should have received copies of the GNU General Public License and the 31GNU Lesser General Public License along with the GNU MP Library. If not, 32see https://www.gnu.org/licenses/. */ 33 34#include "gmp-impl.h" 35#include "longlong.h" 36 37 38/* Set q to a/d, expecting d to be from a GCD and therefore usually small. 39 40 The distribution of GCDs of random numbers can be found in Knuth volume 2 41 section 4.5.2 theorem D. 42 43 GCD chance 44 1 60.8% 45 2^k 20.2% (1<=k<32) 46 3*2^k 9.0% (1<=k<32) 47 other 10.1% 48 49 Only the low limb is examined for optimizations, since GCDs bigger than 50 2^32 (or 2^64) will occur very infrequently. 51 52 Future: This could change to an mpn_divexact_gcd, possibly partly 53 inlined, if/when the relevant mpq functions change to an mpn based 54 implementation. */ 55 56 57#if GMP_NUMB_BITS % 2 == 0 58static void 59mpz_divexact_by3 (mpz_ptr q, mpz_srcptr a) 60{ 61 mp_size_t size = SIZ(a); 62 mp_size_t abs_size = ABS(size); 63 mp_ptr qp; 64 65 qp = MPZ_REALLOC (q, abs_size); 66 67 mpn_bdiv_dbm1 (qp, PTR(a), abs_size, GMP_NUMB_MASK / 3); 68 69 abs_size -= (qp[abs_size-1] == 0); 70 SIZ(q) = (size>0 ? abs_size : -abs_size); 71} 72#endif 73 74#if GMP_NUMB_BITS % 4 == 0 75static void 76mpz_divexact_by5 (mpz_ptr q, mpz_srcptr a) 77{ 78 mp_size_t size = SIZ(a); 79 mp_size_t abs_size = ABS(size); 80 mp_ptr qp; 81 82 qp = MPZ_REALLOC (q, abs_size); 83 84 mpn_bdiv_dbm1 (qp, PTR(a), abs_size, GMP_NUMB_MASK / 5); 85 86 abs_size -= (qp[abs_size-1] == 0); 87 SIZ(q) = (size>0 ? abs_size : -abs_size); 88} 89#endif 90 91static void 92mpz_divexact_limb (mpz_ptr q, mpz_srcptr a, mp_limb_t d) 93{ 94 mp_size_t size = SIZ(a); 95 mp_size_t abs_size = ABS(size); 96 mp_ptr qp; 97 98 qp = MPZ_REALLOC (q, abs_size); 99 100 MPN_DIVREM_OR_DIVEXACT_1 (qp, PTR(a), abs_size, d); 101 102 abs_size -= (qp[abs_size-1] == 0); 103 SIZ(q) = (size > 0 ? abs_size : -abs_size); 104} 105 106void 107mpz_divexact_gcd (mpz_ptr q, mpz_srcptr a, mpz_srcptr d) 108{ 109 ASSERT (mpz_sgn (d) > 0); 110 111 if (SIZ(a) == 0) 112 { 113 SIZ(q) = 0; 114 return; 115 } 116 117 if (SIZ(d) == 1) 118 { 119 mp_limb_t dl = PTR(d)[0]; 120 int twos; 121 122 if ((dl & 1) == 0) 123 { 124 count_trailing_zeros (twos, dl); 125 dl >>= twos; 126 mpz_tdiv_q_2exp (q, a, twos); 127 a = q; 128 } 129 130 if (dl == 1) 131 { 132 if (q != a) 133 mpz_set (q, a); 134 return; 135 } 136#if GMP_NUMB_BITS % 2 == 0 137 if (dl == 3) 138 { 139 mpz_divexact_by3 (q, a); 140 return; 141 } 142#endif 143#if GMP_NUMB_BITS % 4 == 0 144 if (dl == 5) 145 { 146 mpz_divexact_by5 (q, a); 147 return; 148 } 149#endif 150 151 mpz_divexact_limb (q, a, dl); 152 return; 153 } 154 155 mpz_divexact (q, a, d); 156} 157