Annotation of OpenXM_contrib/gmp/mpz/ui_pow_ui.c, Revision 1.1.1.2
1.1 maekawa 1: /* mpz_ui_pow_ui(res, base, exp) -- Set RES to BASE**EXP.
2:
1.1.1.2 ! maekawa 3: Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 Free Software Foundation,
! 4: Inc.
1.1 maekawa 5:
6: This file is part of the GNU MP Library.
7:
8: The GNU MP Library is free software; you can redistribute it and/or modify
1.1.1.2 ! maekawa 9: it under the terms of the GNU Lesser General Public License as published by
! 10: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1 maekawa 11: option) any later version.
12:
13: The GNU MP Library is distributed in the hope that it will be useful, but
14: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.1.1.2 ! maekawa 15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 16: License for more details.
17:
1.1.1.2 ! maekawa 18: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 19: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21: MA 02111-1307, USA. */
22:
23: #include "gmp.h"
24: #include "gmp-impl.h"
25: #include "longlong.h"
26:
1.1.1.2 ! maekawa 27:
! 28: static void mpz_pow2 _PROTO ((mpz_ptr r, mp_limb_t blimb, unsigned long int e, mp_limb_t rl));
! 29:
1.1 maekawa 30: void
31: #if __STDC__
32: mpz_ui_pow_ui (mpz_ptr r, unsigned long int b, unsigned long int e)
33: #else
34: mpz_ui_pow_ui (r, b, e)
35: mpz_ptr r;
36: unsigned long int b;
37: unsigned long int e;
38: #endif
39: {
40: mp_limb_t blimb = b;
1.1.1.2 ! maekawa 41: mp_limb_t rl;
1.1 maekawa 42:
43: if (e == 0)
44: {
1.1.1.2 ! maekawa 45: /* For x^0 we return 1, even if x is 0. */
1.1 maekawa 46: r->_mp_d[0] = 1;
47: r->_mp_size = 1;
48: return;
49: }
50:
1.1.1.2 ! maekawa 51: /* Compute b^e as (b^n)^(e div n) * b^(e mod n), where n is chosen such that
! 52: the latter factor is the largest number small enough to fit in a limb. */
! 53:
! 54: rl = 1;
! 55: while (e != 0 && blimb < ((mp_limb_t) 1 << BITS_PER_MP_LIMB/2))
1.1 maekawa 56: {
1.1.1.2 ! maekawa 57: if ((e & 1) != 0)
! 58: rl = rl * blimb;
! 59: blimb = blimb * blimb;
! 60: e = e >> 1;
1.1 maekawa 61: }
1.1.1.2 ! maekawa 62:
! 63: /* rl is now b^(e mod n). (I.e., the latter factor above.) */
! 64:
! 65: if (e == 0)
1.1 maekawa 66: {
1.1.1.2 ! maekawa 67: r->_mp_d[0] = rl;
! 68: r->_mp_size = rl != 0;
! 69: return;
1.1 maekawa 70: }
71:
1.1.1.2 ! maekawa 72: mpz_pow2 (r, blimb, e, rl);
! 73: }
! 74:
! 75: /* Multi-precision part of expontialization code. */
! 76: static void
! 77: #if __STDC__
! 78: mpz_pow2 (mpz_ptr r, mp_limb_t blimb, unsigned long int e, mp_limb_t rl)
! 79: #else
! 80: mpz_pow2 (r, blimb, e, rl)
! 81: mpz_ptr r;
! 82: mp_limb_t blimb;
! 83: unsigned long int e;
! 84: mp_limb_t rl;
! 85: #endif
! 86: {
! 87: mp_ptr rp, tp;
! 88: mp_size_t ralloc, rsize;
! 89: int cnt, i;
! 90: TMP_DECL (marker);
! 91:
1.1 maekawa 92: TMP_MARK (marker);
93:
1.1.1.2 ! maekawa 94: /* Over-estimate temporary space requirements somewhat. */
! 95: count_leading_zeros (cnt, blimb);
! 96: ralloc = e - cnt * e / BITS_PER_MP_LIMB + 1;
! 97:
! 98: /* The two areas are used to alternatingly hold the input and receive the
! 99: product for mpn_mul. (Needed since mpn_mul_n requires that the product
! 100: is distinct from either input operand.) */
! 101: rp = (mp_ptr) TMP_ALLOC (ralloc * BYTES_PER_MP_LIMB);
! 102: tp = (mp_ptr) TMP_ALLOC (ralloc * BYTES_PER_MP_LIMB);
1.1 maekawa 103:
104: rp[0] = blimb;
105: rsize = 1;
106:
1.1.1.2 ! maekawa 107: count_leading_zeros (cnt, e);
1.1 maekawa 108: for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--)
109: {
110: mpn_mul_n (tp, rp, rp, rsize);
111: rsize = 2 * rsize;
112: rsize -= tp[rsize - 1] == 0;
1.1.1.2 ! maekawa 113: MP_PTR_SWAP (rp, tp);
1.1 maekawa 114:
115: if ((e & ((mp_limb_t) 1 << i)) != 0)
116: {
117: mp_limb_t cy;
1.1.1.2 ! maekawa 118: cy = mpn_mul_1 (rp, rp, rsize, blimb);
! 119: rp[rsize] = cy;
! 120: rsize += cy != 0;
1.1 maekawa 121: }
122: }
123:
1.1.1.2 ! maekawa 124: /* We will need rsize or rsize+1 limbs for the result. */
! 125: if (r->_mp_alloc <= rsize)
! 126: _mpz_realloc (r, rsize + 1);
! 127:
! 128: /* Multiply the two factors (in rp,rsize and rl) and put the final result
! 129: in place. */
! 130: {
! 131: mp_limb_t cy;
! 132: cy = mpn_mul_1 (r->_mp_d, rp, rsize, rl);
! 133: (r->_mp_d)[rsize] = cy;
! 134: rsize += cy != 0;
! 135: }
1.1 maekawa 136:
137: r->_mp_size = rsize;
138: TMP_FREE (marker);
139: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>