xref: /dflybsd-src/contrib/gmp/mpn/generic/remove.c (revision 86d7f5d305c6adaa56ff4582ece9859d73106103)
1*86d7f5d3SJohn Marino /* mpn_remove -- divide out all multiples of odd mpn number from another mpn
2*86d7f5d3SJohn Marino    number.
3*86d7f5d3SJohn Marino 
4*86d7f5d3SJohn Marino    Contributed to the GNU project by Torbjorn Granlund.
5*86d7f5d3SJohn Marino 
6*86d7f5d3SJohn Marino    THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
7*86d7f5d3SJohn Marino    SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
8*86d7f5d3SJohn Marino    GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
9*86d7f5d3SJohn Marino 
10*86d7f5d3SJohn Marino Copyright 2009 Free Software Foundation, Inc.
11*86d7f5d3SJohn Marino 
12*86d7f5d3SJohn Marino This file is part of the GNU MP Library.
13*86d7f5d3SJohn Marino 
14*86d7f5d3SJohn Marino The GNU MP Library is free software; you can redistribute it and/or modify
15*86d7f5d3SJohn Marino it under the terms of the GNU Lesser General Public License as published by
16*86d7f5d3SJohn Marino the Free Software Foundation; either version 3 of the License, or (at your
17*86d7f5d3SJohn Marino option) any later version.
18*86d7f5d3SJohn Marino 
19*86d7f5d3SJohn Marino The GNU MP Library is distributed in the hope that it will be useful, but
20*86d7f5d3SJohn Marino WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
21*86d7f5d3SJohn Marino or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
22*86d7f5d3SJohn Marino License for more details.
23*86d7f5d3SJohn Marino 
24*86d7f5d3SJohn Marino You should have received a copy of the GNU Lesser General Public License
25*86d7f5d3SJohn Marino along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
26*86d7f5d3SJohn Marino 
27*86d7f5d3SJohn Marino #include "gmp.h"
28*86d7f5d3SJohn Marino #include "gmp-impl.h"
29*86d7f5d3SJohn Marino 
30*86d7f5d3SJohn Marino #if GMP_LIMB_BITS > 50
31*86d7f5d3SJohn Marino #define LOG 50
32*86d7f5d3SJohn Marino #else
33*86d7f5d3SJohn Marino #define LOG GMP_LIMB_BITS
34*86d7f5d3SJohn Marino #endif
35*86d7f5d3SJohn Marino 
36*86d7f5d3SJohn Marino 
37*86d7f5d3SJohn Marino /* Input: U = {up,un}, V = {vp,vn} must be odd, cap
38*86d7f5d3SJohn Marino    Ouput  W = {wp,*wn} allocation need is exactly *wn
39*86d7f5d3SJohn Marino 
40*86d7f5d3SJohn Marino    Set W = U / V^k, where k is the largest integer <= cap such that the
41*86d7f5d3SJohn Marino    division yields an integer.
42*86d7f5d3SJohn Marino 
43*86d7f5d3SJohn Marino    FIXME: We currently allow any operand overlap.  This is quite non mpn-ish
44*86d7f5d3SJohn Marino    and might be changed, since it cost significant temporary space.
45*86d7f5d3SJohn Marino    * If we require W to have space for un limbs, we could save qp or qp2 (but
46*86d7f5d3SJohn Marino      we will still need to copy things into wp 50% of the time).
47*86d7f5d3SJohn Marino    * If we allow ourselves to clobber U, we could save the other of qp and qp2.
48*86d7f5d3SJohn Marino */
49*86d7f5d3SJohn Marino 
50*86d7f5d3SJohn 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)51*86d7f5d3SJohn Marino mpn_remove (mp_ptr wp, mp_size_t *wn,
52*86d7f5d3SJohn Marino 	    mp_ptr up, mp_size_t un, mp_ptr vp, mp_size_t vn,
53*86d7f5d3SJohn Marino 	    mp_bitcnt_t cap)
54*86d7f5d3SJohn Marino {
55*86d7f5d3SJohn Marino   mp_ptr    pwpsp[LOG];
56*86d7f5d3SJohn Marino   mp_size_t pwpsn[LOG];
57*86d7f5d3SJohn Marino   mp_size_t npowers;
58*86d7f5d3SJohn Marino   mp_ptr tp, qp, np, pp, qp2, scratch_out;
59*86d7f5d3SJohn Marino   mp_size_t pn, nn, qn, i;
60*86d7f5d3SJohn Marino   mp_bitcnt_t pwr;
61*86d7f5d3SJohn Marino   TMP_DECL;
62*86d7f5d3SJohn Marino 
63*86d7f5d3SJohn Marino   ASSERT (un > 0);
64*86d7f5d3SJohn Marino   ASSERT (vn > 0);
65*86d7f5d3SJohn Marino   ASSERT (vp[0] % 2 != 0);	/* 2-adic division wants odd numbers */
66*86d7f5d3SJohn Marino   ASSERT (vn > 1 || vp[0] > 1);	/* else we would loop indefinitely */
67*86d7f5d3SJohn Marino 
68*86d7f5d3SJohn Marino   TMP_MARK;
69*86d7f5d3SJohn Marino 
70*86d7f5d3SJohn Marino   tp = TMP_ALLOC_LIMBS ((un + vn) / 2); /* remainder */
71*86d7f5d3SJohn Marino   qp = TMP_ALLOC_LIMBS (un);		/* quotient, alternating */
72*86d7f5d3SJohn Marino   qp2 = TMP_ALLOC_LIMBS (un);		/* quotient, alternating */
73*86d7f5d3SJohn Marino   np = TMP_ALLOC_LIMBS (un + LOG);	/* powers of V */
74*86d7f5d3SJohn Marino   pp = vp;
75*86d7f5d3SJohn Marino   pn = vn;
76*86d7f5d3SJohn Marino 
77*86d7f5d3SJohn Marino   /* FIXME: This allocation need indicate a flaw in the current itch mechanism:
78*86d7f5d3SJohn Marino      Which operands not greater than un,un will incur the worst itch?  We need
79*86d7f5d3SJohn Marino      a parallel foo_maxitch set of functions.  */
80*86d7f5d3SJohn Marino   scratch_out = TMP_ALLOC_LIMBS (mpn_bdiv_qr_itch (un, un >> 1));
81*86d7f5d3SJohn Marino 
82*86d7f5d3SJohn Marino   MPN_COPY (qp, up, un);
83*86d7f5d3SJohn Marino   qn = un;
84*86d7f5d3SJohn Marino 
85*86d7f5d3SJohn Marino   npowers = 0;
86*86d7f5d3SJohn Marino   while (qn >= pn)
87*86d7f5d3SJohn Marino     {
88*86d7f5d3SJohn Marino       mpn_bdiv_qr (qp2, tp, qp, qn, pp, pn, scratch_out);
89*86d7f5d3SJohn Marino       if (!mpn_zero_p (tp, pn))
90*86d7f5d3SJohn Marino 	break;			/* could not divide by V^npowers */
91*86d7f5d3SJohn Marino 
92*86d7f5d3SJohn Marino       MP_PTR_SWAP (qp, qp2);
93*86d7f5d3SJohn Marino       qn = qn - pn;
94*86d7f5d3SJohn Marino       qn += qp[qn] != 0;
95*86d7f5d3SJohn Marino 
96*86d7f5d3SJohn Marino       pwpsp[npowers] = pp;
97*86d7f5d3SJohn Marino       pwpsn[npowers] = pn;
98*86d7f5d3SJohn Marino       npowers++;
99*86d7f5d3SJohn Marino 
100*86d7f5d3SJohn Marino       if (((mp_bitcnt_t) 2 << npowers) - 1 > cap)
101*86d7f5d3SJohn Marino 	break;
102*86d7f5d3SJohn Marino 
103*86d7f5d3SJohn Marino       nn = 2 * pn - 1;		/* next power will be at least this many limbs */
104*86d7f5d3SJohn Marino       if (nn > qn)
105*86d7f5d3SJohn Marino 	break;			/* next power would be overlarge */
106*86d7f5d3SJohn Marino 
107*86d7f5d3SJohn Marino       mpn_sqr (np, pp, pn);
108*86d7f5d3SJohn Marino       nn += np[nn] != 0;
109*86d7f5d3SJohn Marino       pp = np;
110*86d7f5d3SJohn Marino       pn = nn;
111*86d7f5d3SJohn Marino       np += nn;
112*86d7f5d3SJohn Marino     }
113*86d7f5d3SJohn Marino 
114*86d7f5d3SJohn Marino   pwr = ((mp_bitcnt_t) 1 << npowers) - 1;
115*86d7f5d3SJohn Marino 
116*86d7f5d3SJohn Marino   for (i = npowers - 1; i >= 0; i--)
117*86d7f5d3SJohn Marino     {
118*86d7f5d3SJohn Marino       pp = pwpsp[i];
119*86d7f5d3SJohn Marino       pn = pwpsn[i];
120*86d7f5d3SJohn Marino       if (qn < pn)
121*86d7f5d3SJohn Marino 	continue;
122*86d7f5d3SJohn Marino 
123*86d7f5d3SJohn Marino       if (pwr + ((mp_bitcnt_t) 1 << i) > cap)
124*86d7f5d3SJohn Marino 	continue;		/* V^i would bring us past cap */
125*86d7f5d3SJohn Marino 
126*86d7f5d3SJohn Marino       mpn_bdiv_qr (qp2, tp, qp, qn, pp, pn, scratch_out);
127*86d7f5d3SJohn Marino       if (!mpn_zero_p (tp, pn))
128*86d7f5d3SJohn Marino 	continue;		/* could not divide by V^i */
129*86d7f5d3SJohn Marino 
130*86d7f5d3SJohn Marino       MP_PTR_SWAP (qp, qp2);
131*86d7f5d3SJohn Marino       qn = qn - pn;
132*86d7f5d3SJohn Marino       qn += qp[qn] != 0;
133*86d7f5d3SJohn Marino 
134*86d7f5d3SJohn Marino       pwr += (mp_bitcnt_t) 1 << i;
135*86d7f5d3SJohn Marino     }
136*86d7f5d3SJohn Marino 
137*86d7f5d3SJohn Marino   MPN_COPY (wp, qp, qn);
138*86d7f5d3SJohn Marino   *wn = qn;
139*86d7f5d3SJohn Marino 
140*86d7f5d3SJohn Marino   TMP_FREE;
141*86d7f5d3SJohn Marino 
142*86d7f5d3SJohn Marino   return pwr;
143*86d7f5d3SJohn Marino }
144