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