remove.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:4k
- /* mpn_remove -- divide out all multiples of odd mpn number from another mpn
- number.
- Contributed to the GNU project by Torbjorn Granlund.
- THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY
- SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
- GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
- Copyright 2009 Free Software Foundation, Inc.
- This file is part of the GNU MP Library.
- The GNU MP Library is free software; you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation; either version 3 of the License, or (at your
- option) any later version.
- The GNU MP Library is distributed in the hope that it will be useful, but
- WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
- License for more details.
- You should have received a copy of the GNU Lesser General Public License
- along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
- #include "gmp.h"
- #include "gmp-impl.h"
- #if GMP_LIMB_BITS > 50
- #define LOG 50
- #else
- #define LOG GMP_LIMB_BITS
- #endif
- /* Input: U = {up,un}, V = {vp,vn} must be odd, cap
- Ouput W = {wp,*wn} allocation need is exactly *wn
- Set W = U / V^k, where k is the largest integer <= cap such that the
- division yields an integer.
- FIXME: We currently allow any operand overlap. This is quite non mpn-ish
- and might be changed, since it cost significant temporary space.
- * If we require W to have space for un limbs, we could save qp or qp2 (but
- we will still need to copy things into wp 50% of the time).
- * If we allow ourselves to clobber U, we could save the other of qp and qp2.
- */
- mp_bitcnt_t
- mpn_remove (mp_ptr wp, mp_size_t *wn,
- mp_ptr up, mp_size_t un, mp_ptr vp, mp_size_t vn,
- mp_bitcnt_t cap)
- {
- mp_ptr pwpsp[LOG];
- mp_size_t pwpsn[LOG];
- mp_size_t npowers;
- mp_ptr tp, qp, np, pp, qp2, scratch_out;
- mp_size_t pn, nn, qn, i;
- mp_bitcnt_t pwr;
- TMP_DECL;
- ASSERT (un > 0);
- ASSERT (vn > 0);
- ASSERT (vp[0] % 2 != 0); /* 2-adic division wants odd numbers */
- ASSERT (vn > 1 || vp[0] > 1); /* else we would loop indefinitely */
- TMP_MARK;
- tp = TMP_ALLOC_LIMBS ((un + vn) / 2); /* remainder */
- qp = TMP_ALLOC_LIMBS (un); /* quotient, alternating */
- qp2 = TMP_ALLOC_LIMBS (un); /* quotient, alternating */
- np = TMP_ALLOC_LIMBS (un + LOG); /* powers of V */
- pp = vp;
- pn = vn;
- /* FIXME: This allocation need indicate a flaw in the current itch mechanism:
- Which operands not greater than un,un will incur the worst itch? We need
- a parallel foo_maxitch set of functions. */
- scratch_out = TMP_ALLOC_LIMBS (mpn_bdiv_qr_itch (un, un >> 1));
- MPN_COPY (qp, up, un);
- qn = un;
- npowers = 0;
- while (qn >= pn)
- {
- mpn_bdiv_qr (qp2, tp, qp, qn, pp, pn, scratch_out);
- if (!mpn_zero_p (tp, pn))
- break; /* could not divide by V^npowers */
- MP_PTR_SWAP (qp, qp2);
- qn = qn - pn;
- qn += qp[qn] != 0;
- pwpsp[npowers] = pp;
- pwpsn[npowers] = pn;
- npowers++;
- if (((mp_bitcnt_t) 2 << npowers) - 1 > cap)
- break;
- nn = 2 * pn - 1; /* next power will be at least this many limbs */
- if (nn > qn)
- break; /* next power would be overlarge */
- mpn_sqr (np, pp, pn);
- nn += np[nn] != 0;
- pp = np;
- pn = nn;
- np += nn;
- }
- pwr = ((mp_bitcnt_t) 1 << npowers) - 1;
- for (i = npowers - 1; i >= 0; i--)
- {
- pp = pwpsp[i];
- pn = pwpsn[i];
- if (qn < pn)
- continue;
- if (pwr + ((mp_bitcnt_t) 1 << i) > cap)
- continue; /* V^i would bring us past cap */
- mpn_bdiv_qr (qp2, tp, qp, qn, pp, pn, scratch_out);
- if (!mpn_zero_p (tp, pn))
- continue; /* could not divide by V^i */
- MP_PTR_SWAP (qp, qp2);
- qn = qn - pn;
- qn += qp[qn] != 0;
- pwr += (mp_bitcnt_t) 1 << i;
- }
- MPN_COPY (wp, qp, qn);
- *wn = qn;
- TMP_FREE;
- return pwr;
- }