Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpn/remove.c
Line
Count
Source (jump to first uncovered line)
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
10
Copyright 2009, 2012-2014, 2017 Free Software Foundation, Inc.
11
12
This file is part of the GNU MP Library.
13
14
The GNU MP Library is free software; you can redistribute it and/or modify
15
it under the terms of either:
16
17
  * the GNU Lesser General Public License as published by the Free
18
    Software Foundation; either version 3 of the License, or (at your
19
    option) any later version.
20
21
or
22
23
  * the GNU General Public License as published by the Free Software
24
    Foundation; either version 2 of the License, or (at your option) any
25
    later version.
26
27
or both in parallel, as here.
28
29
The GNU MP Library is distributed in the hope that it will be useful, but
30
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
31
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
32
for more details.
33
34
You should have received copies of the GNU General Public License and the
35
GNU Lesser General Public License along with the GNU MP Library.  If not,
36
see https://www.gnu.org/licenses/.  */
37
38
#include "gmp-impl.h"
39
40
#if GMP_LIMB_BITS > 50
41
#define LOG 50
42
#else
43
#define LOG GMP_LIMB_BITS
44
#endif
45
46
47
/* Input: U = {up,un}, V = {vp,vn} must be odd, cap
48
   Ouput  W = {wp,*wn} allocation need is exactly *wn
49
50
   Set W = U / V^k, where k is the largest integer <= cap such that the
51
   division yields an integer.
52
53
   FIXME: We currently allow any operand overlap.  This is quite non mpn-ish
54
   and might be changed, since it cost significant temporary space.
55
   * If we require W to have space for un + 1 limbs, we could save qp or qp2
56
     (but we will still need to copy things into wp 50% of the time).
57
   * If we allow ourselves to clobber U, we could save the other of qp and qp2,
58
     and the initial COPY (but also here we would need un + 1 limbs).
59
*/
60
61
/* FIXME: We need to wrap mpn_bdiv_qr due to the itch interface.  This need
62
   indicates a flaw in the current itch mechanism: Which operands not greater
63
   than un,un will incur the worst itch?  We need a parallel foo_maxitch set
64
   of functions.  */
65
static void
66
mpn_bdiv_qr_wrap (mp_ptr qp, mp_ptr rp,
67
      mp_srcptr np, mp_size_t nn,
68
      mp_srcptr dp, mp_size_t dn)
69
1.72k
{
70
1.72k
  mp_ptr scratch_out;
71
1.72k
  TMP_DECL;
72
73
1.72k
  TMP_MARK;
74
1.72k
  scratch_out = TMP_ALLOC_LIMBS (mpn_bdiv_qr_itch (nn, dn));
75
1.72k
  mpn_bdiv_qr (qp, rp, np, nn, dp, dn, scratch_out);
76
77
1.72k
  TMP_FREE;
78
1.72k
}
79
80
mp_bitcnt_t
81
mpn_remove (mp_ptr wp, mp_size_t *wn,
82
      mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn,
83
      mp_bitcnt_t cap)
84
240
{
85
240
  mp_srcptr pwpsp[LOG];
86
240
  mp_size_t pwpsn[LOG];
87
240
  mp_size_t npowers;
88
240
  mp_ptr tp, qp, np, qp2;
89
240
  mp_srcptr pp;
90
240
  mp_size_t pn, nn, qn, i;
91
240
  mp_bitcnt_t pwr;
92
240
  TMP_DECL;
93
94
240
  ASSERT (un > 0);
95
240
  ASSERT (vn > 0);
96
240
  ASSERT (vp[0] % 2 != 0);  /* 2-adic division wants odd numbers */
97
240
  ASSERT (vn > 1 || vp[0] > 1); /* else we would loop indefinitely */
98
99
240
  TMP_MARK;
100
101
240
  TMP_ALLOC_LIMBS_3 (qp, un + 1, /* quotient, alternating */
102
240
         qp2, un + 1, /* quotient, alternating */
103
240
         tp, (un + 1 + vn) / 2); /* remainder */
104
240
  pp = vp;
105
240
  pn = vn;
106
107
240
  MPN_COPY (qp, up, un);
108
240
  qn = un;
109
110
240
  npowers = 0;
111
978
  while (qn >= pn)
112
978
    {
113
978
      qp[qn] = 0;
114
978
      mpn_bdiv_qr_wrap (qp2, tp, qp, qn + 1, pp, pn);
115
978
      if (!mpn_zero_p (tp, pn))
116
978
  {
117
978
    if (mpn_cmp (tp, pp, pn) != 0)
118
224
      break;   /* could not divide by V^npowers */
119
978
  }
120
121
754
      MP_PTR_SWAP (qp, qp2);
122
754
      qn = qn - pn;
123
754
      mpn_neg (qp, qp, qn+1);
124
125
754
      qn += qp[qn] != 0;
126
127
754
      pwpsp[npowers] = pp;
128
754
      pwpsn[npowers] = pn;
129
754
      ++npowers;
130
131
754
      if (((mp_bitcnt_t) 2 << npowers) - 1 > cap)
132
0
  break;
133
134
754
      nn = 2 * pn - 1;    /* next power will be at least this large */
135
754
      if (nn > qn)
136
16
  break;      /* next power would be overlarge */
137
138
738
      if (npowers == 1)    /* Alloc once, but only if it's needed */
139
240
  np = TMP_ALLOC_LIMBS (qn + LOG);  /* powers of V */
140
498
      else
141
498
  np += pn;
142
143
738
      mpn_sqr (np, pp, pn);
144
738
      pn = nn + (np[nn] != 0);
145
738
      pp = np;
146
738
    }
147
148
240
  pwr = ((mp_bitcnt_t) 1 << npowers) - 1;
149
150
994
  for (i = npowers; --i >= 0;)
151
754
    {
152
754
      pn = pwpsn[i];
153
754
      if (qn < pn)
154
5
  continue;
155
156
749
      if (pwr + ((mp_bitcnt_t) 1 << i) > cap)
157
0
  continue;   /* V^i would bring us past cap */
158
159
749
      qp[qn] = 0;
160
749
      mpn_bdiv_qr_wrap (qp2, tp, qp, qn + 1, pwpsp[i], pn);
161
749
      if (!mpn_zero_p (tp, pn))
162
749
  {
163
749
    if (mpn_cmp (tp, pwpsp[i], pn) != 0)
164
346
      continue;   /* could not divide by V^i */
165
749
  }
166
167
403
      MP_PTR_SWAP (qp, qp2);
168
403
      qn = qn - pn;
169
403
      mpn_neg (qp, qp, qn+1);
170
171
403
      qn += qp[qn] != 0;
172
173
403
      pwr += (mp_bitcnt_t) 1 << i;
174
403
    }
175
176
240
  MPN_COPY (wp, qp, qn);
177
240
  *wn = qn;
178
179
240
  TMP_FREE;
180
181
240
  return pwr;
182
240
}