rts/gmp/mpz/ui_pow_ui.c

   1 /* mpz_ui_pow_ui(res, base, exp) -- Set RES to BASE**EXP.
   2
   3 Copyright (C) 1991, 1993, 1994, 1996, 1997, 2000 Free Software Foundation,
   4 Inc.
   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
   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
  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
  15 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16 License for more details.
  17
  18 You should have received a copy of the GNU Lesser General Public License
  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
  27
  28 static void mpz_pow2 _PROTO ((mpz_ptr r, mp_limb_t blimb, unsigned long int e, mp_limb_t rl));
  29
  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;
  41   mp_limb_t rl;
  42
  43   if (e == 0)
  44     {
  45       /* For x^0 we return 1, even if x is 0.  */
  46       r->_mp_d[0] = 1;
  47       r->_mp_size = 1;
  48       return;
  49     }
  50
  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))
  56     {
  57       if ((e & 1) != 0)
  58         rl = rl * blimb;
  59       blimb = blimb * blimb;
  60       e = e >> 1;
  61     }
  62
  63   /* rl is now b^(e mod n).  (I.e., the latter factor above.)  */
  64
  65   if (e == 0)
  66     {
  67       r->_mp_d[0] = rl;
  68       r->_mp_size = rl != 0;
  69       return;
  70     }
  71
  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
  92   TMP_MARK (marker);
  93
  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);
 103
 104   rp[0] = blimb;
 105   rsize = 1;
 106
 107   count_leading_zeros (cnt, e);
 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;
 113       MP_PTR_SWAP (rp, tp);
 114
 115       if ((e & ((mp_limb_t) 1 << i)) != 0)
 116         {
 117           mp_limb_t cy;
 118           cy = mpn_mul_1 (rp, rp, rsize, blimb);
 119           rp[rsize] = cy;
 120           rsize += cy != 0;
 121         }
 122     }
 123
 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   }
 136
 137   r->_mp_size = rsize;
 138   TMP_FREE (marker);
 139 }