Annotation of OpenXM_contrib/gmp/mpz/n_pow_ui.c, Revision 1.1
1.1 ! ohara 1: /* mpz_n_pow_ui -- mpn raised to ulong.
! 2:
! 3: Copyright 2001, 2002 Free Software Foundation, Inc.
! 4:
! 5: This file is part of the GNU MP Library.
! 6:
! 7: The GNU MP Library is free software; you can redistribute it and/or modify
! 8: it under the terms of the GNU Lesser General Public License as published by
! 9: the Free Software Foundation; either version 2.1 of the License, or (at your
! 10: option) any later version.
! 11:
! 12: The GNU MP Library is distributed in the hope that it will be useful, but
! 13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 15: License for more details.
! 16:
! 17: You should have received a copy of the GNU Lesser General Public License
! 18: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 20: MA 02111-1307, USA. */
! 21:
! 22: #include "gmp.h"
! 23: #include "gmp-impl.h"
! 24: #include "longlong.h"
! 25:
! 26:
! 27: /* Change this to "#define TRACE(x) x" for some traces. */
! 28: #define TRACE(x)
! 29:
! 30:
! 31: /* Use this to test the mul_2 code on a CPU without a native version of that
! 32: routine. */
! 33: #if 0
! 34: #define mpn_mul_2 refmpn_mul_2
! 35: #define HAVE_NATIVE_mpn_mul_2 1
! 36: #endif
! 37:
! 38:
! 39: /* mpz_pow_ui and mpz_ui_pow_ui want to share almost all of this code.
! 40: ui_pow_ui doesn't need the mpn_mul based powering loop or the tests on
! 41: bsize==2 or >2, but separating that isn't easy because there's shared
! 42: code both before and after (the size calculations and the powers of 2
! 43: handling).
! 44:
! 45: Alternatives:
! 46:
! 47: It would work to just use the mpn_mul powering loop for 1 and 2 limb
! 48: bases, but the current separate loop allows mul_1 and mul_2 to be done
! 49: in-place, which might help cache locality a bit. If mpn_mul was relaxed
! 50: to allow source==dest when vn==1 or 2 then some pointer twiddling might
! 51: let us get the same effect in one loop.
! 52:
! 53: The initial powering for bsize==1 into blimb or blimb:blimb_low doesn't
! 54: form the biggest possible power of b that fits, only the biggest power of
! 55: 2 power, ie. b^(2^n). It'd be possible to choose a bigger power, perhaps
! 56: using __mp_bases[b].big_base for small b, and thereby get better value
! 57: from mpn_mul_1 or mpn_mul_2 in the bignum powering. It's felt that doing
! 58: so would be more complicated than it's worth, and could well end up being
! 59: a slowdown for small e. For big e on the other hand the algorithm is
! 60: dominated by mpn_sqr_n so there wouldn't much of a saving. The current
! 61: code can be viewed as simply doing the first few steps of the powering in
! 62: a single or double limb where possible.
! 63:
! 64: If r==b, and blow_twos==0, and r must be realloc'ed, then the temporary
! 65: copy made of b is unnecessary. We could just use the old alloc'ed block
! 66: and free it at the end. But arranging this seems like a lot more trouble
! 67: than it's worth. */
! 68:
! 69:
! 70: /* floor(sqrt(GMP_NUMB_MAX)), ie. the biggest value that can be squared in
! 71: a limb without overflowing.
! 72: FIXME: This formula is an underestimate when GMP_NUMB_BITS is odd. */
! 73:
! 74: #define GMP_NUMB_HALFMAX (((mp_limb_t) 1 << GMP_NUMB_BITS/2) - 1)
! 75:
! 76:
! 77: /* The following are for convenience, they update the size and check the
! 78: alloc. */
! 79:
! 80: #define MPN_SQR_N(dst, alloc, src, size) \
! 81: do { \
! 82: ASSERT (2*(size) <= (alloc)); \
! 83: mpn_sqr_n (dst, src, size); \
! 84: (size) *= 2; \
! 85: (size) -= ((dst)[(size)-1] == 0); \
! 86: } while (0)
! 87:
! 88: #define MPN_MUL(dst, alloc, src, size, src2, size2) \
! 89: do { \
! 90: mp_limb_t cy; \
! 91: ASSERT ((size) + (size2) <= (alloc)); \
! 92: cy = mpn_mul (dst, src, size, src2, size2); \
! 93: (size) += (size2) - (cy == 0); \
! 94: } while (0)
! 95:
! 96: #define MPN_MUL_2(ptr, size, alloc, mult) \
! 97: do { \
! 98: mp_limb_t cy; \
! 99: ASSERT ((size)+2 <= (alloc)); \
! 100: cy = mpn_mul_2 (ptr, ptr, size, mult); \
! 101: (size)++; \
! 102: (ptr)[(size)] = cy; \
! 103: (size) += (cy != 0); \
! 104: } while (0)
! 105:
! 106: #define MPN_MUL_1(ptr, size, alloc, limb) \
! 107: do { \
! 108: mp_limb_t cy; \
! 109: ASSERT ((size)+1 <= (alloc)); \
! 110: cy = mpn_mul_1 (ptr, ptr, size, limb); \
! 111: (ptr)[size] = cy; \
! 112: (size) += (cy != 0); \
! 113: } while (0)
! 114:
! 115: #define MPN_LSHIFT(ptr, size, alloc, shift) \
! 116: do { \
! 117: mp_limb_t cy; \
! 118: ASSERT ((size)+1 <= (alloc)); \
! 119: cy = mpn_lshift (ptr, ptr, size, shift); \
! 120: (ptr)[size] = cy; \
! 121: (size) += (cy != 0); \
! 122: } while (0)
! 123:
! 124: #define MPN_RSHIFT_OR_COPY(dst, src, size, shift) \
! 125: do { \
! 126: if ((shift) == 0) \
! 127: MPN_COPY (dst, src, size); \
! 128: else \
! 129: { \
! 130: mpn_rshift (dst, src, size, shift); \
! 131: (size) -= ((dst)[(size)-1] == 0); \
! 132: } \
! 133: } while (0)
! 134:
! 135:
! 136: /* ralloc and talloc are only wanted for ASSERTs, after the initial space
! 137: allocations. Avoid writing values to them in a normal build, to ensure
! 138: the compiler lets them go dead. gcc already figures this out itself
! 139: actually. */
! 140:
! 141: #define SWAP_RP_TP \
! 142: do { \
! 143: MP_PTR_SWAP (rp, tp); \
! 144: ASSERT_CODE (MP_SIZE_T_SWAP (ralloc, talloc)); \
! 145: } while (0)
! 146:
! 147:
! 148: void
! 149: mpz_n_pow_ui (mpz_ptr r, mp_srcptr bp, mp_size_t bsize, unsigned long int e)
! 150: {
! 151: mp_ptr rp;
! 152: mp_size_t rtwos_limbs, ralloc, rsize;
! 153: int rneg, i, cnt, btwos, r_bp_overlap;
! 154: mp_limb_t blimb, rl;
! 155: unsigned long rtwos_bits;
! 156: #if HAVE_NATIVE_mpn_mul_2
! 157: mp_limb_t blimb_low, rl_high;
! 158: #else
! 159: mp_limb_t b_twolimbs[2];
! 160: #endif
! 161: TMP_DECL (marker);
! 162:
! 163: TRACE (printf ("mpz_n_pow_ui rp=0x%lX bp=0x%lX bsize=%ld e=%lu (0x%lX)\n",
! 164: PTR(r), bp, bsize, e, e);
! 165: mpn_trace ("b", bp, bsize));
! 166:
! 167: ASSERT (bsize == 0 || bp[ABS(bsize)-1] != 0);
! 168: ASSERT (MPN_SAME_OR_SEPARATE2_P (PTR(r), ABSIZ(r), bp, bsize));
! 169:
! 170: /* b^0 == 1, including 0^0 == 1 */
! 171: if (e == 0)
! 172: {
! 173: PTR(r)[0] = 1;
! 174: SIZ(r) = 1;
! 175: return;
! 176: }
! 177:
! 178: /* 0^e == 0 apart from 0^0 above */
! 179: if (bsize == 0)
! 180: {
! 181: SIZ(r) = 0;
! 182: return;
! 183: }
! 184:
! 185: /* Sign of the final result. */
! 186: rneg = (bsize < 0 && (e & 1) != 0);
! 187: bsize = ABS (bsize);
! 188: TRACE (printf ("rneg %d\n", rneg));
! 189:
! 190: r_bp_overlap = (PTR(r) == bp);
! 191:
! 192: /* Strip low zero limbs from b. */
! 193: rtwos_limbs = 0;
! 194: for (blimb = *bp; blimb == 0; blimb = *++bp)
! 195: {
! 196: rtwos_limbs += e;
! 197: bsize--; ASSERT (bsize >= 1);
! 198: }
! 199: TRACE (printf ("trailing zero rtwos_limbs=%ld\n", rtwos_limbs));
! 200:
! 201: /* Strip low zero bits from b. */
! 202: count_trailing_zeros (btwos, blimb);
! 203: blimb >>= btwos;
! 204: rtwos_bits = e * btwos;
! 205: rtwos_limbs += rtwos_bits / GMP_NUMB_BITS;
! 206: rtwos_bits %= GMP_NUMB_BITS;
! 207: TRACE (printf ("trailing zero btwos=%d rtwos_limbs=%ld rtwos_bits=%lu\n",
! 208: btwos, rtwos_limbs, rtwos_bits));
! 209:
! 210: TMP_MARK (marker);
! 211:
! 212: rl = 1;
! 213: #if HAVE_NATIVE_mpn_mul_2
! 214: rl_high = 0;
! 215: #endif
! 216:
! 217: if (bsize == 1)
! 218: {
! 219: bsize_1:
! 220: /* Power up as far as possible within blimb. We start here with e!=0,
! 221: but if e is small then we might reach e==0 and the whole b^e in rl.
! 222: Notice this code works when blimb==1 too, reaching e==0. */
! 223:
! 224: while (blimb <= GMP_NUMB_HALFMAX)
! 225: {
! 226: TRACE (printf ("small e=0x%lX blimb=0x%lX rl=0x%lX\n",
! 227: e, blimb, rl));
! 228: ASSERT (e != 0);
! 229: if ((e & 1) != 0)
! 230: rl *= blimb;
! 231: e >>= 1;
! 232: if (e == 0)
! 233: goto got_rl;
! 234: blimb *= blimb;
! 235: }
! 236:
! 237: #if HAVE_NATIVE_mpn_mul_2
! 238: TRACE (printf ("single power, e=0x%lX b=0x%lX rl=0x%lX\n",
! 239: e, blimb, rl));
! 240:
! 241: /* Can power b once more into blimb:blimb_low */
! 242: bsize = 2;
! 243: ASSERT (e != 0);
! 244: if ((e & 1) != 0)
! 245: {
! 246: umul_ppmm (rl_high, rl, rl, blimb << GMP_NAIL_BITS);
! 247: rl >>= GMP_NAIL_BITS;
! 248: }
! 249: e >>= 1;
! 250: umul_ppmm (blimb, blimb_low, blimb, blimb << GMP_NAIL_BITS);
! 251: blimb_low >>= GMP_NAIL_BITS;
! 252:
! 253: got_rl:
! 254: TRACE (printf ("double power e=0x%lX blimb=0x%lX:0x%lX rl=0x%lX:%lX\n",
! 255: e, blimb, blimb_low, rl_high, rl));
! 256:
! 257: /* Combine left-over rtwos_bits into rl_high:rl to be handled by the
! 258: final mul_1 or mul_2 rather than a separate lshift.
! 259: - rl_high:rl mustn't be 1 (since then there's no final mul)
! 260: - rl_high mustn't overflow
! 261: - rl_high mustn't change to non-zero, since mul_1+lshift is
! 262: probably faster than mul_2 (FIXME: is this true?) */
! 263:
! 264: if (rtwos_bits != 0
! 265: && ! (rl_high == 0 && rl == 1)
! 266: && (rl_high >> (GMP_NUMB_BITS-rtwos_bits)) == 0)
! 267: {
! 268: mp_limb_t new_rl_high = (rl_high << rtwos_bits)
! 269: | (rl >> (GMP_NUMB_BITS-rtwos_bits));
! 270: if (! (rl_high == 0 && new_rl_high != 0))
! 271: {
! 272: rl_high = new_rl_high;
! 273: rl <<= rtwos_bits;
! 274: rtwos_bits = 0;
! 275: TRACE (printf ("merged rtwos_bits, rl=0x%lX:%lX\n",
! 276: rl_high, rl));
! 277: }
! 278: }
! 279: #else
! 280: got_rl:
! 281: TRACE (printf ("small power e=0x%lX blimb=0x%lX rl=0x%lX\n",
! 282: e, blimb, rl));
! 283:
! 284: /* Combine left-over rtwos_bits into rl to be handled by the final
! 285: mul_1 rather than a separate lshift.
! 286: - rl mustn't be 1 (since then there's no final mul)
! 287: - rl mustn't overflow */
! 288:
! 289: if (rtwos_bits != 0
! 290: && rl != 1
! 291: && (rl >> (GMP_NUMB_BITS-rtwos_bits)) == 0)
! 292: {
! 293: rl <<= rtwos_bits;
! 294: rtwos_bits = 0;
! 295: TRACE (printf ("merged rtwos_bits, rl=0x%lX\n", rl));
! 296: }
! 297: #endif
! 298: }
! 299: else if (bsize == 2)
! 300: {
! 301: mp_limb_t bsecond = bp[1];
! 302: if (btwos != 0)
! 303: blimb |= (bsecond << (GMP_NUMB_BITS - btwos)) & GMP_NUMB_MASK;
! 304: bsecond >>= btwos;
! 305: if (bsecond == 0)
! 306: {
! 307: /* Two limbs became one after rshift. */
! 308: bsize = 1;
! 309: goto bsize_1;
! 310: }
! 311:
! 312: TRACE (printf ("bsize==2 using b=0x%lX:%lX", bsecond, blimb));
! 313: #if HAVE_NATIVE_mpn_mul_2
! 314: blimb_low = blimb;
! 315: #else
! 316: bp = b_twolimbs;
! 317: b_twolimbs[0] = blimb;
! 318: b_twolimbs[1] = bsecond;
! 319: #endif
! 320: blimb = bsecond;
! 321: }
! 322: else
! 323: {
! 324: if (r_bp_overlap || btwos != 0)
! 325: {
! 326: mp_ptr tp = TMP_ALLOC_LIMBS (bsize);
! 327: MPN_RSHIFT_OR_COPY (tp, bp, bsize, btwos);
! 328: bp = tp;
! 329: TRACE (printf ("rshift or copy bp,bsize, new bsize=%ld\n", bsize));
! 330: }
! 331: #if HAVE_NATIVE_mpn_mul_2
! 332: /* in case 3 limbs rshift to 2 and hence use the mul_2 loop below */
! 333: blimb_low = bp[0];
! 334: #endif
! 335: blimb = bp[bsize-1];
! 336:
! 337: TRACE (printf ("big bsize=%ld ", bsize);
! 338: mpn_trace ("b", bp, bsize));
! 339: }
! 340:
! 341: /* At this point blimb is the most significant limb of the base to use.
! 342:
! 343: Each factor of b takes (bsize*BPML-cnt) bits and there's e of them; +1
! 344: limb to round up the division; +1 for multiplies all using an extra
! 345: limb over the true size; +2 for rl at the end; +1 for lshift at the
! 346: end.
! 347:
! 348: The size calculation here is reasonably accurate. The base is at least
! 349: half a limb, so in 32 bits the worst case is 2^16+1 treated as 17 bits
! 350: when it will power up as just over 16, an overestimate of 17/16 =
! 351: 6.25%. For a 64-bit limb it's half that.
! 352:
! 353: If e==0 then blimb won't be anything useful (though it will be
! 354: non-zero), but that doesn't matter since we just end up with ralloc==5,
! 355: and that's fine for 2 limbs of rl and 1 of lshift. */
! 356:
! 357: ASSERT (blimb != 0);
! 358: count_leading_zeros (cnt, blimb);
! 359: ralloc = (bsize*GMP_NUMB_BITS - cnt + GMP_NAIL_BITS) * e / GMP_NUMB_BITS + 5;
! 360: TRACE (printf ("ralloc %ld, from bsize=%ld blimb=0x%lX cnt=%d\n",
! 361: ralloc, bsize, blimb, cnt));
! 362: MPZ_REALLOC (r, ralloc + rtwos_limbs);
! 363: rp = PTR(r);
! 364:
! 365: /* Low zero limbs resulting from powers of 2. */
! 366: MPN_ZERO (rp, rtwos_limbs);
! 367: rp += rtwos_limbs;
! 368:
! 369: if (e == 0)
! 370: {
! 371: /* Any e==0 other than via bsize==1 or bsize==2 is covered at the
! 372: start. */
! 373: rp[0] = rl;
! 374: rsize = 1;
! 375: #if HAVE_NATIVE_mpn_mul_2
! 376: rp[1] = rl_high;
! 377: rsize += (rl_high != 0);
! 378: #endif
! 379: ASSERT (rp[rsize-1] != 0);
! 380: }
! 381: else
! 382: {
! 383: mp_ptr tp;
! 384: mp_size_t talloc;
! 385:
! 386: /* In the mpn_mul_1 or mpn_mul_2 loops or in the mpn_mul loop when the
! 387: low bit of e is zero, tp only has to hold the second last power
! 388: step, which is half the size of the final result. There's no need
! 389: to round up the divide by 2, since ralloc includes a +2 for rl
! 390: which not needed by tp. In the mpn_mul loop when the low bit of e
! 391: is 1, tp must hold nearly the full result, so just size it the same
! 392: as rp. */
! 393:
! 394: talloc = ralloc;
! 395: #if HAVE_NATIVE_mpn_mul_2
! 396: if (bsize <= 2 || (e & 1) == 0)
! 397: talloc /= 2;
! 398: #else
! 399: if (bsize <= 1 || (e & 1) == 0)
! 400: talloc /= 2;
! 401: #endif
! 402: TRACE (printf ("talloc %ld\n", talloc));
! 403: tp = TMP_ALLOC_LIMBS (talloc);
! 404:
! 405: /* Go from high to low over the bits of e, starting with i pointing at
! 406: the bit below the highest 1 (which will mean i==-1 if e==1). */
! 407: count_leading_zeros (cnt, e);
! 408: i = GMP_LIMB_BITS - cnt - 2;
! 409:
! 410: #if HAVE_NATIVE_mpn_mul_2
! 411: if (bsize <= 2)
! 412: {
! 413: mp_limb_t mult[2];
! 414:
! 415: /* Any bsize==1 will have been powered above to be two limbs. */
! 416: ASSERT (bsize == 2);
! 417: ASSERT (blimb != 0);
! 418:
! 419: /* Arrange the final result ends up in r, not in the temp space */
! 420: if ((i & 1) == 0)
! 421: SWAP_RP_TP;
! 422:
! 423: rp[0] = blimb_low;
! 424: rp[1] = blimb;
! 425: rsize = 2;
! 426:
! 427: mult[0] = blimb_low;
! 428: mult[1] = blimb;
! 429:
! 430: for ( ; i >= 0; i--)
! 431: {
! 432: TRACE (printf ("mul_2 loop i=%d e=0x%lX, rsize=%ld ralloc=%ld talloc=%ld\n",
! 433: i, e, rsize, ralloc, talloc);
! 434: mpn_trace ("r", rp, rsize));
! 435:
! 436: MPN_SQR_N (tp, talloc, rp, rsize);
! 437: SWAP_RP_TP;
! 438: if ((e & (1L << i)) != 0)
! 439: MPN_MUL_2 (rp, rsize, ralloc, mult);
! 440: }
! 441:
! 442: TRACE (mpn_trace ("mul_2 before rl, r", rp, rsize));
! 443: if (rl_high != 0)
! 444: {
! 445: mult[0] = rl;
! 446: mult[1] = rl_high;
! 447: MPN_MUL_2 (rp, rsize, ralloc, mult);
! 448: }
! 449: else if (rl != 1)
! 450: MPN_MUL_1 (rp, rsize, ralloc, rl);
! 451: }
! 452: #else
! 453: if (bsize == 1)
! 454: {
! 455: /* Arrange the final result ends up in r, not in the temp space */
! 456: if ((i & 1) == 0)
! 457: SWAP_RP_TP;
! 458:
! 459: rp[0] = blimb;
! 460: rsize = 1;
! 461:
! 462: for ( ; i >= 0; i--)
! 463: {
! 464: TRACE (printf ("mul_1 loop i=%d e=0x%lX, rsize=%ld ralloc=%ld talloc=%ld\n",
! 465: i, e, rsize, ralloc, talloc);
! 466: mpn_trace ("r", rp, rsize));
! 467:
! 468: MPN_SQR_N (tp, talloc, rp, rsize);
! 469: SWAP_RP_TP;
! 470: if ((e & (1L << i)) != 0)
! 471: MPN_MUL_1 (rp, rsize, ralloc, blimb);
! 472: }
! 473:
! 474: TRACE (mpn_trace ("mul_1 before rl, r", rp, rsize));
! 475: if (rl != 1)
! 476: MPN_MUL_1 (rp, rsize, ralloc, rl);
! 477: }
! 478: #endif
! 479: else
! 480: {
! 481: int parity;
! 482:
! 483: /* Arrange the final result ends up in r, not in the temp space */
! 484: ULONG_PARITY (parity, e);
! 485: if (((parity ^ i) & 1) != 0)
! 486: SWAP_RP_TP;
! 487:
! 488: MPN_COPY (rp, bp, bsize);
! 489: rsize = bsize;
! 490:
! 491: for ( ; i >= 0; i--)
! 492: {
! 493: TRACE (printf ("mul loop i=%d e=0x%lX, rsize=%ld ralloc=%ld talloc=%ld\n",
! 494: i, e, rsize, ralloc, talloc);
! 495: mpn_trace ("r", rp, rsize));
! 496:
! 497: MPN_SQR_N (tp, talloc, rp, rsize);
! 498: SWAP_RP_TP;
! 499: if ((e & (1L << i)) != 0)
! 500: {
! 501: MPN_MUL (tp, talloc, rp, rsize, bp, bsize);
! 502: SWAP_RP_TP;
! 503: }
! 504: }
! 505: }
! 506: }
! 507:
! 508: ASSERT (rp == PTR(r) + rtwos_limbs);
! 509: TRACE (mpn_trace ("end loop r", rp, rsize));
! 510: TMP_FREE (marker);
! 511:
! 512: /* Apply any partial limb factors of 2. */
! 513: if (rtwos_bits != 0)
! 514: {
! 515: MPN_LSHIFT (rp, rsize, ralloc, (unsigned) rtwos_bits);
! 516: TRACE (mpn_trace ("lshift r", rp, rsize));
! 517: }
! 518:
! 519: rsize += rtwos_limbs;
! 520: SIZ(r) = (rneg ? -rsize : rsize);
! 521: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to n_pow_ui.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpz