Home | History | Annotate | Line # | Download | only in mpz 
      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-impl.h"
     32 
     33 mp_bitcnt_t
     34 mpz_remove (mpz_ptr dest, mpz_srcptr src, mpz_srcptr f)
     35 {
     36   mp_bitcnt_t pwr;
     37   mp_srcptr fp;
     38   mp_size_t sn, fn, afn;
     39   mp_limb_t fp0;
     40 
     41   sn = SIZ (src);
     42   fn = SIZ (f);
     43   fp = PTR (f);
     44   afn = ABS (fn);
     45   fp0 = fp[0];
     46 
     47   if (UNLIKELY ((afn <= (fp0 == 1)) /* mpz_cmpabs_ui (f, 1) <= 0 */
     48 		| (sn == 0)))
     49     {
     50       /*  f = 0 or f = +- 1 or src = 0 */
     51       if (afn == 0)
     52 	DIVIDE_BY_ZERO;
     53       mpz_set (dest, src);
     54       return 0;
     55     }
     56 
     57   if ((fp0 & 1) != 0)
     58     { /* f is odd */
     59       mp_ptr dp;
     60       mp_size_t dn;
     61 
     62       dn = ABS (sn);
     63       dp = MPZ_REALLOC (dest, dn);
     64 
     65       pwr = mpn_remove (dp, &dn, PTR(src), dn, PTR(f), afn, ~(mp_bitcnt_t) 0);
     66 
     67       SIZ (dest) = ((pwr & (fn < 0)) ^ (sn < 0)) ? -dn : dn;
     68     }
     69   else if (afn == (fp0 == 2))
     70     { /* mpz_cmpabs_ui (f, 2) == 0 */
     71       pwr = mpz_scan1 (src, 0);
     72       mpz_div_2exp (dest, src, pwr);
     73       if (pwr & (fn < 0)) /*((pwr % 2 == 1) && (SIZ (f) < 0))*/
     74 	mpz_neg (dest, dest);
     75     }
     76   else
     77     { /* f != +-2 */
     78       mpz_t x, rem;
     79 
     80       mpz_init (rem);
     81       mpz_init (x);
     82 
     83       pwr = 0;
     84       mpz_tdiv_qr (x, rem, src, f);
     85       if (SIZ (rem) == 0)
     86 	{
     87 	  mpz_t fpow[GMP_LIMB_BITS];		/* Really MP_SIZE_T_BITS */
     88 	  int p;
     89 
     90 #if WANT_ORIGINAL_DEST
     91 	  mp_ptr dp;
     92 	  dp = PTR (dest);
     93 #endif
     94       /* We could perhaps compute mpz_scan1(src,0)/mpz_scan1(f,0).  It is an
     95 	 upper bound of the result we're seeking.  We could also shift down the
     96 	 operands so that they become odd, to make intermediate values
     97 	 smaller.  */
     98 	  mpz_init_set (fpow[0], f);
     99 	  mpz_swap (dest, x);
    100 
    101 	  p = 1;
    102       /* Divide by f, f^2 ... f^(2^k) until we get a remainder for f^(2^k).  */
    103 	  while (ABSIZ (dest) >= 2 * ABSIZ (fpow[p - 1]) - 1)
    104 	    {
    105 	      mpz_init (fpow[p]);
    106 	      mpz_mul (fpow[p], fpow[p - 1], fpow[p - 1]);
    107 	      mpz_tdiv_qr (x, rem, dest, fpow[p]);
    108 	      if (SIZ (rem) != 0) {
    109 		mpz_clear (fpow[p]);
    110 		break;
    111 	      }
    112 	      mpz_swap (dest, x);
    113 	      p++;
    114 	    }
    115 
    116 	  pwr = ((mp_bitcnt_t)1 << p) - 1;
    117 
    118       /* Divide by f^(2^(k-1)), f^(2^(k-2)), ..., f for all divisors that give
    119 	 a zero remainder.  */
    120 	  while (--p >= 0)
    121 	    {
    122 	      mpz_tdiv_qr (x, rem, dest, fpow[p]);
    123 	      if (SIZ (rem) == 0)
    124 		{
    125 		  pwr += (mp_bitcnt_t)1 << p;
    126 		  mpz_swap (dest, x);
    127 		}
    128 	      mpz_clear (fpow[p]);
    129 	    }
    130 
    131 #if WANT_ORIGINAL_DEST
    132 	  if (PTR (x) == dp) {
    133 	    mpz_swap (dest, x);
    134 	    mpz_set (dest, x);
    135 	  }
    136 #endif
    137 	}
    138       else
    139 	mpz_set (dest, src);
    140 
    141       mpz_clear (x);
    142       mpz_clear (rem);
    143     }
    144 
    145   return pwr;
    146 }
    147