1 /* mpn_remove -- divide out all multiples of odd mpn number from another mpn
4 Contributed to the GNU project by Torbjorn Granlund.
6 THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY
7 SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
8 GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
10 Copyright 2009 Free Software Foundation, Inc.
12 This file is part of the GNU MP Library.
14 The GNU MP Library is free software; you can redistribute it and/or modify
15 it under the terms of the GNU Lesser General Public License as published by
16 the Free Software Foundation; either version 3 of the License, or (at your
17 option) any later version.
19 The GNU MP Library is distributed in the hope that it will be useful, but
20 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
21 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
22 License for more details.
24 You should have received a copy of the GNU Lesser General Public License
25 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
30 #if GMP_LIMB_BITS > 50
33 #define LOG GMP_LIMB_BITS
37 /* Input: U = {up,un}, V = {vp,vn} must be odd, cap
38 Ouput W = {wp,*wn} allocation need is exactly *wn
40 Set W = U / V^k, where k is the largest integer <= cap such that the
41 division yields an integer.
43 FIXME: We currently allow any operand overlap. This is quite non mpn-ish
44 and might be changed, since it cost significant temporary space.
45 * If we require W to have space for un limbs, we could save qp or qp2 (but
46 we will still need to copy things into wp 50% of the time).
47 * If we allow ourselves to clobber U, we could save the other of qp and qp2.
51 mpn_remove (mp_ptr wp, mp_size_t *wn,
52 mp_ptr up, mp_size_t un, mp_ptr vp, mp_size_t vn,
58 mp_ptr tp, qp, np, pp, qp2, scratch_out;
59 mp_size_t pn, nn, qn, i;
65 ASSERT (vp[0] % 2 != 0); /* 2-adic division wants odd numbers */
66 ASSERT (vn > 1 || vp[0] > 1); /* else we would loop indefinitely */
70 tp = TMP_ALLOC_LIMBS ((un + vn) / 2); /* remainder */
71 qp = TMP_ALLOC_LIMBS (un); /* quotient, alternating */
72 qp2 = TMP_ALLOC_LIMBS (un); /* quotient, alternating */
73 np = TMP_ALLOC_LIMBS (un + LOG); /* powers of V */
77 /* FIXME: This allocation need indicate a flaw in the current itch mechanism:
78 Which operands not greater than un,un will incur the worst itch? We need
79 a parallel foo_maxitch set of functions. */
80 scratch_out = TMP_ALLOC_LIMBS (mpn_bdiv_qr_itch (un, un >> 1));
82 MPN_COPY (qp, up, un);
88 mpn_bdiv_qr (qp2, tp, qp, qn, pp, pn, scratch_out);
89 if (!mpn_zero_p (tp, pn))
90 break; /* could not divide by V^npowers */
92 MP_PTR_SWAP (qp, qp2);
100 if (((mp_bitcnt_t) 2 << npowers) - 1 > cap)
103 nn = 2 * pn - 1; /* next power will be at least this many limbs */
105 break; /* next power would be overlarge */
107 mpn_sqr (np, pp, pn);
114 pwr = ((mp_bitcnt_t) 1 << npowers) - 1;
116 for (i = npowers - 1; i >= 0; i--)
123 if (pwr + ((mp_bitcnt_t) 1 << i) > cap)
124 continue; /* V^i would bring us past cap */
126 mpn_bdiv_qr (qp2, tp, qp, qn, pp, pn, scratch_out);
127 if (!mpn_zero_p (tp, pn))
128 continue; /* could not divide by V^i */
130 MP_PTR_SWAP (qp, qp2);
134 pwr += (mp_bitcnt_t) 1 << i;
137 MPN_COPY (wp, qp, qn);