Coverage Report

Created: 2024-11-21 06:47

/src/libgmp/mpn/dcpi1_div_qr.c
Line
Count
Source (jump to first uncovered line)
1
/* mpn_dcpi1_div_qr_n -- recursive divide-and-conquer division for arbitrary
2
   size operands.
3
4
   Contributed to the GNU project by Torbjorn Granlund.
5
6
   THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
7
   SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
8
   GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
9
10
Copyright 2006, 2007, 2009 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
#include "longlong.h"
40
41
42
mp_limb_t
43
mpn_dcpi1_div_qr_n (mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t n,
44
        gmp_pi1_t *dinv, mp_ptr tp)
45
1.94k
{
46
1.94k
  mp_size_t lo, hi;
47
1.94k
  mp_limb_t cy, qh, ql;
48
49
1.94k
  lo = n >> 1;      /* floor(n/2) */
50
1.94k
  hi = n - lo;      /* ceil(n/2) */
51
52
1.94k
  if (BELOW_THRESHOLD (hi, DC_DIV_QR_THRESHOLD))
53
1.94k
    qh = mpn_sbpi1_div_qr (qp + lo, np + 2 * lo, 2 * hi, dp + lo, hi, dinv->inv32);
54
0
  else
55
0
    qh = mpn_dcpi1_div_qr_n (qp + lo, np + 2 * lo, dp + lo, hi, dinv, tp);
56
57
1.94k
  mpn_mul (tp, qp + lo, hi, dp, lo);
58
59
1.94k
  cy = mpn_sub_n (np + lo, np + lo, tp, n);
60
1.94k
  if (qh != 0)
61
35
    cy += mpn_sub_n (np + n, np + n, dp, lo);
62
63
2.49k
  while (cy != 0)
64
551
    {
65
551
      qh -= mpn_sub_1 (qp + lo, qp + lo, hi, 1);
66
551
      cy -= mpn_add_n (np + lo, np + lo, dp, n);
67
551
    }
68
69
1.94k
  if (BELOW_THRESHOLD (lo, DC_DIV_QR_THRESHOLD))
70
1.94k
    ql = mpn_sbpi1_div_qr (qp, np + hi, 2 * lo, dp + hi, lo, dinv->inv32);
71
0
  else
72
0
    ql = mpn_dcpi1_div_qr_n (qp, np + hi, dp + hi, lo, dinv, tp);
73
74
1.94k
  mpn_mul (tp, dp, hi, qp, lo);
75
76
1.94k
  cy = mpn_sub_n (np, np, tp, n);
77
1.94k
  if (ql != 0)
78
2
    cy += mpn_sub_n (np + lo, np + lo, dp, hi);
79
80
2.68k
  while (cy != 0)
81
742
    {
82
742
      mpn_sub_1 (qp, qp, lo, 1);
83
742
      cy -= mpn_add_n (np, np, dp, n);
84
742
    }
85
86
1.94k
  return qh;
87
1.94k
}
88
89
mp_limb_t
90
mpn_dcpi1_div_qr (mp_ptr qp,
91
      mp_ptr np, mp_size_t nn,
92
      mp_srcptr dp, mp_size_t dn,
93
      gmp_pi1_t *dinv)
94
1.94k
{
95
1.94k
  mp_size_t qn;
96
1.94k
  mp_limb_t qh, cy;
97
1.94k
  mp_ptr tp;
98
1.94k
  TMP_DECL;
99
100
1.94k
  TMP_MARK;
101
102
1.94k
  ASSERT (dn >= 6);   /* to adhere to mpn_sbpi1_div_qr's limits */
103
1.94k
  ASSERT (nn - dn >= 3);  /* to adhere to mpn_sbpi1_div_qr's limits */
104
1.94k
  ASSERT (dp[dn-1] & GMP_NUMB_HIGHBIT);
105
106
1.94k
  tp = TMP_ALLOC_LIMBS (dn);
107
108
1.94k
  qn = nn - dn;
109
1.94k
  qp += qn;
110
1.94k
  np += nn;
111
1.94k
  dp += dn;
112
113
1.94k
  if (qn > dn)
114
1.19k
    {
115
      /* Reduce qn mod dn without division, optimizing small operations.  */
116
1.19k
      do
117
1.19k
  qn -= dn;
118
1.19k
      while (qn > dn);
119
120
1.19k
      qp -= qn;     /* point at low limb of next quotient block */
121
1.19k
      np -= qn;     /* point in the middle of partial remainder */
122
123
      /* Perform the typically smaller block first.  */
124
1.19k
      if (qn == 1)
125
496
  {
126
496
    mp_limb_t q, n2, n1, n0, d1, d0;
127
128
    /* Handle qh up front, for simplicity. */
129
496
    qh = mpn_cmp (np - dn + 1, dp - dn, dn) >= 0;
130
496
    if (qh)
131
0
      ASSERT_NOCARRY (mpn_sub_n (np - dn + 1, np - dn + 1, dp - dn, dn));
132
133
    /* A single iteration of schoolbook: One 3/2 division,
134
       followed by the bignum update and adjustment. */
135
496
    n2 = np[0];
136
496
    n1 = np[-1];
137
496
    n0 = np[-2];
138
496
    d1 = dp[-1];
139
496
    d0 = dp[-2];
140
141
496
    ASSERT (n2 < d1 || (n2 == d1 && n1 <= d0));
142
143
496
    if (UNLIKELY (n2 == d1) && n1 == d0)
144
0
      {
145
0
        q = GMP_NUMB_MASK;
146
0
        cy = mpn_submul_1 (np - dn, dp - dn, dn, q);
147
0
        ASSERT (cy == n2);
148
0
      }
149
496
    else
150
496
      {
151
496
        udiv_qr_3by2 (q, n1, n0, n2, n1, n0, d1, d0, dinv->inv32);
152
153
496
        if (dn > 2)
154
496
    {
155
496
      mp_limb_t cy, cy1;
156
496
      cy = mpn_submul_1 (np - dn, dp - dn, dn - 2, q);
157
158
496
      cy1 = n0 < cy;
159
496
      n0 = (n0 - cy) & GMP_NUMB_MASK;
160
496
      cy = n1 < cy1;
161
496
      n1 = (n1 - cy1) & GMP_NUMB_MASK;
162
496
      np[-2] = n0;
163
164
496
      if (UNLIKELY (cy != 0))
165
71
        {
166
71
          n1 += d1 + mpn_add_n (np - dn, np - dn, dp - dn, dn - 1);
167
71
          qh -= (q == 0);
168
71
          q = (q - 1) & GMP_NUMB_MASK;
169
71
        }
170
496
    }
171
0
        else
172
0
    np[-2] = n0;
173
174
496
        np[-1] = n1;
175
496
      }
176
496
    qp[0] = q;
177
496
  }
178
702
      else
179
702
  {
180
    /* Do a 2qn / qn division */
181
702
    if (qn == 2)
182
69
      qh = mpn_divrem_2 (qp, 0L, np - 2, 4, dp - 2); /* FIXME: obsolete function. Use 5/3 division? */
183
633
    else if (BELOW_THRESHOLD (qn, DC_DIV_QR_THRESHOLD))
184
633
      qh = mpn_sbpi1_div_qr (qp, np - qn, 2 * qn, dp - qn, qn, dinv->inv32);
185
0
    else
186
0
      qh = mpn_dcpi1_div_qr_n (qp, np - qn, dp - qn, qn, dinv, tp);
187
188
702
    if (qn != dn)
189
702
      {
190
702
        if (qn > dn - qn)
191
34
    mpn_mul (tp, qp, qn, dp - dn, dn - qn);
192
668
        else
193
668
    mpn_mul (tp, dp - dn, dn - qn, qp, qn);
194
195
702
        cy = mpn_sub_n (np - dn, np - dn, tp, dn);
196
702
        if (qh != 0)
197
0
    cy += mpn_sub_n (np - dn + qn, np - dn + qn, dp - dn, dn - qn);
198
199
892
        while (cy != 0)
200
190
    {
201
190
      qh -= mpn_sub_1 (qp, qp, qn, 1);
202
190
      cy -= mpn_add_n (np - dn, np - dn, dp - dn, dn);
203
190
    }
204
702
      }
205
702
  }
206
207
1.19k
      qn = nn - dn - qn;
208
1.19k
      do
209
1.19k
  {
210
1.19k
    qp -= dn;
211
1.19k
    np -= dn;
212
1.19k
    mpn_dcpi1_div_qr_n (qp, np - dn, dp - dn, dn, dinv, tp);
213
1.19k
    qn -= dn;
214
1.19k
  }
215
1.19k
      while (qn > 0);
216
1.19k
    }
217
745
  else
218
745
    {
219
745
      qp -= qn;     /* point at low limb of next quotient block */
220
745
      np -= qn;     /* point in the middle of partial remainder */
221
222
745
      if (BELOW_THRESHOLD (qn, DC_DIV_QR_THRESHOLD))
223
0
  qh = mpn_sbpi1_div_qr (qp, np - qn, 2 * qn, dp - qn, qn, dinv->inv32);
224
745
      else
225
745
  qh = mpn_dcpi1_div_qr_n (qp, np - qn, dp - qn, qn, dinv, tp);
226
227
745
      if (qn != dn)
228
0
  {
229
0
    if (qn > dn - qn)
230
0
      mpn_mul (tp, qp, qn, dp - dn, dn - qn);
231
0
    else
232
0
      mpn_mul (tp, dp - dn, dn - qn, qp, qn);
233
234
0
    cy = mpn_sub_n (np - dn, np - dn, tp, dn);
235
0
    if (qh != 0)
236
0
      cy += mpn_sub_n (np - dn + qn, np - dn + qn, dp - dn, dn - qn);
237
238
0
    while (cy != 0)
239
0
      {
240
0
        qh -= mpn_sub_1 (qp, qp, qn, 1);
241
0
        cy -= mpn_add_n (np - dn, np - dn, dp - dn, dn);
242
0
      }
243
0
  }
244
745
    }
245
246
1.94k
  TMP_FREE;
247
1.94k
  return qh;
248
1.94k
}