1/* mpn_remove -- divide out all multiples of odd mpn number from another mpn 2 number. 3 4 Contributed to the GNU project by Torbjorn Granlund. 5 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. 9 10Copyright 2009 Free Software Foundation, Inc. 11 12This file is part of the GNU MP Library. 13 14The GNU MP Library is free software; you can redistribute it and/or modify 15it under the terms of the GNU Lesser General Public License as published by 16the Free Software Foundation; either version 3 of the License, or (at your 17option) any later version. 18 19The GNU MP Library is distributed in the hope that it will be useful, but 20WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 21or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 22License for more details. 23 24You should have received a copy of the GNU Lesser General Public License 25along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 26 27#include "gmp.h" 28#include "gmp-impl.h" 29 30#if GMP_LIMB_BITS > 50 31#define LOG 50 32#else 33#define LOG GMP_LIMB_BITS 34#endif 35 36 37/* Input: U = {up,un}, V = {vp,vn} must be odd, cap 38 Ouput W = {wp,*wn} allocation need is exactly *wn 39 40 Set W = U / V^k, where k is the largest integer <= cap such that the 41 division yields an integer. 42 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. 48*/ 49 50mp_bitcnt_t 51mpn_remove (mp_ptr wp, mp_size_t *wn, 52 mp_ptr up, mp_size_t un, mp_ptr vp, mp_size_t vn, 53 mp_bitcnt_t cap) 54{ 55 mp_ptr pwpsp[LOG]; 56 mp_size_t pwpsn[LOG]; 57 mp_size_t npowers; 58 mp_ptr tp, qp, np, pp, qp2, scratch_out; 59 mp_size_t pn, nn, qn, i; 60 mp_bitcnt_t pwr; 61 TMP_DECL; 62 63 ASSERT (un > 0); 64 ASSERT (vn > 0); 65 ASSERT (vp[0] % 2 != 0); /* 2-adic division wants odd numbers */ 66 ASSERT (vn > 1 || vp[0] > 1); /* else we would loop indefinitely */ 67 68 TMP_MARK; 69 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 */ 74 pp = vp; 75 pn = vn; 76 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)); 81 82 MPN_COPY (qp, up, un); 83 qn = un; 84 85 npowers = 0; 86 while (qn >= pn) 87 { 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 */ 91 92 MP_PTR_SWAP (qp, qp2); 93 qn = qn - pn; 94 qn += qp[qn] != 0; 95 96 pwpsp[npowers] = pp; 97 pwpsn[npowers] = pn; 98 npowers++; 99 100 if (((mp_bitcnt_t) 2 << npowers) - 1 > cap) 101 break; 102 103 nn = 2 * pn - 1; /* next power will be at least this many limbs */ 104 if (nn > qn) 105 break; /* next power would be overlarge */ 106 107 mpn_sqr (np, pp, pn); 108 nn += np[nn] != 0; 109 pp = np; 110 pn = nn; 111 np += nn; 112 } 113 114 pwr = ((mp_bitcnt_t) 1 << npowers) - 1; 115 116 for (i = npowers - 1; i >= 0; i--) 117 { 118 pp = pwpsp[i]; 119 pn = pwpsn[i]; 120 if (qn < pn) 121 continue; 122 123 if (pwr + ((mp_bitcnt_t) 1 << i) > cap) 124 continue; /* V^i would bring us past cap */ 125 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 */ 129 130 MP_PTR_SWAP (qp, qp2); 131 qn = qn - pn; 132 qn += qp[qn] != 0; 133 134 pwr += (mp_bitcnt_t) 1 << i; 135 } 136 137 MPN_COPY (wp, qp, qn); 138 *wn = qn; 139 140 TMP_FREE; 141 142 return pwr; 143} 144