Coverage Report

Created: 2024-11-21 06:47

/src/libgmp/mpn/tdiv_qr.c
Line
Count
Source (jump to first uncovered line)
1
/* mpn_tdiv_qr -- Divide the numerator (np,nn) by the denominator (dp,dn) and
2
   write the nn-dn+1 quotient limbs at qp and the dn remainder limbs at rp.  If
3
   qxn is non-zero, generate that many fraction limbs and append them after the
4
   other quotient limbs, and update the remainder accordingly.  The input
5
   operands are unaffected.
6
7
   Preconditions:
8
   1. The most significant limb of the divisor must be non-zero.
9
   2. nn >= dn, even if qxn is non-zero.  (??? relax this ???)
10
11
   The time complexity of this is O(qn*qn+M(dn,qn)), where M(m,n) is the time
12
   complexity of multiplication.
13
14
Copyright 1997, 2000-2002, 2005, 2009, 2015 Free Software Foundation, Inc.
15
16
This file is part of the GNU MP Library.
17
18
The GNU MP Library is free software; you can redistribute it and/or modify
19
it under the terms of either:
20
21
  * the GNU Lesser General Public License as published by the Free
22
    Software Foundation; either version 3 of the License, or (at your
23
    option) any later version.
24
25
or
26
27
  * the GNU General Public License as published by the Free Software
28
    Foundation; either version 2 of the License, or (at your option) any
29
    later version.
30
31
or both in parallel, as here.
32
33
The GNU MP Library is distributed in the hope that it will be useful, but
34
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
35
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
36
for more details.
37
38
You should have received copies of the GNU General Public License and the
39
GNU Lesser General Public License along with the GNU MP Library.  If not,
40
see https://www.gnu.org/licenses/.  */
41
42
#include "gmp-impl.h"
43
#include "longlong.h"
44
45
46
void
47
mpn_tdiv_qr (mp_ptr qp, mp_ptr rp, mp_size_t qxn,
48
       mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
49
1.49M
{
50
1.49M
  ASSERT_ALWAYS (qxn == 0);
51
52
1.49M
  ASSERT (nn >= 0);
53
1.49M
  ASSERT (dn >= 0);
54
1.49M
  ASSERT (dn == 0 || dp[dn - 1] != 0);
55
1.49M
  ASSERT (! MPN_OVERLAP_P (qp, nn - dn + 1 + qxn, np, nn));
56
1.49M
  ASSERT (! MPN_OVERLAP_P (qp, nn - dn + 1 + qxn, dp, dn));
57
58
1.49M
  switch (dn)
59
1.49M
    {
60
0
    case 0:
61
0
      DIVIDE_BY_ZERO;
62
63
2.70k
    case 1:
64
2.70k
      {
65
2.70k
  rp[0] = mpn_divrem_1 (qp, (mp_size_t) 0, np, nn, dp[0]);
66
2.70k
  return;
67
0
      }
68
69
504
    case 2:
70
504
      {
71
504
  mp_ptr n2p;
72
504
  mp_limb_t qhl, cy;
73
504
  TMP_DECL;
74
504
  TMP_MARK;
75
504
  if ((dp[1] & GMP_NUMB_HIGHBIT) == 0)
76
487
    {
77
487
      int cnt;
78
487
      mp_limb_t d2p[2];
79
487
      count_leading_zeros (cnt, dp[1]);
80
487
      cnt -= GMP_NAIL_BITS;
81
487
      d2p[1] = (dp[1] << cnt) | (dp[0] >> (GMP_NUMB_BITS - cnt));
82
487
      d2p[0] = (dp[0] << cnt) & GMP_NUMB_MASK;
83
487
      n2p = TMP_ALLOC_LIMBS (nn + 1);
84
487
      cy = mpn_lshift (n2p, np, nn, cnt);
85
487
      n2p[nn] = cy;
86
487
      qhl = mpn_divrem_2 (qp, 0L, n2p, nn + (cy != 0), d2p);
87
487
      if (cy == 0)
88
110
        qp[nn - 2] = qhl; /* always store nn-2+1 quotient limbs */
89
487
      rp[0] = (n2p[0] >> cnt)
90
487
        | ((n2p[1] << (GMP_NUMB_BITS - cnt)) & GMP_NUMB_MASK);
91
487
      rp[1] = (n2p[1] >> cnt);
92
487
    }
93
17
  else
94
17
    {
95
17
      n2p = TMP_ALLOC_LIMBS (nn);
96
17
      MPN_COPY (n2p, np, nn);
97
17
      qhl = mpn_divrem_2 (qp, 0L, n2p, nn, dp);
98
17
      qp[nn - 2] = qhl; /* always store nn-2+1 quotient limbs */
99
17
      rp[0] = n2p[0];
100
17
      rp[1] = n2p[1];
101
17
    }
102
504
  TMP_FREE;
103
504
  return;
104
0
      }
105
106
1.49M
    default:
107
1.49M
      {
108
1.49M
  int adjust;
109
1.49M
  gmp_pi1_t dinv;
110
1.49M
  TMP_DECL;
111
1.49M
  TMP_MARK;
112
1.49M
  adjust = np[nn - 1] >= dp[dn - 1];  /* conservative tests for quotient size */
113
1.49M
  if (nn + adjust >= 2 * dn)
114
1.47M
    {
115
1.47M
      mp_ptr n2p, d2p;
116
1.47M
      mp_limb_t cy;
117
1.47M
      int cnt;
118
119
1.47M
      qp[nn - dn] = 0;        /* zero high quotient limb */
120
1.47M
      if ((dp[dn - 1] & GMP_NUMB_HIGHBIT) == 0) /* normalize divisor */
121
1.47M
        {
122
1.47M
    count_leading_zeros (cnt, dp[dn - 1]);
123
1.47M
    cnt -= GMP_NAIL_BITS;
124
1.47M
    d2p = TMP_ALLOC_LIMBS (dn);
125
1.47M
    mpn_lshift (d2p, dp, dn, cnt);
126
1.47M
    n2p = TMP_ALLOC_LIMBS (nn + 1);
127
1.47M
    cy = mpn_lshift (n2p, np, nn, cnt);
128
1.47M
    n2p[nn] = cy;
129
1.47M
    nn += adjust;
130
1.47M
        }
131
430
      else
132
430
        {
133
430
    cnt = 0;
134
430
    d2p = (mp_ptr) dp;
135
430
    n2p = TMP_ALLOC_LIMBS (nn + 1);
136
430
    MPN_COPY (n2p, np, nn);
137
430
    n2p[nn] = 0;
138
430
    nn += adjust;
139
430
        }
140
141
1.47M
      invert_pi1 (dinv, d2p[dn - 1], d2p[dn - 2]);
142
1.47M
      if (BELOW_THRESHOLD (dn, DC_DIV_QR_THRESHOLD))
143
1.47M
        mpn_sbpi1_div_qr (qp, n2p, nn, d2p, dn, dinv.inv32);
144
1.34k
      else if (BELOW_THRESHOLD (dn, MUPI_DIV_QR_THRESHOLD) ||   /* fast condition */
145
1.34k
         BELOW_THRESHOLD (nn, 2 * MU_DIV_QR_THRESHOLD) || /* fast condition */
146
1.34k
         (double) (2 * (MU_DIV_QR_THRESHOLD - MUPI_DIV_QR_THRESHOLD)) * dn /* slow... */
147
0
         + (double) MUPI_DIV_QR_THRESHOLD * nn > (double) dn * nn)    /* ...condition */
148
1.34k
        mpn_dcpi1_div_qr (qp, n2p, nn, d2p, dn, &dinv);
149
0
      else
150
0
        {
151
0
    mp_size_t itch = mpn_mu_div_qr_itch (nn, dn, 0);
152
0
    mp_ptr scratch = TMP_ALLOC_LIMBS (itch);
153
0
    mpn_mu_div_qr (qp, rp, n2p, nn, d2p, dn, scratch);
154
0
    n2p = rp;
155
0
        }
156
157
1.47M
      if (cnt != 0)
158
1.47M
        mpn_rshift (rp, n2p, dn, cnt);
159
430
      else
160
430
        MPN_COPY (rp, n2p, dn);
161
1.47M
      TMP_FREE;
162
1.47M
      return;
163
1.47M
    }
164
165
  /* When we come here, the numerator/partial remainder is less
166
     than twice the size of the denominator.  */
167
168
16.1k
    {
169
      /* Problem:
170
171
         Divide a numerator N with nn limbs by a denominator D with dn
172
         limbs forming a quotient of qn=nn-dn+1 limbs.  When qn is small
173
         compared to dn, conventional division algorithms perform poorly.
174
         We want an algorithm that has an expected running time that is
175
         dependent only on qn.
176
177
         Algorithm (very informally stated):
178
179
         1) Divide the 2 x qn most significant limbs from the numerator
180
      by the qn most significant limbs from the denominator.  Call
181
      the result qest.  This is either the correct quotient, but
182
      might be 1 or 2 too large.  Compute the remainder from the
183
      division.  (This step is implemented by an mpn_divrem call.)
184
185
         2) Is the most significant limb from the remainder < p, where p
186
      is the product of the most significant limb from the quotient
187
      and the next(d)?  (Next(d) denotes the next ignored limb from
188
      the denominator.)  If it is, decrement qest, and adjust the
189
      remainder accordingly.
190
191
         3) Is the remainder >= qest?  If it is, qest is the desired
192
      quotient.  The algorithm terminates.
193
194
         4) Subtract qest x next(d) from the remainder.  If there is
195
      borrow out, decrement qest, and adjust the remainder
196
      accordingly.
197
198
         5) Skip one word from the denominator (i.e., let next(d) denote
199
      the next less significant limb.  */
200
201
16.1k
      mp_size_t qn;
202
16.1k
      mp_ptr n2p, d2p;
203
16.1k
      mp_ptr tp;
204
16.1k
      mp_limb_t cy;
205
16.1k
      mp_size_t in, rn;
206
16.1k
      mp_limb_t quotient_too_large;
207
16.1k
      unsigned int cnt;
208
209
16.1k
      qn = nn - dn;
210
16.1k
      qp[qn] = 0;       /* zero high quotient limb */
211
16.1k
      qn += adjust;     /* qn cannot become bigger */
212
213
16.1k
      if (qn == 0)
214
516
        {
215
516
    MPN_COPY (rp, np, dn);
216
516
    TMP_FREE;
217
516
    return;
218
516
        }
219
220
15.6k
      in = dn - qn;   /* (at least partially) ignored # of limbs in ops */
221
      /* Normalize denominator by shifting it to the left such that its
222
         most significant bit is set.  Then shift the numerator the same
223
         amount, to mathematically preserve quotient.  */
224
15.6k
      if ((dp[dn - 1] & GMP_NUMB_HIGHBIT) == 0)
225
14.8k
        {
226
14.8k
    count_leading_zeros (cnt, dp[dn - 1]);
227
14.8k
    cnt -= GMP_NAIL_BITS;
228
229
14.8k
    d2p = TMP_ALLOC_LIMBS (qn);
230
14.8k
    mpn_lshift (d2p, dp + in, qn, cnt);
231
14.8k
    d2p[0] |= dp[in - 1] >> (GMP_NUMB_BITS - cnt);
232
233
14.8k
    n2p = TMP_ALLOC_LIMBS (2 * qn + 1);
234
14.8k
    cy = mpn_lshift (n2p, np + nn - 2 * qn, 2 * qn, cnt);
235
14.8k
    if (adjust)
236
5.99k
      {
237
5.99k
        n2p[2 * qn] = cy;
238
5.99k
        n2p++;
239
5.99k
      }
240
8.83k
    else
241
8.83k
      {
242
8.83k
        n2p[0] |= np[nn - 2 * qn - 1] >> (GMP_NUMB_BITS - cnt);
243
8.83k
      }
244
14.8k
        }
245
781
      else
246
781
        {
247
781
    cnt = 0;
248
781
    d2p = (mp_ptr) dp + in;
249
250
781
    n2p = TMP_ALLOC_LIMBS (2 * qn + 1);
251
781
    MPN_COPY (n2p, np + nn - 2 * qn, 2 * qn);
252
781
    if (adjust)
253
78
      {
254
78
        n2p[2 * qn] = 0;
255
78
        n2p++;
256
78
      }
257
781
        }
258
259
      /* Get an approximate quotient using the extracted operands.  */
260
15.6k
      if (qn == 1)
261
1.98k
        {
262
1.98k
    mp_limb_t q0, r0;
263
1.98k
    udiv_qrnnd (q0, r0, n2p[1], n2p[0] << GMP_NAIL_BITS, d2p[0] << GMP_NAIL_BITS);
264
1.98k
    n2p[0] = r0 >> GMP_NAIL_BITS;
265
1.98k
    qp[0] = q0;
266
1.98k
        }
267
13.6k
      else if (qn == 2)
268
1.38k
        mpn_divrem_2 (qp, 0L, n2p, 4L, d2p); /* FIXME: obsolete function */
269
12.2k
      else
270
12.2k
        {
271
12.2k
    invert_pi1 (dinv, d2p[qn - 1], d2p[qn - 2]);
272
12.2k
    if (BELOW_THRESHOLD (qn, DC_DIV_QR_THRESHOLD))
273
11.6k
      mpn_sbpi1_div_qr (qp, n2p, 2 * qn, d2p, qn, dinv.inv32);
274
594
    else if (BELOW_THRESHOLD (qn, MU_DIV_QR_THRESHOLD))
275
594
      mpn_dcpi1_div_qr (qp, n2p, 2 * qn, d2p, qn, &dinv);
276
0
    else
277
0
      {
278
0
        mp_size_t itch = mpn_mu_div_qr_itch (2 * qn, qn, 0);
279
0
        mp_ptr scratch = TMP_ALLOC_LIMBS (itch);
280
0
        mp_ptr r2p = rp;
281
0
        if (np == r2p) /* If N and R share space, put ... */
282
0
          r2p += nn - qn; /* intermediate remainder at N's upper end. */
283
0
        mpn_mu_div_qr (qp, r2p, n2p, 2 * qn, d2p, qn, scratch);
284
0
        MPN_COPY (n2p, r2p, qn);
285
0
      }
286
12.2k
        }
287
288
15.6k
      rn = qn;
289
      /* Multiply the first ignored divisor limb by the most significant
290
         quotient limb.  If that product is > the partial remainder's
291
         most significant limb, we know the quotient is too large.  This
292
         test quickly catches most cases where the quotient is too large;
293
         it catches all cases where the quotient is 2 too large.  */
294
15.6k
      {
295
15.6k
        mp_limb_t dl, x;
296
15.6k
        mp_limb_t h, dummy;
297
298
15.6k
        if (in - 2 < 0)
299
4.59k
    dl = 0;
300
11.0k
        else
301
11.0k
    dl = dp[in - 2];
302
303
15.6k
#if GMP_NAIL_BITS == 0
304
15.6k
        x = (dp[in - 1] << cnt) | ((dl >> 1) >> ((~cnt) % GMP_LIMB_BITS));
305
#else
306
        x = (dp[in - 1] << cnt) & GMP_NUMB_MASK;
307
        if (cnt != 0)
308
    x |= dl >> (GMP_NUMB_BITS - cnt);
309
#endif
310
15.6k
        umul_ppmm (h, dummy, x, qp[qn - 1] << GMP_NAIL_BITS);
311
312
15.6k
        if (n2p[qn - 1] < h)
313
502
    {
314
502
      mp_limb_t cy;
315
316
502
      mpn_decr_u (qp, (mp_limb_t) 1);
317
502
      cy = mpn_add_n (n2p, n2p, d2p, qn);
318
502
      if (cy)
319
143
        {
320
          /* The partial remainder is safely large.  */
321
143
          n2p[qn] = cy;
322
143
          ++rn;
323
143
        }
324
502
    }
325
15.6k
      }
326
327
15.6k
      quotient_too_large = 0;
328
15.6k
      if (cnt != 0)
329
14.8k
        {
330
14.8k
    mp_limb_t cy1, cy2;
331
332
    /* Append partially used numerator limb to partial remainder.  */
333
14.8k
    cy1 = mpn_lshift (n2p, n2p, rn, GMP_NUMB_BITS - cnt);
334
14.8k
    n2p[0] |= np[in - 1] & (GMP_NUMB_MASK >> cnt);
335
336
    /* Update partial remainder with partially used divisor limb.  */
337
14.8k
    cy2 = mpn_submul_1 (n2p, qp, qn, dp[in - 1] & (GMP_NUMB_MASK >> cnt));
338
14.8k
    if (qn != rn)
339
133
      {
340
133
        ASSERT_ALWAYS (n2p[qn] >= cy2);
341
133
        n2p[qn] -= cy2;
342
133
      }
343
14.6k
    else
344
14.6k
      {
345
14.6k
        n2p[qn] = cy1 - cy2; /* & GMP_NUMB_MASK; */
346
347
14.6k
        quotient_too_large = (cy1 < cy2);
348
14.6k
        ++rn;
349
14.6k
      }
350
14.8k
    --in;
351
14.8k
        }
352
      /* True: partial remainder now is neutral, i.e., it is not shifted up.  */
353
354
15.6k
      tp = TMP_ALLOC_LIMBS (dn);
355
356
15.6k
      if (in < qn)
357
7.37k
        {
358
7.37k
    if (in == 0)
359
4.55k
      {
360
4.55k
        MPN_COPY (rp, n2p, rn);
361
4.55k
        ASSERT_ALWAYS (rn == dn);
362
4.55k
        goto foo;
363
4.55k
      }
364
2.82k
    mpn_mul (tp, qp, qn, dp, in);
365
2.82k
        }
366
8.23k
      else
367
8.23k
        mpn_mul (tp, dp, in, qp, qn);
368
369
11.0k
      cy = mpn_sub (n2p, n2p, rn, tp + in, qn);
370
11.0k
      MPN_COPY (rp + in, n2p, dn - in);
371
11.0k
      quotient_too_large |= cy;
372
11.0k
      cy = mpn_sub_n (rp, np, tp, in);
373
11.0k
      cy = mpn_sub_1 (rp + in, rp + in, rn, cy);
374
11.0k
      quotient_too_large |= cy;
375
15.6k
    foo:
376
15.6k
      if (quotient_too_large)
377
931
        {
378
931
    mpn_decr_u (qp, (mp_limb_t) 1);
379
931
    mpn_add_n (rp, rp, dp, dn);
380
931
        }
381
15.6k
    }
382
15.6k
  TMP_FREE;
383
15.6k
  return;
384
11.0k
      }
385
1.49M
    }
386
1.49M
}