Coverage Report

Created: 2023-02-22 06:39

/src/gmp-6.2.1/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
3.06k
{
46
3.06k
  mp_size_t lo, hi;
47
3.06k
  mp_limb_t cy, qh, ql;
48
49
3.06k
  lo = n >> 1;      /* floor(n/2) */
50
3.06k
  hi = n - lo;      /* ceil(n/2) */
51
52
3.06k
  if (BELOW_THRESHOLD (hi, DC_DIV_QR_THRESHOLD))
53
2.39k
    qh = mpn_sbpi1_div_qr (qp + lo, np + 2 * lo, 2 * hi, dp + lo, hi, dinv->inv32);
54
670
  else
55
670
    qh = mpn_dcpi1_div_qr_n (qp + lo, np + 2 * lo, dp + lo, hi, dinv, tp);
56
57
3.06k
  mpn_mul (tp, qp + lo, hi, dp, lo);
58
59
3.06k
  cy = mpn_sub_n (np + lo, np + lo, tp, n);
60
3.06k
  if (qh != 0)
61
0
    cy += mpn_sub_n (np + n, np + n, dp, lo);
62
63
4.14k
  while (cy != 0)
64
1.08k
    {
65
1.08k
      qh -= mpn_sub_1 (qp + lo, qp + lo, hi, 1);
66
1.08k
      cy -= mpn_add_n (np + lo, np + lo, dp, n);
67
1.08k
    }
68
69
3.06k
  if (BELOW_THRESHOLD (lo, DC_DIV_QR_THRESHOLD))
70
2.39k
    ql = mpn_sbpi1_div_qr (qp, np + hi, 2 * lo, dp + hi, lo, dinv->inv32);
71
669
  else
72
669
    ql = mpn_dcpi1_div_qr_n (qp, np + hi, dp + hi, lo, dinv, tp);
73
74
3.06k
  mpn_mul (tp, dp, hi, qp, lo);
75
76
3.06k
  cy = mpn_sub_n (np, np, tp, n);
77
3.06k
  if (ql != 0)
78
1
    cy += mpn_sub_n (np + lo, np + lo, dp, hi);
79
80
3.89k
  while (cy != 0)
81
834
    {
82
834
      mpn_sub_1 (qp, qp, lo, 1);
83
834
      cy -= mpn_add_n (np, np, dp, n);
84
834
    }
85
86
3.06k
  return qh;
87
3.06k
}
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.50k
{
95
1.50k
  mp_size_t qn;
96
1.50k
  mp_limb_t qh, cy;
97
1.50k
  mp_ptr tp;
98
1.50k
  TMP_DECL;
99
100
1.50k
  TMP_MARK;
101
102
1.50k
  ASSERT (dn >= 6);   /* to adhere to mpn_sbpi1_div_qr's limits */
103
1.50k
  ASSERT (nn - dn >= 3);  /* to adhere to mpn_sbpi1_div_qr's limits */
104
1.50k
  ASSERT (dp[dn-1] & GMP_NUMB_HIGHBIT);
105
106
1.50k
  tp = TMP_ALLOC_LIMBS (dn);
107
108
1.50k
  qn = nn - dn;
109
1.50k
  qp += qn;
110
1.50k
  np += nn;
111
1.50k
  dp += dn;
112
113
1.50k
  if (qn > dn)
114
1.30k
    {
115
      /* Reduce qn mod dn without division, optimizing small operations.  */
116
1.30k
      do
117
1.33k
  qn -= dn;
118
1.33k
      while (qn > dn);
119
120
1.30k
      qp -= qn;     /* point at low limb of next quotient block */
121
1.30k
      np -= qn;     /* point in the middle of partial remainder */
122
123
      /* Perform the typically smaller block first.  */
124
1.30k
      if (qn == 1)
125
90
  {
126
90
    mp_limb_t q, n2, n1, n0, d1, d0;
127
128
    /* Handle qh up front, for simplicity. */
129
90
    qh = mpn_cmp (np - dn + 1, dp - dn, dn) >= 0;
130
90
    if (qh)
131
90
      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
90
    n2 = np[0];
136
90
    n1 = np[-1];
137
90
    n0 = np[-2];
138
90
    d1 = dp[-1];
139
90
    d0 = dp[-2];
140
141
90
    ASSERT (n2 < d1 || (n2 == d1 && n1 <= d0));
142
143
90
    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
90
    else
150
90
      {
151
90
        udiv_qr_3by2 (q, n1, n0, n2, n1, n0, d1, d0, dinv->inv32);
152
153
90
        if (dn > 2)
154
90
    {
155
90
      mp_limb_t cy, cy1;
156
90
      cy = mpn_submul_1 (np - dn, dp - dn, dn - 2, q);
157
158
90
      cy1 = n0 < cy;
159
90
      n0 = (n0 - cy) & GMP_NUMB_MASK;
160
90
      cy = n1 < cy1;
161
90
      n1 = (n1 - cy1) & GMP_NUMB_MASK;
162
90
      np[-2] = n0;
163
164
90
      if (UNLIKELY (cy != 0))
165
0
        {
166
0
          n1 += d1 + mpn_add_n (np - dn, np - dn, dp - dn, dn - 1);
167
0
          qh -= (q == 0);
168
0
          q = (q - 1) & GMP_NUMB_MASK;
169
0
        }
170
90
    }
171
0
        else
172
0
    np[-2] = n0;
173
174
90
        np[-1] = n1;
175
90
      }
176
90
    qp[0] = q;
177
90
  }
178
1.21k
      else
179
1.21k
  {
180
    /* Do a 2qn / qn division */
181
1.21k
    if (qn == 2)
182
15
      qh = mpn_divrem_2 (qp, 0L, np - 2, 4, dp - 2); /* FIXME: obsolete function. Use 5/3 division? */
183
1.19k
    else if (BELOW_THRESHOLD (qn, DC_DIV_QR_THRESHOLD))
184
1.00k
      qh = mpn_sbpi1_div_qr (qp, np - qn, 2 * qn, dp - qn, qn, dinv->inv32);
185
193
    else
186
193
      qh = mpn_dcpi1_div_qr_n (qp, np - qn, dp - qn, qn, dinv, tp);
187
188
1.21k
    if (qn != dn)
189
1.21k
      {
190
1.21k
        if (qn > dn - qn)
191
64
    mpn_mul (tp, qp, qn, dp - dn, dn - qn);
192
1.14k
        else
193
1.14k
    mpn_mul (tp, dp - dn, dn - qn, qp, qn);
194
195
1.21k
        cy = mpn_sub_n (np - dn, np - dn, tp, dn);
196
1.21k
        if (qh != 0)
197
0
    cy += mpn_sub_n (np - dn + qn, np - dn + qn, dp - dn, dn - qn);
198
199
1.25k
        while (cy != 0)
200
48
    {
201
48
      qh -= mpn_sub_1 (qp, qp, qn, 1);
202
48
      cy -= mpn_add_n (np - dn, np - dn, dp - dn, dn);
203
48
    }
204
1.21k
      }
205
1.21k
  }
206
207
1.30k
      qn = nn - dn - qn;
208
1.30k
      do
209
1.33k
  {
210
1.33k
    qp -= dn;
211
1.33k
    np -= dn;
212
1.33k
    mpn_dcpi1_div_qr_n (qp, np - dn, dp - dn, dn, dinv, tp);
213
1.33k
    qn -= dn;
214
1.33k
  }
215
1.33k
      while (qn > 0);
216
1.30k
    }
217
199
  else
218
199
    {
219
199
      qp -= qn;     /* point at low limb of next quotient block */
220
199
      np -= qn;     /* point in the middle of partial remainder */
221
222
199
      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
199
      else
225
199
  qh = mpn_dcpi1_div_qr_n (qp, np - qn, dp - qn, qn, dinv, tp);
226
227
199
      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
199
    }
245
246
1.50k
  TMP_FREE;
247
1.50k
  return qh;
248
1.50k
}