Return to pow_ui.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpz

Annotation of OpenXM_contrib/gmp/mpz/pow_ui.c, Revision 1.1.1.1

1.1       maekawa     1: /* mpz_pow_ui(res, base, exp) -- Set RES to BASE**EXP.
                      2:
                      3: Copyright (C) 1991, 1993, 1994, 1996 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 Library General Public License as published by
                      9: the Free Software Foundation; either version 2 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 Library General Public
                     15: License for more details.
                     16:
                     17: You should have received a copy of the GNU Library 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: #ifdef BERKELEY_MP
                     23: #include "mp.h"
                     24: #endif
                     25: #include "gmp.h"
                     26: #include "gmp-impl.h"
                     27: #include "longlong.h"
                     28:
                     29: #ifndef BERKELEY_MP
                     30: void
                     31: #if __STDC__
                     32: mpz_pow_ui (mpz_ptr r, mpz_srcptr b, unsigned long int e)
                     33: #else
                     34: mpz_pow_ui (r, b, e)
                     35:      mpz_ptr r;
                     36:      mpz_srcptr b;
                     37:      unsigned long int e;
                     38: #endif
                     39: #else /* BERKELEY_MP */
                     40: void
                     41: #if __STDC__
                     42: rpow (const MINT *b, signed short int e, MINT *r)
                     43: #else
                     44: rpow (b, e, r)
                     45:      const MINT *b;
                     46:      signed short int e;
                     47:      MINT *r;
                     48: #endif
                     49: #endif /* BERKELEY_MP */
                     50: {
                     51:   mp_ptr rp, bp, tp, xp;
                     52:   mp_size_t rsize, bsize;
                     53:   int cnt, i;
                     54:   mp_limb_t blimb;
                     55:   TMP_DECL (marker);
                     56:
                     57:   bsize = ABS (b->_mp_size);
                     58:
                     59:   /* Single out cases that give result == 0 or 1.  These tests are here
                     60:      to simplify the general code below, not to optimize.  */
                     61:   if (e == 0)
                     62:     {
                     63:       r->_mp_d[0] = 1;
                     64:       r->_mp_size = 1;
                     65:       return;
                     66:     }
                     67:   if (bsize == 0
                     68: #ifdef BERKELEY_MP
                     69:       || e < 0
                     70: #endif
                     71:       )
                     72:     {
                     73:       r->_mp_size = 0;
                     74:       return;
                     75:     }
                     76:
                     77:   bp = b->_mp_d;
                     78:
                     79:   blimb = bp[bsize - 1];
                     80:   if (bsize == 1 && blimb < 0x100)
                     81:     {
                     82:       /* Estimate space requirements accurately.  Using the code from the
                     83:         `else' path would over-estimate space requirements wildly.   */
                     84:       float lb = __mp_bases[blimb].chars_per_bit_exactly;
                     85:       rsize = 2 + ((mp_size_t) (e / lb) / BITS_PER_MP_LIMB);
                     86:     }
                     87:   else
                     88:     {
                     89:       /* Over-estimate space requirements somewhat.  */
                     90:       count_leading_zeros (cnt, blimb);
                     91:       rsize = bsize * e - cnt * e / BITS_PER_MP_LIMB + 1;
                     92:     }
                     93:
                     94:   TMP_MARK (marker);
                     95:
                     96:   /* The two areas are used to alternatingly hold the input and recieve the
                     97:      product for mpn_mul.  (This scheme is used to fulfill the requirements
                     98:      of mpn_mul; that the product space may not be the same as any of the
                     99:      input operands.)  */
                    100:   rp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
                    101:   tp = (mp_ptr) TMP_ALLOC (rsize * BYTES_PER_MP_LIMB);
                    102:
                    103:   MPN_COPY (rp, bp, bsize);
                    104:   rsize = bsize;
                    105:   count_leading_zeros (cnt, e);
                    106:
                    107:   for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--)
                    108:     {
                    109:       mpn_mul_n (tp, rp, rp, rsize);
                    110:       rsize = 2 * rsize;
                    111:       rsize -= tp[rsize - 1] == 0;
                    112:       xp = tp; tp = rp; rp = xp;
                    113:
                    114:       if ((e & ((mp_limb_t) 1 << i)) != 0)
                    115:        {
                    116:          rsize = rsize + bsize - (mpn_mul (tp, rp, rsize, bp, bsize) == 0);
                    117:          xp = tp; tp = rp; rp = xp;
                    118:        }
                    119:     }
                    120:
                    121:   /* Now then we know the exact space requirements, reallocate if
                    122:      necessary.  */
                    123:   if (r->_mp_alloc < rsize)
                    124:     _mpz_realloc (r, rsize);
                    125:
                    126:   MPN_COPY (r->_mp_d, rp, rsize);
                    127:   r->_mp_size = (e & 1) == 0 || b->_mp_size >= 0 ? rsize : -rsize;
                    128:   TMP_FREE (marker);
                    129: }