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

     1  /* mpz_remove -- divide out a factor and return its multiplicity.
     2  
     3  Copyright 1998-2002, 2012 Free Software Foundation, Inc.
     4  
     5  This file is part of the GNU MP Library.
     6  
     7  The GNU MP Library is free software; you can redistribute it and/or modify
     8  it under the terms of either:
     9  
    10    * the GNU Lesser General Public License as published by the Free
    11      Software Foundation; either version 3 of the License, or (at your
    12      option) any later version.
    13  
    14  or
    15  
    16    * the GNU General Public License as published by the Free Software
    17      Foundation; either version 2 of the License, or (at your option) any
    18      later version.
    19  
    20  or both in parallel, as here.
    21  
    22  The GNU MP Library is distributed in the hope that it will be useful, but
    23  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    24  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    25  for more details.
    26  
    27  You should have received copies of the GNU General Public License and the
    28  GNU Lesser General Public License along with the GNU MP Library.  If not,
    29  see https://www.gnu.org/licenses/.  */
    30  
    31  #include "gmp.h"
    32  #include "gmp-impl.h"
    33  
    34  mp_bitcnt_t
    35  mpz_remove (mpz_ptr dest, mpz_srcptr src, mpz_srcptr f)
    36  {
    37    mp_bitcnt_t pwr;
    38    mp_srcptr fp;
    39    mp_size_t sn, fn, afn;
    40    mp_limb_t fp0;
    41  
    42    sn = SIZ (src);
    43    fn = SIZ (f);
    44    fp = PTR (f);
    45    afn = ABS (fn);
    46    fp0 = fp[0];
    47  
    48    if (UNLIKELY ((afn <= (fp0 == 1)) /* mpz_cmpabs_ui (f, 1) <= 0 */
    49  		| (sn == 0)))
    50      {
    51        /*  f = 0 or f = +- 1 or src = 0 */
    52        if (afn == 0)
    53  	DIVIDE_BY_ZERO;
    54        mpz_set (dest, src);
    55        return 0;
    56      }
    57  
    58    if ((fp0 & 1) != 0)
    59      { /* f is odd */
    60        mp_ptr dp;
    61        mp_size_t dn;
    62  
    63        dn = ABS (sn);
    64        dp = MPZ_REALLOC (dest, dn);
    65  
    66        pwr = mpn_remove (dp, &dn, PTR(src), dn, PTR(f), afn, ~(mp_bitcnt_t) 0);
    67  
    68        SIZ (dest) = ((pwr & (fn < 0)) ^ (sn < 0)) ? -dn : dn;
    69      }
    70    else if (afn == (fp0 == 2))
    71      { /* mpz_cmpabs_ui (f, 2) == 0 */
    72        pwr = mpz_scan1 (src, 0);
    73        mpz_div_2exp (dest, src, pwr);
    74        if (pwr & (fn < 0)) /*((pwr % 2 == 1) && (SIZ (f) < 0))*/
    75  	mpz_neg (dest, dest);
    76      }
    77    else
    78      { /* f != +-2 */
    79        mpz_t x, rem;
    80  
    81        mpz_init (rem);
    82        mpz_init (x);
    83  
    84        pwr = 0;
    85        mpz_tdiv_qr (x, rem, src, f);
    86        if (SIZ (rem) == 0)
    87  	{
    88  	  mpz_t fpow[GMP_LIMB_BITS];		/* Really MP_SIZE_T_BITS */
    89  	  int p;
    90  
    91  #if WANT_ORIGINAL_DEST
    92  	  mp_ptr dp;
    93  	  dp = PTR (dest);
    94  #endif
    95        /* We could perhaps compute mpz_scan1(src,0)/mpz_scan1(f,0).  It is an
    96  	 upper bound of the result we're seeking.  We could also shift down the
    97  	 operands so that they become odd, to make intermediate values
    98  	 smaller.  */
    99  	  mpz_init_set (fpow[0], f);
   100  	  mpz_swap (dest, x);
   101  
   102  	  p = 1;
   103        /* Divide by f, f^2 ... f^(2^k) until we get a remainder for f^(2^k).  */
   104  	  while (ABSIZ (dest) >= 2 * ABSIZ (fpow[p - 1]) - 1)
   105  	    {
   106  	      mpz_init (fpow[p]);
   107  	      mpz_mul (fpow[p], fpow[p - 1], fpow[p - 1]);
   108  	      mpz_tdiv_qr (x, rem, dest, fpow[p]);
   109  	      if (SIZ (rem) != 0) {
   110  		mpz_clear (fpow[p]);
   111  		break;
   112  	      }
   113  	      mpz_swap (dest, x);
   114  	      p++;
   115  	    }
   116  
   117  	  pwr = ((mp_bitcnt_t)1 << p) - 1;
   118  
   119        /* Divide by f^(2^(k-1)), f^(2^(k-2)), ..., f for all divisors that give
   120  	 a zero remainder.  */
   121  	  while (--p >= 0)
   122  	    {
   123  	      mpz_tdiv_qr (x, rem, dest, fpow[p]);
   124  	      if (SIZ (rem) == 0)
   125  		{
   126  		  pwr += (mp_bitcnt_t)1 << p;
   127  		  mpz_swap (dest, x);
   128  		}
   129  	      mpz_clear (fpow[p]);
   130  	    }
   131  
   132  #if WANT_ORIGINAL_DEST
   133  	  if (PTR (x) == dp) {
   134  	    mpz_swap (dest, x);
   135  	    mpz_set (dest, x);
   136  	  }
   137  #endif
   138  	}
   139        else
   140  	mpz_set (dest, src);
   141  
   142        mpz_clear (x);
   143        mpz_clear (rem);
   144      }
   145  
   146    return pwr;
   147  }