Return to pow_ui.c CVS log

Up to [local] / OpenXM_contrib / gmp / mpfr

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

1.1       ohara       1: /* mpfr_pow_ui-- compute the power of a floating-point
                      2:                                   by a machine integer
                      3:
                      4: Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
                      5:
                      6: This file is part of the MPFR Library.
                      7:
                      8: The MPFR 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 MPFR 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 MPFR 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 "mpfr.h"
                     25: #include "mpfr-impl.h"
                     26:
                     27: /* sets x to y^n, and return 0 if exact, non-zero otherwise */
                     28: int
                     29: mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd)
                     30: {
                     31:   long int i, err;
                     32:   unsigned long m;
                     33:   mpfr_t res;
                     34:   mp_prec_t prec;
                     35:   int inexact;
                     36:   mp_rnd_t rnd1;
                     37:
                     38:   if (MPFR_IS_NAN(y))
                     39:     {
                     40:       MPFR_SET_NAN(x);
                     41:       MPFR_RET_NAN;
                     42:     }
                     43:
                     44:   MPFR_CLEAR_NAN(x);
                     45:
                     46:   if (n == 0) /* x^0 = 1 for any x */
                     47:     {
                     48:       /* The return mpfr_set_ui is important as 1 isn't necessarily
                     49:          in the exponent range. */
                     50:       return mpfr_set_ui (x, 1, rnd);
                     51:     }
                     52:
                     53:   if (MPFR_IS_INF(y))
                     54:     {
                     55:       /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */
                     56:       if ((MPFR_SIGN(y) < 0) && (n % 2 == 1))
                     57:         MPFR_SET_NEG(x);
                     58:       else
                     59:         MPFR_SET_POS(x);
                     60:       MPFR_SET_INF(x);
                     61:       MPFR_RET(0);
                     62:     }
                     63:
                     64:   MPFR_CLEAR_INF(x);
                     65:
                     66:   mpfr_init (res);
                     67:
                     68:   prec = MPFR_PREC(x);
                     69:
                     70:   rnd1 = (MPFR_SIGN(y) > 0) ? GMP_RNDU : GMP_RNDD; /* away */
                     71:
                     72:   do
                     73:     {
                     74:       prec += 3;
                     75:       for (m=n, i=0; m; i++, m>>=1, prec++);
                     76:       mpfr_set_prec (res, prec);
                     77:       inexact = mpfr_set (res, y, rnd1);
                     78:       err = 1;
                     79:       /* now 2^(i-1) <= n < 2^i */
                     80:       for (i-=2; i>=0; i--)
                     81:        {
                     82:          if (mpfr_mul (res, res, res, GMP_RNDU))
                     83:            inexact = 1;
                     84:          err++;
                     85:          if (n & (1 << i))
                     86:            if (mpfr_mul (res, res, y, rnd1))
                     87:              inexact = 1;
                     88:        }
                     89:       err = prec - err;
                     90:       if (err < 0)
                     91:        err = 0;
                     92:     }
                     93:   while (inexact && (mpfr_can_round (res, err,
                     94:          (MPFR_SIGN(res) > 0) ? GMP_RNDU : GMP_RNDD, rnd, MPFR_PREC(x)) == 0));
                     95:
                     96:   if (mpfr_set (x, res, rnd))
                     97:     inexact = 1;
                     98:
                     99:   mpfr_clear (res);
                    100:
                    101:   return inexact;
                    102: }
                    103: