1 /* Reference mpn functions, designed to be simple, portable and independent
2 of the normal gmp code. Speed isn't a consideration.
3
4 Copyright 1996-2009, 2011-2014 Free Software Foundation, Inc.
5
6 This file is part of the GNU MP Library test suite.
7
8 The GNU MP Library test suite is free software; you can redistribute it
9 and/or modify it under the terms of the GNU General Public License as
10 published by the Free Software Foundation; either version 3 of the License,
11 or (at your option) any later version.
12
13 The GNU MP Library test suite is distributed in the hope that it will be
14 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
16 Public License for more details.
17
18 You should have received a copy of the GNU General Public License along with
19 the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
20
21
22 /* Most routines have assertions representing what the mpn routines are
23 supposed to accept. Many of these reference routines do sensible things
24 outside these ranges (eg. for size==0), but the assertions are present to
25 pick up bad parameters passed here that are about to be passed the same
26 to a real mpn routine being compared. */
27
28 /* always do assertion checking */
29 #define WANT_ASSERT 1
30
31 #include <stdio.h> /* for NULL */
32 #include <stdlib.h> /* for malloc */
33
34 #include "gmp-impl.h"
35 #include "longlong.h"
36
37 #include "tests.h"
38
39
40
41 /* Return non-zero if regions {xp,xsize} and {yp,ysize} overlap, with sizes
42 in bytes. */
43 int
44 byte_overlap_p (const void *v_xp, mp_size_t xsize,
45 const void *v_yp, mp_size_t ysize)
46 {
47 const char *xp = (const char *) v_xp;
48 const char *yp = (const char *) v_yp;
49
50 ASSERT (xsize >= 0);
51 ASSERT (ysize >= 0);
52
53 /* no wraparounds */
54 ASSERT (xp+xsize >= xp);
55 ASSERT (yp+ysize >= yp);
56
57 if (xp + xsize <= yp)
58 return 0;
59
60 if (yp + ysize <= xp)
61 return 0;
62
63 return 1;
64 }
65
66 /* Return non-zero if limb regions {xp,xsize} and {yp,ysize} overlap. */
67 int
68 refmpn_overlap_p (mp_srcptr xp, mp_size_t xsize, mp_srcptr yp, mp_size_t ysize)
69 {
70 return byte_overlap_p (xp, xsize * GMP_LIMB_BYTES,
71 yp, ysize * GMP_LIMB_BYTES);
72 }
73
74 /* Check overlap for a routine defined to work low to high. */
75 int
76 refmpn_overlap_low_to_high_p (mp_srcptr dst, mp_srcptr src, mp_size_t size)
77 {
78 return (dst <= src || ! refmpn_overlap_p (dst, size, src, size));
79 }
80
81 /* Check overlap for a routine defined to work high to low. */
82 int
83 refmpn_overlap_high_to_low_p (mp_srcptr dst, mp_srcptr src, mp_size_t size)
84 {
85 return (dst >= src || ! refmpn_overlap_p (dst, size, src, size));
86 }
87
88 /* Check overlap for a standard routine requiring equal or separate. */
89 int
90 refmpn_overlap_fullonly_p (mp_srcptr dst, mp_srcptr src, mp_size_t size)
91 {
92 return (dst == src || ! refmpn_overlap_p (dst, size, src, size));
93 }
94 int
95 refmpn_overlap_fullonly_two_p (mp_srcptr dst, mp_srcptr src1, mp_srcptr src2,
96 mp_size_t size)
97 {
98 return (refmpn_overlap_fullonly_p (dst, src1, size)
99 && refmpn_overlap_fullonly_p (dst, src2, size));
100 }
101
102
103 mp_ptr
104 refmpn_malloc_limbs (mp_size_t size)
105 {
106 mp_ptr p;
107 ASSERT (size >= 0);
108 if (size == 0)
109 size = 1;
110 p = (mp_ptr) malloc ((size_t) (size * GMP_LIMB_BYTES));
111 ASSERT (p != NULL);
112 return p;
113 }
114
115 /* Free limbs allocated by refmpn_malloc_limbs. NOTE: Can't free
116 * memory allocated by refmpn_malloc_limbs_aligned. */
117 void
118 refmpn_free_limbs (mp_ptr p)
119 {
120 free (p);
121 }
122
123 mp_ptr
124 refmpn_memdup_limbs (mp_srcptr ptr, mp_size_t size)
125 {
126 mp_ptr p;
127 p = refmpn_malloc_limbs (size);
128 refmpn_copyi (p, ptr, size);
129 return p;
130 }
131
132 /* malloc n limbs on a multiple of m bytes boundary */
133 mp_ptr
134 refmpn_malloc_limbs_aligned (mp_size_t n, size_t m)
135 {
136 return (mp_ptr) align_pointer (refmpn_malloc_limbs (n + m-1), m);
137 }
138
139
140 void
141 refmpn_fill (mp_ptr ptr, mp_size_t size, mp_limb_t value)
142 {
143 mp_size_t i;
144 ASSERT (size >= 0);
145 for (i = 0; i < size; i++)
146 ptr[i] = value;
147 }
148
149 void
150 refmpn_zero (mp_ptr ptr, mp_size_t size)
151 {
152 refmpn_fill (ptr, size, CNST_LIMB(0));
153 }
154
155 void
156 refmpn_zero_extend (mp_ptr ptr, mp_size_t oldsize, mp_size_t newsize)
157 {
158 ASSERT (newsize >= oldsize);
159 refmpn_zero (ptr+oldsize, newsize-oldsize);
160 }
161
162 int
163 refmpn_zero_p (mp_srcptr ptr, mp_size_t size)
164 {
165 mp_size_t i;
166 for (i = 0; i < size; i++)
167 if (ptr[i] != 0)
168 return 0;
169 return 1;
170 }
171
172 mp_size_t
173 refmpn_normalize (mp_srcptr ptr, mp_size_t size)
174 {
175 ASSERT (size >= 0);
176 while (size > 0 && ptr[size-1] == 0)
177 size--;
178 return size;
179 }
180
181 /* the highest one bit in x */
182 mp_limb_t
183 refmpn_msbone (mp_limb_t x)
184 {
185 mp_limb_t n = (mp_limb_t) 1 << (GMP_LIMB_BITS-1);
186
187 while (n != 0)
188 {
189 if (x & n)
190 break;
191 n >>= 1;
192 }
193 return n;
194 }
195
196 /* a mask of the highest one bit plus and all bits below */
197 mp_limb_t
198 refmpn_msbone_mask (mp_limb_t x)
199 {
200 if (x == 0)
201 return 0;
202
203 return (refmpn_msbone (x) << 1) - 1;
204 }
205
206 /* How many digits in the given base will fit in a limb.
207 Notice that the product b is allowed to be equal to the limit
208 2^GMP_NUMB_BITS, this ensures the result for base==2 will be
209 GMP_NUMB_BITS (and similarly other powers of 2). */
210 int
211 refmpn_chars_per_limb (int base)
212 {
213 mp_limb_t limit[2], b[2];
214 int chars_per_limb;
215
216 ASSERT (base >= 2);
217
218 limit[0] = 0; /* limit = 2^GMP_NUMB_BITS */
219 limit[1] = 1;
220 b[0] = 1; /* b = 1 */
221 b[1] = 0;
222
223 chars_per_limb = 0;
224 for (;;)
225 {
226 if (refmpn_mul_1 (b, b, (mp_size_t) 2, (mp_limb_t) base))
227 break;
228 if (refmpn_cmp (b, limit, (mp_size_t) 2) > 0)
229 break;
230 chars_per_limb++;
231 }
232 return chars_per_limb;
233 }
234
235 /* The biggest value base**n which fits in GMP_NUMB_BITS. */
236 mp_limb_t
237 refmpn_big_base (int base)
238 {
239 int chars_per_limb = refmpn_chars_per_limb (base);
240 int i;
241 mp_limb_t bb;
242
243 ASSERT (base >= 2);
244 bb = 1;
245 for (i = 0; i < chars_per_limb; i++)
246 bb *= base;
247 return bb;
248 }
249
250
251 void
252 refmpn_setbit (mp_ptr ptr, unsigned long bit)
253 {
254 ptr[bit/GMP_NUMB_BITS] |= CNST_LIMB(1) << (bit%GMP_NUMB_BITS);
255 }
256
257 void
258 refmpn_clrbit (mp_ptr ptr, unsigned long bit)
259 {
260 ptr[bit/GMP_NUMB_BITS] &= ~ (CNST_LIMB(1) << (bit%GMP_NUMB_BITS));
261 }
262
263 #define REFMPN_TSTBIT(ptr,bit) \
264 (((ptr)[(bit)/GMP_NUMB_BITS] & (CNST_LIMB(1) << ((bit)%GMP_NUMB_BITS))) != 0)
265
266 int
267 refmpn_tstbit (mp_srcptr ptr, unsigned long bit)
268 {
269 return REFMPN_TSTBIT (ptr, bit);
270 }
271
272 unsigned long
273 refmpn_scan0 (mp_srcptr ptr, unsigned long bit)
274 {
275 while (REFMPN_TSTBIT (ptr, bit) != 0)
276 bit++;
277 return bit;
278 }
279
280 unsigned long
281 refmpn_scan1 (mp_srcptr ptr, unsigned long bit)
282 {
283 while (REFMPN_TSTBIT (ptr, bit) == 0)
284 bit++;
285 return bit;
286 }
287
288 void
289 refmpn_copy (mp_ptr rp, mp_srcptr sp, mp_size_t size)
290 {
291 ASSERT (refmpn_overlap_fullonly_p (rp, sp, size));
292 refmpn_copyi (rp, sp, size);
293 }
294
295 void
296 refmpn_copyi (mp_ptr rp, mp_srcptr sp, mp_size_t size)
297 {
298 mp_size_t i;
299
300 ASSERT (refmpn_overlap_low_to_high_p (rp, sp, size));
301 ASSERT (size >= 0);
302
303 for (i = 0; i < size; i++)
304 rp[i] = sp[i];
305 }
306
307 void
308 refmpn_copyd (mp_ptr rp, mp_srcptr sp, mp_size_t size)
309 {
310 mp_size_t i;
311
312 ASSERT (refmpn_overlap_high_to_low_p (rp, sp, size));
313 ASSERT (size >= 0);
314
315 for (i = size-1; i >= 0; i--)
316 rp[i] = sp[i];
317 }
318
319 /* Copy {xp,xsize} to {wp,wsize}. If x is shorter, then pad w with low
320 zeros to wsize. If x is longer, then copy just the high wsize limbs. */
321 void
322 refmpn_copy_extend (mp_ptr wp, mp_size_t wsize, mp_srcptr xp, mp_size_t xsize)
323 {
324 ASSERT (wsize >= 0);
325 ASSERT (xsize >= 0);
326
327 /* high part of x if x bigger than w */
328 if (xsize > wsize)
329 {
330 xp += xsize - wsize;
331 xsize = wsize;
332 }
333
334 refmpn_copy (wp + wsize-xsize, xp, xsize);
335 refmpn_zero (wp, wsize-xsize);
336 }
337
338 int
339 refmpn_cmp (mp_srcptr xp, mp_srcptr yp, mp_size_t size)
340 {
341 mp_size_t i;
342
343 ASSERT (size >= 1);
344 ASSERT_MPN (xp, size);
345 ASSERT_MPN (yp, size);
346
347 for (i = size-1; i >= 0; i--)
348 {
349 if (xp[i] > yp[i]) return 1;
350 if (xp[i] < yp[i]) return -1;
351 }
352 return 0;
353 }
354
355 int
356 refmpn_cmp_allowzero (mp_srcptr xp, mp_srcptr yp, mp_size_t size)
357 {
358 if (size == 0)
359 return 0;
360 else
361 return refmpn_cmp (xp, yp, size);
362 }
363
364 int
365 refmpn_cmp_twosizes (mp_srcptr xp, mp_size_t xsize,
366 mp_srcptr yp, mp_size_t ysize)
367 {
368 int opp, cmp;
369
370 ASSERT_MPN (xp, xsize);
371 ASSERT_MPN (yp, ysize);
372
373 opp = (xsize < ysize);
374 if (opp)
375 MPN_SRCPTR_SWAP (xp,xsize, yp,ysize);
376
377 if (! refmpn_zero_p (xp+ysize, xsize-ysize))
378 cmp = 1;
379 else
380 cmp = refmpn_cmp (xp, yp, ysize);
381
382 return (opp ? -cmp : cmp);
383 }
384
385 int
386 refmpn_equal_anynail (mp_srcptr xp, mp_srcptr yp, mp_size_t size)
387 {
388 mp_size_t i;
389 ASSERT (size >= 0);
390
391 for (i = 0; i < size; i++)
392 if (xp[i] != yp[i])
393 return 0;
394 return 1;
395 }
396
397
398 #define LOGOPS(operation) \
399 { \
400 mp_size_t i; \
401 \
402 ASSERT (refmpn_overlap_fullonly_two_p (rp, s1p, s2p, size)); \
403 ASSERT (size >= 1); \
404 ASSERT_MPN (s1p, size); \
405 ASSERT_MPN (s2p, size); \
406 \
407 for (i = 0; i < size; i++) \
408 rp[i] = operation; \
409 }
410
411 void
412 refmpn_and_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
413 {
414 LOGOPS (s1p[i] & s2p[i]);
415 }
416 void
417 refmpn_andn_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
418 {
419 LOGOPS (s1p[i] & ~s2p[i]);
420 }
421 void
422 refmpn_nand_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
423 {
424 LOGOPS ((s1p[i] & s2p[i]) ^ GMP_NUMB_MASK);
425 }
426 void
427 refmpn_ior_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
428 {
429 LOGOPS (s1p[i] | s2p[i]);
430 }
431 void
432 refmpn_iorn_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
433 {
434 LOGOPS (s1p[i] | (s2p[i] ^ GMP_NUMB_MASK));
435 }
436 void
437 refmpn_nior_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
438 {
439 LOGOPS ((s1p[i] | s2p[i]) ^ GMP_NUMB_MASK);
440 }
441 void
442 refmpn_xor_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
443 {
444 LOGOPS (s1p[i] ^ s2p[i]);
445 }
446 void
447 refmpn_xnor_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
448 {
449 LOGOPS ((s1p[i] ^ s2p[i]) ^ GMP_NUMB_MASK);
450 }
451
452
453 /* set *dh,*dl to mh:ml - sh:sl, in full limbs */
454 void
455 refmpn_sub_ddmmss (mp_limb_t *dh, mp_limb_t *dl,
456 mp_limb_t mh, mp_limb_t ml, mp_limb_t sh, mp_limb_t sl)
457 {
458 *dl = ml - sl;
459 *dh = mh - sh - (ml < sl);
460 }
461
462
463 /* set *w to x+y, return 0 or 1 carry */
464 mp_limb_t
465 ref_addc_limb (mp_limb_t *w, mp_limb_t x, mp_limb_t y)
466 {
467 mp_limb_t sum, cy;
468
469 ASSERT_LIMB (x);
470 ASSERT_LIMB (y);
471
472 sum = x + y;
473 #if GMP_NAIL_BITS == 0
474 *w = sum;
475 cy = (sum < x);
476 #else
477 *w = sum & GMP_NUMB_MASK;
478 cy = (sum >> GMP_NUMB_BITS);
479 #endif
480 return cy;
481 }
482
483 /* set *w to x-y, return 0 or 1 borrow */
484 mp_limb_t
485 ref_subc_limb (mp_limb_t *w, mp_limb_t x, mp_limb_t y)
486 {
487 mp_limb_t diff, cy;
488
489 ASSERT_LIMB (x);
490 ASSERT_LIMB (y);
491
492 diff = x - y;
493 #if GMP_NAIL_BITS == 0
494 *w = diff;
495 cy = (diff > x);
496 #else
497 *w = diff & GMP_NUMB_MASK;
498 cy = (diff >> GMP_NUMB_BITS) & 1;
499 #endif
500 return cy;
501 }
502
503 /* set *w to x+y+c (where c is 0 or 1), return 0 or 1 carry */
504 mp_limb_t
505 adc (mp_limb_t *w, mp_limb_t x, mp_limb_t y, mp_limb_t c)
506 {
507 mp_limb_t r;
508
509 ASSERT_LIMB (x);
510 ASSERT_LIMB (y);
511 ASSERT (c == 0 || c == 1);
512
513 r = ref_addc_limb (w, x, y);
514 return r + ref_addc_limb (w, *w, c);
515 }
516
517 /* set *w to x-y-c (where c is 0 or 1), return 0 or 1 borrow */
518 mp_limb_t
519 sbb (mp_limb_t *w, mp_limb_t x, mp_limb_t y, mp_limb_t c)
520 {
521 mp_limb_t r;
522
523 ASSERT_LIMB (x);
524 ASSERT_LIMB (y);
525 ASSERT (c == 0 || c == 1);
526
527 r = ref_subc_limb (w, x, y);
528 return r + ref_subc_limb (w, *w, c);
529 }
530
531
532 #define AORS_1(operation) \
533 { \
534 mp_size_t i; \
535 \
536 ASSERT (refmpn_overlap_fullonly_p (rp, sp, size)); \
537 ASSERT (size >= 1); \
538 ASSERT_MPN (sp, size); \
539 ASSERT_LIMB (n); \
540 \
541 for (i = 0; i < size; i++) \
542 n = operation (&rp[i], sp[i], n); \
543 return n; \
544 }
545
546 mp_limb_t
547 refmpn_add_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t n)
548 {
549 AORS_1 (ref_addc_limb);
550 }
551 mp_limb_t
552 refmpn_sub_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t n)
553 {
554 AORS_1 (ref_subc_limb);
555 }
556
557 #define AORS_NC(operation) \
558 { \
559 mp_size_t i; \
560 \
561 ASSERT (refmpn_overlap_fullonly_two_p (rp, s1p, s2p, size)); \
562 ASSERT (carry == 0 || carry == 1); \
563 ASSERT (size >= 1); \
564 ASSERT_MPN (s1p, size); \
565 ASSERT_MPN (s2p, size); \
566 \
567 for (i = 0; i < size; i++) \
568 carry = operation (&rp[i], s1p[i], s2p[i], carry); \
569 return carry; \
570 }
571
572 mp_limb_t
573 refmpn_add_nc (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size,
574 mp_limb_t carry)
575 {
576 AORS_NC (adc);
577 }
578 mp_limb_t
579 refmpn_sub_nc (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size,
580 mp_limb_t carry)
581 {
582 AORS_NC (sbb);
583 }
584
585
586 mp_limb_t
587 refmpn_add_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
588 {
589 return refmpn_add_nc (rp, s1p, s2p, size, CNST_LIMB(0));
590 }
591 mp_limb_t
592 refmpn_sub_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
593 {
594 return refmpn_sub_nc (rp, s1p, s2p, size, CNST_LIMB(0));
595 }
596
597 mp_limb_t
598 refmpn_cnd_add_n (mp_limb_t cnd, mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
599 {
600 if (cnd != 0)
601 return refmpn_add_n (rp, s1p, s2p, size);
602 else
603 {
604 refmpn_copyi (rp, s1p, size);
605 return 0;
606 }
607 }
608 mp_limb_t
609 refmpn_cnd_sub_n (mp_limb_t cnd, mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
610 {
611 if (cnd != 0)
612 return refmpn_sub_n (rp, s1p, s2p, size);
613 else
614 {
615 refmpn_copyi (rp, s1p, size);
616 return 0;
617 }
618 }
619
620
621 #define AORS_ERR1_N(operation) \
622 { \
623 mp_size_t i; \
624 mp_limb_t carry2; \
625 \
626 ASSERT (refmpn_overlap_fullonly_p (rp, s1p, size)); \
627 ASSERT (refmpn_overlap_fullonly_p (rp, s2p, size)); \
628 ASSERT (! refmpn_overlap_p (rp, size, yp, size)); \
629 ASSERT (! refmpn_overlap_p (ep, 2, s1p, size)); \
630 ASSERT (! refmpn_overlap_p (ep, 2, s2p, size)); \
631 ASSERT (! refmpn_overlap_p (ep, 2, yp, size)); \
632 ASSERT (! refmpn_overlap_p (ep, 2, rp, size)); \
633 \
634 ASSERT (carry == 0 || carry == 1); \
635 ASSERT (size >= 1); \
636 ASSERT_MPN (s1p, size); \
637 ASSERT_MPN (s2p, size); \
638 ASSERT_MPN (yp, size); \
639 \
640 ep[0] = ep[1] = CNST_LIMB(0); \
641 \
642 for (i = 0; i < size; i++) \
643 { \
644 carry = operation (&rp[i], s1p[i], s2p[i], carry); \
645 if (carry == 1) \
646 { \
647 carry2 = ref_addc_limb (&ep[0], ep[0], yp[size - 1 - i]); \
648 carry2 = ref_addc_limb (&ep[1], ep[1], carry2); \
649 ASSERT (carry2 == 0); \
650 } \
651 } \
652 return carry; \
653 }
654
655 mp_limb_t
656 refmpn_add_err1_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p,
657 mp_ptr ep, mp_srcptr yp,
658 mp_size_t size, mp_limb_t carry)
659 {
660 AORS_ERR1_N (adc);
661 }
662 mp_limb_t
663 refmpn_sub_err1_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p,
664 mp_ptr ep, mp_srcptr yp,
665 mp_size_t size, mp_limb_t carry)
666 {
667 AORS_ERR1_N (sbb);
668 }
669
670
671 #define AORS_ERR2_N(operation) \
672 { \
673 mp_size_t i; \
674 mp_limb_t carry2; \
675 \
676 ASSERT (refmpn_overlap_fullonly_p (rp, s1p, size)); \
677 ASSERT (refmpn_overlap_fullonly_p (rp, s2p, size)); \
678 ASSERT (! refmpn_overlap_p (rp, size, y1p, size)); \
679 ASSERT (! refmpn_overlap_p (rp, size, y2p, size)); \
680 ASSERT (! refmpn_overlap_p (ep, 4, s1p, size)); \
681 ASSERT (! refmpn_overlap_p (ep, 4, s2p, size)); \
682 ASSERT (! refmpn_overlap_p (ep, 4, y1p, size)); \
683 ASSERT (! refmpn_overlap_p (ep, 4, y2p, size)); \
684 ASSERT (! refmpn_overlap_p (ep, 4, rp, size)); \
685 \
686 ASSERT (carry == 0 || carry == 1); \
687 ASSERT (size >= 1); \
688 ASSERT_MPN (s1p, size); \
689 ASSERT_MPN (s2p, size); \
690 ASSERT_MPN (y1p, size); \
691 ASSERT_MPN (y2p, size); \
692 \
693 ep[0] = ep[1] = CNST_LIMB(0); \
694 ep[2] = ep[3] = CNST_LIMB(0); \
695 \
696 for (i = 0; i < size; i++) \
697 { \
698 carry = operation (&rp[i], s1p[i], s2p[i], carry); \
699 if (carry == 1) \
700 { \
701 carry2 = ref_addc_limb (&ep[0], ep[0], y1p[size - 1 - i]); \
702 carry2 = ref_addc_limb (&ep[1], ep[1], carry2); \
703 ASSERT (carry2 == 0); \
704 carry2 = ref_addc_limb (&ep[2], ep[2], y2p[size - 1 - i]); \
705 carry2 = ref_addc_limb (&ep[3], ep[3], carry2); \
706 ASSERT (carry2 == 0); \
707 } \
708 } \
709 return carry; \
710 }
711
712 mp_limb_t
713 refmpn_add_err2_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p,
714 mp_ptr ep, mp_srcptr y1p, mp_srcptr y2p,
715 mp_size_t size, mp_limb_t carry)
716 {
717 AORS_ERR2_N (adc);
718 }
719 mp_limb_t
720 refmpn_sub_err2_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p,
721 mp_ptr ep, mp_srcptr y1p, mp_srcptr y2p,
722 mp_size_t size, mp_limb_t carry)
723 {
724 AORS_ERR2_N (sbb);
725 }
726
727
728 #define AORS_ERR3_N(operation) \
729 { \
730 mp_size_t i; \
731 mp_limb_t carry2; \
732 \
733 ASSERT (refmpn_overlap_fullonly_p (rp, s1p, size)); \
734 ASSERT (refmpn_overlap_fullonly_p (rp, s2p, size)); \
735 ASSERT (! refmpn_overlap_p (rp, size, y1p, size)); \
736 ASSERT (! refmpn_overlap_p (rp, size, y2p, size)); \
737 ASSERT (! refmpn_overlap_p (rp, size, y3p, size)); \
738 ASSERT (! refmpn_overlap_p (ep, 6, s1p, size)); \
739 ASSERT (! refmpn_overlap_p (ep, 6, s2p, size)); \
740 ASSERT (! refmpn_overlap_p (ep, 6, y1p, size)); \
741 ASSERT (! refmpn_overlap_p (ep, 6, y2p, size)); \
742 ASSERT (! refmpn_overlap_p (ep, 6, y3p, size)); \
743 ASSERT (! refmpn_overlap_p (ep, 6, rp, size)); \
744 \
745 ASSERT (carry == 0 || carry == 1); \
746 ASSERT (size >= 1); \
747 ASSERT_MPN (s1p, size); \
748 ASSERT_MPN (s2p, size); \
749 ASSERT_MPN (y1p, size); \
750 ASSERT_MPN (y2p, size); \
751 ASSERT_MPN (y3p, size); \
752 \
753 ep[0] = ep[1] = CNST_LIMB(0); \
754 ep[2] = ep[3] = CNST_LIMB(0); \
755 ep[4] = ep[5] = CNST_LIMB(0); \
756 \
757 for (i = 0; i < size; i++) \
758 { \
759 carry = operation (&rp[i], s1p[i], s2p[i], carry); \
760 if (carry == 1) \
761 { \
762 carry2 = ref_addc_limb (&ep[0], ep[0], y1p[size - 1 - i]); \
763 carry2 = ref_addc_limb (&ep[1], ep[1], carry2); \
764 ASSERT (carry2 == 0); \
765 carry2 = ref_addc_limb (&ep[2], ep[2], y2p[size - 1 - i]); \
766 carry2 = ref_addc_limb (&ep[3], ep[3], carry2); \
767 ASSERT (carry2 == 0); \
768 carry2 = ref_addc_limb (&ep[4], ep[4], y3p[size - 1 - i]); \
769 carry2 = ref_addc_limb (&ep[5], ep[5], carry2); \
770 ASSERT (carry2 == 0); \
771 } \
772 } \
773 return carry; \
774 }
775
776 mp_limb_t
777 refmpn_add_err3_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p,
778 mp_ptr ep, mp_srcptr y1p, mp_srcptr y2p, mp_srcptr y3p,
779 mp_size_t size, mp_limb_t carry)
780 {
781 AORS_ERR3_N (adc);
782 }
783 mp_limb_t
784 refmpn_sub_err3_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p,
785 mp_ptr ep, mp_srcptr y1p, mp_srcptr y2p, mp_srcptr y3p,
786 mp_size_t size, mp_limb_t carry)
787 {
788 AORS_ERR3_N (sbb);
789 }
790
791
792 mp_limb_t
793 refmpn_addlsh_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp,
794 mp_size_t n, unsigned int s)
795 {
796 mp_limb_t cy;
797 mp_ptr tp;
798
799 ASSERT (refmpn_overlap_fullonly_two_p (rp, up, vp, n));
800 ASSERT (n >= 1);
801 ASSERT (0 < s && s < GMP_NUMB_BITS);
802 ASSERT_MPN (up, n);
803 ASSERT_MPN (vp, n);
804
805 tp = refmpn_malloc_limbs (n);
806 cy = refmpn_lshift (tp, vp, n, s);
807 cy += refmpn_add_n (rp, up, tp, n);
808 free (tp);
809 return cy;
810 }
811 mp_limb_t
812 refmpn_addlsh1_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n)
813 {
814 return refmpn_addlsh_n (rp, up, vp, n, 1);
815 }
816 mp_limb_t
817 refmpn_addlsh2_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n)
818 {
819 return refmpn_addlsh_n (rp, up, vp, n, 2);
820 }
821 mp_limb_t
822 refmpn_addlsh_n_ip1 (mp_ptr rp, mp_srcptr vp, mp_size_t n, unsigned int s)
823 {
824 return refmpn_addlsh_n (rp, rp, vp, n, s);
825 }
826 mp_limb_t
827 refmpn_addlsh1_n_ip1 (mp_ptr rp, mp_srcptr vp, mp_size_t n)
828 {
829 return refmpn_addlsh_n (rp, rp, vp, n, 1);
830 }
831 mp_limb_t
832 refmpn_addlsh2_n_ip1 (mp_ptr rp, mp_srcptr vp, mp_size_t n)
833 {
834 return refmpn_addlsh_n (rp, rp, vp, n, 2);
835 }
836 mp_limb_t
837 refmpn_addlsh_n_ip2 (mp_ptr rp, mp_srcptr vp, mp_size_t n, unsigned int s)
838 {
839 return refmpn_addlsh_n (rp, vp, rp, n, s);
840 }
841 mp_limb_t
842 refmpn_addlsh1_n_ip2 (mp_ptr rp, mp_srcptr vp, mp_size_t n)
843 {
844 return refmpn_addlsh_n (rp, vp, rp, n, 1);
845 }
846 mp_limb_t
847 refmpn_addlsh2_n_ip2 (mp_ptr rp, mp_srcptr vp, mp_size_t n)
848 {
849 return refmpn_addlsh_n (rp, vp, rp, n, 2);
850 }
851 mp_limb_t
852 refmpn_addlsh_nc (mp_ptr rp, mp_srcptr up, mp_srcptr vp,
853 mp_size_t n, unsigned int s, mp_limb_t carry)
854 {
855 mp_limb_t cy;
856
857 ASSERT (carry <= (CNST_LIMB(1) << s));
858
859 cy = refmpn_addlsh_n (rp, up, vp, n, s);
860 cy += refmpn_add_1 (rp, rp, n, carry);
861 return cy;
862 }
863 mp_limb_t
864 refmpn_addlsh1_nc (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n, mp_limb_t carry)
865 {
866 return refmpn_addlsh_nc (rp, up, vp, n, 1, carry);
867 }
868 mp_limb_t
869 refmpn_addlsh2_nc (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n, mp_limb_t carry)
870 {
871 return refmpn_addlsh_nc (rp, up, vp, n, 2, carry);
872 }
873
874 mp_limb_t
875 refmpn_sublsh_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp,
876 mp_size_t n, unsigned int s)
877 {
878 mp_limb_t cy;
879 mp_ptr tp;
880
881 ASSERT (refmpn_overlap_fullonly_two_p (rp, up, vp, n));
882 ASSERT (n >= 1);
883 ASSERT (0 < s && s < GMP_NUMB_BITS);
884 ASSERT_MPN (up, n);
885 ASSERT_MPN (vp, n);
886
887 tp = refmpn_malloc_limbs (n);
888 cy = mpn_lshift (tp, vp, n, s);
889 cy += mpn_sub_n (rp, up, tp, n);
890 free (tp);
891 return cy;
892 }
893 mp_limb_t
894 refmpn_sublsh1_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n)
895 {
896 return refmpn_sublsh_n (rp, up, vp, n, 1);
897 }
898 mp_limb_t
899 refmpn_sublsh2_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n)
900 {
901 return refmpn_sublsh_n (rp, up, vp, n, 2);
902 }
903 mp_limb_t
904 refmpn_sublsh_n_ip1 (mp_ptr rp, mp_srcptr vp, mp_size_t n, unsigned int s)
905 {
906 return refmpn_sublsh_n (rp, rp, vp, n, s);
907 }
908 mp_limb_t
909 refmpn_sublsh1_n_ip1 (mp_ptr rp, mp_srcptr vp, mp_size_t n)
910 {
911 return refmpn_sublsh_n (rp, rp, vp, n, 1);
912 }
913 mp_limb_t
914 refmpn_sublsh2_n_ip1 (mp_ptr rp, mp_srcptr vp, mp_size_t n)
915 {
916 return refmpn_sublsh_n (rp, rp, vp, n, 2);
917 }
918 mp_limb_t
919 refmpn_sublsh_n_ip2 (mp_ptr rp, mp_srcptr vp, mp_size_t n, unsigned int s)
920 {
921 return refmpn_sublsh_n (rp, vp, rp, n, s);
922 }
923 mp_limb_t
924 refmpn_sublsh1_n_ip2 (mp_ptr rp, mp_srcptr vp, mp_size_t n)
925 {
926 return refmpn_sublsh_n (rp, vp, rp, n, 1);
927 }
928 mp_limb_t
929 refmpn_sublsh2_n_ip2 (mp_ptr rp, mp_srcptr vp, mp_size_t n)
930 {
931 return refmpn_sublsh_n (rp, vp, rp, n, 2);
932 }
933 mp_limb_t
934 refmpn_sublsh_nc (mp_ptr rp, mp_srcptr up, mp_srcptr vp,
935 mp_size_t n, unsigned int s, mp_limb_t carry)
936 {
937 mp_limb_t cy;
938
939 ASSERT (carry <= (CNST_LIMB(1) << s));
940
941 cy = refmpn_sublsh_n (rp, up, vp, n, s);
942 cy += refmpn_sub_1 (rp, rp, n, carry);
943 return cy;
944 }
945 mp_limb_t
946 refmpn_sublsh1_nc (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n, mp_limb_t carry)
947 {
948 return refmpn_sublsh_nc (rp, up, vp, n, 1, carry);
949 }
950 mp_limb_t
951 refmpn_sublsh2_nc (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n, mp_limb_t carry)
952 {
953 return refmpn_sublsh_nc (rp, up, vp, n, 2, carry);
954 }
955
956 mp_limb_signed_t
957 refmpn_rsblsh_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp,
958 mp_size_t n, unsigned int s)
959 {
960 mp_limb_signed_t cy;
961 mp_ptr tp;
962
963 ASSERT (refmpn_overlap_fullonly_two_p (rp, up, vp, n));
964 ASSERT (n >= 1);
965 ASSERT (0 < s && s < GMP_NUMB_BITS);
966 ASSERT_MPN (up, n);
967 ASSERT_MPN (vp, n);
968
969 tp = refmpn_malloc_limbs (n);
970 cy = mpn_lshift (tp, vp, n, s);
971 cy -= mpn_sub_n (rp, tp, up, n);
972 free (tp);
973 return cy;
974 }
975 mp_limb_signed_t
976 refmpn_rsblsh1_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n)
977 {
978 return refmpn_rsblsh_n (rp, up, vp, n, 1);
979 }
980 mp_limb_signed_t
981 refmpn_rsblsh2_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n)
982 {
983 return refmpn_rsblsh_n (rp, up, vp, n, 2);
984 }
985 mp_limb_signed_t
986 refmpn_rsblsh_nc (mp_ptr rp, mp_srcptr up, mp_srcptr vp,
987 mp_size_t n, unsigned int s, mp_limb_signed_t carry)
988 {
989 mp_limb_signed_t cy;
990
991 ASSERT (carry == -1 || (carry >> s) == 0);
992
993 cy = refmpn_rsblsh_n (rp, up, vp, n, s);
994 if (carry > 0)
995 cy += refmpn_add_1 (rp, rp, n, carry);
996 else
997 cy -= refmpn_sub_1 (rp, rp, n, -carry);
998 return cy;
999 }
1000 mp_limb_signed_t
1001 refmpn_rsblsh1_nc (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n, mp_limb_signed_t carry)
1002 {
1003 return refmpn_rsblsh_nc (rp, up, vp, n, 1, carry);
1004 }
1005 mp_limb_signed_t
1006 refmpn_rsblsh2_nc (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n, mp_limb_signed_t carry)
1007 {
1008 return refmpn_rsblsh_nc (rp, up, vp, n, 2, carry);
1009 }
1010
1011 mp_limb_t
1012 refmpn_rsh1add_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n)
1013 {
1014 mp_limb_t cya, cys;
1015
1016 ASSERT (refmpn_overlap_fullonly_two_p (rp, up, vp, n));
1017 ASSERT (n >= 1);
1018 ASSERT_MPN (up, n);
1019 ASSERT_MPN (vp, n);
1020
1021 cya = mpn_add_n (rp, up, vp, n);
1022 cys = mpn_rshift (rp, rp, n, 1) >> (GMP_NUMB_BITS - 1);
1023 rp[n - 1] |= cya << (GMP_NUMB_BITS - 1);
1024 return cys;
1025 }
1026 mp_limb_t
1027 refmpn_rsh1sub_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n)
1028 {
1029 mp_limb_t cya, cys;
1030
1031 ASSERT (refmpn_overlap_fullonly_two_p (rp, up, vp, n));
1032 ASSERT (n >= 1);
1033 ASSERT_MPN (up, n);
1034 ASSERT_MPN (vp, n);
1035
1036 cya = mpn_sub_n (rp, up, vp, n);
1037 cys = mpn_rshift (rp, rp, n, 1) >> (GMP_NUMB_BITS - 1);
1038 rp[n - 1] |= cya << (GMP_NUMB_BITS - 1);
1039 return cys;
1040 }
1041
1042 /* Twos complement, return borrow. */
1043 mp_limb_t
1044 refmpn_neg (mp_ptr dst, mp_srcptr src, mp_size_t size)
1045 {
1046 mp_ptr zeros;
1047 mp_limb_t ret;
1048
1049 ASSERT (size >= 1);
1050
1051 zeros = refmpn_malloc_limbs (size);
1052 refmpn_fill (zeros, size, CNST_LIMB(0));
1053 ret = refmpn_sub_n (dst, zeros, src, size);
1054 free (zeros);
1055 return ret;
1056 }
1057
1058
1059 #define AORS(aors_n, aors_1) \
1060 { \
1061 mp_limb_t c; \
1062 ASSERT (s1size >= s2size); \
1063 ASSERT (s2size >= 1); \
1064 c = aors_n (rp, s1p, s2p, s2size); \
1065 if (s1size-s2size != 0) \
1066 c = aors_1 (rp+s2size, s1p+s2size, s1size-s2size, c); \
1067 return c; \
1068 }
1069 mp_limb_t
1070 refmpn_add (mp_ptr rp,
1071 mp_srcptr s1p, mp_size_t s1size,
1072 mp_srcptr s2p, mp_size_t s2size)
1073 {
1074 AORS (refmpn_add_n, refmpn_add_1);
1075 }
1076 mp_limb_t
1077 refmpn_sub (mp_ptr rp,
1078 mp_srcptr s1p, mp_size_t s1size,
1079 mp_srcptr s2p, mp_size_t s2size)
1080 {
1081 AORS (refmpn_sub_n, refmpn_sub_1);
1082 }
1083
1084
1085 #define SHIFTHIGH(x) ((x) << GMP_LIMB_BITS/2)
1086 #define SHIFTLOW(x) ((x) >> GMP_LIMB_BITS/2)
1087
1088 #define LOWMASK (((mp_limb_t) 1 << GMP_LIMB_BITS/2)-1)
1089 #define HIGHMASK SHIFTHIGH(LOWMASK)
1090
1091 #define LOWPART(x) ((x) & LOWMASK)
1092 #define HIGHPART(x) SHIFTLOW((x) & HIGHMASK)
1093
1094 /* Set return:*lo to x*y, using full limbs not nails. */
1095 mp_limb_t
1096 refmpn_umul_ppmm (mp_limb_t *lo, mp_limb_t x, mp_limb_t y)
1097 {
1098 mp_limb_t hi, s;
1099
1100 *lo = LOWPART(x) * LOWPART(y);
1101 hi = HIGHPART(x) * HIGHPART(y);
1102
1103 s = LOWPART(x) * HIGHPART(y);
1104 hi += HIGHPART(s);
1105 s = SHIFTHIGH(LOWPART(s));
1106 *lo += s;
1107 hi += (*lo < s);
1108
1109 s = HIGHPART(x) * LOWPART(y);
1110 hi += HIGHPART(s);
1111 s = SHIFTHIGH(LOWPART(s));
1112 *lo += s;
1113 hi += (*lo < s);
1114
1115 return hi;
1116 }
1117
1118 mp_limb_t
1119 refmpn_umul_ppmm_r (mp_limb_t x, mp_limb_t y, mp_limb_t *lo)
1120 {
1121 return refmpn_umul_ppmm (lo, x, y);
1122 }
1123
1124 mp_limb_t
1125 refmpn_mul_1c (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t multiplier,
1126 mp_limb_t carry)
1127 {
1128 mp_size_t i;
1129 mp_limb_t hi, lo;
1130
1131 ASSERT (refmpn_overlap_low_to_high_p (rp, sp, size));
1132 ASSERT (size >= 1);
1133 ASSERT_MPN (sp, size);
1134 ASSERT_LIMB (multiplier);
1135 ASSERT_LIMB (carry);
1136
1137 multiplier <<= GMP_NAIL_BITS;
1138 for (i = 0; i < size; i++)
1139 {
1140 hi = refmpn_umul_ppmm (&lo, sp[i], multiplier);
1141 lo >>= GMP_NAIL_BITS;
1142 ASSERT_NOCARRY (ref_addc_limb (&hi, hi, ref_addc_limb (&lo, lo, carry)));
1143 rp[i] = lo;
1144 carry = hi;
1145 }
1146 return carry;
1147 }
1148
1149 mp_limb_t
1150 refmpn_mul_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t multiplier)
1151 {
1152 return refmpn_mul_1c (rp, sp, size, multiplier, CNST_LIMB(0));
1153 }
1154
1155
1156 mp_limb_t
1157 refmpn_mul_N (mp_ptr dst, mp_srcptr src, mp_size_t size,
1158 mp_srcptr mult, mp_size_t msize)
1159 {
1160 mp_ptr src_copy;
1161 mp_limb_t ret;
1162 mp_size_t i;
1163
1164 ASSERT (refmpn_overlap_fullonly_p (dst, src, size));
1165 ASSERT (! refmpn_overlap_p (dst, size+msize-1, mult, msize));
1166 ASSERT (size >= msize);
1167 ASSERT_MPN (mult, msize);
1168
1169 /* in case dst==src */
1170 src_copy = refmpn_malloc_limbs (size);
1171 refmpn_copyi (src_copy, src, size);
1172 src = src_copy;
1173
1174 dst[size] = refmpn_mul_1 (dst, src, size, mult[0]);
1175 for (i = 1; i < msize-1; i++)
1176 dst[size+i] = refmpn_addmul_1 (dst+i, src, size, mult[i]);
1177 ret = refmpn_addmul_1 (dst+i, src, size, mult[i]);
1178
1179 free (src_copy);
1180 return ret;
1181 }
1182
1183 mp_limb_t
1184 refmpn_mul_2 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
1185 {
1186 return refmpn_mul_N (rp, sp, size, mult, (mp_size_t) 2);
1187 }
1188 mp_limb_t
1189 refmpn_mul_3 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
1190 {
1191 return refmpn_mul_N (rp, sp, size, mult, (mp_size_t) 3);
1192 }
1193 mp_limb_t
1194 refmpn_mul_4 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
1195 {
1196 return refmpn_mul_N (rp, sp, size, mult, (mp_size_t) 4);
1197 }
1198 mp_limb_t
1199 refmpn_mul_5 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
1200 {
1201 return refmpn_mul_N (rp, sp, size, mult, (mp_size_t) 5);
1202 }
1203 mp_limb_t
1204 refmpn_mul_6 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
1205 {
1206 return refmpn_mul_N (rp, sp, size, mult, (mp_size_t) 6);
1207 }
1208
1209 #define AORSMUL_1C(operation_n) \
1210 { \
1211 mp_ptr p; \
1212 mp_limb_t ret; \
1213 \
1214 ASSERT (refmpn_overlap_fullonly_p (rp, sp, size)); \
1215 \
1216 p = refmpn_malloc_limbs (size); \
1217 ret = refmpn_mul_1c (p, sp, size, multiplier, carry); \
1218 ret += operation_n (rp, rp, p, size); \
1219 \
1220 free (p); \
1221 return ret; \
1222 }
1223
1224 mp_limb_t
1225 refmpn_addmul_1c (mp_ptr rp, mp_srcptr sp, mp_size_t size,
1226 mp_limb_t multiplier, mp_limb_t carry)
1227 {
1228 AORSMUL_1C (refmpn_add_n);
1229 }
1230 mp_limb_t
1231 refmpn_submul_1c (mp_ptr rp, mp_srcptr sp, mp_size_t size,
1232 mp_limb_t multiplier, mp_limb_t carry)
1233 {
1234 AORSMUL_1C (refmpn_sub_n);
1235 }
1236
1237
1238 mp_limb_t
1239 refmpn_addmul_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t multiplier)
1240 {
1241 return refmpn_addmul_1c (rp, sp, size, multiplier, CNST_LIMB(0));
1242 }
1243 mp_limb_t
1244 refmpn_submul_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t multiplier)
1245 {
1246 return refmpn_submul_1c (rp, sp, size, multiplier, CNST_LIMB(0));
1247 }
1248
1249
1250 mp_limb_t
1251 refmpn_addmul_N (mp_ptr dst, mp_srcptr src, mp_size_t size,
1252 mp_srcptr mult, mp_size_t msize)
1253 {
1254 mp_ptr src_copy;
1255 mp_limb_t ret;
1256 mp_size_t i;
1257
1258 ASSERT (dst == src || ! refmpn_overlap_p (dst, size+msize-1, src, size));
1259 ASSERT (! refmpn_overlap_p (dst, size+msize-1, mult, msize));
1260 ASSERT (size >= msize);
1261 ASSERT_MPN (mult, msize);
1262
1263 /* in case dst==src */
1264 src_copy = refmpn_malloc_limbs (size);
1265 refmpn_copyi (src_copy, src, size);
1266 src = src_copy;
1267
1268 for (i = 0; i < msize-1; i++)
1269 dst[size+i] = refmpn_addmul_1 (dst+i, src, size, mult[i]);
1270 ret = refmpn_addmul_1 (dst+i, src, size, mult[i]);
1271
1272 free (src_copy);
1273 return ret;
1274 }
1275
1276 mp_limb_t
1277 refmpn_addmul_2 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
1278 {
1279 return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 2);
1280 }
1281 mp_limb_t
1282 refmpn_addmul_3 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
1283 {
1284 return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 3);
1285 }
1286 mp_limb_t
1287 refmpn_addmul_4 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
1288 {
1289 return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 4);
1290 }
1291 mp_limb_t
1292 refmpn_addmul_5 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
1293 {
1294 return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 5);
1295 }
1296 mp_limb_t
1297 refmpn_addmul_6 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
1298 {
1299 return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 6);
1300 }
1301 mp_limb_t
1302 refmpn_addmul_7 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
1303 {
1304 return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 7);
1305 }
1306 mp_limb_t
1307 refmpn_addmul_8 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_srcptr mult)
1308 {
1309 return refmpn_addmul_N (rp, sp, size, mult, (mp_size_t) 8);
1310 }
1311
1312 mp_limb_t
1313 refmpn_add_n_sub_nc (mp_ptr r1p, mp_ptr r2p,
1314 mp_srcptr s1p, mp_srcptr s2p, mp_size_t size,
1315 mp_limb_t carry)
1316 {
1317 mp_ptr p;
1318 mp_limb_t acy, scy;
1319
1320 /* Destinations can't overlap. */
1321 ASSERT (! refmpn_overlap_p (r1p, size, r2p, size));
1322 ASSERT (refmpn_overlap_fullonly_two_p (r1p, s1p, s2p, size));
1323 ASSERT (refmpn_overlap_fullonly_two_p (r2p, s1p, s2p, size));
1324 ASSERT (size >= 1);
1325
1326 /* in case r1p==s1p or r1p==s2p */
1327 p = refmpn_malloc_limbs (size);
1328
1329 acy = refmpn_add_nc (p, s1p, s2p, size, carry >> 1);
1330 scy = refmpn_sub_nc (r2p, s1p, s2p, size, carry & 1);
1331 refmpn_copyi (r1p, p, size);
1332
1333 free (p);
1334 return 2 * acy + scy;
1335 }
1336
1337 mp_limb_t
1338 refmpn_add_n_sub_n (mp_ptr r1p, mp_ptr r2p,
1339 mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
1340 {
1341 return refmpn_add_n_sub_nc (r1p, r2p, s1p, s2p, size, CNST_LIMB(0));
1342 }
1343
1344
1345 /* Right shift hi,lo and return the low limb of the result.
1346 Note a shift by GMP_LIMB_BITS isn't assumed to work (doesn't on x86). */
1347 mp_limb_t
1348 rshift_make (mp_limb_t hi, mp_limb_t lo, unsigned shift)
1349 {
1350 ASSERT (shift < GMP_NUMB_BITS);
1351 if (shift == 0)
1352 return lo;
1353 else
1354 return ((hi << (GMP_NUMB_BITS-shift)) | (lo >> shift)) & GMP_NUMB_MASK;
1355 }
1356
1357 /* Left shift hi,lo and return the high limb of the result.
1358 Note a shift by GMP_LIMB_BITS isn't assumed to work (doesn't on x86). */
1359 mp_limb_t
1360 lshift_make (mp_limb_t hi, mp_limb_t lo, unsigned shift)
1361 {
1362 ASSERT (shift < GMP_NUMB_BITS);
1363 if (shift == 0)
1364 return hi;
1365 else
1366 return ((hi << shift) | (lo >> (GMP_NUMB_BITS-shift))) & GMP_NUMB_MASK;
1367 }
1368
1369
1370 mp_limb_t
1371 refmpn_rshift (mp_ptr rp, mp_srcptr sp, mp_size_t size, unsigned shift)
1372 {
1373 mp_limb_t ret;
1374 mp_size_t i;
1375
1376 ASSERT (refmpn_overlap_low_to_high_p (rp, sp, size));
1377 ASSERT (size >= 1);
1378 ASSERT (shift >= 1 && shift < GMP_NUMB_BITS);
1379 ASSERT_MPN (sp, size);
1380
1381 ret = rshift_make (sp[0], CNST_LIMB(0), shift);
1382
1383 for (i = 0; i < size-1; i++)
1384 rp[i] = rshift_make (sp[i+1], sp[i], shift);
1385
1386 rp[i] = rshift_make (CNST_LIMB(0), sp[i], shift);
1387 return ret;
1388 }
1389
1390 mp_limb_t
1391 refmpn_lshift (mp_ptr rp, mp_srcptr sp, mp_size_t size, unsigned shift)
1392 {
1393 mp_limb_t ret;
1394 mp_size_t i;
1395
1396 ASSERT (refmpn_overlap_high_to_low_p (rp, sp, size));
1397 ASSERT (size >= 1);
1398 ASSERT (shift >= 1 && shift < GMP_NUMB_BITS);
1399 ASSERT_MPN (sp, size);
1400
1401 ret = lshift_make (CNST_LIMB(0), sp[size-1], shift);
1402
1403 for (i = size-2; i >= 0; i--)
1404 rp[i+1] = lshift_make (sp[i+1], sp[i], shift);
1405
1406 rp[i+1] = lshift_make (sp[i+1], CNST_LIMB(0), shift);
1407 return ret;
1408 }
1409
1410 void
1411 refmpn_com (mp_ptr rp, mp_srcptr sp, mp_size_t size)
1412 {
1413 mp_size_t i;
1414
1415 /* We work downwards since mpn_lshiftc needs that. */
1416 ASSERT (refmpn_overlap_high_to_low_p (rp, sp, size));
1417
1418 for (i = size - 1; i >= 0; i--)
1419 rp[i] = (~sp[i]) & GMP_NUMB_MASK;
1420 }
1421
1422 mp_limb_t
1423 refmpn_lshiftc (mp_ptr rp, mp_srcptr sp, mp_size_t size, unsigned shift)
1424 {
1425 mp_limb_t res;
1426
1427 /* No asserts here, refmpn_lshift will assert what we need. */
1428
1429 res = refmpn_lshift (rp, sp, size, shift);
1430 refmpn_com (rp, rp, size);
1431 return res;
1432 }
1433
1434 /* accepting shift==0 and doing a plain copyi or copyd in that case */
1435 mp_limb_t
1436 refmpn_rshift_or_copy (mp_ptr rp, mp_srcptr sp, mp_size_t size, unsigned shift)
1437 {
1438 if (shift == 0)
1439 {
1440 refmpn_copyi (rp, sp, size);
1441 return 0;
1442 }
1443 else
1444 {
1445 return refmpn_rshift (rp, sp, size, shift);
1446 }
1447 }
1448 mp_limb_t
1449 refmpn_lshift_or_copy (mp_ptr rp, mp_srcptr sp, mp_size_t size, unsigned shift)
1450 {
1451 if (shift == 0)
1452 {
1453 refmpn_copyd (rp, sp, size);
1454 return 0;
1455 }
1456 else
1457 {
1458 return refmpn_lshift (rp, sp, size, shift);
1459 }
1460 }
1461
1462 /* accepting size==0 too */
1463 mp_limb_t
1464 refmpn_rshift_or_copy_any (mp_ptr rp, mp_srcptr sp, mp_size_t size,
1465 unsigned shift)
1466 {
1467 return (size == 0 ? 0 : refmpn_rshift_or_copy (rp, sp, size, shift));
1468 }
1469 mp_limb_t
1470 refmpn_lshift_or_copy_any (mp_ptr rp, mp_srcptr sp, mp_size_t size,
1471 unsigned shift)
1472 {
1473 return (size == 0 ? 0 : refmpn_lshift_or_copy (rp, sp, size, shift));
1474 }
1475
1476 /* Divide h,l by d, return quotient, store remainder to *rp.
1477 Operates on full limbs, not nails.
1478 Must have h < d.
1479 __udiv_qrnnd_c isn't simple, and it's a bit slow, but it works. */
1480 mp_limb_t
1481 refmpn_udiv_qrnnd (mp_limb_t *rp, mp_limb_t h, mp_limb_t l, mp_limb_t d)
1482 {
1483 mp_limb_t q, r;
1484 int n;
1485
1486 ASSERT (d != 0);
1487 ASSERT (h < d);
1488
1489 #if 0
1490 udiv_qrnnd (q, r, h, l, d);
1491 *rp = r;
1492 return q;
1493 #endif
1494
1495 n = refmpn_count_leading_zeros (d);
1496 d <<= n;
1497
1498 if (n != 0)
1499 {
1500 h = (h << n) | (l >> (GMP_LIMB_BITS - n));
1501 l <<= n;
1502 }
1503
1504 __udiv_qrnnd_c (q, r, h, l, d);
1505 r >>= n;
1506 *rp = r;
1507 return q;
1508 }
1509
1510 mp_limb_t
1511 refmpn_udiv_qrnnd_r (mp_limb_t h, mp_limb_t l, mp_limb_t d, mp_limb_t *rp)
1512 {
1513 return refmpn_udiv_qrnnd (rp, h, l, d);
1514 }
1515
1516 /* This little subroutine avoids some bad code generation from i386 gcc 3.0
1517 -fPIC -O2 -fomit-frame-pointer (%ebp being used uninitialized). */
1518 static mp_limb_t
1519 refmpn_divmod_1c_workaround (mp_ptr rp, mp_srcptr sp, mp_size_t size,
1520 mp_limb_t divisor, mp_limb_t carry)
1521 {
1522 mp_size_t i;
1523 mp_limb_t rem[1];
1524 for (i = size-1; i >= 0; i--)
1525 {
1526 rp[i] = refmpn_udiv_qrnnd (rem, carry,
1527 sp[i] << GMP_NAIL_BITS,
1528 divisor << GMP_NAIL_BITS);
1529 carry = *rem >> GMP_NAIL_BITS;
1530 }
1531 return carry;
1532 }
1533
1534 mp_limb_t
1535 refmpn_divmod_1c (mp_ptr rp, mp_srcptr sp, mp_size_t size,
1536 mp_limb_t divisor, mp_limb_t carry)
1537 {
1538 mp_ptr sp_orig;
1539 mp_ptr prod;
1540 mp_limb_t carry_orig;
1541
1542 ASSERT (refmpn_overlap_fullonly_p (rp, sp, size));
1543 ASSERT (size >= 0);
1544 ASSERT (carry < divisor);
1545 ASSERT_MPN (sp, size);
1546 ASSERT_LIMB (divisor);
1547 ASSERT_LIMB (carry);
1548
1549 if (size == 0)
1550 return carry;
1551
1552 sp_orig = refmpn_memdup_limbs (sp, size);
1553 prod = refmpn_malloc_limbs (size);
1554 carry_orig = carry;
1555
1556 carry = refmpn_divmod_1c_workaround (rp, sp, size, divisor, carry);
1557
1558 /* check by multiplying back */
1559 #if 0
1560 printf ("size=%ld divisor=0x%lX carry=0x%lX remainder=0x%lX\n",
1561 size, divisor, carry_orig, carry);
1562 mpn_trace("s",sp_copy,size);
1563 mpn_trace("r",rp,size);
1564 printf ("mul_1c %lX\n", refmpn_mul_1c (prod, rp, size, divisor, carry));
1565 mpn_trace("p",prod,size);
1566 #endif
1567 ASSERT (refmpn_mul_1c (prod, rp, size, divisor, carry) == carry_orig);
1568 ASSERT (refmpn_cmp (prod, sp_orig, size) == 0);
1569 free (sp_orig);
1570 free (prod);
1571
1572 return carry;
1573 }
1574
1575 mp_limb_t
1576 refmpn_divmod_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t divisor)
1577 {
1578 return refmpn_divmod_1c (rp, sp, size, divisor, CNST_LIMB(0));
1579 }
1580
1581
1582 mp_limb_t
1583 refmpn_mod_1c (mp_srcptr sp, mp_size_t size, mp_limb_t divisor,
1584 mp_limb_t carry)
1585 {
1586 mp_ptr p = refmpn_malloc_limbs (size);
1587 carry = refmpn_divmod_1c (p, sp, size, divisor, carry);
1588 free (p);
1589 return carry;
1590 }
1591
1592 mp_limb_t
1593 refmpn_mod_1 (mp_srcptr sp, mp_size_t size, mp_limb_t divisor)
1594 {
1595 return refmpn_mod_1c (sp, size, divisor, CNST_LIMB(0));
1596 }
1597
1598 mp_limb_t
1599 refmpn_preinv_mod_1 (mp_srcptr sp, mp_size_t size, mp_limb_t divisor,
1600 mp_limb_t inverse)
1601 {
1602 ASSERT (divisor & GMP_NUMB_HIGHBIT);
1603 ASSERT (inverse == refmpn_invert_limb (divisor));
1604 return refmpn_mod_1 (sp, size, divisor);
1605 }
1606
1607 /* This implementation will be rather slow, but has the advantage of being
1608 in a different style than the libgmp versions. */
1609 mp_limb_t
1610 refmpn_mod_34lsub1 (mp_srcptr p, mp_size_t n)
1611 {
1612 ASSERT ((GMP_NUMB_BITS % 4) == 0);
1613 return mpn_mod_1 (p, n, (CNST_LIMB(1) << (3 * GMP_NUMB_BITS / 4)) - 1);
1614 }
1615
1616
1617 mp_limb_t
1618 refmpn_divrem_1c (mp_ptr rp, mp_size_t xsize,
1619 mp_srcptr sp, mp_size_t size, mp_limb_t divisor,
1620 mp_limb_t carry)
1621 {
1622 mp_ptr z;
1623
1624 z = refmpn_malloc_limbs (xsize);
1625 refmpn_fill (z, xsize, CNST_LIMB(0));
1626
1627 carry = refmpn_divmod_1c (rp+xsize, sp, size, divisor, carry);
1628 carry = refmpn_divmod_1c (rp, z, xsize, divisor, carry);
1629
1630 free (z);
1631 return carry;
1632 }
1633
1634 mp_limb_t
1635 refmpn_divrem_1 (mp_ptr rp, mp_size_t xsize,
1636 mp_srcptr sp, mp_size_t size, mp_limb_t divisor)
1637 {
1638 return refmpn_divrem_1c (rp, xsize, sp, size, divisor, CNST_LIMB(0));
1639 }
1640
1641 mp_limb_t
1642 refmpn_preinv_divrem_1 (mp_ptr rp, mp_size_t xsize,
1643 mp_srcptr sp, mp_size_t size,
1644 mp_limb_t divisor, mp_limb_t inverse, unsigned shift)
1645 {
1646 ASSERT (size >= 0);
1647 ASSERT (shift == refmpn_count_leading_zeros (divisor));
1648 ASSERT (inverse == refmpn_invert_limb (divisor << shift));
1649
1650 return refmpn_divrem_1 (rp, xsize, sp, size, divisor);
1651 }
1652
1653 mp_limb_t
1654 refmpn_divrem_2 (mp_ptr qp, mp_size_t qxn,
1655 mp_ptr np, mp_size_t nn,
1656 mp_srcptr dp)
1657 {
1658 mp_ptr tp;
1659 mp_limb_t qh;
1660
1661 tp = refmpn_malloc_limbs (nn + qxn);
1662 refmpn_zero (tp, qxn);
1663 refmpn_copyi (tp + qxn, np, nn);
1664 qh = refmpn_sb_div_qr (qp, tp, nn + qxn, dp, 2);
1665 refmpn_copyi (np, tp, 2);
1666 free (tp);
1667 return qh;
1668 }
1669
1670 /* Inverse is floor((b*(b-d)-1) / d), per division by invariant integers
1671 paper, figure 8.1 m', where b=2^GMP_LIMB_BITS. Note that -d-1 < d
1672 since d has the high bit set. */
1673
1674 mp_limb_t
1675 refmpn_invert_limb (mp_limb_t d)
1676 {
1677 mp_limb_t r;
1678 ASSERT (d & GMP_LIMB_HIGHBIT);
1679 return refmpn_udiv_qrnnd (&r, -d-1, MP_LIMB_T_MAX, d);
1680 }
1681
1682 void
1683 refmpn_invert (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_ptr scratch)
1684 {
1685 mp_ptr qp, tp;
1686 TMP_DECL;
1687 TMP_MARK;
1688
1689 tp = TMP_ALLOC_LIMBS (2 * n);
1690 qp = TMP_ALLOC_LIMBS (n + 1);
1691
1692 MPN_ZERO (tp, 2 * n); mpn_sub_1 (tp, tp, 2 * n, 1);
1693
1694 refmpn_tdiv_qr (qp, rp, 0, tp, 2 * n, up, n);
1695 refmpn_copyi (rp, qp, n);
1696
1697 TMP_FREE;
1698 }
1699
1700 void
1701 refmpn_binvert (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_ptr scratch)
1702 {
1703 mp_ptr tp;
1704 mp_limb_t binv;
1705 TMP_DECL;
1706 TMP_MARK;
1707
1708 /* We use the library mpn_sbpi1_bdiv_q here, which isn't kosher in testing
1709 code. To make up for it, we check that the inverse is correct using a
1710 multiply. */
1711
1712 tp = TMP_ALLOC_LIMBS (2 * n);
1713
1714 MPN_ZERO (tp, n);
1715 tp[0] = 1;
1716 binvert_limb (binv, up[0]);
1717 mpn_sbpi1_bdiv_q (rp, tp, n, up, n, -binv);
1718
1719 refmpn_mul_n (tp, rp, up, n);
1720 ASSERT_ALWAYS (tp[0] == 1 && mpn_zero_p (tp + 1, n - 1));
1721
1722 TMP_FREE;
1723 }
1724
1725 /* The aim is to produce a dst quotient and return a remainder c, satisfying
1726 c*b^n + src-i == 3*dst, where i is the incoming carry.
1727
1728 Some value c==0, c==1 or c==2 will satisfy, so just try each.
1729
1730 If GMP_NUMB_BITS is even then 2^GMP_NUMB_BITS==1mod3 and a non-zero
1731 remainder from the first division attempt determines the correct
1732 remainder (3-c), but don't bother with that, since we can't guarantee
1733 anything about GMP_NUMB_BITS when using nails.
1734
1735 If the initial src-i produces a borrow then refmpn_sub_1 leaves a twos
1736 complement negative, ie. b^n+a-i, and the calculation produces c1
1737 satisfying c1*b^n + b^n+src-i == 3*dst, from which clearly c=c1+1. This
1738 means it's enough to just add any borrow back at the end.
1739
1740 A borrow only occurs when a==0 or a==1, and, by the same reasoning as in
1741 mpn/generic/diveby3.c, the c1 that results in those cases will only be 0
1742 or 1 respectively, so with 1 added the final return value is still in the
1743 prescribed range 0 to 2. */
1744
1745 mp_limb_t
1746 refmpn_divexact_by3c (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t carry)
1747 {
1748 mp_ptr spcopy;
1749 mp_limb_t c, cs;
1750
1751 ASSERT (refmpn_overlap_fullonly_p (rp, sp, size));
1752 ASSERT (size >= 1);
1753 ASSERT (carry <= 2);
1754 ASSERT_MPN (sp, size);
1755
1756 spcopy = refmpn_malloc_limbs (size);
1757 cs = refmpn_sub_1 (spcopy, sp, size, carry);
1758
1759 for (c = 0; c <= 2; c++)
1760 if (refmpn_divmod_1c (rp, spcopy, size, CNST_LIMB(3), c) == 0)
1761 goto done;
1762 ASSERT_FAIL (no value of c satisfies);
1763
1764 done:
1765 c += cs;
1766 ASSERT (c <= 2);
1767
1768 free (spcopy);
1769 return c;
1770 }
1771
1772 mp_limb_t
1773 refmpn_divexact_by3 (mp_ptr rp, mp_srcptr sp, mp_size_t size)
1774 {
1775 return refmpn_divexact_by3c (rp, sp, size, CNST_LIMB(0));
1776 }
1777
1778
1779 /* The same as mpn/generic/mul_basecase.c, but using refmpn functions. */
1780 void
1781 refmpn_mul_basecase (mp_ptr prodp,
1782 mp_srcptr up, mp_size_t usize,
1783 mp_srcptr vp, mp_size_t vsize)
1784 {
1785 mp_size_t i;
1786
1787 ASSERT (! refmpn_overlap_p (prodp, usize+vsize, up, usize));
1788 ASSERT (! refmpn_overlap_p (prodp, usize+vsize, vp, vsize));
1789 ASSERT (usize >= vsize);
1790 ASSERT (vsize >= 1);
1791 ASSERT_MPN (up, usize);
1792 ASSERT_MPN (vp, vsize);
1793
1794 prodp[usize] = refmpn_mul_1 (prodp, up, usize, vp[0]);
1795 for (i = 1; i < vsize; i++)
1796 prodp[usize+i] = refmpn_addmul_1 (prodp+i, up, usize, vp[i]);
1797 }
1798
1799
1800 /* The same as mpn/generic/mulmid_basecase.c, but using refmpn functions. */
1801 void
1802 refmpn_mulmid_basecase (mp_ptr rp,
1803 mp_srcptr up, mp_size_t un,
1804 mp_srcptr vp, mp_size_t vn)
1805 {
1806 mp_limb_t cy;
1807 mp_size_t i;
1808
1809 ASSERT (un >= vn);
1810 ASSERT (vn >= 1);
1811 ASSERT (! refmpn_overlap_p (rp, un - vn + 3, up, un));
1812 ASSERT (! refmpn_overlap_p (rp, un - vn + 3, vp, vn));
1813 ASSERT_MPN (up, un);
1814 ASSERT_MPN (vp, vn);
1815
1816 rp[un - vn + 1] = refmpn_mul_1 (rp, up + vn - 1, un - vn + 1, vp[0]);
1817 rp[un - vn + 2] = CNST_LIMB (0);
1818 for (i = 1; i < vn; i++)
1819 {
1820 cy = refmpn_addmul_1 (rp, up + vn - i - 1, un - vn + 1, vp[i]);
1821 cy = ref_addc_limb (&rp[un - vn + 1], rp[un - vn + 1], cy);
1822 cy = ref_addc_limb (&rp[un - vn + 2], rp[un - vn + 2], cy);
1823 ASSERT (cy == 0);
1824 }
1825 }
1826
1827 void
1828 refmpn_toom42_mulmid (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n,
1829 mp_ptr scratch)
1830 {
1831 refmpn_mulmid_basecase (rp, up, 2*n - 1, vp, n);
1832 }
1833
1834 void
1835 refmpn_mulmid_n (mp_ptr rp, mp_srcptr up, mp_srcptr vp, mp_size_t n)
1836 {
1837 /* FIXME: this could be made faster by using refmpn_mul and then subtracting
1838 off products near the middle product region boundary */
1839 refmpn_mulmid_basecase (rp, up, 2*n - 1, vp, n);
1840 }
1841
1842 void
1843 refmpn_mulmid (mp_ptr rp, mp_srcptr up, mp_size_t un,
1844 mp_srcptr vp, mp_size_t vn)
1845 {
1846 /* FIXME: this could be made faster by using refmpn_mul and then subtracting
1847 off products near the middle product region boundary */
1848 refmpn_mulmid_basecase (rp, up, un, vp, vn);
1849 }
1850
1851
1852
1853 #define TOOM3_THRESHOLD (MAX (MUL_TOOM33_THRESHOLD, SQR_TOOM3_THRESHOLD))
1854 #define TOOM4_THRESHOLD (MAX (MUL_TOOM44_THRESHOLD, SQR_TOOM4_THRESHOLD))
1855 #define TOOM6_THRESHOLD (MAX (MUL_TOOM6H_THRESHOLD, SQR_TOOM6_THRESHOLD))
1856 #if WANT_FFT
1857 #define FFT_THRESHOLD (MAX (MUL_FFT_THRESHOLD, SQR_FFT_THRESHOLD))
1858 #else
1859 #define FFT_THRESHOLD MP_SIZE_T_MAX /* don't use toom44 here */
1860 #endif
1861
1862 void
1863 refmpn_mul (mp_ptr wp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
1864 {
1865 mp_ptr tp, rp;
1866 mp_size_t tn;
1867
1868 if (vn < TOOM3_THRESHOLD)
1869 {
1870 /* In the mpn_mul_basecase and toom2 ranges, use our own mul_basecase. */
1871 if (vn != 0)
1872 refmpn_mul_basecase (wp, up, un, vp, vn);
1873 else
1874 MPN_ZERO (wp, un);
1875 return;
1876 }
1877
1878 MPN_ZERO (wp, vn);
1879 rp = refmpn_malloc_limbs (2 * vn);
1880
1881 if (vn < TOOM4_THRESHOLD)
1882 tn = mpn_toom22_mul_itch (vn, vn);
1883 else if (vn < TOOM6_THRESHOLD)
1884 tn = mpn_toom33_mul_itch (vn, vn);
1885 else if (vn < FFT_THRESHOLD)
1886 tn = mpn_toom44_mul_itch (vn, vn);
1887 else
1888 tn = mpn_toom6h_mul_itch (vn, vn);
1889 tp = refmpn_malloc_limbs (tn);
1890
1891 while (un >= vn)
1892 {
1893 if (vn < TOOM4_THRESHOLD)
1894 /* In the toom3 range, use mpn_toom22_mul. */
1895 mpn_toom22_mul (rp, up, vn, vp, vn, tp);
1896 else if (vn < TOOM6_THRESHOLD)
1897 /* In the toom4 range, use mpn_toom33_mul. */
1898 mpn_toom33_mul (rp, up, vn, vp, vn, tp);
1899 else if (vn < FFT_THRESHOLD)
1900 /* In the toom6 range, use mpn_toom44_mul. */
1901 mpn_toom44_mul (rp, up, vn, vp, vn, tp);
1902 else
1903 /* For the largest operands, use mpn_toom6h_mul. */
1904 mpn_toom6h_mul (rp, up, vn, vp, vn, tp);
1905
1906 ASSERT_NOCARRY (refmpn_add (wp, rp, 2 * vn, wp, vn));
1907 wp += vn;
1908
1909 up += vn;
1910 un -= vn;
1911 }
1912
1913 free (tp);
1914
1915 if (un != 0)
1916 {
1917 refmpn_mul (rp, vp, vn, up, un);
1918 ASSERT_NOCARRY (refmpn_add (wp, rp, un + vn, wp, vn));
1919 }
1920 free (rp);
1921 }
1922
1923 void
1924 refmpn_mul_n (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
1925 {
1926 refmpn_mul (prodp, up, size, vp, size);
1927 }
1928
1929 void
1930 refmpn_mullo_n (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
1931 {
1932 mp_ptr tp = refmpn_malloc_limbs (2*size);
1933 refmpn_mul (tp, up, size, vp, size);
1934 refmpn_copyi (prodp, tp, size);
1935 free (tp);
1936 }
1937
1938 void
1939 refmpn_sqr (mp_ptr dst, mp_srcptr src, mp_size_t size)
1940 {
1941 refmpn_mul (dst, src, size, src, size);
1942 }
1943
1944 void
1945 refmpn_sqrlo (mp_ptr dst, mp_srcptr src, mp_size_t size)
1946 {
1947 refmpn_mullo_n (dst, src, src, size);
1948 }
1949
1950 /* Allowing usize<vsize, usize==0 or vsize==0. */
1951 void
1952 refmpn_mul_any (mp_ptr prodp,
1953 mp_srcptr up, mp_size_t usize,
1954 mp_srcptr vp, mp_size_t vsize)
1955 {
1956 ASSERT (! refmpn_overlap_p (prodp, usize+vsize, up, usize));
1957 ASSERT (! refmpn_overlap_p (prodp, usize+vsize, vp, vsize));
1958 ASSERT (usize >= 0);
1959 ASSERT (vsize >= 0);
1960 ASSERT_MPN (up, usize);
1961 ASSERT_MPN (vp, vsize);
1962
1963 if (usize == 0)
1964 {
1965 refmpn_fill (prodp, vsize, CNST_LIMB(0));
1966 return;
1967 }
1968
1969 if (vsize == 0)
1970 {
1971 refmpn_fill (prodp, usize, CNST_LIMB(0));
1972 return;
1973 }
1974
1975 if (usize >= vsize)
1976 refmpn_mul (prodp, up, usize, vp, vsize);
1977 else
1978 refmpn_mul (prodp, vp, vsize, up, usize);
1979 }
1980
1981
1982 mp_limb_t
1983 refmpn_gcd_11 (mp_limb_t x, mp_limb_t y)
1984 {
1985 /* The non-ref function also requires input operands to be odd, but
1986 below refmpn_gcd_1 doesn't guarantee that. */
1987 ASSERT (x > 0);
1988 ASSERT (y > 0);
1989 do
1990 {
1991 while ((x & 1) == 0) x >>= 1;
1992 while ((y & 1) == 0) y >>= 1;
1993
1994 if (x < y)
1995 MP_LIMB_T_SWAP (x, y);
1996
1997 x -= y;
1998 }
1999 while (x != 0);
2000
2001 return y;
2002 }
2003
2004 mp_double_limb_t
2005 refmpn_gcd_22 (mp_limb_t x1, mp_limb_t x0, mp_limb_t y1, mp_limb_t y0)
2006 {
2007 ASSERT ((x0 & 1) != 0);
2008 ASSERT ((y0 & 1) != 0);
2009 mp_double_limb_t g;
2010 mp_limb_t cy;
2011
2012 do
2013 {
2014 while ((x0 & 1) == 0)
2015 {
2016 x0 = (x1 << (GMP_NUMB_BITS - 1)) | (x0 >> 1);
2017 x1 >>= 1;
2018 }
2019 while ((y0 & 1) == 0)
2020 {
2021 y0 = (y1 << (GMP_NUMB_BITS - 1)) | (y0 >> 1);
2022 y1 >>= 1;
2023 }
2024
2025
2026 if (x1 < y1 || (x1 == y1 && x0 < y0))
2027 {
2028 mp_limb_t t;
2029 t = x1; x1 = y1; y1 = t;
2030 t = x0; x0 = y0; y0 = t;
2031 }
2032
2033 cy = (x0 < y0);
2034 x0 -= y0;
2035 x1 -= y1 + cy;
2036 }
2037 while ((x1 | x0) != 0);
2038
2039 g.d1 = y1;
2040 g.d0 = y0;
2041 return g;
2042 }
2043
2044 mp_limb_t
2045 refmpn_gcd_1 (mp_srcptr xp, mp_size_t xsize, mp_limb_t y)
2046 {
2047 mp_limb_t x;
2048 int twos;
2049
2050 ASSERT (y != 0);
2051 ASSERT (! refmpn_zero_p (xp, xsize));
2052 ASSERT_MPN (xp, xsize);
2053 ASSERT_LIMB (y);
2054
2055 x = refmpn_mod_1 (xp, xsize, y);
2056 if (x == 0)
2057 return y;
2058
2059 twos = 0;
2060 while ((x & 1) == 0 && (y & 1) == 0)
2061 {
2062 x >>= 1;
2063 y >>= 1;
2064 twos++;
2065 }
2066
2067 return refmpn_gcd_11 (x, y) << twos;
2068 }
2069
2070
2071 /* Based on the full limb x, not nails. */
2072 unsigned
2073 refmpn_count_leading_zeros (mp_limb_t x)
2074 {
2075 unsigned n = 0;
2076
2077 ASSERT (x != 0);
2078
2079 while ((x & GMP_LIMB_HIGHBIT) == 0)
2080 {
2081 x <<= 1;
2082 n++;
2083 }
2084 return n;
2085 }
2086
2087 /* Full limbs allowed, not limited to nails. */
2088 unsigned
2089 refmpn_count_trailing_zeros (mp_limb_t x)
2090 {
2091 unsigned n = 0;
2092
2093 ASSERT (x != 0);
2094 ASSERT_LIMB (x);
2095
2096 while ((x & 1) == 0)
2097 {
2098 x >>= 1;
2099 n++;
2100 }
2101 return n;
2102 }
2103
2104 /* Strip factors of two (low zero bits) from {p,size} by right shifting.
2105 The return value is the number of twos stripped. */
2106 mp_size_t
2107 refmpn_strip_twos (mp_ptr p, mp_size_t size)
2108 {
2109 mp_size_t limbs;
2110 unsigned shift;
2111
2112 ASSERT (size >= 1);
2113 ASSERT (! refmpn_zero_p (p, size));
2114 ASSERT_MPN (p, size);
2115
2116 for (limbs = 0; p[0] == 0; limbs++)
2117 {
2118 refmpn_copyi (p, p+1, size-1);
2119 p[size-1] = 0;
2120 }
2121
2122 shift = refmpn_count_trailing_zeros (p[0]);
2123 if (shift)
2124 refmpn_rshift (p, p, size, shift);
2125
2126 return limbs*GMP_NUMB_BITS + shift;
2127 }
2128
2129 mp_limb_t
2130 refmpn_gcd (mp_ptr gp, mp_ptr xp, mp_size_t xsize, mp_ptr yp, mp_size_t ysize)
2131 {
2132 int cmp;
2133
2134 ASSERT (ysize >= 1);
2135 ASSERT (xsize >= ysize);
2136 ASSERT ((xp[0] & 1) != 0);
2137 ASSERT ((yp[0] & 1) != 0);
2138 /* ASSERT (xp[xsize-1] != 0); */ /* don't think x needs to be odd */
2139 ASSERT (yp[ysize-1] != 0);
2140 ASSERT (refmpn_overlap_fullonly_p (gp, xp, xsize));
2141 ASSERT (refmpn_overlap_fullonly_p (gp, yp, ysize));
2142 ASSERT (! refmpn_overlap_p (xp, xsize, yp, ysize));
2143 if (xsize == ysize)
2144 ASSERT (refmpn_msbone (xp[xsize-1]) >= refmpn_msbone (yp[ysize-1]));
2145 ASSERT_MPN (xp, xsize);
2146 ASSERT_MPN (yp, ysize);
2147
2148 refmpn_strip_twos (xp, xsize);
2149 MPN_NORMALIZE (xp, xsize);
2150 MPN_NORMALIZE (yp, ysize);
2151
2152 for (;;)
2153 {
2154 cmp = refmpn_cmp_twosizes (xp, xsize, yp, ysize);
2155 if (cmp == 0)
2156 break;
2157 if (cmp < 0)
2158 MPN_PTR_SWAP (xp,xsize, yp,ysize);
2159
2160 ASSERT_NOCARRY (refmpn_sub (xp, xp, xsize, yp, ysize));
2161
2162 refmpn_strip_twos (xp, xsize);
2163 MPN_NORMALIZE (xp, xsize);
2164 }
2165
2166 refmpn_copyi (gp, xp, xsize);
2167 return xsize;
2168 }
2169
2170 unsigned long
2171 ref_popc_limb (mp_limb_t src)
2172 {
2173 unsigned long count;
2174 int i;
2175
2176 count = 0;
2177 for (i = 0; i < GMP_LIMB_BITS; i++)
2178 {
2179 count += (src & 1);
2180 src >>= 1;
2181 }
2182 return count;
2183 }
2184
2185 unsigned long
2186 refmpn_popcount (mp_srcptr sp, mp_size_t size)
2187 {
2188 unsigned long count = 0;
2189 mp_size_t i;
2190
2191 ASSERT (size >= 0);
2192 ASSERT_MPN (sp, size);
2193
2194 for (i = 0; i < size; i++)
2195 count += ref_popc_limb (sp[i]);
2196 return count;
2197 }
2198
2199 unsigned long
2200 refmpn_hamdist (mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
2201 {
2202 mp_ptr d;
2203 unsigned long count;
2204
2205 ASSERT (size >= 0);
2206 ASSERT_MPN (s1p, size);
2207 ASSERT_MPN (s2p, size);
2208
2209 if (size == 0)
2210 return 0;
2211
2212 d = refmpn_malloc_limbs (size);
2213 refmpn_xor_n (d, s1p, s2p, size);
2214 count = refmpn_popcount (d, size);
2215 free (d);
2216 return count;
2217 }
2218
2219
2220 /* set r to a%d */
2221 void
2222 refmpn_mod2 (mp_limb_t r[2], const mp_limb_t a[2], const mp_limb_t d[2])
2223 {
2224 mp_limb_t D[2];
2225 int n;
2226
2227 ASSERT (! refmpn_overlap_p (r, (mp_size_t) 2, d, (mp_size_t) 2));
2228 ASSERT_MPN (a, 2);
2229 ASSERT_MPN (d, 2);
2230
2231 D[1] = d[1], D[0] = d[0];
2232 r[1] = a[1], r[0] = a[0];
2233 n = 0;
2234
2235 for (;;)
2236 {
2237 if (D[1] & GMP_NUMB_HIGHBIT)
2238 break;
2239 if (refmpn_cmp (r, D, (mp_size_t) 2) <= 0)
2240 break;
2241 refmpn_lshift (D, D, (mp_size_t) 2, 1);
2242 n++;
2243 ASSERT (n <= GMP_NUMB_BITS);
2244 }
2245
2246 while (n >= 0)
2247 {
2248 if (refmpn_cmp (r, D, (mp_size_t) 2) >= 0)
2249 ASSERT_NOCARRY (refmpn_sub_n (r, r, D, (mp_size_t) 2));
2250 refmpn_rshift (D, D, (mp_size_t) 2, 1);
2251 n--;
2252 }
2253
2254 ASSERT (refmpn_cmp (r, d, (mp_size_t) 2) < 0);
2255 }
2256
2257
2258
2259 /* Similar to the old mpn/generic/sb_divrem_mn.c, but somewhat simplified, in
2260 particular the trial quotient is allowed to be 2 too big. */
2261 mp_limb_t
2262 refmpn_sb_div_qr (mp_ptr qp,
2263 mp_ptr np, mp_size_t nsize,
2264 mp_srcptr dp, mp_size_t dsize)
2265 {
2266 mp_limb_t retval = 0;
2267 mp_size_t i;
2268 mp_limb_t d1 = dp[dsize-1];
2269 mp_ptr np_orig = refmpn_memdup_limbs (np, nsize);
2270
2271 ASSERT (nsize >= dsize);
2272 /* ASSERT (dsize > 2); */
2273 ASSERT (dsize >= 2);
2274 ASSERT (dp[dsize-1] & GMP_NUMB_HIGHBIT);
2275 ASSERT (! refmpn_overlap_p (qp, nsize-dsize, np, nsize) || qp+dsize >= np);
2276 ASSERT_MPN (np, nsize);
2277 ASSERT_MPN (dp, dsize);
2278
2279 i = nsize-dsize;
2280 if (refmpn_cmp (np+i, dp, dsize) >= 0)
2281 {
2282 ASSERT_NOCARRY (refmpn_sub_n (np+i, np+i, dp, dsize));
2283 retval = 1;
2284 }
2285
2286 for (i--; i >= 0; i--)
2287 {
2288 mp_limb_t n0 = np[i+dsize];
2289 mp_limb_t n1 = np[i+dsize-1];
2290 mp_limb_t q, dummy_r;
2291
2292 ASSERT (n0 <= d1);
2293 if (n0 == d1)
2294 q = GMP_NUMB_MAX;
2295 else
2296 q = refmpn_udiv_qrnnd (&dummy_r, n0, n1 << GMP_NAIL_BITS,
2297 d1 << GMP_NAIL_BITS);
2298
2299 n0 -= refmpn_submul_1 (np+i, dp, dsize, q);
2300 ASSERT (n0 == 0 || n0 == MP_LIMB_T_MAX);
2301 if (n0)
2302 {
2303 q--;
2304 if (! refmpn_add_n (np+i, np+i, dp, dsize))
2305 {
2306 q--;
2307 ASSERT_CARRY (refmpn_add_n (np+i, np+i, dp, dsize));
2308 }
2309 }
2310 np[i+dsize] = 0;
2311
2312 qp[i] = q;
2313 }
2314
2315 /* remainder < divisor */
2316 #if 0 /* ASSERT triggers gcc 4.2.1 bug */
2317 ASSERT (refmpn_cmp (np, dp, dsize) < 0);
2318 #endif
2319
2320 /* multiply back to original */
2321 {
2322 mp_ptr mp = refmpn_malloc_limbs (nsize);
2323
2324 refmpn_mul_any (mp, qp, nsize-dsize, dp, dsize);
2325 if (retval)
2326 ASSERT_NOCARRY (refmpn_add_n (mp+nsize-dsize,mp+nsize-dsize, dp, dsize));
2327 ASSERT_NOCARRY (refmpn_add (mp, mp, nsize, np, dsize));
2328 ASSERT (refmpn_cmp (mp, np_orig, nsize) == 0);
2329
2330 free (mp);
2331 }
2332
2333 free (np_orig);
2334 return retval;
2335 }
2336
2337 /* Similar to the old mpn/generic/sb_divrem_mn.c, but somewhat simplified, in
2338 particular the trial quotient is allowed to be 2 too big. */
2339 void
2340 refmpn_tdiv_qr (mp_ptr qp, mp_ptr rp, mp_size_t qxn,
2341 mp_ptr np, mp_size_t nsize,
2342 mp_srcptr dp, mp_size_t dsize)
2343 {
2344 ASSERT (qxn == 0);
2345 ASSERT_MPN (np, nsize);
2346 ASSERT_MPN (dp, dsize);
2347 ASSERT (dsize > 0);
2348 ASSERT (dp[dsize-1] != 0);
2349
2350 if (dsize == 1)
2351 {
2352 rp[0] = refmpn_divmod_1 (qp, np, nsize, dp[0]);
2353 return;
2354 }
2355 else
2356 {
2357 mp_ptr n2p = refmpn_malloc_limbs (nsize+1);
2358 mp_ptr d2p = refmpn_malloc_limbs (dsize);
2359 int norm = refmpn_count_leading_zeros (dp[dsize-1]) - GMP_NAIL_BITS;
2360
2361 n2p[nsize] = refmpn_lshift_or_copy (n2p, np, nsize, norm);
2362 ASSERT_NOCARRY (refmpn_lshift_or_copy (d2p, dp, dsize, norm));
2363
2364 refmpn_sb_div_qr (qp, n2p, nsize+1, d2p, dsize);
2365 refmpn_rshift_or_copy (rp, n2p, dsize, norm);
2366
2367 /* ASSERT (refmpn_zero_p (tp+dsize, nsize-dsize)); */
2368 free (n2p);
2369 free (d2p);
2370 }
2371 }
2372
2373 mp_limb_t
2374 refmpn_redc_1 (mp_ptr rp, mp_ptr up, mp_srcptr mp, mp_size_t n, mp_limb_t invm)
2375 {
2376 mp_size_t j;
2377 mp_limb_t cy;
2378
2379 ASSERT_MPN (up, 2*n);
2380 /* ASSERT about directed overlap rp, up */
2381 /* ASSERT about overlap rp, mp */
2382 /* ASSERT about overlap up, mp */
2383
2384 for (j = n - 1; j >= 0; j--)
2385 {
2386 up[0] = refmpn_addmul_1 (up, mp, n, (up[0] * invm) & GMP_NUMB_MASK);
2387 up++;
2388 }
2389 cy = mpn_add_n (rp, up, up - n, n);
2390 return cy;
2391 }
2392
2393 size_t
2394 refmpn_get_str (unsigned char *dst, int base, mp_ptr src, mp_size_t size)
2395 {
2396 unsigned char *d;
2397 size_t dsize;
2398
2399 ASSERT (size >= 0);
2400 ASSERT (base >= 2);
2401 ASSERT (base < numberof (mp_bases));
2402 ASSERT (size == 0 || src[size-1] != 0);
2403 ASSERT_MPN (src, size);
2404
2405 MPN_SIZEINBASE (dsize, src, size, base);
2406 ASSERT (dsize >= 1);
2407 ASSERT (! byte_overlap_p (dst, (mp_size_t) dsize, src, size * GMP_LIMB_BYTES));
2408
2409 if (size == 0)
2410 {
2411 dst[0] = 0;
2412 return 1;
2413 }
2414
2415 /* don't clobber input for power of 2 bases */
2416 if (POW2_P (base))
2417 src = refmpn_memdup_limbs (src, size);
2418
2419 d = dst + dsize;
2420 do
2421 {
2422 d--;
2423 ASSERT (d >= dst);
2424 *d = refmpn_divrem_1 (src, (mp_size_t) 0, src, size, (mp_limb_t) base);
2425 size -= (src[size-1] == 0);
2426 }
2427 while (size != 0);
2428
2429 /* Move result back and decrement dsize if we didn't generate
2430 the maximum possible digits. */
2431 if (d != dst)
2432 {
2433 size_t i;
2434 dsize -= d - dst;
2435 for (i = 0; i < dsize; i++)
2436 dst[i] = d[i];
2437 }
2438
2439 if (POW2_P (base))
2440 free (src);
2441
2442 return dsize;
2443 }
2444
2445
2446 mp_limb_t
2447 ref_bswap_limb (mp_limb_t src)
2448 {
2449 mp_limb_t dst;
2450 int i;
2451
2452 dst = 0;
2453 for (i = 0; i < GMP_LIMB_BYTES; i++)
2454 {
2455 dst = (dst << 8) + (src & 0xFF);
2456 src >>= 8;
2457 }
2458 return dst;
2459 }
2460
2461
2462 /* These random functions are mostly for transitional purposes while adding
2463 nail support, since they're independent of the normal mpn routines. They
2464 can probably be removed when those normal routines are reliable, though
2465 perhaps something independent would still be useful at times. */
2466
2467 #if GMP_LIMB_BITS == 32
2468 #define RAND_A CNST_LIMB(0x29CF535)
2469 #endif
2470 #if GMP_LIMB_BITS == 64
2471 #define RAND_A CNST_LIMB(0xBAECD515DAF0B49D)
2472 #endif
2473
2474 mp_limb_t refmpn_random_seed;
2475
2476 mp_limb_t
2477 refmpn_random_half (void)
2478 {
2479 refmpn_random_seed = refmpn_random_seed * RAND_A + 1;
2480 return (refmpn_random_seed >> GMP_LIMB_BITS/2);
2481 }
2482
2483 mp_limb_t
2484 refmpn_random_limb (void)
2485 {
2486 return ((refmpn_random_half () << (GMP_LIMB_BITS/2))
2487 | refmpn_random_half ()) & GMP_NUMB_MASK;
2488 }
2489
2490 void
2491 refmpn_random (mp_ptr ptr, mp_size_t size)
2492 {
2493 mp_size_t i;
2494 if (GMP_NAIL_BITS == 0)
2495 {
2496 mpn_random (ptr, size);
2497 return;
2498 }
2499
2500 for (i = 0; i < size; i++)
2501 ptr[i] = refmpn_random_limb ();
2502 }
2503
2504 void
2505 refmpn_random2 (mp_ptr ptr, mp_size_t size)
2506 {
2507 mp_size_t i;
2508 mp_limb_t bit, mask, limb;
2509 int run;
2510
2511 if (GMP_NAIL_BITS == 0)
2512 {
2513 mpn_random2 (ptr, size);
2514 return;
2515 }
2516
2517 #define RUN_MODULUS 32
2518
2519 /* start with ones at a random pos in the high limb */
2520 bit = CNST_LIMB(1) << (refmpn_random_half () % GMP_NUMB_BITS);
2521 mask = 0;
2522 run = 0;
2523
2524 for (i = size-1; i >= 0; i--)
2525 {
2526 limb = 0;
2527 do
2528 {
2529 if (run == 0)
2530 {
2531 run = (refmpn_random_half () % RUN_MODULUS) + 1;
2532 mask = ~mask;
2533 }
2534
2535 limb |= (bit & mask);
2536 bit >>= 1;
2537 run--;
2538 }
2539 while (bit != 0);
2540
2541 ptr[i] = limb;
2542 bit = GMP_NUMB_HIGHBIT;
2543 }
2544 }
2545
2546 /* This is a simple bitwise algorithm working high to low across "s" and
2547 testing each time whether setting the bit would make s^2 exceed n. */
2548 mp_size_t
2549 refmpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr np, mp_size_t nsize)
2550 {
2551 mp_ptr tp, dp;
2552 mp_size_t ssize, talloc, tsize, dsize, ret, ilimbs;
2553 unsigned ibit;
2554 long i;
2555 mp_limb_t c;
2556
2557 ASSERT (nsize >= 0);
2558
2559 /* If n==0, then s=0 and r=0. */
2560 if (nsize == 0)
2561 return 0;
2562
2563 ASSERT (np[nsize - 1] != 0);
2564 ASSERT (rp == NULL || MPN_SAME_OR_SEPARATE_P (np, rp, nsize));
2565 ASSERT (rp == NULL || ! MPN_OVERLAP_P (sp, (nsize + 1) / 2, rp, nsize));
2566 ASSERT (! MPN_OVERLAP_P (sp, (nsize + 1) / 2, np, nsize));
2567
2568 /* root */
2569 ssize = (nsize+1)/2;
2570 refmpn_zero (sp, ssize);
2571
2572 /* the remainder so far */
2573 dp = refmpn_memdup_limbs (np, nsize);
2574 dsize = nsize;
2575
2576 /* temporary */
2577 talloc = 2*ssize + 1;
2578 tp = refmpn_malloc_limbs (talloc);
2579
2580 for (i = GMP_NUMB_BITS * ssize - 1; i >= 0; i--)
2581 {
2582 /* t = 2*s*2^i + 2^(2*i), being the amount s^2 will increase by if 2^i
2583 is added to it */
2584
2585 ilimbs = (i+1) / GMP_NUMB_BITS;
2586 ibit = (i+1) % GMP_NUMB_BITS;
2587 refmpn_zero (tp, ilimbs);
2588 c = refmpn_lshift_or_copy (tp+ilimbs, sp, ssize, ibit);
2589 tsize = ilimbs + ssize;
2590 tp[tsize] = c;
2591 tsize += (c != 0);
2592
2593 ilimbs = (2*i) / GMP_NUMB_BITS;
2594 ibit = (2*i) % GMP_NUMB_BITS;
2595 if (ilimbs + 1 > tsize)
2596 {
2597 refmpn_zero_extend (tp, tsize, ilimbs + 1);
2598 tsize = ilimbs + 1;
2599 }
2600 c = refmpn_add_1 (tp+ilimbs, tp+ilimbs, tsize-ilimbs,
2601 CNST_LIMB(1) << ibit);
2602 ASSERT (tsize < talloc);
2603 tp[tsize] = c;
2604 tsize += (c != 0);
2605
2606 if (refmpn_cmp_twosizes (dp, dsize, tp, tsize) >= 0)
2607 {
2608 /* set this bit in s and subtract from the remainder */
2609 refmpn_setbit (sp, i);
2610
2611 ASSERT_NOCARRY (refmpn_sub_n (dp, dp, tp, dsize));
2612 dsize = refmpn_normalize (dp, dsize);
2613 }
2614 }
2615
2616 if (rp == NULL)
2617 {
2618 ret = ! refmpn_zero_p (dp, dsize);
2619 }
2620 else
2621 {
2622 ASSERT (dsize == 0 || dp[dsize-1] != 0);
2623 refmpn_copy (rp, dp, dsize);
2624 ret = dsize;
2625 }
2626
2627 free (dp);
2628 free (tp);
2629 return ret;
2630 }