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 the GNU Lesser General Public License as published by
17 the Free Software Foundation; either version 3 of the License, or (at your
18 option) any later version.
19
20 The GNU MP Library is distributed in the hope that it will be useful, but
21 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
22 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
23 License for more details.
24
25 You should have received a copy of the GNU Lesser General Public License
26 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
27
28
29 #include "gmp.h"
30 #include "gmp-impl.h"
31 #include "longlong.h"
32
33 mp_limb_t
mpn_sbpi1_divappr_q(mp_ptr qp,mp_ptr np,mp_size_t nn,mp_srcptr dp,mp_size_t dn,mp_limb_t dinv)34 mpn_sbpi1_divappr_q (mp_ptr qp,
35 mp_ptr np, mp_size_t nn,
36 mp_srcptr dp, mp_size_t dn,
37 mp_limb_t dinv)
38 {
39 mp_limb_t qh;
40 mp_size_t qn, i;
41 mp_limb_t n1, n0;
42 mp_limb_t d1, d0;
43 mp_limb_t cy, cy1;
44 mp_limb_t q;
45 mp_limb_t flag;
46
47 ASSERT (dn > 2);
48 ASSERT (nn >= dn);
49 ASSERT ((dp[dn-1] & GMP_NUMB_HIGHBIT) != 0);
50
51 np += nn;
52
53 qn = nn - dn;
54 if (qn + 1 < dn)
55 {
56 dp += dn - (qn + 1);
57 dn = qn + 1;
58 }
59
60 qh = mpn_cmp (np - dn, dp, dn) >= 0;
61 if (qh != 0)
62 mpn_sub_n (np - dn, np - dn, dp, dn);
63
64 qp += qn;
65
66 dn -= 2; /* offset dn by 2 for main division loops,
67 saving two iterations in mpn_submul_1. */
68 d1 = dp[dn + 1];
69 d0 = dp[dn + 0];
70
71 np -= 2;
72
73 n1 = np[1];
74
75 for (i = qn - (dn + 2); i >= 0; i--)
76 {
77 np--;
78 if (UNLIKELY (n1 == d1) && np[1] == d0)
79 {
80 q = GMP_NUMB_MASK;
81 mpn_submul_1 (np - dn, dp, dn + 2, q);
82 n1 = np[1]; /* update n1, last loop's value will now be invalid */
83 }
84 else
85 {
86 udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
87
88 cy = mpn_submul_1 (np - dn, dp, dn, q);
89
90 cy1 = n0 < cy;
91 n0 = (n0 - cy) & GMP_NUMB_MASK;
92 cy = n1 < cy1;
93 n1 -= cy1;
94 np[0] = n0;
95
96 if (UNLIKELY (cy != 0))
97 {
98 n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1);
99 q--;
100 }
101 }
102
103 *--qp = q;
104 }
105
106 flag = ~CNST_LIMB(0);
107
108 if (dn >= 0)
109 {
110 for (i = dn; i > 0; i--)
111 {
112 np--;
113 if (UNLIKELY (n1 >= (d1 & flag)))
114 {
115 q = GMP_NUMB_MASK;
116 cy = mpn_submul_1 (np - dn, dp, dn + 2, q);
117
118 if (UNLIKELY (n1 != cy))
119 {
120 if (n1 < (cy & flag))
121 {
122 q--;
123 mpn_add_n (np - dn, np - dn, dp, dn + 2);
124 }
125 else
126 flag = 0;
127 }
128 n1 = np[1];
129 }
130 else
131 {
132 udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
133
134 cy = mpn_submul_1 (np - dn, dp, dn, q);
135
136 cy1 = n0 < cy;
137 n0 = (n0 - cy) & GMP_NUMB_MASK;
138 cy = n1 < cy1;
139 n1 -= cy1;
140 np[0] = n0;
141
142 if (UNLIKELY (cy != 0))
143 {
144 n1 += d1 + mpn_add_n (np - dn, np - dn, dp, dn + 1);
145 q--;
146 }
147 }
148
149 *--qp = q;
150
151 /* Truncate operands. */
152 dn--;
153 dp++;
154 }
155
156 np--;
157 if (UNLIKELY (n1 >= (d1 & flag)))
158 {
159 q = GMP_NUMB_MASK;
160 cy = mpn_submul_1 (np, dp, 2, q);
161
162 if (UNLIKELY (n1 != cy))
163 {
164 if (n1 < (cy & flag))
165 {
166 q--;
167 add_ssaaaa (np[1], np[0], np[1], np[0], dp[1], dp[0]);
168 }
169 else
170 flag = 0;
171 }
172 n1 = np[1];
173 }
174 else
175 {
176 udiv_qr_3by2 (q, n1, n0, n1, np[1], np[0], d1, d0, dinv);
177
178 np[1] = n1;
179 np[0] = n0;
180 }
181
182 *--qp = q;
183 }
184
185 ASSERT_ALWAYS (np[1] == n1);
186
187 return qh;
188 }
189