--- /dev/null
+/* mpz_n_pow_ui -- mpn raised to ulong.
+
+ THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST
+ CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
+ FUTURE GNU MP RELEASES.
+
+Copyright 2001, 2002, 2005 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 3 of the License, or (at your
+option) any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+
+/* Change this to "#define TRACE(x) x" for some traces. */
+#define TRACE(x)
+
+
+/* Use this to test the mul_2 code on a CPU without a native version of that
+ routine. */
+#if 0
+#define mpn_mul_2 refmpn_mul_2
+#define HAVE_NATIVE_mpn_mul_2 1
+#endif
+
+
+/* mpz_pow_ui and mpz_ui_pow_ui want to share almost all of this code.
+ ui_pow_ui doesn't need the mpn_mul based powering loop or the tests on
+ bsize==2 or >2, but separating that isn't easy because there's shared
+ code both before and after (the size calculations and the powers of 2
+ handling).
+
+ Alternatives:
+
+ It would work to just use the mpn_mul powering loop for 1 and 2 limb
+ bases, but the current separate loop allows mul_1 and mul_2 to be done
+ in-place, which might help cache locality a bit. If mpn_mul was relaxed
+ to allow source==dest when vn==1 or 2 then some pointer twiddling might
+ let us get the same effect in one loop.
+
+ The initial powering for bsize==1 into blimb or blimb:blimb_low doesn't
+ form the biggest possible power of b that fits, only the biggest power of
+ 2 power, ie. b^(2^n). It'd be possible to choose a bigger power, perhaps
+ using __mp_bases[b].big_base for small b, and thereby get better value
+ from mpn_mul_1 or mpn_mul_2 in the bignum powering. It's felt that doing
+ so would be more complicated than it's worth, and could well end up being
+ a slowdown for small e. For big e on the other hand the algorithm is
+ dominated by mpn_sqr_n so there wouldn't much of a saving. The current
+ code can be viewed as simply doing the first few steps of the powering in
+ a single or double limb where possible.
+
+ If r==b, and blow_twos==0, and r must be realloc'ed, then the temporary
+ copy made of b is unnecessary. We could just use the old alloc'ed block
+ and free it at the end. But arranging this seems like a lot more trouble
+ than it's worth. */
+
+
+/* floor(sqrt(GMP_NUMB_MAX)), ie. the biggest value that can be squared in
+ a limb without overflowing.
+ FIXME: This formula is an underestimate when GMP_NUMB_BITS is odd. */
+
+#define GMP_NUMB_HALFMAX (((mp_limb_t) 1 << GMP_NUMB_BITS/2) - 1)
+
+
+/* The following are for convenience, they update the size and check the
+ alloc. */
+
+#define MPN_SQR_N(dst, alloc, src, size) \
+ do { \
+ ASSERT (2*(size) <= (alloc)); \
+ mpn_sqr_n (dst, src, size); \
+ (size) *= 2; \
+ (size) -= ((dst)[(size)-1] == 0); \
+ } while (0)
+
+#define MPN_MUL(dst, alloc, src, size, src2, size2) \
+ do { \
+ mp_limb_t cy; \
+ ASSERT ((size) + (size2) <= (alloc)); \
+ cy = mpn_mul (dst, src, size, src2, size2); \
+ (size) += (size2) - (cy == 0); \
+ } while (0)
+
+#define MPN_MUL_2(ptr, size, alloc, mult) \
+ do { \
+ mp_limb_t cy; \
+ ASSERT ((size)+2 <= (alloc)); \
+ cy = mpn_mul_2 (ptr, ptr, size, mult); \
+ (size)++; \
+ (ptr)[(size)] = cy; \
+ (size) += (cy != 0); \
+ } while (0)
+
+#define MPN_MUL_1(ptr, size, alloc, limb) \
+ do { \
+ mp_limb_t cy; \
+ ASSERT ((size)+1 <= (alloc)); \
+ cy = mpn_mul_1 (ptr, ptr, size, limb); \
+ (ptr)[size] = cy; \
+ (size) += (cy != 0); \
+ } while (0)
+
+#define MPN_LSHIFT(ptr, size, alloc, shift) \
+ do { \
+ mp_limb_t cy; \
+ ASSERT ((size)+1 <= (alloc)); \
+ cy = mpn_lshift (ptr, ptr, size, shift); \
+ (ptr)[size] = cy; \
+ (size) += (cy != 0); \
+ } while (0)
+
+#define MPN_RSHIFT_OR_COPY(dst, src, size, shift) \
+ do { \
+ if ((shift) == 0) \
+ MPN_COPY (dst, src, size); \
+ else \
+ { \
+ mpn_rshift (dst, src, size, shift); \
+ (size) -= ((dst)[(size)-1] == 0); \
+ } \
+ } while (0)
+
+
+/* ralloc and talloc are only wanted for ASSERTs, after the initial space
+ allocations. Avoid writing values to them in a normal build, to ensure
+ the compiler lets them go dead. gcc already figures this out itself
+ actually. */
+
+#define SWAP_RP_TP \
+ do { \
+ MP_PTR_SWAP (rp, tp); \
+ ASSERT_CODE (MP_SIZE_T_SWAP (ralloc, talloc)); \
+ } while (0)
+
+
+void
+mpz_n_pow_ui (mpz_ptr r, mp_srcptr bp, mp_size_t bsize, unsigned long int e)
+{
+ mp_ptr rp;
+ mp_size_t rtwos_limbs, ralloc, rsize;
+ int rneg, i, cnt, btwos, r_bp_overlap;
+ mp_limb_t blimb, rl;
+ unsigned long rtwos_bits;
+#if HAVE_NATIVE_mpn_mul_2
+ mp_limb_t blimb_low, rl_high;
+#else
+ mp_limb_t b_twolimbs[2];
+#endif
+ TMP_DECL;
+
+ TRACE (printf ("mpz_n_pow_ui rp=0x%lX bp=0x%lX bsize=%ld e=%lu (0x%lX)\n",
+ PTR(r), bp, bsize, e, e);
+ mpn_trace ("b", bp, bsize));
+
+ ASSERT (bsize == 0 || bp[ABS(bsize)-1] != 0);
+ ASSERT (MPN_SAME_OR_SEPARATE2_P (PTR(r), ABSIZ(r), bp, bsize));
+
+ /* b^0 == 1, including 0^0 == 1 */
+ if (e == 0)
+ {
+ PTR(r)[0] = 1;
+ SIZ(r) = 1;
+ return;
+ }
+
+ /* 0^e == 0 apart from 0^0 above */
+ if (bsize == 0)
+ {
+ SIZ(r) = 0;
+ return;
+ }
+
+ /* Sign of the final result. */
+ rneg = (bsize < 0 && (e & 1) != 0);
+ bsize = ABS (bsize);
+ TRACE (printf ("rneg %d\n", rneg));
+
+ r_bp_overlap = (PTR(r) == bp);
+
+ /* Strip low zero limbs from b. */
+ rtwos_limbs = 0;
+ for (blimb = *bp; blimb == 0; blimb = *++bp)
+ {
+ rtwos_limbs += e;
+ bsize--; ASSERT (bsize >= 1);
+ }
+ TRACE (printf ("trailing zero rtwos_limbs=%ld\n", rtwos_limbs));
+
+ /* Strip low zero bits from b. */
+ count_trailing_zeros (btwos, blimb);
+ blimb >>= btwos;
+ rtwos_bits = e * btwos;
+ rtwos_limbs += rtwos_bits / GMP_NUMB_BITS;
+ rtwos_bits %= GMP_NUMB_BITS;
+ TRACE (printf ("trailing zero btwos=%d rtwos_limbs=%ld rtwos_bits=%lu\n",
+ btwos, rtwos_limbs, rtwos_bits));
+
+ TMP_MARK;
+
+ rl = 1;
+#if HAVE_NATIVE_mpn_mul_2
+ rl_high = 0;
+#endif
+
+ if (bsize == 1)
+ {
+ bsize_1:
+ /* Power up as far as possible within blimb. We start here with e!=0,
+ but if e is small then we might reach e==0 and the whole b^e in rl.
+ Notice this code works when blimb==1 too, reaching e==0. */
+
+ while (blimb <= GMP_NUMB_HALFMAX)
+ {
+ TRACE (printf ("small e=0x%lX blimb=0x%lX rl=0x%lX\n",
+ e, blimb, rl));
+ ASSERT (e != 0);
+ if ((e & 1) != 0)
+ rl *= blimb;
+ e >>= 1;
+ if (e == 0)
+ goto got_rl;
+ blimb *= blimb;
+ }
+
+#if HAVE_NATIVE_mpn_mul_2
+ TRACE (printf ("single power, e=0x%lX b=0x%lX rl=0x%lX\n",
+ e, blimb, rl));
+
+ /* Can power b once more into blimb:blimb_low */
+ bsize = 2;
+ ASSERT (e != 0);
+ if ((e & 1) != 0)
+ {
+ umul_ppmm (rl_high, rl, rl, blimb << GMP_NAIL_BITS);
+ rl >>= GMP_NAIL_BITS;
+ }
+ e >>= 1;
+ umul_ppmm (blimb, blimb_low, blimb, blimb << GMP_NAIL_BITS);
+ blimb_low >>= GMP_NAIL_BITS;
+
+ got_rl:
+ TRACE (printf ("double power e=0x%lX blimb=0x%lX:0x%lX rl=0x%lX:%lX\n",
+ e, blimb, blimb_low, rl_high, rl));
+
+ /* Combine left-over rtwos_bits into rl_high:rl to be handled by the
+ final mul_1 or mul_2 rather than a separate lshift.
+ - rl_high:rl mustn't be 1 (since then there's no final mul)
+ - rl_high mustn't overflow
+ - rl_high mustn't change to non-zero, since mul_1+lshift is
+ probably faster than mul_2 (FIXME: is this true?) */
+
+ if (rtwos_bits != 0
+ && ! (rl_high == 0 && rl == 1)
+ && (rl_high >> (GMP_NUMB_BITS-rtwos_bits)) == 0)
+ {
+ mp_limb_t new_rl_high = (rl_high << rtwos_bits)
+ | (rl >> (GMP_NUMB_BITS-rtwos_bits));
+ if (! (rl_high == 0 && new_rl_high != 0))
+ {
+ rl_high = new_rl_high;
+ rl <<= rtwos_bits;
+ rtwos_bits = 0;
+ TRACE (printf ("merged rtwos_bits, rl=0x%lX:%lX\n",
+ rl_high, rl));
+ }
+ }
+#else
+ got_rl:
+ TRACE (printf ("small power e=0x%lX blimb=0x%lX rl=0x%lX\n",
+ e, blimb, rl));
+
+ /* Combine left-over rtwos_bits into rl to be handled by the final
+ mul_1 rather than a separate lshift.
+ - rl mustn't be 1 (since then there's no final mul)
+ - rl mustn't overflow */
+
+ if (rtwos_bits != 0
+ && rl != 1
+ && (rl >> (GMP_NUMB_BITS-rtwos_bits)) == 0)
+ {
+ rl <<= rtwos_bits;
+ rtwos_bits = 0;
+ TRACE (printf ("merged rtwos_bits, rl=0x%lX\n", rl));
+ }
+#endif
+ }
+ else if (bsize == 2)
+ {
+ mp_limb_t bsecond = bp[1];
+ if (btwos != 0)
+ blimb |= (bsecond << (GMP_NUMB_BITS - btwos)) & GMP_NUMB_MASK;
+ bsecond >>= btwos;
+ if (bsecond == 0)
+ {
+ /* Two limbs became one after rshift. */
+ bsize = 1;
+ goto bsize_1;
+ }
+
+ TRACE (printf ("bsize==2 using b=0x%lX:%lX", bsecond, blimb));
+#if HAVE_NATIVE_mpn_mul_2
+ blimb_low = blimb;
+#else
+ bp = b_twolimbs;
+ b_twolimbs[0] = blimb;
+ b_twolimbs[1] = bsecond;
+#endif
+ blimb = bsecond;
+ }
+ else
+ {
+ if (r_bp_overlap || btwos != 0)
+ {
+ mp_ptr tp = TMP_ALLOC_LIMBS (bsize);
+ MPN_RSHIFT_OR_COPY (tp, bp, bsize, btwos);
+ bp = tp;
+ TRACE (printf ("rshift or copy bp,bsize, new bsize=%ld\n", bsize));
+ }
+#if HAVE_NATIVE_mpn_mul_2
+ /* in case 3 limbs rshift to 2 and hence use the mul_2 loop below */
+ blimb_low = bp[0];
+#endif
+ blimb = bp[bsize-1];
+
+ TRACE (printf ("big bsize=%ld ", bsize);
+ mpn_trace ("b", bp, bsize));
+ }
+
+ /* At this point blimb is the most significant limb of the base to use.
+
+ Each factor of b takes (bsize*BPML-cnt) bits and there's e of them; +1
+ limb to round up the division; +1 for multiplies all using an extra
+ limb over the true size; +2 for rl at the end; +1 for lshift at the
+ end.
+
+ The size calculation here is reasonably accurate. The base is at least
+ half a limb, so in 32 bits the worst case is 2^16+1 treated as 17 bits
+ when it will power up as just over 16, an overestimate of 17/16 =
+ 6.25%. For a 64-bit limb it's half that.
+
+ If e==0 then blimb won't be anything useful (though it will be
+ non-zero), but that doesn't matter since we just end up with ralloc==5,
+ and that's fine for 2 limbs of rl and 1 of lshift. */
+
+ ASSERT (blimb != 0);
+ count_leading_zeros (cnt, blimb);
+ ralloc = (bsize*GMP_NUMB_BITS - cnt + GMP_NAIL_BITS) * e / GMP_NUMB_BITS + 5;
+ TRACE (printf ("ralloc %ld, from bsize=%ld blimb=0x%lX cnt=%d\n",
+ ralloc, bsize, blimb, cnt));
+ MPZ_REALLOC (r, ralloc + rtwos_limbs);
+ rp = PTR(r);
+
+ /* Low zero limbs resulting from powers of 2. */
+ MPN_ZERO (rp, rtwos_limbs);
+ rp += rtwos_limbs;
+
+ if (e == 0)
+ {
+ /* Any e==0 other than via bsize==1 or bsize==2 is covered at the
+ start. */
+ rp[0] = rl;
+ rsize = 1;
+#if HAVE_NATIVE_mpn_mul_2
+ rp[1] = rl_high;
+ rsize += (rl_high != 0);
+#endif
+ ASSERT (rp[rsize-1] != 0);
+ }
+ else
+ {
+ mp_ptr tp;
+ mp_size_t talloc;
+
+ /* In the mpn_mul_1 or mpn_mul_2 loops or in the mpn_mul loop when the
+ low bit of e is zero, tp only has to hold the second last power
+ step, which is half the size of the final result. There's no need
+ to round up the divide by 2, since ralloc includes a +2 for rl
+ which not needed by tp. In the mpn_mul loop when the low bit of e
+ is 1, tp must hold nearly the full result, so just size it the same
+ as rp. */
+
+ talloc = ralloc;
+#if HAVE_NATIVE_mpn_mul_2
+ if (bsize <= 2 || (e & 1) == 0)
+ talloc /= 2;
+#else
+ if (bsize <= 1 || (e & 1) == 0)
+ talloc /= 2;
+#endif
+ TRACE (printf ("talloc %ld\n", talloc));
+ tp = TMP_ALLOC_LIMBS (talloc);
+
+ /* Go from high to low over the bits of e, starting with i pointing at
+ the bit below the highest 1 (which will mean i==-1 if e==1). */
+ count_leading_zeros (cnt, e);
+ i = GMP_LIMB_BITS - cnt - 2;
+
+#if HAVE_NATIVE_mpn_mul_2
+ if (bsize <= 2)
+ {
+ mp_limb_t mult[2];
+
+ /* Any bsize==1 will have been powered above to be two limbs. */
+ ASSERT (bsize == 2);
+ ASSERT (blimb != 0);
+
+ /* Arrange the final result ends up in r, not in the temp space */
+ if ((i & 1) == 0)
+ SWAP_RP_TP;
+
+ rp[0] = blimb_low;
+ rp[1] = blimb;
+ rsize = 2;
+
+ mult[0] = blimb_low;
+ mult[1] = blimb;
+
+ for ( ; i >= 0; i--)
+ {
+ TRACE (printf ("mul_2 loop i=%d e=0x%lX, rsize=%ld ralloc=%ld talloc=%ld\n",
+ i, e, rsize, ralloc, talloc);
+ mpn_trace ("r", rp, rsize));
+
+ MPN_SQR_N (tp, talloc, rp, rsize);
+ SWAP_RP_TP;
+ if ((e & (1L << i)) != 0)
+ MPN_MUL_2 (rp, rsize, ralloc, mult);
+ }
+
+ TRACE (mpn_trace ("mul_2 before rl, r", rp, rsize));
+ if (rl_high != 0)
+ {
+ mult[0] = rl;
+ mult[1] = rl_high;
+ MPN_MUL_2 (rp, rsize, ralloc, mult);
+ }
+ else if (rl != 1)
+ MPN_MUL_1 (rp, rsize, ralloc, rl);
+ }
+#else
+ if (bsize == 1)
+ {
+ /* Arrange the final result ends up in r, not in the temp space */
+ if ((i & 1) == 0)
+ SWAP_RP_TP;
+
+ rp[0] = blimb;
+ rsize = 1;
+
+ for ( ; i >= 0; i--)
+ {
+ TRACE (printf ("mul_1 loop i=%d e=0x%lX, rsize=%ld ralloc=%ld talloc=%ld\n",
+ i, e, rsize, ralloc, talloc);
+ mpn_trace ("r", rp, rsize));
+
+ MPN_SQR_N (tp, talloc, rp, rsize);
+ SWAP_RP_TP;
+ if ((e & (1L << i)) != 0)
+ MPN_MUL_1 (rp, rsize, ralloc, blimb);
+ }
+
+ TRACE (mpn_trace ("mul_1 before rl, r", rp, rsize));
+ if (rl != 1)
+ MPN_MUL_1 (rp, rsize, ralloc, rl);
+ }
+#endif
+ else
+ {
+ int parity;
+
+ /* Arrange the final result ends up in r, not in the temp space */
+ ULONG_PARITY (parity, e);
+ if (((parity ^ i) & 1) != 0)
+ SWAP_RP_TP;
+
+ MPN_COPY (rp, bp, bsize);
+ rsize = bsize;
+
+ for ( ; i >= 0; i--)
+ {
+ TRACE (printf ("mul loop i=%d e=0x%lX, rsize=%ld ralloc=%ld talloc=%ld\n",
+ i, e, rsize, ralloc, talloc);
+ mpn_trace ("r", rp, rsize));
+
+ MPN_SQR_N (tp, talloc, rp, rsize);
+ SWAP_RP_TP;
+ if ((e & (1L << i)) != 0)
+ {
+ MPN_MUL (tp, talloc, rp, rsize, bp, bsize);
+ SWAP_RP_TP;
+ }
+ }
+ }
+ }
+
+ ASSERT (rp == PTR(r) + rtwos_limbs);
+ TRACE (mpn_trace ("end loop r", rp, rsize));
+ TMP_FREE;
+
+ /* Apply any partial limb factors of 2. */
+ if (rtwos_bits != 0)
+ {
+ MPN_LSHIFT (rp, rsize, ralloc, (unsigned) rtwos_bits);
+ TRACE (mpn_trace ("lshift r", rp, rsize));
+ }
+
+ rsize += rtwos_limbs;
+ SIZ(r) = (rneg ? -rsize : rsize);
+}
projects / msp430-gcc.git / blobdiff