2 * IBM Accurate Mathematical Library
3 * written by International Business Machines Corp.
4 * Copyright (C) 2001 Free Software Foundation
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 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, write to the Free Software
18 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
20 /***************************************************************************/
21 /* MODULE_NAME:uexp.c */
26 /* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h */
27 /* mpa.c mpexp.x slowexp.c */
29 /* An ultimate exp routine. Given an IEEE double machine number x */
30 /* it computes the correctly rounded (to nearest) value of e^x */
31 /* Assumption: Machine arithmetic operations are performed in */
32 /* round to nearest mode of IEEE 754 standard. */
34 /***************************************************************************/
41 #include "math_private.h"
43 double __slowexp(double);
45 /***************************************************************************/
46 /* An ultimate exp routine. Given an IEEE double machine number x */
47 /* it computes the correctly rounded (to nearest) value of e^x */
48 /***************************************************************************/
49 double __ieee754_exp(double x) {
50 double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
51 mynumber junk1, junk2, binexp = {{0,0}};
58 m = junk1.i[HIGH_HALF];
61 if (n > smallint && n < bigint) {
63 y = x*log2e.x + three51.x;
64 bexp = y - three51.x; /* multiply the result by 2**bexp */
68 eps = bexp*ln_two2.x; /* x = bexp*ln(2) + t - eps */
69 t = x - bexp*ln_two1.x;
72 base = y - three33.x; /* t rounded to a multiple of 2**-18 */
74 del = (t - base) - eps; /* x = bexp*ln(2) + base + del */
75 eps = del + del*del*(p3.x*del + p2.x);
77 binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;
79 i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
80 j = (junk2.i[LOW_HALF]&511)<<1;
82 al = coar.x[i]*fine.x[j];
83 bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
85 rem=(bet + bet*eps)+al*eps;
87 cor = (al - res) + rem;
88 if (res == (res+cor*err_0)) return res*binexp.x;
89 else return __slowexp(x); /*if error is over bound */
92 if (n <= smallint) return 1.0;
95 if (n > infint) return(zero/zero); /* x is NaN, return invalid */
96 if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
97 /* x is finite, cause either overflow or underflow */
98 if (junk1.i[LOW_HALF] != 0) return (zero/zero); /* x is NaN */
99 return ((x>0)?inf.x:zero ); /* |x| = inf; return either inf or 0 */
102 y = x*log2e.x + three51.x;
103 bexp = y - three51.x;
105 eps = bexp*ln_two2.x;
106 t = x - bexp*ln_two1.x;
108 base = y - three33.x;
110 del = (t - base) - eps;
111 eps = del + del*del*(p3.x*del + p2.x);
112 i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
113 j = (junk2.i[LOW_HALF]&511)<<1;
114 al = coar.x[i]*fine.x[j];
115 bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
116 rem=(bet + bet*eps)+al*eps;
118 cor = (al - res) + rem;
120 ex=junk1.i[LOW_HALF];
121 if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
123 binexp.i[HIGH_HALF] = (1023+ex)<<20;
124 if (res == (res+cor*err_0)) return res*binexp.x;
125 else return __slowexp(x); /*if error is over bound */
128 binexp.i[HIGH_HALF] = (1023-ex)<<20;
131 eps=1.0000000001+err_0*binexp.x;
133 y = ((1.0-t)+res)+cor;
136 if (res == (res + eps*cor))
137 { binexp.i[HIGH_HALF] = 0x00100000;
138 return (res-1.0)*binexp.x;
140 else return __slowexp(x); /* if error is over bound */
143 binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
144 if (res == (res+cor*err_0)) return res*binexp.x*t256.x;
145 else return __slowexp(x);
149 /************************************************************************/
150 /* Compute e^(x+xx)(Double-Length number) .The routine also receive */
151 /* bound of error of previous calculation .If after computing exp */
152 /* error bigger than allows routine return non positive number */
153 /*else return e^(x + xx) (always positive ) */
154 /************************************************************************/
156 double __exp1(double x, double xx, double error) {
157 double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
158 mynumber junk1, junk2, binexp = {{0,0}};
165 m = junk1.i[HIGH_HALF];
166 n = m&hugeint; /* no sign */
168 if (n > smallint && n < bigint) {
169 y = x*log2e.x + three51.x;
170 bexp = y - three51.x; /* multiply the result by 2**bexp */
174 eps = bexp*ln_two2.x; /* x = bexp*ln(2) + t - eps */
175 t = x - bexp*ln_two1.x;
178 base = y - three33.x; /* t rounded to a multiple of 2**-18 */
180 del = (t - base) + (xx-eps); /* x = bexp*ln(2) + base + del */
181 eps = del + del*del*(p3.x*del + p2.x);
183 binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;
185 i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
186 j = (junk2.i[LOW_HALF]&511)<<1;
188 al = coar.x[i]*fine.x[j];
189 bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
191 rem=(bet + bet*eps)+al*eps;
193 cor = (al - res) + rem;
194 if (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
198 if (n <= smallint) return 1.0; /* if x->0 e^x=1 */
201 if (n > infint) return(zero/zero); /* x is NaN, return invalid */
202 if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
203 /* x is finite, cause either overflow or underflow */
204 if (junk1.i[LOW_HALF] != 0) return (zero/zero); /* x is NaN */
205 return ((x>0)?inf.x:zero ); /* |x| = inf; return either inf or 0 */
208 y = x*log2e.x + three51.x;
209 bexp = y - three51.x;
211 eps = bexp*ln_two2.x;
212 t = x - bexp*ln_two1.x;
214 base = y - three33.x;
216 del = (t - base) + (xx-eps);
217 eps = del + del*del*(p3.x*del + p2.x);
218 i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
219 j = (junk2.i[LOW_HALF]&511)<<1;
220 al = coar.x[i]*fine.x[j];
221 bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
222 rem=(bet + bet*eps)+al*eps;
224 cor = (al - res) + rem;
226 ex=junk1.i[LOW_HALF];
227 if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
229 binexp.i[HIGH_HALF] = (1023+ex)<<20;
230 if (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
234 binexp.i[HIGH_HALF] = (1023-ex)<<20;
237 eps=1.00000000001+(error+err_1)*binexp.x;
239 y = ((1.0-t)+res)+cor;
242 if (res == (res + eps*cor))
243 {binexp.i[HIGH_HALF] = 0x00100000; return (res-1.0)*binexp.x;}
247 binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
248 if (res == (res+cor*(1.0+error+err_1)))
249 return res*binexp.x*t256.x;