Annotation of OpenXM_contrib/gmp/tests/refmpn.c, Revision 1.1.1.1
1.1 ohara 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, 1997, 1998, 1999, 2000, 2001, 2002 Free Software Foundation,
5: Inc.
6:
7: This file is part of the GNU MP Library.
8:
9: The GNU MP Library is free software; you can redistribute it and/or modify
10: it under the terms of the GNU Lesser General Public License as published by
11: the Free Software Foundation; either version 2.1 of the License, or (at your
12: option) any later version.
13:
14: The GNU MP Library is distributed in the hope that it will be useful, but
15: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17: License for more details.
18:
19: You should have received a copy of the GNU Lesser General Public License
20: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
21: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
22: MA 02111-1307, USA. */
23:
24:
25: /* Most routines have assertions representing what the mpn routines are
26: supposed to accept. Many of these reference routines do sensible things
27: outside these ranges (eg. for size==0), but the assertions are present to
28: pick up bad parameters passed here that are about to be passed the same
29: to a real mpn routine being compared. */
30:
31: /* always do assertion checking */
32: #define WANT_ASSERT 1
33:
34: #include <stdio.h> /* for NULL */
35: #include <stdlib.h> /* for malloc */
36:
37: #include "gmp.h"
38: #include "gmp-impl.h"
39: #include "longlong.h"
40:
41: #include "tests.h"
42:
43:
44:
45: /* Return non-zero if regions {xp,xsize} and {yp,ysize} overlap, with sizes
46: in bytes. */
47: int
48: byte_overlap_p (const void *v_xp, mp_size_t xsize,
49: const void *v_yp, mp_size_t ysize)
50: {
51: const char *xp = v_xp;
52: const char *yp = v_yp;
53:
54: ASSERT (xsize >= 0);
55: ASSERT (ysize >= 0);
56:
57: /* no wraparounds */
58: ASSERT (xp+xsize >= xp);
59: ASSERT (yp+ysize >= yp);
60:
61: if (xp + xsize <= yp)
62: return 0;
63:
64: if (yp + ysize <= xp)
65: return 0;
66:
67: return 1;
68: }
69:
70: /* Return non-zero if limb regions {xp,xsize} and {yp,ysize} overlap. */
71: int
72: refmpn_overlap_p (mp_srcptr xp, mp_size_t xsize, mp_srcptr yp, mp_size_t ysize)
73: {
74: return byte_overlap_p (xp, xsize * BYTES_PER_MP_LIMB,
75: yp, ysize * BYTES_PER_MP_LIMB);
76: }
77:
78: /* Check overlap for a routine defined to work low to high. */
79: int
80: refmpn_overlap_low_to_high_p (mp_srcptr dst, mp_srcptr src, mp_size_t size)
81: {
82: return (dst <= src || ! refmpn_overlap_p (dst, size, src, size));
83: }
84:
85: /* Check overlap for a routine defined to work high to low. */
86: int
87: refmpn_overlap_high_to_low_p (mp_srcptr dst, mp_srcptr src, mp_size_t size)
88: {
89: return (dst >= src || ! refmpn_overlap_p (dst, size, src, size));
90: }
91:
92: /* Check overlap for a standard routine requiring equal or separate. */
93: int
94: refmpn_overlap_fullonly_p (mp_srcptr dst, mp_srcptr src, mp_size_t size)
95: {
96: return (dst == src || ! refmpn_overlap_p (dst, size, src, size));
97: }
98: int
99: refmpn_overlap_fullonly_two_p (mp_srcptr dst, mp_srcptr src1, mp_srcptr src2,
100: mp_size_t size)
101: {
102: return (refmpn_overlap_fullonly_p (dst, src1, size)
103: && refmpn_overlap_fullonly_p (dst, src2, size));
104: }
105:
106:
107: mp_ptr
108: refmpn_malloc_limbs (mp_size_t size)
109: {
110: mp_ptr p;
111: ASSERT (size >= 0);
112: if (size == 0)
113: size = 1;
114: p = (mp_ptr) malloc (size * BYTES_PER_MP_LIMB);
115: ASSERT (p != NULL);
116: return p;
117: }
118:
119: mp_ptr
120: refmpn_memdup_limbs (mp_srcptr ptr, mp_size_t size)
121: {
122: mp_ptr p;
123: p = refmpn_malloc_limbs (size);
124: refmpn_copyi (p, ptr, size);
125: return p;
126: }
127:
128: /* malloc n limbs on a multiple of m bytes boundary */
129: mp_ptr
130: refmpn_malloc_limbs_aligned (size_t n, size_t m)
131: {
132: return (mp_ptr) align_pointer (refmpn_malloc_limbs (n + m-1), m);
133: }
134:
135:
136: void
137: refmpn_fill (mp_ptr ptr, mp_size_t size, mp_limb_t value)
138: {
139: mp_size_t i;
140: ASSERT (size >= 0);
141: for (i = 0; i < size; i++)
142: ptr[i] = value;
143: }
144:
145: void
146: refmpn_zero (mp_ptr ptr, mp_size_t size)
147: {
148: refmpn_fill (ptr, size, CNST_LIMB(0));
149: }
150:
151: int
152: refmpn_zero_p (mp_srcptr ptr, mp_size_t size)
153: {
154: mp_size_t i;
155: for (i = 0; i < size; i++)
156: if (ptr[i] != 0)
157: return 0;
158: return 1;
159: }
160:
161: /* the highest one bit in x */
162: mp_limb_t
163: refmpn_msbone (mp_limb_t x)
164: {
165: mp_limb_t n = (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1);
166:
167: while (n != 0)
168: {
169: if (x & n)
170: break;
171: n >>= 1;
172: }
173: return n;
174: }
175:
176: /* a mask of the highest one bit plus and all bits below */
177: mp_limb_t
178: refmpn_msbone_mask (mp_limb_t x)
179: {
180: if (x == 0)
181: return 0;
182:
183: return (refmpn_msbone (x) << 1) - 1;
184: }
185:
186:
187: void
188: refmpn_setbit (mp_ptr ptr, unsigned long bit)
189: {
190: ptr[bit/GMP_NUMB_BITS] |= CNST_LIMB(1) << (bit%GMP_NUMB_BITS);
191: }
192:
193: void
194: refmpn_clrbit (mp_ptr ptr, unsigned long bit)
195: {
196: ptr[bit/GMP_NUMB_BITS] &= ~ (CNST_LIMB(1) << (bit%GMP_NUMB_BITS));
197: }
198:
199: #define REFMPN_TSTBIT(ptr,bit) \
200: (((ptr)[(bit)/GMP_NUMB_BITS] & (CNST_LIMB(1) << ((bit)%GMP_NUMB_BITS))) != 0)
201:
202: int
203: refmpn_tstbit (mp_srcptr ptr, unsigned long bit)
204: {
205: return REFMPN_TSTBIT (ptr, bit);
206: }
207:
208: unsigned long
209: refmpn_scan0 (mp_srcptr ptr, unsigned long bit)
210: {
211: while (REFMPN_TSTBIT (ptr, bit) != 0)
212: bit++;
213: return bit;
214: }
215:
216: unsigned long
217: refmpn_scan1 (mp_srcptr ptr, unsigned long bit)
218: {
219: while (REFMPN_TSTBIT (ptr, bit) == 0)
220: bit++;
221: return bit;
222: }
223:
224: void
225: refmpn_copy (mp_ptr rp, mp_srcptr sp, mp_size_t size)
226: {
227: ASSERT (refmpn_overlap_fullonly_p (rp, sp, size));
228: refmpn_copyi (rp, sp, size);
229: }
230:
231: void
232: refmpn_copyi (mp_ptr rp, mp_srcptr sp, mp_size_t size)
233: {
234: mp_size_t i;
235:
236: ASSERT (refmpn_overlap_low_to_high_p (rp, sp, size));
237: ASSERT (size >= 0);
238:
239: for (i = 0; i < size; i++)
240: rp[i] = sp[i];
241: }
242:
243: void
244: refmpn_copyd (mp_ptr rp, mp_srcptr sp, mp_size_t size)
245: {
246: mp_size_t i;
247:
248: ASSERT (refmpn_overlap_high_to_low_p (rp, sp, size));
249: ASSERT (size >= 0);
250:
251: for (i = size-1; i >= 0; i--)
252: rp[i] = sp[i];
253: }
254:
255: void
256: refmpn_com_n (mp_ptr rp, mp_srcptr sp, mp_size_t size)
257: {
258: mp_size_t i;
259:
260: ASSERT (refmpn_overlap_fullonly_p (rp, sp, size));
261: ASSERT (size >= 1);
262: ASSERT_MPN (sp, size);
263:
264: for (i = 0; i < size; i++)
265: rp[i] = sp[i] ^ GMP_NUMB_MASK;
266: }
267:
268:
269: int
270: refmpn_cmp (mp_srcptr xp, mp_srcptr yp, mp_size_t size)
271: {
272: mp_size_t i;
273:
274: ASSERT (size >= 1);
275: ASSERT_MPN (xp, size);
276: ASSERT_MPN (yp, size);
277:
278: for (i = size-1; i >= 0; i--)
279: {
280: if (xp[i] > yp[i]) return 1;
281: if (xp[i] < yp[i]) return -1;
282: }
283: return 0;
284: }
285:
286: int
287: refmpn_cmp_allowzero (mp_srcptr xp, mp_srcptr yp, mp_size_t size)
288: {
289: if (size == 0)
290: return 0;
291: else
292: return refmpn_cmp (xp, yp, size);
293: }
294:
295: int
296: refmpn_cmp_twosizes (mp_srcptr xp, mp_size_t xsize,
297: mp_srcptr yp, mp_size_t ysize)
298: {
299: int opp, cmp;
300:
301: ASSERT_MPN (xp, xsize);
302: ASSERT_MPN (yp, ysize);
303:
304: opp = (xsize < ysize);
305: if (opp)
306: MPN_SRCPTR_SWAP (xp,xsize, yp,ysize);
307:
308: if (! refmpn_zero_p (xp+ysize, xsize-ysize))
309: cmp = 1;
310: else
311: cmp = refmpn_cmp (xp, yp, ysize);
312:
313: return (opp ? -cmp : cmp);
314: }
315:
316: int
317: refmpn_equal_anynail (mp_srcptr xp, mp_srcptr yp, mp_size_t size)
318: {
319: mp_size_t i;
320: ASSERT (size >= 0);
321:
322: for (i = 0; i < size; i++)
323: if (xp[i] != yp[i])
324: return 0;
325: return 1;
326: }
327:
328:
329: #define LOGOPS(operation) \
330: { \
331: mp_size_t i; \
332: \
333: ASSERT (refmpn_overlap_fullonly_two_p (rp, s1p, s2p, size)); \
334: ASSERT (size >= 1); \
335: ASSERT_MPN (s1p, size); \
336: ASSERT_MPN (s2p, size); \
337: \
338: for (i = 0; i < size; i++) \
339: rp[i] = operation; \
340: }
341:
342: void
343: refmpn_and_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
344: {
345: LOGOPS (s1p[i] & s2p[i]);
346: }
347: void
348: refmpn_andn_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
349: {
350: LOGOPS (s1p[i] & ~s2p[i]);
351: }
352: void
353: refmpn_nand_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
354: {
355: LOGOPS ((s1p[i] & s2p[i]) ^ GMP_NUMB_MASK);
356: }
357: void
358: refmpn_ior_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
359: {
360: LOGOPS (s1p[i] | s2p[i]);
361: }
362: void
363: refmpn_iorn_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
364: {
365: LOGOPS (s1p[i] | (s2p[i] ^ GMP_NUMB_MASK));
366: }
367: void
368: refmpn_nior_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
369: {
370: LOGOPS ((s1p[i] | s2p[i]) ^ GMP_NUMB_MASK);
371: }
372: void
373: refmpn_xor_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
374: {
375: LOGOPS (s1p[i] ^ s2p[i]);
376: }
377: void
378: refmpn_xnor_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
379: {
380: LOGOPS ((s1p[i] ^ s2p[i]) ^ GMP_NUMB_MASK);
381: }
382:
383: /* set *w to x+y, return 0 or 1 carry */
384: mp_limb_t
385: add (mp_limb_t *w, mp_limb_t x, mp_limb_t y)
386: {
387: mp_limb_t sum, cy;
388:
389: ASSERT_LIMB (x);
390: ASSERT_LIMB (y);
391:
392: sum = x + y;
393: #if GMP_NAIL_BITS == 0
394: *w = sum;
395: cy = (sum < x);
396: #else
397: *w = sum & GMP_NUMB_MASK;
398: cy = (sum >> GMP_NUMB_BITS);
399: #endif
400: return cy;
401: }
402:
403: /* set *w to x-y, return 0 or 1 borrow */
404: mp_limb_t
405: sub (mp_limb_t *w, mp_limb_t x, mp_limb_t y)
406: {
407: mp_limb_t diff, cy;
408:
409: ASSERT_LIMB (x);
410: ASSERT_LIMB (y);
411:
412: diff = x - y;
413: #if GMP_NAIL_BITS == 0
414: *w = diff;
415: cy = (diff > x);
416: #else
417: *w = diff & GMP_NUMB_MASK;
418: cy = (diff >> GMP_NUMB_BITS) & 1;
419: #endif
420: return cy;
421: }
422:
423: /* set *w to x+y+c (where c is 0 or 1), return 0 or 1 carry */
424: mp_limb_t
425: adc (mp_limb_t *w, mp_limb_t x, mp_limb_t y, mp_limb_t c)
426: {
427: mp_limb_t r;
428:
429: ASSERT_LIMB (x);
430: ASSERT_LIMB (y);
431: ASSERT (c == 0 || c == 1);
432:
433: r = add (w, x, y);
434: return r + add (w, *w, c);
435: }
436:
437: /* set *w to x-y-c (where c is 0 or 1), return 0 or 1 borrow */
438: mp_limb_t
439: sbb (mp_limb_t *w, mp_limb_t x, mp_limb_t y, mp_limb_t c)
440: {
441: mp_limb_t r;
442:
443: ASSERT_LIMB (x);
444: ASSERT_LIMB (y);
445: ASSERT (c == 0 || c == 1);
446:
447: r = sub (w, x, y);
448: return r + sub (w, *w, c);
449: }
450:
451:
452: #define AORS_1(operation) \
453: { \
454: mp_limb_t i; \
455: \
456: ASSERT (refmpn_overlap_fullonly_p (rp, sp, size)); \
457: ASSERT (size >= 1); \
458: ASSERT_MPN (sp, size); \
459: ASSERT_LIMB (n); \
460: \
461: for (i = 0; i < size; i++) \
462: n = operation (&rp[i], sp[i], n); \
463: return n; \
464: }
465:
466: mp_limb_t
467: refmpn_add_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t n)
468: {
469: AORS_1 (add);
470: }
471: mp_limb_t
472: refmpn_sub_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t n)
473: {
474: AORS_1 (sub);
475: }
476:
477: #define AORS_NC(operation) \
478: { \
479: mp_size_t i; \
480: \
481: ASSERT (refmpn_overlap_fullonly_two_p (rp, s1p, s2p, size)); \
482: ASSERT (carry == 0 || carry == 1); \
483: ASSERT (size >= 1); \
484: ASSERT_MPN (s1p, size); \
485: ASSERT_MPN (s2p, size); \
486: \
487: for (i = 0; i < size; i++) \
488: carry = operation (&rp[i], s1p[i], s2p[i], carry); \
489: return carry; \
490: }
491:
492: mp_limb_t
493: refmpn_add_nc (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size,
494: mp_limb_t carry)
495: {
496: AORS_NC (adc);
497: }
498: mp_limb_t
499: refmpn_sub_nc (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size,
500: mp_limb_t carry)
501: {
502: AORS_NC (sbb);
503: }
504:
505:
506: mp_limb_t
507: refmpn_add_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
508: {
509: return refmpn_add_nc (rp, s1p, s2p, size, CNST_LIMB(0));
510: }
511: mp_limb_t
512: refmpn_sub_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
513: {
514: return refmpn_sub_nc (rp, s1p, s2p, size, CNST_LIMB(0));
515: }
516:
517:
518: /* Twos complement, return borrow. */
519: mp_limb_t
520: refmpn_neg_n (mp_ptr dst, mp_srcptr src, mp_size_t size)
521: {
522: mp_ptr zeros;
523: mp_limb_t ret;
524:
525: ASSERT (size >= 1);
526:
527: zeros = refmpn_malloc_limbs (size);
528: refmpn_fill (zeros, size, CNST_LIMB(0));
529: ret = refmpn_sub_n (dst, zeros, src, size);
530: free (zeros);
531: return ret;
532: }
533:
534:
535: #define AORS(aors_n, aors_1) \
536: { \
537: mp_limb_t c; \
538: ASSERT (s1size >= s2size); \
539: ASSERT (s2size >= 1); \
540: c = aors_n (rp, s1p, s2p, s2size); \
541: if (s1size-s2size != 0) \
542: c = aors_1 (rp+s2size, s1p+s2size, s1size-s2size, c); \
543: return c; \
544: }
545: mp_limb_t
546: refmpn_add (mp_ptr rp,
547: mp_srcptr s1p, mp_size_t s1size,
548: mp_srcptr s2p, mp_size_t s2size)
549: {
550: AORS (refmpn_add_n, refmpn_add_1);
551: }
552: mp_limb_t
553: refmpn_sub (mp_ptr rp,
554: mp_srcptr s1p, mp_size_t s1size,
555: mp_srcptr s2p, mp_size_t s2size)
556: {
557: AORS (refmpn_sub_n, refmpn_sub_1);
558: }
559:
560:
561: #define SHIFTHIGH(x) ((x) << BITS_PER_MP_LIMB/2)
562: #define SHIFTLOW(x) ((x) >> BITS_PER_MP_LIMB/2)
563:
564: #define LOWMASK (((mp_limb_t) 1 << BITS_PER_MP_LIMB/2)-1)
565: #define HIGHMASK SHIFTHIGH(LOWMASK)
566:
567: #define LOWPART(x) ((x) & LOWMASK)
568: #define HIGHPART(x) SHIFTLOW((x) & HIGHMASK)
569:
570: /* Set *hi,*lo to x*y, using full limbs not nails. */
571: mp_limb_t
572: refmpn_umul_ppmm (mp_limb_t *lo, mp_limb_t x, mp_limb_t y)
573: {
574: mp_limb_t hi, s;
575:
576: *lo = LOWPART(x) * LOWPART(y);
577: hi = HIGHPART(x) * HIGHPART(y);
578:
579: s = LOWPART(x) * HIGHPART(y);
580: hi += HIGHPART(s);
581: s = SHIFTHIGH(LOWPART(s));
582: *lo += s;
583: hi += (*lo < s);
584:
585: s = HIGHPART(x) * LOWPART(y);
586: hi += HIGHPART(s);
587: s = SHIFTHIGH(LOWPART(s));
588: *lo += s;
589: hi += (*lo < s);
590:
591: return hi;
592: }
593:
594: mp_limb_t
595: refmpn_umul_ppmm_r (mp_limb_t x, mp_limb_t y, mp_limb_t *lo)
596: {
597: return refmpn_umul_ppmm (lo, x, y);
598: }
599:
600: mp_limb_t
601: refmpn_mul_1c (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t multiplier,
602: mp_limb_t carry)
603: {
604: mp_size_t i;
605: mp_limb_t hi, lo;
606:
607: ASSERT (refmpn_overlap_low_to_high_p (rp, sp, size));
608: ASSERT (size >= 1);
609: ASSERT_MPN (sp, size);
610: ASSERT_LIMB (multiplier);
611: ASSERT_LIMB (carry);
612:
613: multiplier <<= GMP_NAIL_BITS;
614: for (i = 0; i < size; i++)
615: {
616: hi = refmpn_umul_ppmm (&lo, sp[i], multiplier);
617: lo >>= GMP_NAIL_BITS;
618: ASSERT_NOCARRY (add (&hi, hi, add (&lo, lo, carry)));
619: rp[i] = lo;
620: carry = hi;
621: }
622: return carry;
623: }
624:
625: mp_limb_t
626: refmpn_mul_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t multiplier)
627: {
628: return refmpn_mul_1c (rp, sp, size, multiplier, CNST_LIMB(0));
629: }
630:
631:
632: mp_limb_t
633: refmpn_mul_2 (mp_ptr dst, mp_srcptr src, mp_size_t size, mp_srcptr mult)
634: {
635: mp_ptr src_copy;
636: mp_limb_t c;
637:
638: ASSERT (refmpn_overlap_fullonly_p (dst, src, size));
639: ASSERT (! refmpn_overlap_p (dst, size+1, mult, 2));
640: ASSERT (size >= 1);
641: ASSERT_MPN (mult, 2);
642:
643: /* in case dst==src */
644: src_copy = refmpn_malloc_limbs (size);
645: refmpn_copyi (src_copy, src, size);
646: src = src_copy;
647:
648: dst[size] = refmpn_mul_1 (dst, src, size, mult[0]);
649: c = refmpn_addmul_1 (dst+1, src, size, mult[1]);
650: free (src_copy);
651: return c;
652: }
653:
654:
655: #define AORSMUL_1C(operation_n) \
656: { \
657: mp_ptr p; \
658: mp_limb_t ret; \
659: \
660: ASSERT (refmpn_overlap_fullonly_p (rp, sp, size)); \
661: \
662: p = refmpn_malloc_limbs (size); \
663: ret = refmpn_mul_1c (p, sp, size, multiplier, carry); \
664: ret += operation_n (rp, rp, p, size); \
665: \
666: free (p); \
667: return ret; \
668: }
669:
670: mp_limb_t
671: refmpn_addmul_1c (mp_ptr rp, mp_srcptr sp, mp_size_t size,
672: mp_limb_t multiplier, mp_limb_t carry)
673: {
674: AORSMUL_1C (refmpn_add_n);
675: }
676: mp_limb_t
677: refmpn_submul_1c (mp_ptr rp, mp_srcptr sp, mp_size_t size,
678: mp_limb_t multiplier, mp_limb_t carry)
679: {
680: AORSMUL_1C (refmpn_sub_n);
681: }
682:
683:
684: mp_limb_t
685: refmpn_addmul_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t multiplier)
686: {
687: return refmpn_addmul_1c (rp, sp, size, multiplier, CNST_LIMB(0));
688: }
689: mp_limb_t
690: refmpn_submul_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t multiplier)
691: {
692: return refmpn_submul_1c (rp, sp, size, multiplier, CNST_LIMB(0));
693: }
694:
695:
696: mp_limb_t
697: refmpn_addsub_nc (mp_ptr r1p, mp_ptr r2p,
698: mp_srcptr s1p, mp_srcptr s2p, mp_size_t size,
699: mp_limb_t carry)
700: {
701: mp_ptr p;
702: mp_limb_t acy, scy;
703:
704: /* Destinations can't overlap. */
705: ASSERT (! refmpn_overlap_p (r1p, size, r2p, size));
706: ASSERT (refmpn_overlap_fullonly_two_p (r1p, s1p, s2p, size));
707: ASSERT (refmpn_overlap_fullonly_two_p (r2p, s1p, s2p, size));
708: ASSERT (size >= 1);
709:
710: /* in case r1p==s1p or r1p==s2p */
711: p = refmpn_malloc_limbs (size);
712:
713: acy = refmpn_add_nc (p, s1p, s2p, size, carry >> 1);
714: scy = refmpn_sub_nc (r2p, s1p, s2p, size, carry & 1);
715: refmpn_copyi (r1p, p, size);
716:
717: free (p);
718: return 2 * acy + scy;
719: }
720:
721: mp_limb_t
722: refmpn_addsub_n (mp_ptr r1p, mp_ptr r2p,
723: mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
724: {
725: return refmpn_addsub_nc (r1p, r2p, s1p, s2p, size, CNST_LIMB(0));
726: }
727:
728:
729: /* Right shift hi,lo and return the low limb of the result.
730: Note a shift by BITS_PER_MP_LIMB isn't assumed to work (doesn't on x86). */
731: mp_limb_t
732: rshift_make (mp_limb_t hi, mp_limb_t lo, unsigned shift)
733: {
734: ASSERT (shift < GMP_NUMB_BITS);
735: if (shift == 0)
736: return lo;
737: else
738: return ((hi << (GMP_NUMB_BITS-shift)) | (lo >> shift)) & GMP_NUMB_MASK;
739: }
740:
741: /* Left shift hi,lo and return the high limb of the result.
742: Note a shift by BITS_PER_MP_LIMB isn't assumed to work (doesn't on x86). */
743: mp_limb_t
744: lshift_make (mp_limb_t hi, mp_limb_t lo, unsigned shift)
745: {
746: ASSERT (shift < GMP_NUMB_BITS);
747: if (shift == 0)
748: return hi;
749: else
750: return ((hi << shift) | (lo >> (GMP_NUMB_BITS-shift))) & GMP_NUMB_MASK;
751: }
752:
753:
754: mp_limb_t
755: refmpn_rshift (mp_ptr rp, mp_srcptr sp, mp_size_t size, unsigned shift)
756: {
757: mp_limb_t ret;
758: mp_size_t i;
759:
760: ASSERT (refmpn_overlap_low_to_high_p (rp, sp, size));
761: ASSERT (size >= 1);
762: ASSERT (shift >= 1 && shift < GMP_NUMB_BITS);
763: ASSERT_MPN (sp, size);
764:
765: ret = rshift_make (sp[0], CNST_LIMB(0), shift);
766:
767: for (i = 0; i < size-1; i++)
768: rp[i] = rshift_make (sp[i+1], sp[i], shift);
769:
770: rp[i] = rshift_make (CNST_LIMB(0), sp[i], shift);
771: return ret;
772: }
773:
774: mp_limb_t
775: refmpn_lshift (mp_ptr rp, mp_srcptr sp, mp_size_t size, unsigned shift)
776: {
777: mp_limb_t ret;
778: mp_size_t i;
779:
780: ASSERT (refmpn_overlap_high_to_low_p (rp, sp, size));
781: ASSERT (size >= 1);
782: ASSERT (shift >= 1 && shift < GMP_NUMB_BITS);
783: ASSERT_MPN (sp, size);
784:
785: ret = lshift_make (CNST_LIMB(0), sp[size-1], shift);
786:
787: for (i = size-2; i >= 0; i--)
788: rp[i+1] = lshift_make (sp[i+1], sp[i], shift);
789:
790: rp[i+1] = lshift_make (sp[i+1], CNST_LIMB(0), shift);
791: return ret;
792: }
793:
794:
795: mp_limb_t
796: refmpn_rshift_or_copy (mp_ptr rp, mp_srcptr sp, mp_size_t size, unsigned shift)
797: {
798: if (shift == 0)
799: {
800: refmpn_copyi (rp, sp, size);
801: return 0;
802: }
803: else
804: {
805: return refmpn_rshift (rp, sp, size, shift);
806: }
807: }
808:
809: mp_limb_t
810: refmpn_lshift_or_copy (mp_ptr rp, mp_srcptr sp, mp_size_t size, unsigned shift)
811: {
812: if (shift == 0)
813: {
814: refmpn_copyd (rp, sp, size);
815: return 0;
816: }
817: else
818: {
819: return refmpn_lshift (rp, sp, size, shift);
820: }
821: }
822:
823:
824: /* Divide h,l by d, return quotient, store remainder to *rp.
825: Operates on full limbs, not nails.
826: Must have h < d.
827: __udiv_qrnnd_c isn't simple, and it's a bit slow, but it works. */
828: mp_limb_t
829: refmpn_udiv_qrnnd (mp_limb_t *rp, mp_limb_t h, mp_limb_t l, mp_limb_t d)
830: {
831: mp_limb_t q, r;
832: int n;
833:
834: ASSERT (d != 0);
835: ASSERT (h < d);
836:
837: #if 0
838: udiv_qrnnd (q, r, h, l, d);
839: *rp = r;
840: return q;
841: #endif
842:
843: n = refmpn_count_leading_zeros (d);
844: d <<= n;
845:
846: if (n != 0)
847: {
848: h = (h << n) | (l >> (GMP_LIMB_BITS - n));
849: l <<= n;
850: }
851:
852: __udiv_qrnnd_c (q, r, h, l, d);
853: r >>= n;
854: *rp = r;
855: return q;
856: }
857:
858: mp_limb_t
859: refmpn_udiv_qrnnd_r (mp_limb_t h, mp_limb_t l, mp_limb_t d, mp_limb_t *rp)
860: {
861: return refmpn_udiv_qrnnd (rp, h, l, d);
862: }
863:
864: /* This little subroutine avoids some bad code generation from i386 gcc 3.0
865: -fPIC -O2 -fomit-frame-pointer (%ebp being used uninitialized). */
866: static mp_limb_t
867: refmpn_divmod_1c_workaround (mp_ptr rp, mp_srcptr sp, mp_size_t size,
868: mp_limb_t divisor, mp_limb_t carry)
869: {
870: mp_size_t i;
871: for (i = size-1; i >= 0; i--)
872: {
873: rp[i] = refmpn_udiv_qrnnd (&carry, carry,
874: sp[i] << GMP_NAIL_BITS,
875: divisor << GMP_NAIL_BITS);
876: carry >>= GMP_NAIL_BITS;
877: }
878: return carry;
879: }
880:
881: mp_limb_t
882: refmpn_divmod_1c (mp_ptr rp, mp_srcptr sp, mp_size_t size,
883: mp_limb_t divisor, mp_limb_t carry)
884: {
885: mp_ptr sp_orig;
886: mp_ptr prod;
887: mp_limb_t carry_orig;
888:
889: ASSERT (refmpn_overlap_fullonly_p (rp, sp, size));
890: ASSERT (size >= 0);
891: ASSERT (carry < divisor);
892: ASSERT_MPN (sp, size);
893: ASSERT_LIMB (divisor);
894: ASSERT_LIMB (carry);
895:
896: if (size == 0)
897: return carry;
898:
899: sp_orig = refmpn_memdup_limbs (sp, size);
900: prod = refmpn_malloc_limbs (size);
901: carry_orig = carry;
902:
903: carry = refmpn_divmod_1c_workaround (rp, sp, size, divisor, carry);
904:
905: /* check by multiplying back */
906: #if 0
907: printf ("size=%ld divisor=0x%lX carry=0x%lX remainder=0x%lX\n",
908: size, divisor, carry_orig, carry);
909: mpn_trace("s",sp_copy,size);
910: mpn_trace("r",rp,size);
911: printf ("mul_1c %lX\n", refmpn_mul_1c (prod, rp, size, divisor, carry));
912: mpn_trace("p",prod,size);
913: #endif
914: ASSERT (refmpn_mul_1c (prod, rp, size, divisor, carry) == carry_orig);
915: ASSERT (refmpn_cmp (prod, sp_orig, size) == 0);
916: free (sp_orig);
917: free (prod);
918:
919: return carry;
920: }
921:
922: mp_limb_t
923: refmpn_divmod_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t divisor)
924: {
925: return refmpn_divmod_1c (rp, sp, size, divisor, CNST_LIMB(0));
926: }
927:
928:
929: mp_limb_t
930: refmpn_mod_1c (mp_srcptr sp, mp_size_t size, mp_limb_t divisor,
931: mp_limb_t carry)
932: {
933: mp_ptr p = refmpn_malloc_limbs (size);
934: carry = refmpn_divmod_1c (p, sp, size, divisor, carry);
935: free (p);
936: return carry;
937: }
938:
939: mp_limb_t
940: refmpn_mod_1 (mp_srcptr sp, mp_size_t size, mp_limb_t divisor)
941: {
942: return refmpn_mod_1c (sp, size, divisor, CNST_LIMB(0));
943: }
944:
945: mp_limb_t
946: refmpn_preinv_mod_1 (mp_srcptr sp, mp_size_t size, mp_limb_t divisor,
947: mp_limb_t inverse)
948: {
949: ASSERT (divisor & GMP_NUMB_HIGHBIT);
950: ASSERT (inverse == refmpn_invert_limb (divisor));
951: return refmpn_mod_1 (sp, size, divisor);
952: }
953:
954: /* This implementation will be rather slow, but has the advantage of being
955: in a different style than the libgmp versions. */
956: mp_limb_t
957: refmpn_mod_34lsub1 (mp_srcptr p, mp_size_t n)
958: {
959: ASSERT ((GMP_NUMB_BITS % 4) == 0);
960: return mpn_mod_1 (p, n, (CNST_LIMB(1) << (3 * GMP_NUMB_BITS / 4)) - 1);
961: }
962:
963:
964: mp_limb_t
965: refmpn_divrem_1c (mp_ptr rp, mp_size_t xsize,
966: mp_srcptr sp, mp_size_t size, mp_limb_t divisor,
967: mp_limb_t carry)
968: {
969: mp_ptr z;
970:
971: z = refmpn_malloc_limbs (xsize);
972: refmpn_fill (z, xsize, CNST_LIMB(0));
973:
974: carry = refmpn_divmod_1c (rp+xsize, sp, size, divisor, carry);
975: carry = refmpn_divmod_1c (rp, z, xsize, divisor, carry);
976:
977: free (z);
978: return carry;
979: }
980:
981: mp_limb_t
982: refmpn_divrem_1 (mp_ptr rp, mp_size_t xsize,
983: mp_srcptr sp, mp_size_t size, mp_limb_t divisor)
984: {
985: return refmpn_divrem_1c (rp, xsize, sp, size, divisor, CNST_LIMB(0));
986: }
987:
988: mp_limb_t
989: refmpn_preinv_divrem_1 (mp_ptr rp, mp_size_t xsize,
990: mp_srcptr sp, mp_size_t size,
991: mp_limb_t divisor, mp_limb_t inverse, unsigned shift)
992: {
993: ASSERT (size >= 0);
994: ASSERT (shift == refmpn_count_leading_zeros (divisor));
995: ASSERT (inverse == refmpn_invert_limb (divisor << shift));
996:
997: return refmpn_divrem_1 (rp, xsize, sp, size, divisor);
998: }
999:
1000:
1001: /* Inverse is floor((b*(b-d)-1) / d), per division by invariant integers
1002: paper, figure 8.1 m', where b=2^BITS_PER_MP_LIMB. Note that -d-1 < d
1003: since d has the high bit set. */
1004:
1005: mp_limb_t
1006: refmpn_invert_limb (mp_limb_t d)
1007: {
1008: mp_limb_t r;
1009: ASSERT (d & GMP_LIMB_HIGHBIT);
1010: return refmpn_udiv_qrnnd (&r, -d-1, -1, d);
1011: }
1012:
1013:
1014: /* The aim is to produce a dst quotient and return a remainder c, satisfying
1015: c*b^n + src-i == 3*dst, where i is the incoming carry.
1016:
1017: Some value c==0, c==1 or c==2 will satisfy, so just try each.
1018:
1019: If GMP_NUMB_BITS is even then 2^GMP_NUMB_BITS==1mod3 and a non-zero
1020: remainder from the first division attempt determines the correct
1021: remainder (3-c), but don't bother with that, since we can't guarantee
1022: anything about GMP_NUMB_BITS when using nails.
1023:
1024: If the initial src-i produces a borrow then refmpn_sub_1 leaves a twos
1025: complement negative, ie. b^n+a-i, and the calculation produces c1
1026: satisfying c1*b^n + b^n+src-i == 3*dst, from which clearly c=c1+1. This
1027: means it's enough to just add any borrow back at the end.
1028:
1029: A borrow only occurs when a==0 or a==1, and, by the same reasoning as in
1030: mpn/generic/diveby3.c, the c1 that results in those cases will only be 0
1031: or 1 respectively, so with 1 added the final return value is still in the
1032: prescribed range 0 to 2. */
1033:
1034: mp_limb_t
1035: refmpn_divexact_by3c (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t carry)
1036: {
1037: mp_ptr spcopy;
1038: mp_limb_t c, cs;
1039:
1040: ASSERT (refmpn_overlap_fullonly_p (rp, sp, size));
1041: ASSERT (size >= 1);
1042: ASSERT (carry <= 2);
1043: ASSERT_MPN (sp, size);
1044:
1045: spcopy = refmpn_malloc_limbs (size);
1046: cs = refmpn_sub_1 (spcopy, sp, size, carry);
1047:
1048: for (c = 0; c <= 2; c++)
1049: if (refmpn_divmod_1c (rp, spcopy, size, 3, c) == 0)
1050: goto done;
1051: ASSERT_FAIL (no value of c satisfies);
1052:
1053: done:
1054: c += cs;
1055: ASSERT (c <= 2);
1056:
1057: free (spcopy);
1058: return c;
1059: }
1060:
1061: mp_limb_t
1062: refmpn_divexact_by3 (mp_ptr rp, mp_srcptr sp, mp_size_t size)
1063: {
1064: return refmpn_divexact_by3c (rp, sp, size, CNST_LIMB(0));
1065: }
1066:
1067:
1068: /* The same as mpn/generic/mul_basecase.c, but using refmpn functions. */
1069: void
1070: refmpn_mul_basecase (mp_ptr prodp,
1071: mp_srcptr up, mp_size_t usize,
1072: mp_srcptr vp, mp_size_t vsize)
1073: {
1074: mp_size_t i;
1075:
1076: ASSERT (! refmpn_overlap_p (prodp, usize+vsize, up, usize));
1077: ASSERT (! refmpn_overlap_p (prodp, usize+vsize, vp, vsize));
1078: ASSERT (usize >= vsize);
1079: ASSERT (vsize >= 1);
1080: ASSERT_MPN (up, usize);
1081: ASSERT_MPN (vp, vsize);
1082:
1083: prodp[usize] = refmpn_mul_1 (prodp, up, usize, vp[0]);
1084: for (i = 1; i < vsize; i++)
1085: prodp[usize+i] = refmpn_addmul_1 (prodp+i, up, usize, vp[i]);
1086: }
1087:
1088: void
1089: refmpn_mul_n (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
1090: {
1091: refmpn_mul_basecase (prodp, up, size, vp, size);
1092: }
1093:
1094: void
1095: refmpn_sqr (mp_ptr dst, mp_srcptr src, mp_size_t size)
1096: {
1097: refmpn_mul_basecase (dst, src, size, src, size);
1098: }
1099:
1100: /* Allowing usize<vsize, usize==0 or vsize==0. */
1101: void
1102: refmpn_mul_any (mp_ptr prodp,
1103: mp_srcptr up, mp_size_t usize,
1104: mp_srcptr vp, mp_size_t vsize)
1105: {
1106: ASSERT (! refmpn_overlap_p (prodp, usize+vsize, up, usize));
1107: ASSERT (! refmpn_overlap_p (prodp, usize+vsize, vp, vsize));
1108: ASSERT (usize >= 0);
1109: ASSERT (vsize >= 0);
1110: ASSERT_MPN (up, usize);
1111: ASSERT_MPN (vp, vsize);
1112:
1113: if (usize == 0)
1114: {
1115: refmpn_fill (prodp, vsize, CNST_LIMB(0));
1116: return;
1117: }
1118:
1119: if (vsize == 0)
1120: {
1121: refmpn_fill (prodp, usize, CNST_LIMB(0));
1122: return;
1123: }
1124:
1125: if (usize >= vsize)
1126: refmpn_mul_basecase (prodp, up, usize, vp, vsize);
1127: else
1128: refmpn_mul_basecase (prodp, vp, vsize, up, usize);
1129: }
1130:
1131:
1132: mp_limb_t
1133: refmpn_gcd_1 (mp_srcptr xp, mp_size_t xsize, mp_limb_t y)
1134: {
1135: mp_limb_t x;
1136: int twos;
1137:
1138: ASSERT (y != 0);
1139: ASSERT (! refmpn_zero_p (xp, xsize));
1140: ASSERT_MPN (xp, xsize);
1141: ASSERT_LIMB (y);
1142:
1143: x = refmpn_mod_1 (xp, xsize, y);
1144: if (x == 0)
1145: return y;
1146:
1147: twos = 0;
1148: while ((x & 1) == 0 && (y & 1) == 0)
1149: {
1150: x >>= 1;
1151: y >>= 1;
1152: twos++;
1153: }
1154:
1155: for (;;)
1156: {
1157: while ((x & 1) == 0) x >>= 1;
1158: while ((y & 1) == 0) y >>= 1;
1159:
1160: if (x < y)
1161: MP_LIMB_T_SWAP (x, y);
1162:
1163: x -= y;
1164: if (x == 0)
1165: break;
1166: }
1167:
1168: return y << twos;
1169: }
1170:
1171:
1172: /* Based on the full limb x, not nails. */
1173: unsigned
1174: refmpn_count_leading_zeros (mp_limb_t x)
1175: {
1176: unsigned n = 0;
1177:
1178: ASSERT (x != 0);
1179:
1180: while ((x & GMP_LIMB_HIGHBIT) == 0)
1181: {
1182: x <<= 1;
1183: n++;
1184: }
1185: return n;
1186: }
1187:
1188: /* Full limbs allowed, not limited to nails. */
1189: unsigned
1190: refmpn_count_trailing_zeros (mp_limb_t x)
1191: {
1192: unsigned n = 0;
1193:
1194: ASSERT (x != 0);
1195: ASSERT_LIMB (x);
1196:
1197: while ((x & 1) == 0)
1198: {
1199: x >>= 1;
1200: n++;
1201: }
1202: return n;
1203: }
1204:
1205: /* Strip factors of two (low zero bits) from {p,size} by right shifting.
1206: The return value is the number of twos stripped. */
1207: mp_size_t
1208: refmpn_strip_twos (mp_ptr p, mp_size_t size)
1209: {
1210: mp_size_t limbs;
1211: unsigned shift;
1212:
1213: ASSERT (size >= 1);
1214: ASSERT (! refmpn_zero_p (p, size));
1215: ASSERT_MPN (p, size);
1216:
1217: for (limbs = 0; p[0] == 0; limbs++)
1218: {
1219: refmpn_copyi (p, p+1, size-1);
1220: p[size-1] = 0;
1221: }
1222:
1223: shift = refmpn_count_trailing_zeros (p[0]);
1224: if (shift)
1225: refmpn_rshift (p, p, size, shift);
1226:
1227: return limbs*GMP_NUMB_BITS + shift;
1228: }
1229:
1230: mp_limb_t
1231: refmpn_gcd (mp_ptr gp, mp_ptr xp, mp_size_t xsize, mp_ptr yp, mp_size_t ysize)
1232: {
1233: int cmp;
1234:
1235: ASSERT (ysize >= 1);
1236: ASSERT (xsize >= ysize);
1237: ASSERT ((xp[0] & 1) != 0);
1238: ASSERT ((yp[0] & 1) != 0);
1239: /* ASSERT (xp[xsize-1] != 0); */ /* don't think x needs to be odd */
1240: ASSERT (yp[ysize-1] != 0);
1241: ASSERT (refmpn_overlap_fullonly_p (gp, xp, xsize));
1242: ASSERT (refmpn_overlap_fullonly_p (gp, yp, ysize));
1243: ASSERT (! refmpn_overlap_p (xp, xsize, yp, ysize));
1244: if (xsize == ysize)
1245: ASSERT (refmpn_msbone (xp[xsize-1]) >= refmpn_msbone (yp[ysize-1]));
1246: ASSERT_MPN (xp, xsize);
1247: ASSERT_MPN (yp, ysize);
1248:
1249: refmpn_strip_twos (xp, xsize);
1250: MPN_NORMALIZE (xp, xsize);
1251: MPN_NORMALIZE (yp, ysize);
1252:
1253: for (;;)
1254: {
1255: cmp = refmpn_cmp_twosizes (xp, xsize, yp, ysize);
1256: if (cmp == 0)
1257: break;
1258: if (cmp < 0)
1259: MPN_PTR_SWAP (xp,xsize, yp,ysize);
1260:
1261: ASSERT_NOCARRY (refmpn_sub (xp, xp, xsize, yp, ysize));
1262:
1263: refmpn_strip_twos (xp, xsize);
1264: MPN_NORMALIZE (xp, xsize);
1265: }
1266:
1267: refmpn_copyi (gp, xp, xsize);
1268: return xsize;
1269: }
1270:
1271:
1272: unsigned long
1273: refmpn_popcount (mp_srcptr sp, mp_size_t size)
1274: {
1275: unsigned long count = 0;
1276: mp_size_t i;
1277: int j;
1278: mp_limb_t l;
1279:
1280: ASSERT (size >= 0);
1281: ASSERT_MPN (sp, size);
1282:
1283: for (i = 0; i < size; i++)
1284: {
1285: l = sp[i];
1286: for (j = 0; j < GMP_NUMB_BITS; j++)
1287: {
1288: count += (l & 1);
1289: l >>= 1;
1290: }
1291: }
1292: return count;
1293: }
1294:
1295: unsigned long
1296: refmpn_hamdist (mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
1297: {
1298: mp_ptr d;
1299: unsigned long count;
1300:
1301: ASSERT (size >= 0);
1302: ASSERT_MPN (s1p, size);
1303: ASSERT_MPN (s2p, size);
1304:
1305: if (size == 0)
1306: return 0;
1307:
1308: d = refmpn_malloc_limbs (size);
1309: refmpn_xor_n (d, s1p, s2p, size);
1310: count = refmpn_popcount (d, size);
1311: free (d);
1312: return count;
1313: }
1314:
1315:
1316: /* set r to a%d */
1317: void
1318: refmpn_mod2 (mp_limb_t r[2], const mp_limb_t a[2], const mp_limb_t d[2])
1319: {
1320: mp_limb_t D[2];
1321: int n;
1322:
1323: ASSERT (! refmpn_overlap_p (r, 2, a, 2));
1324: ASSERT (! refmpn_overlap_p (r, 2, d, 2));
1325: ASSERT_MPN (a, 2);
1326: ASSERT_MPN (d, 2);
1327:
1328: D[1] = d[1], D[0] = d[0];
1329: r[1] = a[1], r[0] = a[0];
1330: n = 0;
1331:
1332: for (;;)
1333: {
1334: if (D[1] & GMP_NUMB_HIGHBIT)
1335: break;
1336: if (refmpn_cmp (r, D, 2) <= 0)
1337: break;
1338: refmpn_lshift (D, D, 2, 1);
1339: n++;
1340: ASSERT (n <= GMP_NUMB_BITS);
1341: }
1342:
1343: while (n >= 0)
1344: {
1345: if (refmpn_cmp (r, D, 2) >= 0)
1346: ASSERT_NOCARRY (refmpn_sub_n (r, r, D, 2));
1347: refmpn_rshift (D, D, 2, 1);
1348: n--;
1349: }
1350:
1351: ASSERT (refmpn_cmp (r, d, 2) < 0);
1352: }
1353:
1354:
1355: /* Find n where 0<n<2^GMP_NUMB_BITS and n==c mod m */
1356: mp_limb_t
1357: refmpn_gcd_finda (const mp_limb_t c[2])
1358: {
1359: mp_limb_t n1[2], n2[2];
1360:
1361: ASSERT (c[0] != 0);
1362: ASSERT (c[1] != 0);
1363: ASSERT_MPN (c, 2);
1364:
1365: n1[0] = c[0];
1366: n1[1] = c[1];
1367:
1368: n2[0] = -n1[0];
1369: n2[1] = ~n1[1];
1370:
1371: while (n2[1] != 0)
1372: {
1373: refmpn_mod2 (n1, n1, n2);
1374:
1375: MP_LIMB_T_SWAP (n1[0], n2[0]);
1376: MP_LIMB_T_SWAP (n1[1], n2[1]);
1377: }
1378:
1379: return n2[0];
1380: }
1381:
1382:
1383: /* Similar to mpn/generic/sb_divrem_mn.c, but somewhat simplified, in
1384: particular the trial quotient is allowed to be 2 too big. */
1385: mp_limb_t
1386: refmpn_sb_divrem_mn (mp_ptr qp,
1387: mp_ptr np, mp_size_t nsize,
1388: mp_srcptr dp, mp_size_t dsize)
1389: {
1390: mp_limb_t retval = 0;
1391: mp_size_t i;
1392: mp_limb_t d1 = dp[dsize-1];
1393: mp_ptr np_orig = refmpn_memdup_limbs (np, nsize);
1394:
1395: ASSERT (nsize >= dsize);
1396: /* ASSERT (dsize > 2); */
1397: ASSERT (dsize >= 2);
1398: ASSERT (dp[dsize-1] & GMP_LIMB_HIGHBIT);
1399: ASSERT (! refmpn_overlap_p (qp, nsize-dsize, np, nsize) || qp+dsize >= np);
1400: ASSERT_MPN (np, nsize);
1401: ASSERT_MPN (dp, dsize);
1402:
1403: i = nsize-dsize;
1404: if (refmpn_cmp (np+i, dp, dsize) >= 0)
1405: {
1406: ASSERT_NOCARRY (refmpn_sub_n (np+i, np+i, dp, dsize));
1407: retval = 1;
1408: }
1409:
1410: for (i--; i >= 0; i--)
1411: {
1412: mp_limb_t n0 = np[i+dsize];
1413: mp_limb_t n1 = np[i+dsize-1];
1414: mp_limb_t q, dummy_r;
1415:
1416: ASSERT (n0 <= d1);
1417: if (n0 == d1)
1418: q = MP_LIMB_T_MAX;
1419: else
1420: q = refmpn_udiv_qrnnd (&dummy_r, n0, n1, d1);
1421:
1422: n0 -= refmpn_submul_1 (np+i, dp, dsize, q);
1423: ASSERT (n0 == 0 || n0 == MP_LIMB_T_MAX);
1424: if (n0)
1425: {
1426: q--;
1427: if (! refmpn_add_n (np+i, np+i, dp, dsize))
1428: {
1429: q--;
1430: ASSERT (refmpn_add_n (np+i, np+i, dp, dsize) != 0);
1431: }
1432: }
1433: np[i+dsize] = 0;
1434:
1435: qp[i] = q;
1436: }
1437:
1438: /* remainder < divisor */
1439: ASSERT (refmpn_cmp (np, dp, dsize) < 0);
1440:
1441: /* multiply back to original */
1442: {
1443: mp_ptr mp = refmpn_malloc_limbs (nsize);
1444:
1445: refmpn_mul_any (mp, qp, nsize-dsize, dp, dsize);
1446: if (retval)
1447: ASSERT_NOCARRY (refmpn_add_n (mp+nsize-dsize,mp+nsize-dsize, dp, dsize));
1448: ASSERT_NOCARRY (refmpn_add (mp, mp, nsize, np, dsize));
1449: ASSERT (refmpn_cmp (mp, np_orig, nsize) == 0);
1450:
1451: free (mp);
1452: }
1453:
1454: free (np_orig);
1455: return retval;
1456: }
1457:
1458: /* Similar to mpn/generic/sb_divrem_mn.c, but somewhat simplified, in
1459: particular the trial quotient is allowed to be 2 too big. */
1460: void
1461: refmpn_tdiv_qr (mp_ptr qp, mp_ptr rp, mp_size_t qxn,
1462: mp_ptr np, mp_size_t nsize,
1463: mp_srcptr dp, mp_size_t dsize)
1464: {
1465: ASSERT (qxn == 0);
1466: ASSERT_MPN (np, nsize);
1467: ASSERT_MPN (dp, dsize);
1468:
1469: if (dsize == 1)
1470: {
1471: rp[0] = refmpn_divmod_1 (qp, np, nsize, dp[0]);
1472: return;
1473: }
1474: else
1475: {
1476: mp_ptr n2p = refmpn_malloc_limbs (nsize+1);
1477: mp_ptr d2p = refmpn_malloc_limbs (dsize);
1478: int norm = refmpn_count_leading_zeros (dp[dsize-1]);
1479:
1480: n2p[nsize] = refmpn_lshift_or_copy (n2p, np, nsize, norm);
1481: ASSERT_NOCARRY (refmpn_lshift_or_copy (d2p, dp, dsize, norm));
1482:
1483: refmpn_sb_divrem_mn (qp, n2p, nsize+1, d2p, dsize);
1484: refmpn_rshift_or_copy (rp, n2p, dsize, norm);
1485:
1486: /* ASSERT (refmpn_zero_p (tp+dsize, nsize-dsize)); */
1487: free (n2p);
1488: free (d2p);
1489: }
1490: }
1491:
1492: size_t
1493: refmpn_get_str (unsigned char *dst, int base, mp_ptr src, mp_size_t size)
1494: {
1495: unsigned char *d;
1496: size_t dsize;
1497:
1498: ASSERT (size >= 0);
1499: ASSERT (base >= 2);
1500: ASSERT (base < numberof (__mp_bases));
1501: ASSERT (size == 0 || src[size-1] != 0);
1502: ASSERT_MPN (src, size);
1503:
1504: MPN_SIZEINBASE (dsize, src, size, base);
1505: ASSERT (dsize >= 1);
1506: ASSERT (! byte_overlap_p (dst, dsize, src, size * BYTES_PER_MP_LIMB));
1507:
1508: if (size == 0)
1509: {
1510: dst[0] = 0;
1511: return 1;
1512: }
1513:
1514: /* don't clobber input for power of 2 bases */
1515: if (POW2_P (base))
1516: src = refmpn_memdup_limbs (src, size);
1517:
1518: d = dst + dsize;
1519: do
1520: {
1521: d--;
1522: ASSERT (d >= dst);
1523: *d = refmpn_divrem_1 (src, 0, src, size, (mp_limb_t) base);
1524: size -= (src[size-1] == 0);
1525: }
1526: while (size != 0);
1527:
1528: /* at most one leading zero */
1529: ASSERT (d == dst || d == dst+1);
1530: if (d != dst)
1531: *dst = 0;
1532:
1533: if (POW2_P (base))
1534: free (src);
1535:
1536: return dsize;
1537: }
1538:
1539:
1540: mp_limb_t
1541: refmpn_bswap_limb (mp_limb_t src)
1542: {
1543: mp_limb_t dst;
1544: int i;
1545:
1546: dst = 0;
1547: for (i = 0; i < BYTES_PER_MP_LIMB; i++)
1548: {
1549: dst = (dst << 8) + (src & 0xFF);
1550: src >>= 8;
1551: }
1552: return dst;
1553: }
1554:
1555:
1556: /* These random functions are mostly for transitional purposes while adding
1557: nail support, since they're independent of the normal mpn routines. They
1558: can probably be removed when those normal routines are reliable, though
1559: perhaps something independent would still be useful at times. */
1560:
1561: #if BITS_PER_MP_LIMB == 32
1562: #define RAND_A CNST_LIMB(0x29CF535)
1563: #endif
1564: #if BITS_PER_MP_LIMB == 64
1565: #define RAND_A CNST_LIMB(0xBAECD515DAF0B49D)
1566: #endif
1567:
1568: mp_limb_t refmpn_random_seed;
1569:
1570: mp_limb_t
1571: refmpn_random_half (void)
1572: {
1573: refmpn_random_seed = refmpn_random_seed * RAND_A + 1;
1574: return (refmpn_random_seed >> BITS_PER_MP_LIMB/2);
1575: }
1576:
1577: mp_limb_t
1578: refmpn_random_limb (void)
1579: {
1580: return ((refmpn_random_half () << (BITS_PER_MP_LIMB/2))
1581: | refmpn_random_half ()) & GMP_NUMB_MASK;
1582: }
1583:
1584: void
1585: refmpn_random (mp_ptr ptr, mp_size_t size)
1586: {
1587: mp_size_t i;
1588: if (GMP_NAIL_BITS == 0)
1589: {
1590: mpn_random (ptr, size);
1591: return;
1592: }
1593:
1594: for (i = 0; i < size; i++)
1595: ptr[i] = refmpn_random_limb ();
1596: }
1597:
1598: void
1599: refmpn_random2 (mp_ptr ptr, mp_size_t size)
1600: {
1601: mp_size_t i;
1602: mp_limb_t bit, mask, limb;
1603: int run;
1604:
1605: if (GMP_NAIL_BITS == 0)
1606: {
1607: mpn_random2 (ptr, size);
1608: return;
1609: }
1610:
1611: #define RUN_MODULUS 32
1612:
1613: /* start with ones at a random pos in the high limb */
1614: bit = CNST_LIMB(1) << (refmpn_random_half () % GMP_NUMB_BITS);
1615: mask = 0;
1616: run = 0;
1617:
1618: for (i = size-1; i >= 0; i--)
1619: {
1620: limb = 0;
1621: do
1622: {
1623: if (run == 0)
1624: {
1625: run = (refmpn_random_half () % RUN_MODULUS) + 1;
1626: mask = ~mask;
1627: }
1628:
1629: limb |= (bit & mask);
1630: bit >>= 1;
1631: run--;
1632: }
1633: while (bit != 0);
1634:
1635: ptr[i] = limb;
1636: bit = GMP_NUMB_HIGHBIT;
1637: }
1638: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to refmpn.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / tests