]>
Commit | Line | Data |
---|---|---|
04067002 UD |
1 | /* |
2 | * IBM Accurate Mathematical Library | |
3 | * written by International Business Machines Corp. | |
f4cf5f2d | 4 | * Copyright (C) 2001-2012 Free Software Foundation, Inc. |
04067002 UD |
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 | |
59ba27a6 | 17 | * along with this program; if not, see <http://www.gnu.org/licenses/>. |
04067002 UD |
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 | /* recompute. */ | |
32 | /*************************************************************************/ | |
33 | ||
34 | #include "mpa.h" | |
1ed0291c | 35 | #include <math_private.h> |
04067002 UD |
36 | |
37 | void __mpexp (mp_no * x, mp_no * y, int p); | |
38 | void __mplog (mp_no * x, mp_no * y, int p); | |
39 | double ulog (double); | |
40 | double __halfulp (double x, double y); | |
41 | ||
42 | double | |
43 | __slowpow (double x, double y, double z) | |
44 | { | |
45 | double res, res1; | |
46 | long double ldw, ldz, ldpp; | |
47 | static const long double ldeps = 0x4.0p-96; | |
48 | ||
49 | res = __halfulp (x, y); /* halfulp() returns -10 or x^y */ | |
50 | if (res >= 0) | |
51 | return res; /* if result was really computed by halfulp */ | |
52 | /* else, if result was not really computed by halfulp */ | |
53 | ||
54 | /* Compute pow as long double, 106 bits */ | |
55 | ldz = __ieee754_logl ((long double) x); | |
56 | ldw = (long double) y *ldz; | |
57 | ldpp = __ieee754_expl (ldw); | |
58 | res = (double) (ldpp + ldeps); | |
59 | res1 = (double) (ldpp - ldeps); | |
60 | ||
61 | if (res != res1) /* if result still not accurate enough */ | |
62 | { /* use mpa for higher persision. */ | |
63 | mp_no mpx, mpy, mpz, mpw, mpp, mpr, mpr1; | |
64 | static const mp_no eps = { -3, {1.0, 4.0} }; | |
65 | int p; | |
66 | ||
67 | p = 10; /* p=precision 240 bits */ | |
68 | __dbl_mp (x, &mpx, p); | |
69 | __dbl_mp (y, &mpy, p); | |
70 | __dbl_mp (z, &mpz, p); | |
71 | __mplog (&mpx, &mpz, p); /* log(x) = z */ | |
72 | __mul (&mpy, &mpz, &mpw, p); /* y * z =w */ | |
73 | __mpexp (&mpw, &mpp, p); /* e^w =pp */ | |
74 | __add (&mpp, &eps, &mpr, p); /* pp+eps =r */ | |
75 | __mp_dbl (&mpr, &res, p); | |
76 | __sub (&mpp, &eps, &mpr1, p); /* pp -eps =r1 */ | |
77 | __mp_dbl (&mpr1, &res1, p); /* converting into double precision */ | |
78 | if (res == res1) | |
79 | return res; | |
80 | ||
81 | /* if we get here result wasn't calculated exactly, continue for | |
82 | more exact calculation using 768 bits. */ | |
83 | p = 32; | |
84 | __dbl_mp (x, &mpx, p); | |
85 | __dbl_mp (y, &mpy, p); | |
86 | __dbl_mp (z, &mpz, p); | |
87 | __mplog (&mpx, &mpz, p); /* log(c)=z */ | |
88 | __mul (&mpy, &mpz, &mpw, p); /* y*z =w */ | |
89 | __mpexp (&mpw, &mpp, p); /* e^w=pp */ | |
90 | __mp_dbl (&mpp, &res, p); /* converting into double precision */ | |
91 | } | |
92 | return res; | |
93 | } |