Annotation of OpenXM_contrib/gmp/mpfr/pow.c, Revision 1.1.1.2
1.1.1.2 ! ohara 1: /* mpfr_pow -- power function x^y
1.1 maekawa 2:
1.1.1.2 ! ohara 3: Copyright 2001, 2002 Free Software Foundation, Inc.
1.1 maekawa 4:
5: This file is part of the MPFR Library.
6:
7: The MPFR Library is free software; you can redistribute it and/or modify
1.1.1.2 ! ohara 8: it under the terms of the GNU Lesser General Public License as published by
! 9: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1 maekawa 10: option) any later version.
11:
12: The MPFR Library is distributed in the hope that it will be useful, but
13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.1.1.2 ! ohara 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 15: License for more details.
16:
1.1.1.2 ! ohara 17: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 18: along with the MPFR 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: #include "gmp.h"
1.1.1.2 ! ohara 23: #include "gmp-impl.h"
! 24: #include "longlong.h"
1.1 maekawa 25: #include "mpfr.h"
1.1.1.2 ! ohara 26: #include "mpfr-impl.h"
1.1 maekawa 27:
1.1.1.2 ! ohara 28: static int mpfr_pow_is_exact _PROTO((mpfr_srcptr, mpfr_srcptr));
! 29:
! 30: /* return non zero iff x^y is exact.
! 31: Assumes x and y are ordinary numbers (neither NaN nor Inf),
! 32: and y is not zero.
! 33: */
! 34: int
! 35: mpfr_pow_is_exact (mpfr_srcptr x, mpfr_srcptr y)
1.1 maekawa 36: {
1.1.1.2 ! ohara 37: mp_exp_t d;
! 38: unsigned long i, c;
! 39: mp_limb_t *yp;
1.1 maekawa 40:
1.1.1.2 ! ohara 41: if ((mpfr_sgn (x) < 0) && (mpfr_isinteger (y) == 0))
! 42: return 0;
! 43:
! 44: if (mpfr_sgn (y) < 0)
! 45: return mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_EXP(x) - 1) == 0;
! 46:
! 47: /* compute d such that y = c*2^d with c odd integer */
! 48: d = MPFR_EXP(y) - MPFR_PREC(y);
! 49: /* since y is not zero, necessarily one of the mantissa limbs is not zero,
! 50: thus we can simply loop until we find a non zero limb */
! 51: yp = MPFR_MANT(y);
! 52: for (i = 0; yp[i] == 0; i++, d += BITS_PER_MP_LIMB);
! 53: /* now yp[i] is not zero */
! 54: count_trailing_zeros (c, yp[i]);
! 55: d += c;
! 56:
! 57: if (d < 0)
! 58: {
! 59: mpz_t a;
! 60: mp_exp_t b;
! 61:
! 62: mpz_init (a);
! 63: b = mpfr_get_z_exp (a, x); /* x = a * 2^b */
! 64: c = mpz_scan1 (a, 0);
! 65: mpz_div_2exp (a, a, c);
! 66: b += c;
! 67: /* now a is odd */
! 68: while (d != 0)
! 69: {
! 70: if (mpz_perfect_square_p (a))
! 71: {
! 72: d++;
! 73: mpz_sqrt (a, a);
! 74: }
! 75: else
! 76: {
! 77: mpz_clear (a);
! 78: return 0;
! 79: }
! 80: }
! 81: mpz_clear (a);
! 82: }
! 83:
! 84: return 1;
1.1 maekawa 85: }
86:
1.1.1.2 ! ohara 87: /* The computation of z = pow(x,y) is done by
! 88: z = exp(y * log(x)) = x^y */
! 89: int
! 90: mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
1.1 maekawa 91: {
1.1.1.2 ! ohara 92: int inexact = 0;
! 93:
! 94: if (MPFR_IS_NAN(x) || MPFR_IS_NAN(y))
! 95: {
! 96: MPFR_SET_NAN(z);
! 97: MPFR_RET_NAN;
! 98: }
! 99:
! 100: if (MPFR_IS_INF(y))
! 101: {
! 102: mpfr_t one;
! 103: int cmp;
! 104:
! 105: if (MPFR_SIGN(y) > 0)
! 106: {
! 107: if (MPFR_IS_INF(x))
! 108: {
! 109: MPFR_CLEAR_FLAGS(z);
! 110: if (MPFR_SIGN(x) > 0)
! 111: MPFR_SET_INF(z);
! 112: else
! 113: MPFR_SET_ZERO(z);
! 114: MPFR_SET_POS(z);
! 115: MPFR_RET(0);
! 116: }
! 117: MPFR_CLEAR_FLAGS(z);
! 118: if (MPFR_IS_ZERO(x))
! 119: {
! 120: MPFR_SET_ZERO(z);
! 121: MPFR_SET_POS(z);
! 122: MPFR_RET(0);
! 123: }
! 124: mpfr_init2(one, BITS_PER_MP_LIMB);
! 125: mpfr_set_ui(one, 1, GMP_RNDN);
! 126: cmp = mpfr_cmp_abs(x, one);
! 127: mpfr_clear(one);
! 128: if (cmp > 0)
! 129: {
! 130: MPFR_SET_INF(z);
! 131: MPFR_SET_POS(z);
! 132: MPFR_RET(0);
! 133: }
! 134: else if (cmp < 0)
! 135: {
! 136: MPFR_SET_ZERO(z);
! 137: MPFR_SET_POS(z);
! 138: MPFR_RET(0);
! 139: }
! 140: else
! 141: {
! 142: MPFR_SET_NAN(z);
! 143: MPFR_RET_NAN;
! 144: }
! 145: }
! 146: else
! 147: {
! 148: if (MPFR_IS_INF(x))
! 149: {
! 150: MPFR_CLEAR_FLAGS(z);
! 151: if (MPFR_SIGN(x) > 0)
! 152: MPFR_SET_ZERO(z);
! 153: else
! 154: MPFR_SET_INF(z);
! 155: MPFR_SET_POS(z);
! 156: MPFR_RET(0);
! 157: }
! 158: if (MPFR_IS_ZERO(x))
! 159: {
! 160: MPFR_SET_INF(z);
! 161: MPFR_SET_POS(z);
! 162: MPFR_RET(0);
! 163: }
! 164: mpfr_init2(one, BITS_PER_MP_LIMB);
! 165: mpfr_set_ui(one, 1, GMP_RNDN);
! 166: cmp = mpfr_cmp_abs(x, one);
! 167: mpfr_clear(one);
! 168: MPFR_CLEAR_FLAGS(z);
! 169: if (cmp > 0)
! 170: {
! 171: MPFR_SET_ZERO(z);
! 172: MPFR_SET_POS(z);
! 173: MPFR_RET(0);
! 174: }
! 175: else if (cmp < 0)
! 176: {
! 177: MPFR_SET_INF(z);
! 178: MPFR_SET_POS(z);
! 179: MPFR_RET(0);
! 180: }
! 181: else
! 182: {
! 183: MPFR_SET_NAN(z);
! 184: MPFR_RET_NAN;
! 185: }
! 186: }
! 187: }
! 188:
! 189: if (MPFR_IS_ZERO(y))
! 190: {
! 191: return mpfr_set_ui (z, 1, GMP_RNDN);
! 192: }
! 193:
! 194: if (mpfr_isinteger (y))
! 195: {
! 196: mpz_t zi;
! 197: long int zii;
! 198: int exptol;
! 199:
! 200: mpz_init(zi);
! 201: exptol = mpfr_get_z_exp (zi, y);
! 202:
! 203: if (exptol>0)
! 204: mpz_mul_2exp(zi, zi, exptol);
! 205: else
! 206: mpz_tdiv_q_2exp(zi, zi, (unsigned long int) (-exptol));
! 207:
! 208: zii=mpz_get_ui(zi);
! 209:
! 210: mpz_clear(zi);
! 211: return mpfr_pow_si (z, x, zii, rnd_mode);
! 212: }
! 213:
! 214: if (MPFR_IS_INF(x))
! 215: {
! 216: if (MPFR_SIGN(x) > 0)
! 217: {
! 218: MPFR_CLEAR_FLAGS(z);
! 219: if (MPFR_SIGN(y) > 0)
! 220: MPFR_SET_INF(z);
! 221: else
! 222: MPFR_SET_ZERO(z);
! 223: MPFR_SET_POS(z);
! 224: MPFR_RET(0);
! 225: }
! 226: else
! 227: {
! 228: MPFR_SET_NAN(z);
! 229: MPFR_RET_NAN;
! 230: }
! 231: }
! 232:
! 233: if (MPFR_IS_ZERO(x))
! 234: {
! 235: MPFR_CLEAR_FLAGS(z);
! 236: MPFR_SET_ZERO(z);
! 237: MPFR_SET_SAME_SIGN(z, x);
! 238: MPFR_RET(0);
! 239: }
! 240:
! 241: if (MPFR_SIGN(x) < 0)
! 242: {
! 243: MPFR_SET_NAN(z);
! 244: MPFR_RET_NAN;
! 245: }
! 246:
! 247: MPFR_CLEAR_FLAGS(z);
! 248:
! 249: /* General case */
! 250: {
! 251: /* Declaration of the intermediary variable */
! 252: mpfr_t t, te, ti;
! 253: int loop = 0, ok;
! 254:
! 255: /* Declaration of the size variable */
! 256: mp_prec_t Nx = MPFR_PREC(x); /* Precision of input variable */
! 257: mp_prec_t Ny = MPFR_PREC(y); /* Precision of input variable */
! 258:
! 259: mp_prec_t Nt; /* Precision of the intermediary variable */
! 260: long int err; /* Precision of error */
! 261:
! 262: /* compute the precision of intermediary variable */
! 263: Nt=MAX(Nx,Ny);
! 264: /* the optimal number of bits : see algorithms.ps */
! 265: Nt=Nt+5+_mpfr_ceil_log2(Nt);
1.1 maekawa 266:
1.1.1.2 ! ohara 267: /* initialise of intermediary variable */
! 268: mpfr_init(t);
! 269: mpfr_init(ti);
! 270: mpfr_init(te);
! 271:
! 272: do
! 273: {
! 274:
! 275: loop ++;
! 276:
! 277: /* reactualisation of the precision */
! 278: mpfr_set_prec (t, Nt);
! 279: mpfr_set_prec (ti, Nt);
! 280: mpfr_set_prec (te, Nt);
! 281:
! 282: /* compute exp(y*ln(x)) */
! 283: mpfr_log (ti, x, GMP_RNDU); /* ln(n) */
! 284: mpfr_mul (te, y, ti, GMP_RNDU); /* y*ln(n) */
! 285: mpfr_exp (t, te, GMP_RNDN); /* exp(x*ln(n))*/
! 286:
! 287: /* estimation of the error -- see pow function in algorithms.ps*/
! 288: err = Nt - (MPFR_EXP(te) + 3);
! 289:
! 290: /* actualisation of the precision */
! 291: Nt += 10;
! 292:
! 293: ok = mpfr_can_round (t, err, GMP_RNDN, rnd_mode, Ny);
! 294:
! 295: /* check exact power */
! 296: if (ok == 0 && loop == 1)
! 297: ok = mpfr_pow_is_exact (x, y);
! 298:
! 299: }
! 300: while (err < 0 || ok == 0);
! 301:
! 302: inexact = mpfr_set (z, t, rnd_mode);
! 303:
! 304: mpfr_clear (t);
! 305: mpfr_clear (ti);
! 306: mpfr_clear (te);
! 307: }
! 308: return inexact;
1.1 maekawa 309: }
1.1.1.2 ! ohara 310:
! 311:
! 312:
! 313:
! 314:
! 315:
! 316:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to pow.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr