1 /* mpn_invertappr and helper functions. Compute I such that
2 floor((B^{2n}-1)/U - 1 <= I + B^n <= floor((B^{2n}-1)/U.
3
4 Contributed to the GNU project by Marco Bodrato.
5
6 The algorithm used here was inspired by ApproximateReciprocal from "Modern
7 Computer Arithmetic", by Richard P. Brent and Paul Zimmermann. Special
8 thanks to Paul Zimmermann for his very valuable suggestions on all the
9 theoretical aspects during the work on this code.
10
11 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY
12 SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
13 GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
14
15 Copyright (C) 2007, 2009, 2010 Free Software Foundation, Inc.
16
17 This file is part of the GNU MP Library.
18
19 The GNU MP Library is free software; you can redistribute it and/or modify
20 it under the terms of the GNU Lesser General Public License as published by
21 the Free Software Foundation; either version 3 of the License, or (at your
22 option) any later version.
23
24 The GNU MP Library is distributed in the hope that it will be useful, but
25 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
26 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
27 License for more details.
28
29 You should have received a copy of the GNU Lesser General Public License
30 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
31
32 /* FIXME: Remove NULL and TMP_*, as soon as all the callers properly
33 allocate and pass the scratch to the function. */
34 #include <stdlib.h> /* for NULL */
35
36 #include "gmp.h"
37 #include "gmp-impl.h"
38 #include "longlong.h"
39
40 /* FIXME: The iterative version splits the operand in two slighty unbalanced
41 parts, the use of log_2 (or counting the bits) underestimate the maximum
42 number of iterations. */
43
44 /* This is intended for constant THRESHOLDs only, where the compiler
45 can completely fold the result. */
46 #define LOG2C(n) \
47 (((n) >= 0x1) + ((n) >= 0x2) + ((n) >= 0x4) + ((n) >= 0x8) + \
48 ((n) >= 0x10) + ((n) >= 0x20) + ((n) >= 0x40) + ((n) >= 0x80) + \
49 ((n) >= 0x100) + ((n) >= 0x200) + ((n) >= 0x400) + ((n) >= 0x800) + \
50 ((n) >= 0x1000) + ((n) >= 0x2000) + ((n) >= 0x4000) + ((n) >= 0x8000))
51
52 #if TUNE_PROGRAM_BUILD
53 #define NPOWS \
54 ((sizeof(mp_size_t) > 6 ? 48 : 8*sizeof(mp_size_t)))
55 #define MAYBE_dcpi1_divappr 1
56 #else
57 #define NPOWS \
58 ((sizeof(mp_size_t) > 6 ? 48 : 8*sizeof(mp_size_t)) - LOG2C (INV_NEWTON_THRESHOLD))
59 #define MAYBE_dcpi1_divappr \
60 (INV_NEWTON_THRESHOLD < DC_DIVAPPR_Q_THRESHOLD)
61 #if (INV_NEWTON_THRESHOLD > INV_MULMOD_BNM1_THRESHOLD) && \
62 (INV_APPR_THRESHOLD > INV_MULMOD_BNM1_THRESHOLD)
63 #undef INV_MULMOD_BNM1_THRESHOLD
64 #define INV_MULMOD_BNM1_THRESHOLD 0 /* always when Newton */
65 #endif
66 #endif
67
68 /* All the three functions mpn{,_bc,_ni}_invertappr (ip, dp, n, scratch), take
69 the strictly normalised value {dp,n} (i.e., most significant bit must be set)
70 as an input, and compute {ip,n}: the approximate reciprocal of {dp,n}.
71
72 Let e = mpn*_invertappr (ip, dp, n, scratch) be the returned value; the
73 following conditions are satisfied by the output:
74 0 <= e <= 1;
75 {dp,n}*(B^n+{ip,n}) < B^{2n} <= {dp,n}*(B^n+{ip,n}+1+e) .
76 I.e. e=0 means that the result {ip,n} equals the one given by mpn_invert.
77 e=1 means that the result _may_ be one less than expected.
78
79 The _bc version returns e=1 most of the time.
80 The _ni version should return e=0 most of the time; only about 1% of
81 possible random input should give e=1.
82
83 When the strict result is needed, i.e., e=0 in the relation above:
84 {dp,n}*(B^n+{ip,n}) < B^{2n} <= {dp,n}*(B^n+{ip,n}+1) ;
85 the function mpn_invert (ip, dp, n, scratch) should be used instead. */
86
87 /* Maximum scratch needed by this branch (at tp): 3*n + 2 */
88 static mp_limb_t
mpn_bc_invertappr(mp_ptr ip,mp_srcptr dp,mp_size_t n,mp_ptr tp)89 mpn_bc_invertappr (mp_ptr ip, mp_srcptr dp, mp_size_t n, mp_ptr tp)
90 {
91 mp_ptr xp;
92
93 ASSERT (n > 0);
94 ASSERT (dp[n-1] & GMP_NUMB_HIGHBIT);
95 ASSERT (! MPN_OVERLAP_P (ip, n, dp, n));
96 ASSERT (! MPN_OVERLAP_P (ip, n, tp, mpn_invertappr_itch(n)));
97 ASSERT (! MPN_OVERLAP_P (dp, n, tp, mpn_invertappr_itch(n)));
98
99 /* Compute a base value of r limbs. */
100 if (n == 1)
101 invert_limb (*ip, *dp);
102 else {
103 mp_size_t i;
104 xp = tp + n + 2; /* 2 * n limbs */
105
106 for (i = n - 1; i >= 0; i--)
107 xp[i] = GMP_NUMB_MAX;
108 mpn_com (xp + n, dp, n);
109
110 /* Now xp contains B^2n - {dp,n}*B^n - 1 */
111
112 /* FIXME: if mpn_*pi1_divappr_q handles n==2, use it! */
113 if (n == 2) {
114 mpn_divrem_2 (ip, 0, xp, 4, dp);
115 } else {
116 gmp_pi1_t inv;
117 invert_pi1 (inv, dp[n-1], dp[n-2]);
118 if (! MAYBE_dcpi1_divappr
119 || BELOW_THRESHOLD (n, DC_DIVAPPR_Q_THRESHOLD))
120 mpn_sbpi1_divappr_q (ip, xp, 2 * n, dp, n, inv.inv32);
121 else
122 mpn_dcpi1_divappr_q (ip, xp, 2 * n, dp, n, &inv);
123 MPN_DECR_U(ip, n, 1);
124 return 1;
125 }
126 }
127 return 0;
128 }
129
130 /* mpn_ni_invertappr: computes the approximate reciprocal using Newton's
131 iterations (at least one).
132
133 Inspired by Algorithm "ApproximateReciprocal", published in "Modern Computer
134 Arithmetic" by Richard P. Brent and Paul Zimmermann, algorithm 3.5, page 121
135 in version 0.4 of the book.
136
137 Some adaptations were introduced, to allow product mod B^m-1 and return the
138 value e.
139
140 USE_MUL_N = 1 (default) introduces a correction in such a way that "the
141 value of B^{n+h}-T computed at step 8 cannot exceed B^n-1" (the book reads
142 "2B^n-1"). This correction should not require to modify the proof.
143
144 We use a wrapped product modulo B^m-1. NOTE: is there any normalisation
145 problem for the [0] class? It shouldn't: we compute 2*|A*X_h - B^{n+h}| <
146 B^m-1. We may get [0] if and only if we get AX_h = B^{n+h}. This can
147 happen only if A=B^{n}/2, but this implies X_h = B^{h}*2-1 i.e., AX_h =
148 B^{n+h} - A, then we get into the "negative" branch, where X_h is not
149 incremented (because A < B^n).
150
151 FIXME: the scratch for mulmod_bnm1 does not currently fit in the scratch, it
152 is allocated apart. */
153
154 #define USE_MUL_N 1
155
156 mp_limb_t
mpn_ni_invertappr(mp_ptr ip,mp_srcptr dp,mp_size_t n,mp_ptr scratch)157 mpn_ni_invertappr (mp_ptr ip, mp_srcptr dp, mp_size_t n, mp_ptr scratch)
158 {
159 mp_limb_t cy;
160 mp_ptr xp;
161 mp_size_t rn, mn;
162 mp_size_t sizes[NPOWS], *sizp;
163 mp_ptr tp;
164 TMP_DECL;
165 #define rp scratch
166
167 ASSERT (n > 2);
168 ASSERT (dp[n-1] & GMP_NUMB_HIGHBIT);
169 ASSERT (! MPN_OVERLAP_P (ip, n, dp, n));
170 ASSERT (! MPN_OVERLAP_P (ip, n, scratch, mpn_invertappr_itch(n)));
171 ASSERT (! MPN_OVERLAP_P (dp, n, scratch, mpn_invertappr_itch(n)));
172
173 /* Compute the computation precisions from highest to lowest, leaving the
174 base case size in 'rn'. */
175 sizp = sizes;
176 rn = n;
177 do {
178 *sizp = rn;
179 rn = ((rn) >> 1) + 1;
180 sizp ++;
181 } while (ABOVE_THRESHOLD (rn, INV_NEWTON_THRESHOLD));
182
183 /* We search the inverse of 0.{dp,n}, we compute it as 1.{ip,n} */
184 dp += n;
185 ip += n;
186
187 /* Compute a base value of rn limbs. */
188 mpn_bc_invertappr (ip - rn, dp - rn, rn, scratch);
189
190 TMP_MARK;
191
192 if (ABOVE_THRESHOLD (n, INV_MULMOD_BNM1_THRESHOLD))
193 {
194 mn = mpn_mulmod_bnm1_next_size (n + 1);
195 tp = TMP_ALLOC_LIMBS (mpn_mulmod_bnm1_itch (mn, n, (n >> 1) + 1));
196 }
197 /* Use Newton's iterations to get the desired precision.*/
198
199 /* define rp scratch; 2rn + 1 limbs <= 2(n>>1 + 1) + 1 <= n + 3 limbs */
200 /* Maximum scratch needed by this branch <= 3*n + 2 */
201 xp = scratch + n + 3; /* n + rn limbs */
202 while (1) {
203 mp_limb_t method;
204
205 n = *--sizp;
206 /*
207 v n v
208 +----+--+
209 ^ rn ^
210 */
211
212 /* Compute i_jd . */
213 if (BELOW_THRESHOLD (n, INV_MULMOD_BNM1_THRESHOLD)
214 || ((mn = mpn_mulmod_bnm1_next_size (n + 1)) > (n + rn))) {
215 /* FIXME: We do only need {xp,n+1}*/
216 mpn_mul (xp, dp - n, n, ip - rn, rn);
217 mpn_add_n (xp + rn, xp + rn, dp - n, n - rn + 1);
218 method = 1; /* Remember we used (truncated) product */
219 /* We computed cy.{xp,rn+n} <- 1.{ip,rn} * 0.{dp,n} */
220 } else { /* Use B^n-1 wraparound */
221 mpn_mulmod_bnm1 (xp, mn, dp - n, n, ip - rn, rn, tp);
222 /* We computed {xp,mn} <- {ip,rn} * {dp,n} mod (B^mn-1) */
223 /* We know that 2*|ip*dp + dp*B^rn - B^{rn+n}| < B^mn-1 */
224 /* Add dp*B^rn mod (B^mn-1) */
225 ASSERT (n >= mn - rn);
226 xp[mn] = 1 + mpn_add_n (xp + rn, xp + rn, dp - n, mn - rn);
227 cy = mpn_add_n (xp, xp, dp - (n - (mn - rn)), n - (mn - rn));
228 MPN_INCR_U (xp + n - (mn - rn), mn + 1 - n + (mn - rn), cy);
229 ASSERT (n + rn >= mn);
230 /* Subtract B^{rn+n} */
231 MPN_DECR_U (xp + rn + n - mn, 2*mn + 1 - rn - n, 1);
232 if (xp[mn])
233 MPN_INCR_U (xp, mn, xp[mn] - 1);
234 else
235 MPN_DECR_U (xp, mn, 1);
236 method = 0; /* Remember we are working Mod B^m-1 */
237 }
238
239 if (xp[n] < 2) { /* "positive" residue class */
240 cy = 1;
241 while (xp[n] || mpn_cmp (xp, dp - n, n)>0) {
242 xp[n] -= mpn_sub_n (xp, xp, dp - n, n);
243 cy ++;
244 }
245 MPN_DECR_U(ip - rn, rn, cy);
246 ASSERT (cy <= 4); /* at most 3 cycles for the while above */
247 ASSERT_NOCARRY (mpn_sub_n (xp, dp - n, xp, n));
248 ASSERT (xp[n] == 0);
249 } else { /* "negative" residue class */
250 mpn_com (xp, xp, n + 1);
251 MPN_INCR_U(xp, n + 1, method);
252 ASSERT (xp[n] <= 1);
253 #if USE_MUL_N
254 if (xp[n]) {
255 MPN_INCR_U(ip - rn, rn, 1);
256 ASSERT_CARRY (mpn_sub_n (xp, xp, dp - n, n));
257 }
258 #endif
259 }
260
261 /* Compute x_ju_j. FIXME:We need {rp+rn,rn}, mulhi? */
262 #if USE_MUL_N
263 mpn_mul_n (rp, xp + n - rn, ip - rn, rn);
264 #else
265 rp[2*rn] = 0;
266 mpn_mul (rp, xp + n - rn, rn + xp[n], ip - rn, rn);
267 #endif
268 /* We need _only_ the carry from the next addition */
269 /* Anyway 2rn-n <= 2... we don't need to optimise. */
270 cy = mpn_add_n (rp + rn, rp + rn, xp + n - rn, 2*rn - n);
271 cy = mpn_add_nc (ip - n, rp + 3*rn - n, xp + rn, n - rn, cy);
272 MPN_INCR_U (ip - rn, rn, cy + (1-USE_MUL_N)*(rp[2*rn] + xp[n]));
273 if (sizp == sizes) { /* Get out of the cycle */
274 /* Check for possible carry propagation from below. */
275 cy = rp[3*rn - n - 1] > GMP_NUMB_MAX - 7; /* Be conservative. */
276 /* cy = mpn_add_1 (rp + rn, rp + rn, 2*rn - n, 4); */
277 break;
278 }
279 rn = n;
280 }
281 TMP_FREE;
282
283 return cy;
284 #undef rp
285 }
286
287 mp_limb_t
mpn_invertappr(mp_ptr ip,mp_srcptr dp,mp_size_t n,mp_ptr scratch)288 mpn_invertappr (mp_ptr ip, mp_srcptr dp, mp_size_t n, mp_ptr scratch)
289 {
290 mp_limb_t res;
291 TMP_DECL;
292
293 TMP_MARK;
294
295 if (scratch == NULL)
296 scratch = TMP_ALLOC_LIMBS (mpn_invertappr_itch (n));
297
298 ASSERT (n > 0);
299 ASSERT (dp[n-1] & GMP_NUMB_HIGHBIT);
300 ASSERT (! MPN_OVERLAP_P (ip, n, dp, n));
301 ASSERT (! MPN_OVERLAP_P (ip, n, scratch, mpn_invertappr_itch(n)));
302 ASSERT (! MPN_OVERLAP_P (dp, n, scratch, mpn_invertappr_itch(n)));
303
304 if (BELOW_THRESHOLD (n, INV_NEWTON_THRESHOLD))
305 res = mpn_bc_invertappr (ip, dp, n, scratch);
306 else
307 res = mpn_ni_invertappr (ip, dp, n, scratch);
308
309 TMP_FREE;
310 return res;
311 }
312