xref: /dflybsd-src/contrib/gmp/mpn/generic/remove.c (revision 86d7f5d305c6adaa56ff4582ece9859d73106103)
186d7f5d3SJohn Marino /* mpn_remove -- divide out all multiples of odd mpn number from another mpn
286d7f5d3SJohn Marino    number.
386d7f5d3SJohn Marino 
486d7f5d3SJohn Marino    Contributed to the GNU project by Torbjorn Granlund.
586d7f5d3SJohn Marino 
686d7f5d3SJohn Marino    THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
786d7f5d3SJohn Marino    SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
886d7f5d3SJohn Marino    GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
986d7f5d3SJohn Marino 
1086d7f5d3SJohn Marino Copyright 2009 Free Software Foundation, Inc.
1186d7f5d3SJohn Marino 
1286d7f5d3SJohn Marino This file is part of the GNU MP Library.
1386d7f5d3SJohn Marino 
1486d7f5d3SJohn Marino The GNU MP Library is free software; you can redistribute it and/or modify
1586d7f5d3SJohn Marino it under the terms of the GNU Lesser General Public License as published by
1686d7f5d3SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
1786d7f5d3SJohn Marino option) any later version.
1886d7f5d3SJohn Marino 
1986d7f5d3SJohn Marino The GNU MP Library is distributed in the hope that it will be useful, but
2086d7f5d3SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
2186d7f5d3SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
2286d7f5d3SJohn Marino License for more details.
2386d7f5d3SJohn Marino 
2486d7f5d3SJohn Marino You should have received a copy of the GNU Lesser General Public License
2586d7f5d3SJohn Marino along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
2686d7f5d3SJohn Marino 
2786d7f5d3SJohn Marino #include "gmp.h"
2886d7f5d3SJohn Marino #include "gmp-impl.h"
2986d7f5d3SJohn Marino 
3086d7f5d3SJohn Marino #if GMP_LIMB_BITS > 50
3186d7f5d3SJohn Marino #define LOG 50
3286d7f5d3SJohn Marino #else
3386d7f5d3SJohn Marino #define LOG GMP_LIMB_BITS
3486d7f5d3SJohn Marino #endif
3586d7f5d3SJohn Marino 
3686d7f5d3SJohn Marino 
3786d7f5d3SJohn Marino /* Input: U = {up,un}, V = {vp,vn} must be odd, cap
3886d7f5d3SJohn Marino    Ouput  W = {wp,*wn} allocation need is exactly *wn
3986d7f5d3SJohn Marino 
4086d7f5d3SJohn Marino    Set W = U / V^k, where k is the largest integer <= cap such that the
4186d7f5d3SJohn Marino    division yields an integer.
4286d7f5d3SJohn Marino 
4386d7f5d3SJohn Marino    FIXME: We currently allow any operand overlap.  This is quite non mpn-ish
4486d7f5d3SJohn Marino    and might be changed, since it cost significant temporary space.
4586d7f5d3SJohn Marino    * If we require W to have space for un limbs, we could save qp or qp2 (but
4686d7f5d3SJohn Marino      we will still need to copy things into wp 50% of the time).
4786d7f5d3SJohn Marino    * If we allow ourselves to clobber U, we could save the other of qp and qp2.
4886d7f5d3SJohn Marino */
4986d7f5d3SJohn Marino 
5086d7f5d3SJohn Marino 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)5186d7f5d3SJohn Marino mpn_remove (mp_ptr wp, mp_size_t *wn,
5286d7f5d3SJohn Marino 	    mp_ptr up, mp_size_t un, mp_ptr vp, mp_size_t vn,
5386d7f5d3SJohn Marino 	    mp_bitcnt_t cap)
5486d7f5d3SJohn Marino {
5586d7f5d3SJohn Marino   mp_ptr    pwpsp[LOG];
5686d7f5d3SJohn Marino   mp_size_t pwpsn[LOG];
5786d7f5d3SJohn Marino   mp_size_t npowers;
5886d7f5d3SJohn Marino   mp_ptr tp, qp, np, pp, qp2, scratch_out;
5986d7f5d3SJohn Marino   mp_size_t pn, nn, qn, i;
6086d7f5d3SJohn Marino   mp_bitcnt_t pwr;
6186d7f5d3SJohn Marino   TMP_DECL;
6286d7f5d3SJohn Marino 
6386d7f5d3SJohn Marino   ASSERT (un > 0);
6486d7f5d3SJohn Marino   ASSERT (vn > 0);
6586d7f5d3SJohn Marino   ASSERT (vp[0] % 2 != 0);	/* 2-adic division wants odd numbers */
6686d7f5d3SJohn Marino   ASSERT (vn > 1 || vp[0] > 1);	/* else we would loop indefinitely */
6786d7f5d3SJohn Marino 
6886d7f5d3SJohn Marino   TMP_MARK;
6986d7f5d3SJohn Marino 
7086d7f5d3SJohn Marino   tp = TMP_ALLOC_LIMBS ((un + vn) / 2); /* remainder */
7186d7f5d3SJohn Marino   qp = TMP_ALLOC_LIMBS (un);		/* quotient, alternating */
7286d7f5d3SJohn Marino   qp2 = TMP_ALLOC_LIMBS (un);		/* quotient, alternating */
7386d7f5d3SJohn Marino   np = TMP_ALLOC_LIMBS (un + LOG);	/* powers of V */
7486d7f5d3SJohn Marino   pp = vp;
7586d7f5d3SJohn Marino   pn = vn;
7686d7f5d3SJohn Marino 
7786d7f5d3SJohn Marino   /* FIXME: This allocation need indicate a flaw in the current itch mechanism:
7886d7f5d3SJohn Marino      Which operands not greater than un,un will incur the worst itch?  We need
7986d7f5d3SJohn Marino      a parallel foo_maxitch set of functions.  */
8086d7f5d3SJohn Marino   scratch_out = TMP_ALLOC_LIMBS (mpn_bdiv_qr_itch (un, un >> 1));
8186d7f5d3SJohn Marino 
8286d7f5d3SJohn Marino   MPN_COPY (qp, up, un);
8386d7f5d3SJohn Marino   qn = un;
8486d7f5d3SJohn Marino 
8586d7f5d3SJohn Marino   npowers = 0;
8686d7f5d3SJohn Marino   while (qn >= pn)
8786d7f5d3SJohn Marino     {
8886d7f5d3SJohn Marino       mpn_bdiv_qr (qp2, tp, qp, qn, pp, pn, scratch_out);
8986d7f5d3SJohn Marino       if (!mpn_zero_p (tp, pn))
9086d7f5d3SJohn Marino 	break;			/* could not divide by V^npowers */
9186d7f5d3SJohn Marino 
9286d7f5d3SJohn Marino       MP_PTR_SWAP (qp, qp2);
9386d7f5d3SJohn Marino       qn = qn - pn;
9486d7f5d3SJohn Marino       qn += qp[qn] != 0;
9586d7f5d3SJohn Marino 
9686d7f5d3SJohn Marino       pwpsp[npowers] = pp;
9786d7f5d3SJohn Marino       pwpsn[npowers] = pn;
9886d7f5d3SJohn Marino       npowers++;
9986d7f5d3SJohn Marino 
10086d7f5d3SJohn Marino       if (((mp_bitcnt_t) 2 << npowers) - 1 > cap)
10186d7f5d3SJohn Marino 	break;
10286d7f5d3SJohn Marino 
10386d7f5d3SJohn Marino       nn = 2 * pn - 1;		/* next power will be at least this many limbs */
10486d7f5d3SJohn Marino       if (nn > qn)
10586d7f5d3SJohn Marino 	break;			/* next power would be overlarge */
10686d7f5d3SJohn Marino 
10786d7f5d3SJohn Marino       mpn_sqr (np, pp, pn);
10886d7f5d3SJohn Marino       nn += np[nn] != 0;
10986d7f5d3SJohn Marino       pp = np;
11086d7f5d3SJohn Marino       pn = nn;
11186d7f5d3SJohn Marino       np += nn;
11286d7f5d3SJohn Marino     }
11386d7f5d3SJohn Marino 
11486d7f5d3SJohn Marino   pwr = ((mp_bitcnt_t) 1 << npowers) - 1;
11586d7f5d3SJohn Marino 
11686d7f5d3SJohn Marino   for (i = npowers - 1; i >= 0; i--)
11786d7f5d3SJohn Marino     {
11886d7f5d3SJohn Marino       pp = pwpsp[i];
11986d7f5d3SJohn Marino       pn = pwpsn[i];
12086d7f5d3SJohn Marino       if (qn < pn)
12186d7f5d3SJohn Marino 	continue;
12286d7f5d3SJohn Marino 
12386d7f5d3SJohn Marino       if (pwr + ((mp_bitcnt_t) 1 << i) > cap)
12486d7f5d3SJohn Marino 	continue;		/* V^i would bring us past cap */
12586d7f5d3SJohn Marino 
12686d7f5d3SJohn Marino       mpn_bdiv_qr (qp2, tp, qp, qn, pp, pn, scratch_out);
12786d7f5d3SJohn Marino       if (!mpn_zero_p (tp, pn))
12886d7f5d3SJohn Marino 	continue;		/* could not divide by V^i */
12986d7f5d3SJohn Marino 
13086d7f5d3SJohn Marino       MP_PTR_SWAP (qp, qp2);
13186d7f5d3SJohn Marino       qn = qn - pn;
13286d7f5d3SJohn Marino       qn += qp[qn] != 0;
13386d7f5d3SJohn Marino 
13486d7f5d3SJohn Marino       pwr += (mp_bitcnt_t) 1 << i;
13586d7f5d3SJohn Marino     }
13686d7f5d3SJohn Marino 
13786d7f5d3SJohn Marino   MPN_COPY (wp, qp, qn);
13886d7f5d3SJohn Marino   *wn = qn;
13986d7f5d3SJohn Marino 
14086d7f5d3SJohn Marino   TMP_FREE;
14186d7f5d3SJohn Marino 
14286d7f5d3SJohn Marino   return pwr;
14386d7f5d3SJohn Marino }
144