source/libs/gmp/gmp-src/mpz/divegcd.c

   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 }