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.
4 Contributed to the GNU project by Marco Bodrato.
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.
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.
15 Copyright (C) 2007, 2009, 2010 Free Software Foundation, Inc.
17 This file is part of the GNU MP Library.
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.
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.
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/. */
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 */
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. */
44 /* This is intended for constant THRESHOLDs only, where the compiler
45 can completely fold the result. */
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))
52 #if TUNE_PROGRAM_BUILD
54 ((sizeof(mp_size_t) > 6 ? 48 : 8*sizeof(mp_size_t)))
55 #define MAYBE_dcpi1_divappr 1
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 */
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}.
72 Let e = mpn*_invertappr (ip, dp, n, scratch) be the returned value; the
73 following conditions are satisfied by the output:
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.
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.
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. */
87 /* Maximum scratch needed by this branch (at tp): 3*n + 2 */
89 mpn_bc_invertappr (mp_ptr ip, mp_srcptr dp, mp_size_t n, mp_ptr tp)
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)));
99 /* Compute a base value of r limbs. */
101 invert_limb (*ip, *dp);
104 xp = tp + n + 2; /* 2 * n limbs */
106 for (i = n - 1; i >= 0; i--)
107 xp[i] = GMP_NUMB_MAX;
108 mpn_com (xp + n, dp, n);
110 /* Now xp contains B^2n - {dp,n}*B^n - 1 */
112 /* FIXME: if mpn_*pi1_divappr_q handles n==2, use it! */
114 mpn_divrem_2 (ip, 0, xp, 4, dp);
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);
122 mpn_dcpi1_divappr_q (ip, xp, 2 * n, dp, n, &inv);
123 MPN_DECR_U(ip, n, 1);
130 /* mpn_ni_invertappr: computes the approximate reciprocal using Newton's
131 iterations (at least one).
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.
137 Some adaptations were introduced, to allow product mod B^m-1 and return the
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.
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).
151 FIXME: the scratch for mulmod_bnm1 does not currently fit in the scratch, it
152 is allocated apart. */
157 mpn_ni_invertappr (mp_ptr ip, mp_srcptr dp, mp_size_t n, mp_ptr scratch)
162 mp_size_t sizes[NPOWS], *sizp;
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)));
173 /* Compute the computation precisions from highest to lowest, leaving the
174 base case size in 'rn'. */
179 rn = ((rn) >> 1) + 1;
181 } while (ABOVE_THRESHOLD (rn, INV_NEWTON_THRESHOLD));
183 /* We search the inverse of 0.{dp,n}, we compute it as 1.{ip,n} */
187 /* Compute a base value of rn limbs. */
188 mpn_bc_invertappr (ip - rn, dp - rn, rn, scratch);
192 if (ABOVE_THRESHOLD (n, INV_MULMOD_BNM1_THRESHOLD))
194 mn = mpn_mulmod_bnm1_next_size (n + 1);
195 tp = TMP_ALLOC_LIMBS (mpn_mulmod_bnm1_itch (mn, n, (n >> 1) + 1));
197 /* Use Newton's iterations to get the desired precision.*/
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 */
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);
233 MPN_INCR_U (xp, mn, xp[mn] - 1);
235 MPN_DECR_U (xp, mn, 1);
236 method = 0; /* Remember we are working Mod B^m-1 */
239 if (xp[n] < 2) { /* "positive" residue class */
241 while (xp[n] || mpn_cmp (xp, dp - n, n)>0) {
242 xp[n] -= mpn_sub_n (xp, xp, dp - n, n);
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));
249 } else { /* "negative" residue class */
250 mpn_com (xp, xp, n + 1);
251 MPN_INCR_U(xp, n + 1, method);
255 MPN_INCR_U(ip - rn, rn, 1);
256 ASSERT_CARRY (mpn_sub_n (xp, xp, dp - n, n));
261 /* Compute x_ju_j. FIXME:We need {rp+rn,rn}, mulhi? */
263 mpn_mul_n (rp, xp + n - rn, ip - rn, rn);
266 mpn_mul (rp, xp + n - rn, rn + xp[n], ip - rn, rn);
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); */
288 mpn_invertappr (mp_ptr ip, mp_srcptr dp, mp_size_t n, mp_ptr scratch)
296 scratch = TMP_ALLOC_LIMBS (mpn_invertappr_itch (n));
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)));
304 if (BELOW_THRESHOLD (n, INV_NEWTON_THRESHOLD))
305 res = mpn_bc_invertappr (ip, dp, n, scratch);
307 res = mpn_ni_invertappr (ip, dp, n, scratch);