github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/mpz/divegcd.c (about)

     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  }