1 /* Interpolaton for the algorithm Toom-Cook 8.5-way.
2
3 Contributed to the GNU project by Marco Bodrato.
4
5 THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY
6 SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
7 GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
8
9 Copyright 2009, 2010, 2012 Free Software Foundation, Inc.
10
11 This file is part of the GNU MP Library.
12
13 The GNU MP Library is free software; you can redistribute it and/or modify
14 it under the terms of the GNU Lesser General Public License as published by
15 the Free Software Foundation; either version 3 of the License, or (at your
16 option) any later version.
17
18 The GNU MP Library is distributed in the hope that it will be useful, but
19 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
20 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
21 License for more details.
22
23 You should have received a copy of the GNU Lesser General Public License
24 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
25
26
27 #include "gmp.h"
28 #include "gmp-impl.h"
29
30 #if GMP_NUMB_BITS < 29
31 #error Not implemented: Both sublsh_n(,,,28) should be corrected; r2 and r5 need one more LIMB.
32 #endif
33
34 #if GMP_NUMB_BITS < 28
35 #error Not implemented: divexact_by188513325 and _by182712915 will not work.
36 #endif
37
38
39 #if HAVE_NATIVE_mpn_sublsh_n
40 #define DO_mpn_sublsh_n(dst,src,n,s,ws) mpn_sublsh_n(dst,dst,src,n,s)
41 #else
42 static mp_limb_t
DO_mpn_sublsh_n(mp_ptr dst,mp_srcptr src,mp_size_t n,unsigned int s,mp_ptr ws)43 DO_mpn_sublsh_n(mp_ptr dst, mp_srcptr src, mp_size_t n, unsigned int s, mp_ptr ws)
44 {
45 #if USE_MUL_1 && 0
46 return mpn_submul_1(dst,src,n,CNST_LIMB(1) <<(s));
47 #else
48 mp_limb_t __cy;
49 __cy = mpn_lshift(ws,src,n,s);
50 return __cy + mpn_sub_n(dst,dst,ws,n);
51 #endif
52 }
53 #endif
54
55 #if HAVE_NATIVE_mpn_addlsh_n
56 #define DO_mpn_addlsh_n(dst,src,n,s,ws) mpn_addlsh_n(dst,dst,src,n,s)
57 #else
58 static mp_limb_t
DO_mpn_addlsh_n(mp_ptr dst,mp_srcptr src,mp_size_t n,unsigned int s,mp_ptr ws)59 DO_mpn_addlsh_n(mp_ptr dst, mp_srcptr src, mp_size_t n, unsigned int s, mp_ptr ws)
60 {
61 #if USE_MUL_1 && 0
62 return mpn_addmul_1(dst,src,n,CNST_LIMB(1) <<(s));
63 #else
64 mp_limb_t __cy;
65 __cy = mpn_lshift(ws,src,n,s);
66 return __cy + mpn_add_n(dst,dst,ws,n);
67 #endif
68 }
69 #endif
70
71 #if HAVE_NATIVE_mpn_subrsh
72 #define DO_mpn_subrsh(dst,nd,src,ns,s,ws) mpn_subrsh(dst,nd,src,ns,s)
73 #else
74 /* FIXME: This is not a correct definition, it assumes no carry */
75 #define DO_mpn_subrsh(dst,nd,src,ns,s,ws) \
76 do { \
77 mp_limb_t __cy; \
78 MPN_DECR_U (dst, nd, src[0] >> s); \
79 __cy = DO_mpn_sublsh_n (dst, src + 1, ns - 1, GMP_NUMB_BITS - s, ws); \
80 MPN_DECR_U (dst + ns - 1, nd - ns + 1, __cy); \
81 } while (0)
82 #endif
83
84
85 /* FIXME: tuneup should decide the best variant */
86 #ifndef AORSMUL_FASTER_AORS_AORSLSH
87 #define AORSMUL_FASTER_AORS_AORSLSH 1
88 #endif
89 #ifndef AORSMUL_FASTER_AORS_2AORSLSH
90 #define AORSMUL_FASTER_AORS_2AORSLSH 1
91 #endif
92 #ifndef AORSMUL_FASTER_2AORSLSH
93 #define AORSMUL_FASTER_2AORSLSH 1
94 #endif
95 #ifndef AORSMUL_FASTER_3AORSLSH
96 #define AORSMUL_FASTER_3AORSLSH 1
97 #endif
98
99 #if GMP_NUMB_BITS < 43
100 #define BIT_CORRECTION 1
101 #define CORRECTION_BITS GMP_NUMB_BITS
102 #else
103 #define BIT_CORRECTION 0
104 #define CORRECTION_BITS 0
105 #endif
106
107 #define BINVERT_9 \
108 ((((GMP_NUMB_MAX / 9) << (6 - GMP_NUMB_BITS % 6)) * 8 & GMP_NUMB_MAX) | 0x39)
109
110 #define BINVERT_255 \
111 (GMP_NUMB_MAX - ((GMP_NUMB_MAX / 255) << (8 - GMP_NUMB_BITS % 8)))
112
113 /* FIXME: find some more general expressions for inverses */
114 #if GMP_LIMB_BITS == 32
115 #define BINVERT_2835 (GMP_NUMB_MASK & CNST_LIMB(0x53E3771B))
116 #define BINVERT_42525 (GMP_NUMB_MASK & CNST_LIMB(0x9F314C35))
117 #define BINVERT_182712915 (GMP_NUMB_MASK & CNST_LIMB(0x550659DB))
118 #define BINVERT_188513325 (GMP_NUMB_MASK & CNST_LIMB(0xFBC333A5))
119 #define BINVERT_255x182712915L (GMP_NUMB_MASK & CNST_LIMB(0x6FC4CB25))
120 #define BINVERT_255x188513325L (GMP_NUMB_MASK & CNST_LIMB(0x6864275B))
121 #if GMP_NAIL_BITS == 0
122 #define BINVERT_255x182712915H CNST_LIMB(0x1B649A07)
123 #define BINVERT_255x188513325H CNST_LIMB(0x06DB993A)
124 #else /* GMP_NAIL_BITS != 0 */
125 #define BINVERT_255x182712915H \
126 (GMP_NUMB_MASK & CNST_LIMB((0x1B649A07<<GMP_NAIL_BITS) | (0x6FC4CB25>>GMP_NUMB_BITS)))
127 #define BINVERT_255x188513325H \
128 (GMP_NUMB_MASK & CNST_LIMB((0x06DB993A<<GMP_NAIL_BITS) | (0x6864275B>>GMP_NUMB_BITS)))
129 #endif
130 #else
131 #if GMP_LIMB_BITS == 64
132 #define BINVERT_2835 (GMP_NUMB_MASK & CNST_LIMB(0x938CC70553E3771B))
133 #define BINVERT_42525 (GMP_NUMB_MASK & CNST_LIMB(0xE7B40D449F314C35))
134 #define BINVERT_255x182712915 (GMP_NUMB_MASK & CNST_LIMB(0x1B649A076FC4CB25))
135 #define BINVERT_255x188513325 (GMP_NUMB_MASK & CNST_LIMB(0x06DB993A6864275B))
136 #endif
137 #endif
138
139 #ifndef mpn_divexact_by255
140 #if GMP_NUMB_BITS % 8 == 0
141 #define mpn_divexact_by255(dst,src,size) \
142 (255 & 1 * mpn_bdiv_dbm1 (dst, src, size, __GMP_CAST (mp_limb_t, GMP_NUMB_MASK / 255)))
143 #else
144 #if HAVE_NATIVE_mpn_pi1_bdiv_q_1
145 #define mpn_divexact_by255(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,CNST_LIMB(255),BINVERT_255,0)
146 #else
147 #define mpn_divexact_by255(dst,src,size) mpn_divexact_1(dst,src,size,CNST_LIMB(255))
148 #endif
149 #endif
150 #endif
151
152 #ifndef mpn_divexact_by255x4
153 #if HAVE_NATIVE_mpn_pi1_bdiv_q_1
154 #define mpn_divexact_by255x4(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,CNST_LIMB(255),BINVERT_255,2)
155 #else
156 #define mpn_divexact_by255x4(dst,src,size) mpn_divexact_1(dst,src,size,CNST_LIMB(255)<<2)
157 #endif
158 #endif
159
160 #ifndef mpn_divexact_by9x16
161 #if HAVE_NATIVE_mpn_pi1_bdiv_q_1
162 #define mpn_divexact_by9x16(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,CNST_LIMB(9),BINVERT_9,4)
163 #else
164 #define mpn_divexact_by9x16(dst,src,size) mpn_divexact_1(dst,src,size,CNST_LIMB(9)<<4)
165 #endif
166 #endif
167
168 #ifndef mpn_divexact_by42525x16
169 #if HAVE_NATIVE_mpn_pi1_bdiv_q_1 && defined(BINVERT_42525)
170 #define mpn_divexact_by42525x16(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,CNST_LIMB(42525),BINVERT_42525,4)
171 #else
172 #define mpn_divexact_by42525x16(dst,src,size) mpn_divexact_1(dst,src,size,CNST_LIMB(42525)<<4)
173 #endif
174 #endif
175
176 #ifndef mpn_divexact_by2835x64
177 #if HAVE_NATIVE_mpn_pi1_bdiv_q_1 && defined(BINVERT_2835)
178 #define mpn_divexact_by2835x64(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,CNST_LIMB(2835),BINVERT_2835,6)
179 #else
180 #define mpn_divexact_by2835x64(dst,src,size) mpn_divexact_1(dst,src,size,CNST_LIMB(2835)<<6)
181 #endif
182 #endif
183
184 #ifndef mpn_divexact_by255x182712915
185 #if GMP_NUMB_BITS < 36
186 #if HAVE_NATIVE_mpn_bdiv_q_2_pi2 && defined(BINVERT_255x182712915H)
187 /* FIXME: use mpn_bdiv_q_2_pi2 */
188 #endif
189 #if HAVE_NATIVE_mpn_pi1_bdiv_q_1 && defined(BINVERT_182712915)
190 #define mpn_divexact_by255x182712915(dst,src,size) \
191 do { \
192 mpn_pi1_bdiv_q_1(dst,src,size,CNST_LIMB(182712915),BINVERT_182712915,0); \
193 mpn_divexact_by255(dst,dst,size); \
194 } while(0)
195 #else
196 #define mpn_divexact_by255x182712915(dst,src,size) \
197 do { \
198 mpn_divexact_1(dst,src,size,CNST_LIMB(182712915)); \
199 mpn_divexact_by255(dst,dst,size); \
200 } while(0)
201 #endif
202 #else /* GMP_NUMB_BITS > 35 */
203 #if HAVE_NATIVE_mpn_pi1_bdiv_q_1 && defined(BINVERT_255x182712915)
204 #define mpn_divexact_by255x182712915(dst,src,size) \
205 mpn_pi1_bdiv_q_1(dst,src,size,255*CNST_LIMB(182712915),BINVERT_255x182712915,0)
206 #else
207 #define mpn_divexact_by255x182712915(dst,src,size) mpn_divexact_1(dst,src,size,255*CNST_LIMB(182712915))
208 #endif
209 #endif /* GMP_NUMB_BITS >?< 36 */
210 #endif
211
212 #ifndef mpn_divexact_by255x188513325
213 #if GMP_NUMB_BITS < 36
214 #if HAVE_NATIVE_mpn_bdiv_q_1_pi2 && defined(BINVERT_255x188513325H)
215 /* FIXME: use mpn_bdiv_q_1_pi2 */
216 #endif
217 #if HAVE_NATIVE_mpn_pi1_bdiv_q_1 && defined(BINVERT_188513325)
218 #define mpn_divexact_by255x188513325(dst,src,size) \
219 do { \
220 mpn_pi1_bdiv_q_1(dst,src,size,CNST_LIMB(188513325),BINVERT_188513325,0); \
221 mpn_divexact_by255(dst,dst,size); \
222 } while(0)
223 #else
224 #define mpn_divexact_by255x188513325(dst,src,size) \
225 do { \
226 mpn_divexact_1(dst,src,size,CNST_LIMB(188513325)); \
227 mpn_divexact_by255(dst,dst,size); \
228 } while(0)
229 #endif
230 #else /* GMP_NUMB_BITS > 35 */
231 #if HAVE_NATIVE_mpn_pi1_bdiv_q_1 && defined(BINVERT_255x188513325)
232 #define mpn_divexact_by255x188513325(dst,src,size) \
233 mpn_pi1_bdiv_q_1(dst,src,size,255*CNST_LIMB(188513325),BINVERT_255x188513325,0)
234 #else
235 #define mpn_divexact_by255x188513325(dst,src,size) mpn_divexact_1(dst,src,size,255*CNST_LIMB(188513325))
236 #endif
237 #endif /* GMP_NUMB_BITS >?< 36 */
238 #endif
239
240 /* Interpolation for Toom-8.5 (or Toom-8), using the evaluation
241 points: infinity(8.5 only), +-8, +-4, +-2, +-1, +-1/4, +-1/2,
242 +-1/8, 0. More precisely, we want to compute
243 f(2^(GMP_NUMB_BITS * n)) for a polynomial f of degree 15 (or
244 14), given the 16 (rsp. 15) values:
245
246 r0 = limit at infinity of f(x) / x^7,
247 r1 = f(8),f(-8),
248 r2 = f(4),f(-4),
249 r3 = f(2),f(-2),
250 r4 = f(1),f(-1),
251 r5 = f(1/4),f(-1/4),
252 r6 = f(1/2),f(-1/2),
253 r7 = f(1/8),f(-1/8),
254 r8 = f(0).
255
256 All couples of the form f(n),f(-n) must be already mixed with
257 toom_couple_handling(f(n),...,f(-n),...)
258
259 The result is stored in {pp, spt + 7*n (or 8*n)}.
260 At entry, r8 is stored at {pp, 2n},
261 r6 is stored at {pp + 3n, 3n + 1}.
262 r4 is stored at {pp + 7n, 3n + 1}.
263 r2 is stored at {pp +11n, 3n + 1}.
264 r0 is stored at {pp +15n, spt}.
265
266 The other values are 3n+1 limbs each (with most significant limbs small).
267
268 Negative intermediate results are stored two-complemented.
269 Inputs are destroyed.
270 */
271
272 void
mpn_toom_interpolate_16pts(mp_ptr pp,mp_ptr r1,mp_ptr r3,mp_ptr r5,mp_ptr r7,mp_size_t n,mp_size_t spt,int half,mp_ptr wsi)273 mpn_toom_interpolate_16pts (mp_ptr pp, mp_ptr r1, mp_ptr r3, mp_ptr r5, mp_ptr r7,
274 mp_size_t n, mp_size_t spt, int half, mp_ptr wsi)
275 {
276 mp_limb_t cy;
277 mp_size_t n3;
278 mp_size_t n3p1;
279 n3 = 3 * n;
280 n3p1 = n3 + 1;
281
282 #define r6 (pp + n3) /* 3n+1 */
283 #define r4 (pp + 7 * n) /* 3n+1 */
284 #define r2 (pp +11 * n) /* 3n+1 */
285 #define r0 (pp +15 * n) /* s+t <= 2*n */
286
287 ASSERT( spt <= 2 * n );
288 /******************************* interpolation *****************************/
289 if( half != 0) {
290 cy = mpn_sub_n (r4, r4, r0, spt);
291 MPN_DECR_U (r4 + spt, n3p1 - spt, cy);
292
293 cy = DO_mpn_sublsh_n (r3, r0, spt, 14, wsi);
294 MPN_DECR_U (r3 + spt, n3p1 - spt, cy);
295 DO_mpn_subrsh(r6, n3p1, r0, spt, 2, wsi);
296
297 cy = DO_mpn_sublsh_n (r2, r0, spt, 28, wsi);
298 MPN_DECR_U (r2 + spt, n3p1 - spt, cy);
299 DO_mpn_subrsh(r5, n3p1, r0, spt, 4, wsi);
300
301 cy = DO_mpn_sublsh_n (r1 + BIT_CORRECTION, r0, spt, 42 - CORRECTION_BITS, wsi);
302 #if BIT_CORRECTION
303 cy = mpn_sub_1 (r1 + spt + BIT_CORRECTION, r1 + spt + BIT_CORRECTION,
304 n3p1 - spt - BIT_CORRECTION, cy);
305 ASSERT (BIT_CORRECTION > 0 || cy == 0);
306 /* FIXME: assumes r7[n3p1] is writable (it is if r5 follows). */
307 cy = r7[n3p1];
308 r7[n3p1] = 0x80;
309 #else
310 MPN_DECR_U (r1 + spt + BIT_CORRECTION, n3p1 - spt - BIT_CORRECTION, cy);
311 #endif
312 DO_mpn_subrsh(r7, n3p1 + BIT_CORRECTION, r0, spt, 6, wsi);
313 #if BIT_CORRECTION
314 /* FIXME: assumes r7[n3p1] is writable. */
315 ASSERT ( BIT_CORRECTION > 0 || r7[n3p1] == 0x80 );
316 r7[n3p1] = cy;
317 #endif
318 };
319
320 r5[n3] -= DO_mpn_sublsh_n (r5 + n, pp, 2 * n, 28, wsi);
321 DO_mpn_subrsh(r2 + n, 2 * n + 1, pp, 2 * n, 4, wsi);
322
323 #if HAVE_NATIVE_mpn_add_n_sub_n
324 mpn_add_n_sub_n (r2, r5, r5, r2, n3p1);
325 #else
326 mpn_sub_n (wsi, r5, r2, n3p1); /* can be negative */
327 ASSERT_NOCARRY(mpn_add_n (r2, r2, r5, n3p1));
328 MP_PTR_SWAP(r5, wsi);
329 #endif
330
331 r6[n3] -= DO_mpn_sublsh_n (r6 + n, pp, 2 * n, 14, wsi);
332 DO_mpn_subrsh(r3 + n, 2 * n + 1, pp, 2 * n, 2, wsi);
333
334 #if HAVE_NATIVE_mpn_add_n_sub_n
335 mpn_add_n_sub_n (r3, r6, r6, r3, n3p1);
336 #else
337 ASSERT_NOCARRY(mpn_add_n (wsi, r3, r6, n3p1));
338 mpn_sub_n (r6, r6, r3, n3p1); /* can be negative */
339 MP_PTR_SWAP(r3, wsi);
340 #endif
341
342 cy = DO_mpn_sublsh_n (r7 + n + BIT_CORRECTION, pp, 2 * n, 42 - CORRECTION_BITS, wsi);
343 #if BIT_CORRECTION
344 MPN_DECR_U (r1 + n, 2 * n + 1, pp[0] >> 6);
345 cy = DO_mpn_sublsh_n (r1 + n, pp + 1, 2 * n - 1, GMP_NUMB_BITS - 6, wsi);
346 cy = mpn_sub_1(r1 + 3 * n - 1, r1 + 3 * n - 1, 2, cy);
347 ASSERT ( BIT_CORRECTION > 0 || cy != 0 );
348 #else
349 r7[n3] -= cy;
350 DO_mpn_subrsh(r1 + n, 2 * n + 1, pp, 2 * n, 6, wsi);
351 #endif
352
353 #if HAVE_NATIVE_mpn_add_n_sub_n
354 mpn_add_n_sub_n (r1, r7, r7, r1, n3p1);
355 #else
356 mpn_sub_n (wsi, r7, r1, n3p1); /* can be negative */
357 mpn_add_n (r1, r1, r7, n3p1); /* if BIT_CORRECTION != 0, can give a carry. */
358 MP_PTR_SWAP(r7, wsi);
359 #endif
360
361 r4[n3] -= mpn_sub_n (r4+n, r4+n, pp, 2 * n);
362
363 #if AORSMUL_FASTER_2AORSLSH
364 mpn_submul_1 (r5, r6, n3p1, 1028); /* can be negative */
365 #else
366 DO_mpn_sublsh_n (r5, r6, n3p1, 2, wsi); /* can be negative */
367 DO_mpn_sublsh_n (r5, r6, n3p1,10, wsi); /* can be negative */
368 #endif
369
370 mpn_submul_1 (r7, r5, n3p1, 1300); /* can be negative */
371 #if AORSMUL_FASTER_3AORSLSH
372 mpn_submul_1 (r7, r6, n3p1, 1052688); /* can be negative */
373 #else
374 DO_mpn_sublsh_n (r7, r6, n3p1, 4, wsi); /* can be negative */
375 DO_mpn_sublsh_n (r7, r6, n3p1,12, wsi); /* can be negative */
376 DO_mpn_sublsh_n (r7, r6, n3p1,20, wsi); /* can be negative */
377 #endif
378 mpn_divexact_by255x188513325(r7, r7, n3p1);
379
380 mpn_submul_1 (r5, r7, n3p1, 12567555); /* can be negative */
381 /* A division by 2835x64 followsi. Warning: the operand can be negative! */
382 mpn_divexact_by2835x64(r5, r5, n3p1);
383 if ((r5[n3] & (GMP_NUMB_MAX << (GMP_NUMB_BITS-7))) != 0)
384 r5[n3] |= (GMP_NUMB_MAX << (GMP_NUMB_BITS-6));
385
386 #if AORSMUL_FASTER_AORS_AORSLSH
387 mpn_submul_1 (r6, r7, n3p1, 4095); /* can be negative */
388 #else
389 mpn_add_n (r6, r6, r7, n3p1); /* can give a carry */
390 DO_mpn_sublsh_n (r6, r7, n3p1, 12, wsi); /* can be negative */
391 #endif
392 #if AORSMUL_FASTER_2AORSLSH
393 mpn_addmul_1 (r6, r5, n3p1, 240); /* can be negative */
394 #else
395 DO_mpn_addlsh_n (r6, r5, n3p1, 8, wsi); /* can give a carry */
396 DO_mpn_sublsh_n (r6, r5, n3p1, 4, wsi); /* can be negative */
397 #endif
398 /* A division by 255x4 followsi. Warning: the operand can be negative! */
399 mpn_divexact_by255x4(r6, r6, n3p1);
400 if ((r6[n3] & (GMP_NUMB_MAX << (GMP_NUMB_BITS-3))) != 0)
401 r6[n3] |= (GMP_NUMB_MAX << (GMP_NUMB_BITS-2));
402
403 ASSERT_NOCARRY(DO_mpn_sublsh_n (r3, r4, n3p1, 7, wsi));
404
405 ASSERT_NOCARRY(DO_mpn_sublsh_n (r2, r4, n3p1, 13, wsi));
406 ASSERT_NOCARRY(mpn_submul_1 (r2, r3, n3p1, 400));
407
408 /* If GMP_NUMB_BITS < 42 next operations on r1 can give a carry!*/
409 DO_mpn_sublsh_n (r1, r4, n3p1, 19, wsi);
410 mpn_submul_1 (r1, r2, n3p1, 1428);
411 mpn_submul_1 (r1, r3, n3p1, 112896);
412 mpn_divexact_by255x182712915(r1, r1, n3p1);
413
414 ASSERT_NOCARRY(mpn_submul_1 (r2, r1, n3p1, 15181425));
415 mpn_divexact_by42525x16(r2, r2, n3p1);
416
417 #if AORSMUL_FASTER_AORS_2AORSLSH
418 ASSERT_NOCARRY(mpn_submul_1 (r3, r1, n3p1, 3969));
419 #else
420 ASSERT_NOCARRY(mpn_sub_n (r3, r3, r1, n3p1));
421 ASSERT_NOCARRY(DO_mpn_addlsh_n (r3, r1, n3p1, 7, wsi));
422 ASSERT_NOCARRY(DO_mpn_sublsh_n (r3, r1, n3p1, 12, wsi));
423 #endif
424 ASSERT_NOCARRY(mpn_submul_1 (r3, r2, n3p1, 900));
425 mpn_divexact_by9x16(r3, r3, n3p1);
426
427 ASSERT_NOCARRY(mpn_sub_n (r4, r4, r1, n3p1));
428 ASSERT_NOCARRY(mpn_sub_n (r4, r4, r3, n3p1));
429 ASSERT_NOCARRY(mpn_sub_n (r4, r4, r2, n3p1));
430
431 mpn_add_n (r6, r2, r6, n3p1);
432 ASSERT_NOCARRY(mpn_rshift(r6, r6, n3p1, 1));
433 ASSERT_NOCARRY(mpn_sub_n (r2, r2, r6, n3p1));
434
435 mpn_sub_n (r5, r3, r5, n3p1);
436 ASSERT_NOCARRY(mpn_rshift(r5, r5, n3p1, 1));
437 ASSERT_NOCARRY(mpn_sub_n (r3, r3, r5, n3p1));
438
439 mpn_add_n (r7, r1, r7, n3p1);
440 ASSERT_NOCARRY(mpn_rshift(r7, r7, n3p1, 1));
441 ASSERT_NOCARRY(mpn_sub_n (r1, r1, r7, n3p1));
442
443 /* last interpolation steps... */
444 /* ... could be mixed with recomposition
445 ||H-r7|M-r7|L-r7| ||H-r5|M-r5|L-r5|
446 */
447
448 /***************************** recomposition *******************************/
449 /*
450 pp[] prior to operations:
451 |M r0|L r0|___||H r2|M r2|L r2|___||H r4|M r4|L r4|___||H r6|M r6|L r6|____|H_r8|L r8|pp
452
453 summation scheme for remaining operations:
454 |__16|n_15|n_14|n_13|n_12|n_11|n_10|n__9|n__8|n__7|n__6|n__5|n__4|n__3|n__2|n___|n___|pp
455 |M r0|L r0|___||H r2|M r2|L r2|___||H r4|M r4|L r4|___||H r6|M r6|L r6|____|H_r8|L r8|pp
456 ||H r1|M r1|L r1| ||H r3|M r3|L r3| ||H_r5|M_r5|L_r5| ||H r7|M r7|L r7|
457 */
458
459 cy = mpn_add_n (pp + n, pp + n, r7, n);
460 cy = mpn_add_1 (pp + 2 * n, r7 + n, n, cy);
461 #if HAVE_NATIVE_mpn_add_nc
462 cy = r7[n3] + mpn_add_nc(pp + n3, pp + n3, r7 + 2 * n, n, cy);
463 #else
464 MPN_INCR_U (r7 + 2 * n, n + 1, cy);
465 cy = r7[n3] + mpn_add_n (pp + n3, pp + n3, r7 + 2 * n, n);
466 #endif
467 MPN_INCR_U (pp + 4 * n, 2 * n + 1, cy);
468
469 pp[2 * n3]+= mpn_add_n (pp + 5 * n, pp + 5 * n, r5, n);
470 cy = mpn_add_1 (pp + 2 * n3, r5 + n, n, pp[2 * n3]);
471 #if HAVE_NATIVE_mpn_add_nc
472 cy = r5[n3] + mpn_add_nc(pp + 7 * n, pp + 7 * n, r5 + 2 * n, n, cy);
473 #else
474 MPN_INCR_U (r5 + 2 * n, n + 1, cy);
475 cy = r5[n3] + mpn_add_n (pp + 7 * n, pp + 7 * n, r5 + 2 * n, n);
476 #endif
477 MPN_INCR_U (pp + 8 * n, 2 * n + 1, cy);
478
479 pp[10 * n]+= mpn_add_n (pp + 9 * n, pp + 9 * n, r3, n);
480 cy = mpn_add_1 (pp + 10 * n, r3 + n, n, pp[10 * n]);
481 #if HAVE_NATIVE_mpn_add_nc
482 cy = r3[n3] + mpn_add_nc(pp +11 * n, pp +11 * n, r3 + 2 * n, n, cy);
483 #else
484 MPN_INCR_U (r3 + 2 * n, n + 1, cy);
485 cy = r3[n3] + mpn_add_n (pp +11 * n, pp +11 * n, r3 + 2 * n, n);
486 #endif
487 MPN_INCR_U (pp +12 * n, 2 * n + 1, cy);
488
489 pp[14 * n]+=mpn_add_n (pp +13 * n, pp +13 * n, r1, n);
490 if ( half ) {
491 cy = mpn_add_1 (pp + 14 * n, r1 + n, n, pp[14 * n]);
492 #if HAVE_NATIVE_mpn_add_nc
493 if(LIKELY(spt > n)) {
494 cy = r1[n3] + mpn_add_nc(pp + 15 * n, pp + 15 * n, r1 + 2 * n, n, cy);
495 MPN_INCR_U (pp + 16 * n, spt - n, cy);
496 } else {
497 ASSERT_NOCARRY(mpn_add_nc(pp + 15 * n, pp + 15 * n, r1 + 2 * n, spt, cy));
498 }
499 #else
500 MPN_INCR_U (r1 + 2 * n, n + 1, cy);
501 if(LIKELY(spt > n)) {
502 cy = r1[n3] + mpn_add_n (pp + 15 * n, pp + 15 * n, r1 + 2 * n, n);
503 MPN_INCR_U (pp + 16 * n, spt - n, cy);
504 } else {
505 ASSERT_NOCARRY(mpn_add_n (pp + 15 * n, pp + 15 * n, r1 + 2 * n, spt));
506 }
507 #endif
508 } else {
509 ASSERT_NOCARRY(mpn_add_1 (pp + 14 * n, r1 + n, spt, pp[14 * n]));
510 }
511
512 #undef r0
513 #undef r2
514 #undef r4
515 #undef r6
516 }
517