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