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