ghc/runtime/gmp/mpz_pow_ui.c

   1 /* mpz_pow_ui(res, base, exp) -- Set RES to BASE**EXP.
   2
   3 Copyright (C) 1991 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 General Public License as published by
   9 the Free Software Foundation; either version 2, or (at your option)
  10 any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful,
  13 but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 GNU General Public License for more details.
  16
  17 You should have received a copy of the GNU General Public License
  18 along with the GNU MP Library; see the file COPYING.  If not, write to
  19 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  20
  21 #include "gmp.h"
  22 #include "gmp-impl.h"
  23 #include "longlong.h"
  24
  25 #ifndef BERKELEY_MP
  26 void
  27 #ifdef __STDC__
  28 mpz_pow_ui (MP_INT *r, const MP_INT *b, unsigned long int e)
  29 #else
  30 mpz_pow_ui (r, b, e)
  31      MP_INT *r;
  32      const MP_INT *b;
  33      unsigned long int e;
  34 #endif
  35 #else /* BERKELEY_MP */
  36 void
  37 #ifdef __STDC__
  38 rpow (const MP_INT *b, signed short int e, MP_INT *r)
  39 #else
  40 rpow (b, e, r)
  41      const MP_INT *b;
  42      signed short int e;
  43      MP_INT *r;
  44 #endif
  45 #endif /* BERKELEY_MP */
  46 {
  47   mp_ptr rp, bp, tp, xp;
  48   mp_size rsize, bsize;
  49   int cnt, i;
  50
  51   bsize = ABS (b->size);
  52
  53   /* Single out cases that give result == 0 or 1.  These tests are here
  54      to simplify the general code below, not to optimize.  */
  55   if (bsize == 0
  56 #ifdef BERKELEY_MP
  57       || e < 0
  58 #endif
  59       )
  60     {
  61       r->size = 0;
  62       return;
  63     }
  64   if (e == 0)
  65     {
  66       r->d[0] = 1;
  67       r->size = 1;
  68       return;
  69     }
  70
  71   /* Count the number of leading zero bits of the base's most
  72      significant limb.  */
  73   count_leading_zeros (cnt, b->d[bsize - 1]);
  74
  75   /* Over-estimate space requirements and allocate enough space for the
  76      final result in two temporary areas.  The two areas are used to
  77      alternately hold the input and recieve the product for mpn_mul.
  78      (This scheme is used to fulfill the requirements of mpn_mul; that
  79      the product space may not be the same as any of the input operands.)  */
  80   rsize = bsize * e - cnt * e / BITS_PER_MP_LIMB;
  81
  82   rp = (mp_ptr) alloca (rsize * BYTES_PER_MP_LIMB);
  83   tp = (mp_ptr) alloca (rsize * BYTES_PER_MP_LIMB);
  84   bp = b->d;
  85
  86   MPN_COPY (rp, bp, bsize);
  87   rsize = bsize;
  88   count_leading_zeros (cnt, e);
  89
  90   for (i = BITS_PER_MP_LIMB - cnt - 2; i >= 0; i--)
  91     {
  92       rsize = mpn_mul (tp, rp, rsize, rp, rsize);
  93       xp = tp; tp = rp; rp = xp;
  94
  95       if ((e & ((mp_limb) 1 << i)) != 0)
  96         {
  97           rsize = mpn_mul (tp, rp, rsize, bp, bsize);
  98           xp = tp; tp = rp; rp = xp;
  99         }
 100     }
 101
 102   /* Now then we know the exact space requirements, reallocate if
 103      necessary.  */
 104   if (r->alloc < rsize)
 105     _mpz_realloc (r, rsize);
 106
 107   MPN_COPY (r->d, rp, rsize);
 108   r->size = (e & 1) == 0 || b->size >= 0 ? rsize : -rsize;
 109   alloca (0);
 110 }