Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpn/sbpi1_divappr_q.c
Line
Count
Source (jump to first uncovered line)
1
/* mpn_sbpi1_divappr_q -- Schoolbook division using the Möller-Granlund 3/2
2
   division algorithm, returning approximate quotient.  The quotient returned
3
   is either correct, or one too large.
4
5
   Contributed to the GNU project by Torbjorn Granlund.
6
7
   THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
8
   SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
9
   GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
10
11
Copyright 2007, 2009 Free Software Foundation, Inc.
12
13
This file is part of the GNU MP Library.
14
15
The GNU MP Library is free software; you can redistribute it and/or modify
16
it under the terms of either:
17
18
  * the GNU Lesser General Public License as published by the Free
19
    Software Foundation; either version 3 of the License, or (at your
20
    option) any later version.
21
22
or
23
24
  * the GNU General Public License as published by the Free Software
25
    Foundation; either version 2 of the License, or (at your option) any
26
    later version.
27
28
or both in parallel, as here.
29
30
The GNU MP Library is distributed in the hope that it will be useful, but
31
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
32
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
33
for more details.
34
35
You should have received copies of the GNU General Public License and the
36
GNU Lesser General Public License along with the GNU MP Library.  If not,
37
see https://www.gnu.org/licenses/.  */
38
39
40
#include "gmp-impl.h"
41
#include "longlong.h"
42
43
mp_limb_t
44
mpn_sbpi1_divappr_q (mp_ptr qp,
45
         mp_ptr np, mp_size_t nn,
46
         mp_srcptr dp, mp_size_t dn,
47
         mp_limb_t dinv)
48
250
{
49
250
  mp_limb_t qh;
50
250
  mp_size_t qn, i;
51
250
  mp_limb_t n1, n0;
52
250
  mp_limb_t d1, d0;
53
250
  mp_limb_t cy, cy1;
54
250
  mp_limb_t q;
55
250
  mp_limb_t flag;
56
57
250
  ASSERT (dn > 2);
58
250
  ASSERT (nn >= dn);
59
250
  ASSERT ((dp[dn-1] & GMP_NUMB_HIGHBIT) != 0);
60
61
250
  np += nn;
62
63
250
  qn = nn - dn;
64
250
  if (qn + 1 < dn)
65
0
    {
66
0
      dp += dn - (qn + 1);
67
0
      dn = qn + 1;
68
0
    }
69
70
250
  qh = mpn_cmp (np - dn, dp, dn) >= 0;
71
250
  if (qh != 0)
72
32
    mpn_sub_n (np - dn, np - dn, dp, dn);
73
74
250
  qp += qn;
75
76
250
  dn -= 2;      /* offset dn by 2 for main division loops,
77
           saving two iterations in mpn_submul_1.  */
78
250
  d1 = dp[dn + 1];
79
250
  d0 = dp[dn + 0];
80
81
250
  np -= 2;
82
83
250
  n1 = np[1];
84
85
407
  for (i = qn - (dn + 2); i >= 0; i--)
86
157
    {
87
157
      np--;
88
157
      if (UNLIKELY (n1 == d1) && np[1] == d0)
89
0
  {
90
0
    q = GMP_NUMB_MASK;
91
0
    mpn_submul_1 (np - dn, dp, dn + 2, q);
92
0
    n1 = np[1];   /* update n1, last loop's value will now be invalid */
93
0
  }
94
157
      else
95
157
  {
96
157
    udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
97
98
157
    cy = mpn_submul_1 (np - dn, dp, dn, q);
99
100
157
    cy1 = n0 < cy;
101
157
    n0 = (n0 - cy) & GMP_NUMB_MASK;
102
157
    cy = n1 < cy1;
103
157
    n1 -= cy1;
104
157
    np[0] = n0;
105
106
157
    if (UNLIKELY (cy != 0))
107
1
      {
108
1
        n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1);
109
1
        q--;
110
1
      }
111
157
  }
112
113
157
      *--qp = q;
114
157
    }
115
116
250
  flag = ~CNST_LIMB(0);
117
118
250
  if (dn >= 0)
119
250
    {
120
3.03k
      for (i = dn; i > 0; i--)
121
2.78k
  {
122
2.78k
    np--;
123
2.78k
    if (UNLIKELY (n1 >= (d1 & flag)))
124
19
      {
125
19
        q = GMP_NUMB_MASK;
126
19
        cy = mpn_submul_1 (np - dn, dp, dn + 2, q);
127
128
19
        if (UNLIKELY (n1 != cy))
129
0
    {
130
0
      if (n1 < (cy & flag))
131
0
        {
132
0
          q--;
133
0
          mpn_add_n (np - dn, np - dn, dp, dn + 2);
134
0
        }
135
0
      else
136
0
        flag = 0;
137
0
    }
138
19
        n1 = np[1];
139
19
      }
140
2.76k
    else
141
2.76k
      {
142
2.76k
        udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
143
144
2.76k
        cy = mpn_submul_1 (np - dn, dp, dn, q);
145
146
2.76k
        cy1 = n0 < cy;
147
2.76k
        n0 = (n0 - cy) & GMP_NUMB_MASK;
148
2.76k
        cy = n1 < cy1;
149
2.76k
        n1 -= cy1;
150
2.76k
        np[0] = n0;
151
152
2.76k
        if (UNLIKELY (cy != 0))
153
2
    {
154
2
      n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1);
155
2
      q--;
156
2
    }
157
2.76k
      }
158
159
2.78k
    *--qp = q;
160
161
    /* Truncate operands.  */
162
2.78k
    dn--;
163
2.78k
    dp++;
164
2.78k
  }
165
166
250
      np--;
167
250
      if (UNLIKELY (n1 >= (d1 & flag)))
168
1
  {
169
1
    q = GMP_NUMB_MASK;
170
1
    cy = mpn_submul_1 (np, dp, 2, q);
171
172
1
    if (UNLIKELY (n1 != cy))
173
0
      {
174
0
        if (n1 < (cy & flag))
175
0
    {
176
0
      q--;
177
0
      add_ssaaaa (np[1], np[0], np[1], np[0], dp[1], dp[0]);
178
0
    }
179
0
        else
180
0
    flag = 0;
181
0
      }
182
1
    n1 = np[1];
183
1
  }
184
249
      else
185
249
  {
186
249
    udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
187
188
249
    np[1] = n1;
189
249
    np[0] = n0;
190
249
  }
191
192
250
      *--qp = q;
193
250
    }
194
195
250
  ASSERT_ALWAYS (np[1] == n1);
196
197
250
  return qh;
198
250
}