1 | /* |
2 | * IBM Accurate Mathematical Library |
3 | * written by International Business Machines Corp. |
4 | * Copyright (C) 2001-2018 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 | |
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); |
47 | |
48 | double |
49 | SECTION |
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 | |