Lines Matching defs:n*
1 /* mulmod_bnm1.c -- multiplication mod B^n-1.
25 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
33 mod B^rn - 1, and values are semi-normalised; zero is represented
34 as either 0 or B^n - 1. Needs a scratch of 2rn limbs at tp.
47 * be no overflow when adding in the carry. */
53 semi-normalised representation, computation is mod B^rn + 1. Needs
55 Output is normalised. */
77 * B^rn-1. This should not be a problem if mulmod_bnm1 is used to
78 * combine results and obtain a natural number when one knows in
80 * Moreover it should not be a problem if mulmod_bnm1 is used to
85 * Scratch need: rn + (need for recursive call OR rn + 4). This gives
87 * S(n) <= rn + MAX (rn + 4, S(n/2)) <= 2rn + 4
117 mp_size_t n;
121 n = rn >> 1;
123 /* We need at least an + bn >= n, to be able to fit one of the
129 ASSERT (an + bn > n);
131 /* Compute xm = a*b mod (B^n - 1), xp = a*b mod (B^n + 1)
134 x = -xp * B^n + (B^n + 1) * [ (xp + xm)/2 mod (B^n-1)]
138 #define a1 (ap + n)
140 #define b1 (bp + n)
142 #define xp tp /* 2n + 2 */
143 /* am1 maybe in {xp, n} */
144 /* bm1 maybe in {xp + n, n} */
145 #define sp1 (tp + 2*n + 2)
146 /* ap1 maybe in {sp1, n + 1} */
147 /* bp1 maybe in {sp1 + n + 1, n + 1} */
154 if (LIKELY (an > n))
157 cy = mpn_add (xp, a0, n, a1, an - n);
158 MPN_INCR_U (xp, n, cy);
159 anm = n;
160 if (LIKELY (bn > n))
162 bm1 = xp + n;
163 cy = mpn_add (xp + n, b0, n, b1, bn - n);
164 MPN_INCR_U (xp + n, n, cy);
165 bnm = n;
166 so = xp + 2*n;
170 so = xp + n;
184 mpn_mulmod_bnm1 (rp, n, am1, anm, bm1, bnm, so);
192 if (LIKELY (an > n)) {
194 cy = mpn_sub (sp1, a0, n, a1, an - n);
195 sp1[n] = 0;
196 MPN_INCR_U (sp1, n + 1, cy);
197 anp = n + ap1[n];
203 if (LIKELY (bn > n)) {
204 bp1 = sp1 + n + 1;
205 cy = mpn_sub (sp1 + n + 1, b0, n, b1, bn - n);
206 sp1[2*n+1] = 0;
207 MPN_INCR_U (sp1 + n + 1, n + 1, cy);
208 bnp = n + bp1[n];
214 if (BELOW_THRESHOLD (n, MUL_FFT_MODF_THRESHOLD))
219 k = mpn_fft_best_k (n, 0);
221 while (n & mask) {k--; mask >>=1;};
224 xp[n] = mpn_mul_fft (xp, n, ap1, anp, bp1, bnp, k);
227 ASSERT (anp + bnp <= 2*n+1);
228 ASSERT (anp + bnp > n);
231 anp = anp + bnp - n;
232 ASSERT (anp <= n || xp[2*n]==0);
233 anp-= anp > n;
234 cy = mpn_sub (xp, xp, n, xp + n, anp);
235 xp[n] = 0;
236 MPN_INCR_U (xp, n+1, cy);
239 mpn_bc_mulmod_bnp1 (xp, ap1, bp1, n, xp);
244 xm <- (xp + xm)/2 = (xp + xm)B^n/2 mod (B^n-1)
247 Assumes xp normalised mod (B^n+1).
249 The residue class [0] is represented by [B^n-1]; except when
255 cy = mpn_rsh1add_nc(rp, rp, xp, n, xp[n]); /* B^n = 1 */
258 /* next update of rp[n-1] will set cy = 1 only if rp[n-1]+=hi
259 overflows, i.e. a further increment will not overflow again. */
261 cy = xp[n] + mpn_rsh1add_n(rp, rp, xp, n); /* B^n = 1 */
264 /* cy = 1 only if xp[n] = 1 i.e. {xp,n} = ZERO, this implies that
268 add_ssaaaa(cy, rp[n-1], cy, rp[n-1], 0, hi);
270 cy += (hi & rp[n-1]) >> (GMP_NUMB_BITS-1);
271 rp[n-1] ^= hi;
275 cy = mpn_add_nc(rp, rp, xp, n, xp[n]);
277 cy = xp[n] + mpn_add_n(rp, rp, xp, n); /* xp[n] == 1 implies {xp,n} == ZERO */
280 mpn_rshift(rp, rp, n, 1);
285 ASSERT ((rp[n-1] & GMP_NUMB_HIGHBIT) == 0);
286 rp[n-1] |= hi;
287 /* ... rp[n-1] + cy can not overflow, the following INCR is correct. */
290 /* Next increment can not overflow, read the previous comments about cy. */
291 ASSERT ((cy == 0) || ((rp[n-1] & GMP_NUMB_HIGHBIT) == 0));
292 MPN_INCR_U(rp, n, cy);
295 ([(xp + xm)/2 mod (B^n-1)] - xp ) * B^n
302 reconstruction is zero, not B^{rn} - 1. Which is good,
305 cy = mpn_sub_n (rp + n, rp, xp, an + bn - n);
307 /* FIXME: This subtraction of the high parts is not really
308 necessary, we do it to get the carry out, and for sanity
310 cy = xp[n] + mpn_sub_nc (xp + an + bn - n, rp + an + bn - n,
311 xp + an + bn - n, rn - (an + bn), cy);
313 mpn_zero_p (xp + an + bn - n + 1, rn - 1 - (an + bn)));
315 ASSERT (cy == (xp + an + bn - n)[0]);
319 cy = xp[n] + mpn_sub_n (rp + n, rp, xp, n);
320 /* cy = 1 only if {xp,n+1} is not ZERO, i.e. {rp,n} is not ZERO.
321 DECR will affect _at most_ the lowest n limbs. */
322 MPN_DECR_U (rp, 2*n, cy);
334 mpn_mulmod_bnm1_next_size (mp_size_t n)
336 mp_size_t nh;
338 if (BELOW_THRESHOLD (n, MULMOD_BNM1_THRESHOLD))
339 return n;
340 if (BELOW_THRESHOLD (n, 4 * (MULMOD_BNM1_THRESHOLD - 1) + 1))
341 return (n + (2-1)) & (-2);
342 if (BELOW_THRESHOLD (n, 8 * (MULMOD_BNM1_THRESHOLD - 1) + 1))
343 return (n + (4-1)) & (-4);
345 nh = (n + 1) >> 1;
347 if (BELOW_THRESHOLD (nh, MUL_FFT_MODF_THRESHOLD))
348 return (n + (8-1)) & (-8);
350 return 2 * mpn_fft_next_size (nh, mpn_fft_best_k (nh, 0));