1 /* mpz_remove -- divide out a factor and return its multiplicity.
3 Copyright 1998, 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
7 The GNU MP Library is free software; you can redistribute it and/or modify
8 it under the terms of the GNU Lesser General Public License as published by
9 the Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
12 The GNU MP Library is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 License for more details.
17 You should have received a copy of the GNU Lesser General Public License
18 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
24 mpz_remove (mpz_ptr dest
, mpz_srcptr src
, mpz_srcptr f
)
26 mpz_t fpow
[GMP_LIMB_BITS
]; /* Really MP_SIZE_T_BITS */
31 if (mpz_cmp_ui (f
, 1) <= 0)
41 if (mpz_cmp_ui (f
, 2) == 0)
44 s0
= mpz_scan1 (src
, 0);
45 mpz_div_2exp (dest
, src
, s0
);
49 /* We could perhaps compute mpz_scan1(src,0)/mpz_scan1(f,0). It is an
50 upper bound of the result we're seeking. We could also shift down the
51 operands so that they become odd, to make intermediate values smaller. */
61 /* Divide by f, f^2, ..., f^(2^k) until we get a remainder for f^(2^k). */
64 mpz_tdiv_qr (x
, rem
, dest
, fpow
[p
]);
67 mpz_init (fpow
[p
+ 1]);
68 mpz_mul (fpow
[p
+ 1], fpow
[p
], fpow
[p
]);
76 /* Divide by f^(2^(k-1)), f^(2^(k-2)), ..., f for all divisors that give a
80 mpz_tdiv_qr (x
, rem
, dest
, fpow
[p
]);