2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2015 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 */
22 /* FUNCTION:slowpow */
24 /*FILES NEEDED:mpa.h */
25 /* mpa.c mpexp.c mplog.c halfulp.c */
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 /*************************************************************************/
35 #include <math_private.h>
37 #include <stap-probe.h>
43 void __mpexp (mp_no
*x
, mp_no
*y
, int p
);
44 void __mplog (mp_no
*x
, mp_no
*y
, int p
);
46 double __halfulp (double x
, double y
);
50 __slowpow (double x
, double y
, double z
)
53 mp_no mpx
, mpy
, mpz
, mpw
, mpp
, mpr
, mpr1
;
54 static const mp_no eps
= {-3, {1.0, 4.0}};
57 /* __HALFULP returns -10 or X^Y. */
58 res
= __halfulp (x
, y
);
60 /* Return if the result was computed by __HALFULP. */
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. */
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. */
84 __dbl_mp (x
, &mpx
, p
);
85 __dbl_mp (y
, &mpy
, p
);
86 __dbl_mp (z
, &mpz
, p
);
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
);
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
);
109 /* If we don't, then we repeat using a higher precision. 768 bits of
110 precision ought to be enough for anybody. */
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
);