Annotation of OpenXM_contrib/gmp/mpn/tests/ref.c, Revision 1.1
1.1 ! maekawa 1: /* Reference mpn functions, designed to be simple, portable and independent
! 2: of the normal gmp code. Speed isn't a consideration. */
! 3:
! 4: /*
! 5: Copyright (C) 1996, 1997, 1998, 1999, 2000 Free Software Foundation, 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:
! 26: /* Most routines have assertions representing what the mpn routines are
! 27: supposed to accept. Many of these reference routines do sensible things
! 28: outside these ranges (eg. for size==0), but the assertions are present to
! 29: pick up bad parameters passed here that are about to be passed the same
! 30: to a real mpn routine being compared. */
! 31:
! 32: /* always do assertion checking */
! 33: #define WANT_ASSERT 1
! 34:
! 35: #include <stdio.h> /* for NULL */
! 36: #include <stdlib.h> /* for malloc */
! 37: #include "gmp.h"
! 38: #include "gmp-impl.h"
! 39: #include "longlong.h"
! 40:
! 41: #include "ref.h"
! 42:
! 43:
! 44: /* Check overlap for a routine defined to work low to high. */
! 45: int
! 46: refmpn_overlap_low_to_high_p (mp_srcptr dst, mp_srcptr src, mp_size_t size)
! 47: {
! 48: return (!MPN_OVERLAP_P (dst, size, src, size) || dst <= src);
! 49: }
! 50:
! 51: /* Check overlap for a routine defined to work high to low. */
! 52: int
! 53: refmpn_overlap_high_to_low_p (mp_srcptr dst, mp_srcptr src, mp_size_t size)
! 54: {
! 55: return (!MPN_OVERLAP_P (dst, size, src, size) || dst >= src);
! 56: }
! 57:
! 58: /* Check overlap for a standard routine requiring equal or separate. */
! 59: int
! 60: refmpn_overlap_fullonly_p (mp_srcptr dst, mp_srcptr src, mp_size_t size)
! 61: {
! 62: return (dst == src || !MPN_OVERLAP_P (dst, size, src, size));
! 63: }
! 64: int
! 65: refmpn_overlap_fullonly_two_p (mp_srcptr dst, mp_srcptr src1, mp_srcptr src2,
! 66: mp_size_t size)
! 67: {
! 68: return (refmpn_overlap_fullonly_p (dst, src1, size)
! 69: && refmpn_overlap_fullonly_p (dst, src2, size));
! 70: }
! 71:
! 72:
! 73: mp_ptr
! 74: refmpn_malloc_limbs (mp_size_t size)
! 75: {
! 76: mp_ptr p;
! 77: ASSERT (size >= 0);
! 78: if (size == 0)
! 79: size = 1;
! 80: p = (mp_ptr) malloc (size * BYTES_PER_MP_LIMB);
! 81: ASSERT (p != NULL);
! 82: return p;
! 83: }
! 84:
! 85: mp_ptr
! 86: refmpn_memdup_limbs (mp_srcptr ptr, mp_size_t size)
! 87: {
! 88: mp_ptr p;
! 89: p = refmpn_malloc_limbs (size);
! 90: refmpn_copyi (p, ptr, size);
! 91: return p;
! 92: }
! 93:
! 94: void
! 95: refmpn_fill (mp_ptr ptr, mp_size_t size, mp_limb_t value)
! 96: {
! 97: mp_size_t i;
! 98: ASSERT (size >= 0);
! 99: for (i = 0; i < size; i++)
! 100: ptr[i] = value;
! 101: }
! 102:
! 103: int
! 104: refmpn_zero_p (mp_srcptr ptr, mp_size_t size)
! 105: {
! 106: mp_size_t i;
! 107: for (i = 0; i < size; i++)
! 108: if (ptr[i] != 0)
! 109: return 0;
! 110: return 1;
! 111: }
! 112:
! 113: mp_limb_t
! 114: refmpn_msbone (mp_limb_t x)
! 115: {
! 116: mp_limb_t n = (mp_limb_t) 1 << (BITS_PER_MP_LIMB-1);
! 117:
! 118: while (n != 0)
! 119: {
! 120: if (x & n)
! 121: break;
! 122: n >>= 1;
! 123: }
! 124: return n;
! 125: }
! 126:
! 127: /* a mask of the MSB one bit and all bits below */
! 128: mp_limb_t
! 129: refmpn_msbone_mask (mp_limb_t x)
! 130: {
! 131: if (x == 0)
! 132: return 0;
! 133:
! 134: return (refmpn_msbone (x) << 1) - 1;
! 135: }
! 136:
! 137: void
! 138: refmpn_copyi (mp_ptr rp, mp_srcptr sp, mp_size_t size)
! 139: {
! 140: mp_size_t i;
! 141:
! 142: ASSERT (refmpn_overlap_low_to_high_p (rp, sp, size));
! 143: ASSERT (size >= 0);
! 144:
! 145: for (i = 0; i < size; i++)
! 146: rp[i] = sp[i];
! 147: }
! 148:
! 149: void
! 150: refmpn_copyd (mp_ptr rp, mp_srcptr sp, mp_size_t size)
! 151: {
! 152: mp_size_t i;
! 153:
! 154: ASSERT (refmpn_overlap_high_to_low_p (rp, sp, size));
! 155: ASSERT (size >= 0);
! 156:
! 157: for (i = size-1; i >= 0; i--)
! 158: rp[i] = sp[i];
! 159: }
! 160:
! 161: void
! 162: refmpn_com_n (mp_ptr rp, mp_srcptr sp, mp_size_t size)
! 163: {
! 164: mp_size_t i;
! 165:
! 166: ASSERT (refmpn_overlap_fullonly_p (rp, sp, size));
! 167: ASSERT (size >= 1);
! 168:
! 169: for (i = 0; i < size; i++)
! 170: rp[i] = ~sp[i];
! 171: }
! 172:
! 173:
! 174: int
! 175: refmpn_cmp (mp_srcptr xp, mp_srcptr yp, mp_size_t size)
! 176: {
! 177: mp_size_t i;
! 178:
! 179: ASSERT (size >= 1);
! 180:
! 181: for (i = size-1; i >= 0; i--)
! 182: {
! 183: if (xp[i] > yp[i]) return 1;
! 184: if (xp[i] < yp[i]) return -1;
! 185: }
! 186: return 0;
! 187: }
! 188:
! 189: int
! 190: refmpn_cmp_twosizes (mp_srcptr xp, mp_size_t xsize,
! 191: mp_srcptr yp, mp_size_t ysize)
! 192: {
! 193: int opp, cmp;
! 194:
! 195: opp = (xsize < ysize);
! 196: if (opp)
! 197: MPN_SRCPTR_SWAP (xp,xsize, yp,ysize);
! 198:
! 199: if (! refmpn_zero_p (xp+ysize, xsize-ysize))
! 200: cmp = 1;
! 201: else
! 202: cmp = refmpn_cmp (xp, yp, ysize);
! 203:
! 204: return (opp ? -cmp : cmp);
! 205: }
! 206:
! 207:
! 208: #define LOGOPS(operation) \
! 209: { \
! 210: mp_size_t i; \
! 211: \
! 212: ASSERT (refmpn_overlap_fullonly_two_p (rp, s1p, s2p, size)); \
! 213: ASSERT (size >= 1); \
! 214: \
! 215: for (i = 0; i < size; i++) \
! 216: rp[i] = operation; \
! 217: }
! 218:
! 219: void
! 220: refmpn_and_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
! 221: {
! 222: LOGOPS (s1p[i] & s2p[i]);
! 223: }
! 224: void
! 225: refmpn_andn_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
! 226: {
! 227: LOGOPS (s1p[i] & ~s2p[i]);
! 228: }
! 229: void
! 230: refmpn_nand_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
! 231: {
! 232: LOGOPS (~(s1p[i] & s2p[i]));
! 233: }
! 234: void
! 235: refmpn_ior_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
! 236: {
! 237: LOGOPS (s1p[i] | s2p[i]);
! 238: }
! 239: void
! 240: refmpn_iorn_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
! 241: {
! 242: LOGOPS (s1p[i] | ~s2p[i]);
! 243: }
! 244: void
! 245: refmpn_nior_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
! 246: {
! 247: LOGOPS (~(s1p[i] | s2p[i]));
! 248: }
! 249: void
! 250: refmpn_xor_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
! 251: {
! 252: LOGOPS (s1p[i] ^ s2p[i]);
! 253: }
! 254: void
! 255: refmpn_xnor_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
! 256: {
! 257: LOGOPS (~(s1p[i] ^ s2p[i]));
! 258: }
! 259:
! 260: /* set *w to x+y, return 0 or 1 carry */
! 261: mp_limb_t
! 262: add (mp_limb_t *w, mp_limb_t x, mp_limb_t y)
! 263: {
! 264: *w = x + y;
! 265: return *w < x;
! 266: }
! 267:
! 268: /* set *w to x-y, return 0 or 1 borrow */
! 269: mp_limb_t
! 270: sub (mp_limb_t *w, mp_limb_t x, mp_limb_t y)
! 271: {
! 272: *w = x - y;
! 273: return *w > x;
! 274: }
! 275:
! 276: /* set *w to x+y+c (where c is 0 or 1), return 0 or 1 carry */
! 277: mp_limb_t
! 278: adc (mp_limb_t *w, mp_limb_t x, mp_limb_t y, mp_limb_t c)
! 279: {
! 280: mp_limb_t r;
! 281: ASSERT (c == 0 || c == 1);
! 282: r = add (w, x, y);
! 283: return r + add (w, *w, c);
! 284: }
! 285:
! 286: /* set *w to x-y-c (where c is 0 or 1), return 0 or 1 borrow */
! 287: mp_limb_t
! 288: sbb (mp_limb_t *w, mp_limb_t x, mp_limb_t y, mp_limb_t c)
! 289: {
! 290: mp_limb_t r;
! 291: ASSERT (c == 0 || c == 1);
! 292: r = sub (w, x, y);
! 293: return r + sub (w, *w, c);
! 294: }
! 295:
! 296:
! 297: #define AORS_1(operation) \
! 298: { \
! 299: mp_limb_t i; \
! 300: \
! 301: ASSERT (refmpn_overlap_fullonly_p (rp, sp, size)); \
! 302: ASSERT (size >= 1); \
! 303: \
! 304: for (i = 0; i < size; i++) \
! 305: n = operation (&rp[i], sp[i], n); \
! 306: return n; \
! 307: }
! 308:
! 309: mp_limb_t
! 310: refmpn_add_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t n)
! 311: {
! 312: AORS_1 (add);
! 313: }
! 314: mp_limb_t
! 315: refmpn_sub_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t n)
! 316: {
! 317: AORS_1 (sub);
! 318: }
! 319:
! 320: #define AORS_NC(operation) \
! 321: { \
! 322: mp_size_t i; \
! 323: \
! 324: ASSERT (refmpn_overlap_fullonly_two_p (rp, s1p, s2p, size)); \
! 325: ASSERT (carry == 0 || carry == 1); \
! 326: ASSERT (size >= 1); \
! 327: \
! 328: for (i = 0; i < size; i++) \
! 329: carry = operation (&rp[i], s1p[i], s2p[i], carry); \
! 330: return carry; \
! 331: }
! 332:
! 333: mp_limb_t
! 334: refmpn_add_nc (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size,
! 335: mp_limb_t carry)
! 336: {
! 337: AORS_NC (adc);
! 338: }
! 339: mp_limb_t
! 340: refmpn_sub_nc (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size,
! 341: mp_limb_t carry)
! 342: {
! 343: AORS_NC (sbb);
! 344: }
! 345:
! 346:
! 347: mp_limb_t
! 348: refmpn_add_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
! 349: {
! 350: return refmpn_add_nc (rp, s1p, s2p, size, 0);
! 351: }
! 352: mp_limb_t
! 353: refmpn_sub_n (mp_ptr rp, mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
! 354: {
! 355: return refmpn_sub_nc (rp, s1p, s2p, size, 0);
! 356: }
! 357:
! 358:
! 359: #define AORS(aors_n, aors_1) \
! 360: { \
! 361: mp_limb_t c; \
! 362: ASSERT (s1size >= s2size); \
! 363: ASSERT (s2size >= 1); \
! 364: c = aors_n (rp, s1p, s2p, s2size); \
! 365: if (s1size-s2size != 0) \
! 366: c = aors_1 (rp+s2size, s1p+s2size, s1size-s2size, c); \
! 367: return c; \
! 368: }
! 369: mp_limb_t
! 370: refmpn_add (mp_ptr rp,
! 371: mp_srcptr s1p, mp_size_t s1size,
! 372: mp_srcptr s2p, mp_size_t s2size)
! 373: {
! 374: AORS (refmpn_add_n, refmpn_add_1);
! 375: }
! 376: mp_limb_t
! 377: refmpn_sub (mp_ptr rp,
! 378: mp_srcptr s1p, mp_size_t s1size,
! 379: mp_srcptr s2p, mp_size_t s2size)
! 380: {
! 381: AORS (refmpn_sub_n, refmpn_sub_1);
! 382: }
! 383:
! 384:
! 385: #define SHIFTHIGH(x) ((x) << BITS_PER_MP_LIMB/2)
! 386: #define SHIFTLOW(x) ((x) >> BITS_PER_MP_LIMB/2)
! 387:
! 388: #define LOWMASK (((mp_limb_t) 1 << BITS_PER_MP_LIMB/2)-1)
! 389: #define HIGHMASK SHIFTHIGH(LOWMASK)
! 390:
! 391: #define LOWPART(x) ((x) & LOWMASK)
! 392: #define HIGHPART(x) SHIFTLOW((x) & HIGHMASK)
! 393:
! 394: /* set *hi,*lo to x*y */
! 395: void
! 396: mul (mp_limb_t *hi, mp_limb_t *lo, mp_limb_t x, mp_limb_t y)
! 397: {
! 398: mp_limb_t s;
! 399:
! 400: *lo = LOWPART(x) * LOWPART(y);
! 401: *hi = HIGHPART(x) * HIGHPART(y);
! 402:
! 403: s = LOWPART(x) * HIGHPART(y);
! 404: ASSERT_NOCARRY (add (hi, *hi, add (lo, *lo, SHIFTHIGH(LOWPART(s)))));
! 405: ASSERT_NOCARRY (add (hi, *hi, HIGHPART(s)));
! 406:
! 407: s = HIGHPART(x) * LOWPART(y);
! 408: ASSERT_NOCARRY (add (hi, *hi, add (lo, *lo, SHIFTHIGH(LOWPART(s)))));
! 409: ASSERT_NOCARRY (add (hi, *hi, HIGHPART(s)));
! 410: }
! 411:
! 412: mp_limb_t
! 413: refmpn_mul_1c (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t multiplier,
! 414: mp_limb_t carry)
! 415: {
! 416: mp_size_t i;
! 417: mp_limb_t hi, lo;
! 418:
! 419: ASSERT (refmpn_overlap_fullonly_p (rp, sp, size));
! 420: ASSERT (size >= 1);
! 421:
! 422: for (i = 0; i < size; i++)
! 423: {
! 424: mul (&hi, &lo, sp[i], multiplier);
! 425: ASSERT_NOCARRY (add (&hi, hi, add (&lo, lo, carry)));
! 426: rp[i] = lo;
! 427: carry = hi;
! 428: }
! 429:
! 430: return carry;
! 431: }
! 432:
! 433: mp_limb_t
! 434: refmpn_mul_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t multiplier)
! 435: {
! 436: return refmpn_mul_1c (rp, sp, size, multiplier, 0);
! 437: }
! 438:
! 439:
! 440: #define AORSMUL_1C(operation_n) \
! 441: { \
! 442: mp_ptr p = refmpn_malloc_limbs (size); \
! 443: mp_limb_t ret; \
! 444: \
! 445: ASSERT (refmpn_overlap_fullonly_p (rp, sp, size)); \
! 446: ASSERT (size >= 1); \
! 447: \
! 448: ret = refmpn_mul_1c (p, sp, size, multiplier, carry); \
! 449: ret += operation_n (rp, rp, p, size); \
! 450: \
! 451: free (p); \
! 452: return ret; \
! 453: }
! 454:
! 455: mp_limb_t
! 456: refmpn_addmul_1c (mp_ptr rp, mp_srcptr sp, mp_size_t size,
! 457: mp_limb_t multiplier, mp_limb_t carry)
! 458: {
! 459: AORSMUL_1C (refmpn_add_n);
! 460: }
! 461: mp_limb_t
! 462: refmpn_submul_1c (mp_ptr rp, mp_srcptr sp, mp_size_t size,
! 463: mp_limb_t multiplier, mp_limb_t carry)
! 464: {
! 465: AORSMUL_1C (refmpn_sub_n);
! 466: }
! 467:
! 468:
! 469: mp_limb_t
! 470: refmpn_addmul_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t multiplier)
! 471: {
! 472: return refmpn_addmul_1c (rp, sp, size, multiplier, 0);
! 473: }
! 474: mp_limb_t
! 475: refmpn_submul_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t multiplier)
! 476: {
! 477: return refmpn_submul_1c (rp, sp, size, multiplier, 0);
! 478: }
! 479:
! 480:
! 481: mp_limb_t
! 482: refmpn_addsub_nc (mp_ptr r1p, mp_ptr r2p,
! 483: mp_srcptr s1p, mp_srcptr s2p, mp_size_t size,
! 484: mp_limb_t carry)
! 485: {
! 486: mp_ptr p;
! 487: mp_limb_t acy, scy;
! 488:
! 489: /* Destinations can't overlap. */
! 490: ASSERT (! MPN_OVERLAP_P (r1p, size, r2p, size));
! 491: ASSERT (refmpn_overlap_fullonly_two_p (r1p, s1p, s2p, size));
! 492: ASSERT (refmpn_overlap_fullonly_two_p (r2p, s1p, s2p, size));
! 493: ASSERT (size >= 1);
! 494:
! 495: p = refmpn_malloc_limbs (size);
! 496: acy = refmpn_add_nc (p, s1p, s2p, size, carry >> 1);
! 497: scy = refmpn_sub_nc (r2p, s1p, s2p, size, carry & 1);
! 498: refmpn_copyi (r1p, p, size);
! 499: free (p);
! 500: return 2 * acy + scy;
! 501: }
! 502:
! 503: mp_limb_t
! 504: refmpn_addsub_n (mp_ptr r1p, mp_ptr r2p,
! 505: mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
! 506: {
! 507: return refmpn_addsub_nc (r1p, r2p, s1p, s2p, size, 0);
! 508: }
! 509:
! 510:
! 511: /* Right shift hi,lo and return the low limb of the result.
! 512: Note a shift by BITS_PER_MP_LIMB isn't assumed to work (doesn't on x86). */
! 513: mp_limb_t
! 514: rshift_make (mp_limb_t hi, mp_limb_t lo, unsigned shift)
! 515: {
! 516: ASSERT (shift >= 0 && shift < BITS_PER_MP_LIMB);
! 517: if (shift == 0)
! 518: return lo;
! 519: else
! 520: return (hi << (BITS_PER_MP_LIMB-shift)) | (lo >> shift);
! 521: }
! 522:
! 523: /* Left shift hi,lo and return the high limb of the result.
! 524: Note a shift by BITS_PER_MP_LIMB isn't assumed to work (doesn't on x86). */
! 525: mp_limb_t
! 526: lshift_make (mp_limb_t hi, mp_limb_t lo, unsigned shift)
! 527: {
! 528: ASSERT (shift >= 0 && shift < BITS_PER_MP_LIMB);
! 529: if (shift == 0)
! 530: return hi;
! 531: else
! 532: return (hi << shift) | (lo >> (BITS_PER_MP_LIMB-shift));
! 533: }
! 534:
! 535:
! 536: mp_limb_t
! 537: refmpn_rshift (mp_ptr rp, mp_srcptr sp, mp_size_t size, unsigned shift)
! 538: {
! 539: mp_limb_t ret;
! 540: mp_size_t i;
! 541:
! 542: ASSERT (refmpn_overlap_low_to_high_p (rp, sp, size));
! 543: ASSERT (size >= 1);
! 544: ASSERT (shift >= 1 && shift < BITS_PER_MP_LIMB);
! 545:
! 546: ret = rshift_make (sp[0], 0, shift);
! 547:
! 548: for (i = 0; i < size-1; i++)
! 549: rp[i] = rshift_make (sp[i+1], sp[i], shift);
! 550:
! 551: rp[i] = rshift_make (0, sp[i], shift);
! 552: return ret;
! 553: }
! 554:
! 555: mp_limb_t
! 556: refmpn_lshift (mp_ptr rp, mp_srcptr sp, mp_size_t size, unsigned shift)
! 557: {
! 558: mp_limb_t ret;
! 559: mp_size_t i;
! 560:
! 561: ASSERT (refmpn_overlap_high_to_low_p (rp, sp, size));
! 562: ASSERT (size >= 1);
! 563: ASSERT (shift >= 1 && shift < BITS_PER_MP_LIMB);
! 564:
! 565: ret = lshift_make (0, sp[size-1], shift);
! 566:
! 567: for (i = size-2; i >= 0; i--)
! 568: rp[i+1] = lshift_make (sp[i+1], sp[i], shift);
! 569:
! 570: rp[i+1] = lshift_make (sp[i+1], 0, shift);
! 571: return ret;
! 572: }
! 573:
! 574:
! 575: /* Divide h,l by d, producing a quotient *q and remainder *r.
! 576: Must have h < d.
! 577: __udiv_qrnnd_c isn't simple, and it's a bit slow, but it works. */
! 578: void
! 579: div1 (mp_limb_t *q, mp_limb_t *r, mp_limb_t h, mp_limb_t l, mp_limb_t d)
! 580: {
! 581: int n;
! 582:
! 583: ASSERT (d != 0);
! 584: ASSERT (h < d);
! 585:
! 586: #if 0
! 587: udiv_qrnnd (*q, *r, h, l, d);
! 588: return;
! 589: #endif
! 590:
! 591: for (n = 0; !(d & MP_LIMB_T_HIGHBIT); n++)
! 592: d <<= 1;
! 593:
! 594: h = lshift_make (h, l, n);
! 595: l <<= n;
! 596:
! 597: __udiv_qrnnd_c (*q, *r, h, l, d);
! 598: *r >>= n;
! 599: }
! 600:
! 601: mp_limb_t
! 602: refmpn_divmod_1c (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t divisor,
! 603: mp_limb_t carry)
! 604: {
! 605: mp_ptr sp_orig;
! 606: mp_ptr prod;
! 607: mp_limb_t carry_orig;
! 608: mp_size_t i;
! 609:
! 610: ASSERT (refmpn_overlap_fullonly_p (rp, sp, size));
! 611: ASSERT (size >= 0);
! 612: ASSERT (carry < divisor);
! 613:
! 614: if (size == 0)
! 615: return carry;
! 616:
! 617: sp_orig = refmpn_memdup_limbs (sp, size);
! 618: prod = refmpn_malloc_limbs (size);
! 619: carry_orig = carry;
! 620:
! 621: for (i = size-1; i >= 0; i--)
! 622: div1 (&rp[i], &carry, carry, sp[i], divisor);
! 623:
! 624: /* check by multiplying back */
! 625: #if 0
! 626: printf ("size=%ld divisor=0x%lX carry=0x%lX remainder=0x%lX\n",
! 627: size, divisor, carry_orig, carry);
! 628: mpn_trace("s",sp_copy,size);
! 629: mpn_trace("r",rp,size);
! 630: printf ("mul_1c %lX\n", refmpn_mul_1c (prod, rp, size, divisor, carry));
! 631: mpn_trace("p",prod,size);
! 632: #endif
! 633: ASSERT (refmpn_mul_1c (prod, rp, size, divisor, carry) == carry_orig);
! 634: ASSERT (refmpn_cmp (prod, sp_orig, size) == 0);
! 635: free (sp_orig);
! 636: free (prod);
! 637:
! 638: return carry;
! 639: }
! 640:
! 641: mp_limb_t
! 642: refmpn_divmod_1 (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t divisor)
! 643: {
! 644: return refmpn_divmod_1c (rp, sp, size, divisor, 0);
! 645: }
! 646:
! 647:
! 648: mp_limb_t
! 649: refmpn_mod_1c (mp_srcptr sp, mp_size_t size, mp_limb_t divisor,
! 650: mp_limb_t carry)
! 651: {
! 652: mp_ptr p = refmpn_malloc_limbs (size);
! 653: carry = refmpn_divmod_1c (p, sp, size, divisor, carry);
! 654: free (p);
! 655: return carry;
! 656: }
! 657:
! 658: mp_limb_t
! 659: refmpn_mod_1 (mp_srcptr sp, mp_size_t size, mp_limb_t divisor)
! 660: {
! 661: return refmpn_mod_1c (sp, size, divisor, 0);
! 662: }
! 663:
! 664: mp_limb_t
! 665: refmpn_mod_1_rshift (mp_srcptr sp, mp_size_t size,
! 666: unsigned shift, mp_limb_t divisor)
! 667: {
! 668: mp_limb_t r;
! 669: mp_ptr p = refmpn_malloc_limbs (size);
! 670: refmpn_rshift (p, sp, size, shift);
! 671: r = refmpn_mod_1 (p, size, divisor);
! 672: free (p);
! 673: return r;
! 674: }
! 675:
! 676:
! 677: mp_limb_t
! 678: refmpn_divrem_1c (mp_ptr rp, mp_size_t xsize,
! 679: mp_srcptr sp, mp_size_t size, mp_limb_t divisor,
! 680: mp_limb_t carry)
! 681: {
! 682: mp_ptr z;
! 683:
! 684: z = refmpn_malloc_limbs (xsize);
! 685: refmpn_fill (z, xsize, 0);
! 686:
! 687: carry = refmpn_divmod_1c (rp+xsize, sp, size, divisor, carry);
! 688: carry = refmpn_divmod_1c (rp, z, xsize, divisor, carry);
! 689:
! 690: free (z);
! 691: return carry;
! 692: }
! 693:
! 694: mp_limb_t
! 695: refmpn_divrem_1 (mp_ptr rp, mp_size_t xsize,
! 696: mp_srcptr sp, mp_size_t size, mp_limb_t divisor)
! 697: {
! 698: return refmpn_divrem_1c (rp, xsize, sp, size, divisor, 0);
! 699: }
! 700:
! 701:
! 702: /* As given in the manual, the divexact method gives quotient q and return
! 703: value c satisfying
! 704:
! 705: c*b^n + a-i == 3*q
! 706:
! 707: where a=dividend, i=initial carry, b=2^BITS_PER_MP_LIMB, and n=size.
! 708:
! 709: If a-i is divisible by 3 then c==0 and a plain divmod gives the quotient.
! 710: If (a-i)%3==r then c is a high limb tacked on that will turn r into 0.
! 711: Because 2^BITS_PER_MP_LIMB==1mod3 (so long as BITS_PER_MP_LIMB is even)
! 712: it's enough to set c=3-r, ie. if r=1 then c=2, or if r=2 then c=1.
! 713:
! 714: If a-i produces a borrow then refmpn_sub_1 leaves a twos complement
! 715: negative, ie. b^n+a-i, and the calculation produces c1 satisfying
! 716:
! 717: c1*b^n + b^n+a-i == 3*q
! 718:
! 719: From which clearly c=c1+1, so it's enough to just add any borrow to the
! 720: return value otherwise calculated.
! 721:
! 722: A borrow only occurs when a==0 or a==1, and, by the same reasoning as in
! 723: mpn/generic/diveby3.c, the c1 that results in those cases will only be 0
! 724: or 1 respectively, so with 1 added the final return value is still in the
! 725: prescribed range 0 to 2. */
! 726:
! 727: mp_limb_t
! 728: refmpn_divexact_by3c (mp_ptr rp, mp_srcptr sp, mp_size_t size, mp_limb_t carry)
! 729: {
! 730: mp_ptr spcopy;
! 731: mp_limb_t c, cs;
! 732:
! 733: ASSERT (refmpn_overlap_fullonly_p (rp, sp, size));
! 734: ASSERT (size >= 1);
! 735: ASSERT (carry <= 2);
! 736:
! 737: spcopy = refmpn_memdup_limbs (sp, size);
! 738: cs = refmpn_sub_1 (spcopy, spcopy, size, carry);
! 739:
! 740: c = refmpn_divmod_1 (rp, spcopy, size, 3);
! 741: if (c != 0)
! 742: {
! 743: ASSERT ((BITS_PER_MP_LIMB % 2) == 0);
! 744: c = 3-c;
! 745: ASSERT_NOCARRY (refmpn_divmod_1c (rp, spcopy, size, 3, c));
! 746: }
! 747:
! 748: c += cs;
! 749: ASSERT (c <= 2);
! 750:
! 751: free (spcopy);
! 752: return c;
! 753: }
! 754:
! 755: mp_limb_t
! 756: refmpn_divexact_by3 (mp_ptr rp, mp_srcptr sp, mp_size_t size)
! 757: {
! 758: return refmpn_divexact_by3c (rp, sp, size, 0);
! 759: }
! 760:
! 761:
! 762: /* The same as mpn/generic/mul_basecase.c, but using refmpn functions. */
! 763: void
! 764: refmpn_mul_basecase (mp_ptr prodp,
! 765: mp_srcptr up, mp_size_t usize,
! 766: mp_srcptr vp, mp_size_t vsize)
! 767: {
! 768: mp_size_t i;
! 769:
! 770: ASSERT (! MPN_OVERLAP_P (prodp, usize+vsize, up, usize));
! 771: ASSERT (! MPN_OVERLAP_P (prodp, usize+vsize, vp, vsize));
! 772: ASSERT (usize >= vsize);
! 773: ASSERT (vsize >= 1);
! 774:
! 775: prodp[usize] = refmpn_mul_1 (prodp, up, usize, vp[0]);
! 776: for (i = 1; i < vsize; i++)
! 777: prodp[usize+i] = refmpn_addmul_1 (prodp+i, up, usize, vp[i]);
! 778: }
! 779:
! 780: void
! 781: refmpn_mul_n (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
! 782: {
! 783: refmpn_mul_basecase (prodp, up, size, vp, size);
! 784: }
! 785:
! 786: void
! 787: refmpn_sqr (mp_ptr dst, mp_srcptr src, mp_size_t size)
! 788: {
! 789: refmpn_mul_basecase (dst, src, size, src, size);
! 790: }
! 791:
! 792:
! 793: mp_limb_t
! 794: refmpn_gcd_1 (mp_srcptr xp, mp_size_t xsize, mp_limb_t y)
! 795: {
! 796: mp_limb_t x;
! 797: int twos;
! 798:
! 799: ASSERT (y != 0);
! 800: ASSERT (! refmpn_zero_p (xp, xsize));
! 801:
! 802: x = mpn_mod_1 (xp, xsize, y);
! 803: if (x == 0)
! 804: return y;
! 805:
! 806: twos = 0;
! 807: while ((x & 1) == 0 && (y & 1) == 0)
! 808: {
! 809: x >>= 1;
! 810: y >>= 1;
! 811: twos++;
! 812: }
! 813:
! 814: for (;;)
! 815: {
! 816: while ((x & 1) == 0) x >>= 1;
! 817: while ((y & 1) == 0) y >>= 1;
! 818:
! 819: if (x < y)
! 820: MP_LIMB_T_SWAP (x, y);
! 821:
! 822: x -= y;
! 823: if (x == 0)
! 824: break;
! 825: }
! 826:
! 827: return y << twos;
! 828: }
! 829:
! 830:
! 831: unsigned
! 832: refmpn_count_trailing_zeros (mp_limb_t x)
! 833: {
! 834: unsigned n = 0;
! 835:
! 836: ASSERT (x != 0);
! 837: while ((x & 1) == 0)
! 838: {
! 839: x >>= 1;
! 840: n++;
! 841: }
! 842: return n;
! 843: }
! 844:
! 845: mp_size_t
! 846: refmpn_strip_twos (mp_ptr p, mp_size_t size)
! 847: {
! 848: mp_size_t limbs;
! 849: unsigned shift;
! 850:
! 851: ASSERT (size >= 1);
! 852: ASSERT (! refmpn_zero_p (p, size));
! 853:
! 854: for (limbs = 0; p[0] == 0; limbs++)
! 855: {
! 856: MPN_COPY_INCR (p, p+1, size-1);
! 857: p[size-1] = 0;
! 858: }
! 859:
! 860: shift = refmpn_count_trailing_zeros (p[0]);
! 861: if (shift)
! 862: refmpn_rshift (p, p, size, shift);
! 863:
! 864: return limbs*BITS_PER_MP_LIMB + shift;
! 865: }
! 866:
! 867: mp_limb_t
! 868: refmpn_gcd (mp_ptr gp, mp_ptr xp, mp_size_t xsize, mp_ptr yp, mp_size_t ysize)
! 869: {
! 870: int cmp;
! 871:
! 872: ASSERT (ysize >= 1);
! 873: ASSERT (xsize >= ysize);
! 874: ASSERT ((xp[0] & 1) != 0);
! 875: ASSERT ((yp[0] & 1) != 0);
! 876: ASSERT (xp[xsize-1] != 0);
! 877: ASSERT (yp[ysize-1] != 0);
! 878: ASSERT (refmpn_overlap_fullonly_p (gp, xp, xsize));
! 879: ASSERT (refmpn_overlap_fullonly_p (gp, yp, ysize));
! 880: ASSERT (! MPN_OVERLAP_P (xp, xsize, yp, ysize));
! 881: if (xsize == ysize)
! 882: ASSERT (refmpn_msbone (xp[xsize-1]) >= refmpn_msbone (yp[ysize-1]));
! 883:
! 884: refmpn_strip_twos (xp, xsize);
! 885: MPN_NORMALIZE (xp, xsize);
! 886: MPN_NORMALIZE (yp, ysize);
! 887:
! 888: for (;;)
! 889: {
! 890: cmp = refmpn_cmp_twosizes (xp, xsize, yp, ysize);
! 891: if (cmp == 0)
! 892: break;
! 893: if (cmp < 0)
! 894: MPN_PTR_SWAP (xp,xsize, yp,ysize);
! 895:
! 896: ASSERT_NOCARRY (refmpn_sub (xp, xp, xsize, yp, ysize));
! 897:
! 898: refmpn_strip_twos (xp, xsize);
! 899: MPN_NORMALIZE (xp, xsize);
! 900: }
! 901:
! 902: refmpn_copyi (gp, xp, xsize);
! 903: return xsize;
! 904: }
! 905:
! 906:
! 907: unsigned long
! 908: refmpn_popcount (mp_srcptr sp, mp_size_t size)
! 909: {
! 910: unsigned long count = 0;
! 911: mp_size_t i;
! 912: int j;
! 913: mp_limb_t l;
! 914:
! 915: ASSERT (size >= 0);
! 916: for (i = 0; i < size; i++)
! 917: {
! 918: l = sp[i];
! 919: for (j = 0; j < BITS_PER_MP_LIMB; j++)
! 920: {
! 921: count += (l & 1);
! 922: l >>= 1;
! 923: }
! 924: }
! 925: return count;
! 926: }
! 927:
! 928: unsigned long
! 929: refmpn_hamdist (mp_srcptr s1p, mp_srcptr s2p, mp_size_t size)
! 930: {
! 931: mp_ptr d;
! 932: unsigned long count;
! 933:
! 934: if (size == 0)
! 935: return 0;
! 936:
! 937: d = refmpn_malloc_limbs (size);
! 938: refmpn_xor_n (d, s1p, s2p, size);
! 939: count = refmpn_popcount (d, size);
! 940: free (d);
! 941: return count;
! 942: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ref.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpn / tests