1/*
2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001-2016 Free Software Foundation, Inc.
5 *
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.
10 *
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.
15 *
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/>.
18 */
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/*************************************************************************/
33
34#include "mpa.h"
35#include <math_private.h>
36
37#include <stap-probe.h>
38
39#ifndef SECTION
40# define SECTION
41#endif
42
43void __mpexp (mp_no *x, mp_no *y, int p);
44void __mplog (mp_no *x, mp_no *y, int p);
45double ulog (double);
46double __halfulp (double x, double y);
47
48double
49SECTION
50__slowpow (double x, double y, double z)
51{
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;
56
57 /* __HALFULP returns -10 or X^Y. */
58 res = __halfulp (x, y);
59
60 /* Return if the result was computed by __HALFULP. */
61 if (res >= 0)
62 return res;
63
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;
69
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);
75
76 /* Return the result if it is accurate enough. */
77 if (res == res1)
78 return res;
79#endif
80
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);
87
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);
94
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)
102 {
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;
107 }
108
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);
119
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);
123
124 return res;
125}
126