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 }