Update copyright notices with scripts/update-copyrights
[glibc.git] / sysdeps / ieee754 / dbl-64 / slowpow.c
blobb2e403308e62f65bea3e3fdb15cf7ee294859041
1 /*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2014 Free Software Foundation, Inc.
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; either version 2.1 of the License, or
9 * (at your option) any later version.
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, see <http://www.gnu.org/licenses/>.
19 /*************************************************************************/
20 /* MODULE_NAME:slowpow.c */
21 /* */
22 /* FUNCTION:slowpow */
23 /* */
24 /*FILES NEEDED:mpa.h */
25 /* mpa.c mpexp.c mplog.c halfulp.c */
26 /* */
27 /* Given two IEEE double machine numbers y,x , routine computes the */
28 /* correctly rounded (to nearest) value of x^y. Result calculated by */
29 /* multiplication (in halfulp.c) or if result isn't accurate enough */
30 /* then routine converts x and y into multi-precision doubles and */
31 /* calls to mpexp routine */
32 /*************************************************************************/
34 #include "mpa.h"
35 #include <math_private.h>
37 #include <stap-probe.h>
39 #ifndef SECTION
40 # define SECTION
41 #endif
43 void __mpexp (mp_no *x, mp_no *y, int p);
44 void __mplog (mp_no *x, mp_no *y, int p);
45 double ulog (double);
46 double __halfulp (double x, double y);
48 double
49 SECTION
50 __slowpow (double x, double y, double z)
52 double res, res1;
53 mp_no mpx, mpy, mpz, mpw, mpp, mpr, mpr1;
54 static const mp_no eps = {-3, {1.0, 4.0}};
55 int p;
57 /* __HALFULP returns -10 or X^Y. */
58 res = __halfulp (x, y);
60 /* Return if the result was computed by __HALFULP. */
61 if (res >= 0)
62 return res;
64 /* Compute pow as long double. This is currently only used by powerpc, where
65 one may get 106 bits of accuracy. */
66 #ifdef USE_LONG_DOUBLE_FOR_MP
67 long double ldw, ldz, ldpp;
68 static const long double ldeps = 0x4.0p-96;
70 ldz = __ieee754_logl ((long double) x);
71 ldw = (long double) y *ldz;
72 ldpp = __ieee754_expl (ldw);
73 res = (double) (ldpp + ldeps);
74 res1 = (double) (ldpp - ldeps);
76 /* Return the result if it is accurate enough. */
77 if (res == res1)
78 return res;
79 #endif
81 /* Or else, calculate using multiple precision. P = 10 implies accuracy of
82 240 bits accuracy, since MP_NO has a radix of 2^24. */
83 p = 10;
84 __dbl_mp (x, &mpx, p);
85 __dbl_mp (y, &mpy, p);
86 __dbl_mp (z, &mpz, p);
88 /* z = x ^ y
89 log (z) = y * log (x)
90 z = exp (y * log (x)) */
91 __mplog (&mpx, &mpz, p);
92 __mul (&mpy, &mpz, &mpw, p);
93 __mpexp (&mpw, &mpp, p);
95 /* Add and subtract EPS to ensure that the result remains unchanged, i.e. we
96 have last bit accuracy. */
97 __add (&mpp, &eps, &mpr, p);
98 __mp_dbl (&mpr, &res, p);
99 __sub (&mpp, &eps, &mpr1, p);
100 __mp_dbl (&mpr1, &res1, p);
101 if (res == res1)
103 /* Track how often we get to the slow pow code plus
104 its input/output values. */
105 LIBC_PROBE (slowpow_p10, 4, &x, &y, &z, &res);
106 return res;
109 /* If we don't, then we repeat using a higher precision. 768 bits of
110 precision ought to be enough for anybody. */
111 p = 32;
112 __dbl_mp (x, &mpx, p);
113 __dbl_mp (y, &mpy, p);
114 __dbl_mp (z, &mpz, p);
115 __mplog (&mpx, &mpz, p);
116 __mul (&mpy, &mpz, &mpw, p);
117 __mpexp (&mpw, &mpp, p);
118 __mp_dbl (&mpp, &res, p);
120 /* Track how often we get to the uber-slow pow code plus
121 its input/output values. */
122 LIBC_PROBE (slowpow_p32, 4, &x, &y, &z, &res);
124 return res;