mirror/glibc
2
0
Fork 0
mirror of git://sourceware.org/git/glibc.git synced 2024-11-21 01:12:26 +08:00
Code Issues Packages Projects Releases Wiki Activity
0b20008e04
glibc/sysdeps/ieee754/dbl-64/e_exp.c
Andreas Jaeger 3737a0d70c Do not raise execptions for exp(NaN).
2003-12-28 18:57:35 +00:00

253 lines
8.7 KiB
C
Raw Blame History

/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
* Copyright (C) 2001 Free Software Foundation
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation; either version 2.1 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
/***************************************************************************/
/* MODULE_NAME:uexp.c */
/* */
/* FUNCTION:uexp */
/* exp1 */
/* */
/* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h */
/* mpa.c mpexp.x slowexp.c */
/* */
/* An ultimate exp routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of e^x */
/* Assumption: Machine arithmetic operations are performed in */
/* round to nearest mode of IEEE 754 standard. */
/* */
/***************************************************************************/
#include "endian.h"
#include "uexp.h"
#include "mydefs.h"
#include "MathLib.h"
#include "uexp.tbl"
#include "math_private.h"
double __slowexp(double);
/***************************************************************************/
/* An ultimate exp routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of e^x */
/***************************************************************************/
double __ieee754_exp(double x) {
double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
mynumber junk1, junk2, binexp = {{0,0}};
#if 0
int4 k;
#endif
int4 i,j,m,n,ex;
junk1.x = x;
m = junk1.i[HIGH_HALF];
n = m&hugeint;
if (n > smallint && n < bigint) {
y = x*log2e.x + three51.x;
bexp = y - three51.x; /* multiply the result by 2**bexp */
junk1.x = y;
eps = bexp*ln_two2.x; /* x = bexp*ln(2) + t - eps */
t = x - bexp*ln_two1.x;
y = t + three33.x;
base = y - three33.x; /* t rounded to a multiple of 2**-18 */
junk2.x = y;
del = (t - base) - eps; /* x = bexp*ln(2) + base + del */
eps = del + del*del*(p3.x*del + p2.x);
binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;
i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
j = (junk2.i[LOW_HALF]&511)<<1;
al = coar.x[i]*fine.x[j];
bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
rem=(bet + bet*eps)+al*eps;
res = al + rem;
cor = (al - res) + rem;
if (res == (res+cor*err_0)) return res*binexp.x;
else return __slowexp(x); /*if error is over bound */
}
if (n <= smallint) return 1.0;
if (n >= badint) {
if (n > infint) return(x+x); /* x is NaN */
if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
/* x is finite, cause either overflow or underflow */
if (junk1.i[LOW_HALF] != 0) return (x+x); /* x is NaN */
return ((x>0)?inf.x:zero ); /* |x| = inf; return either inf or 0 */
}
y = x*log2e.x + three51.x;
bexp = y - three51.x;
junk1.x = y;
eps = bexp*ln_two2.x;
t = x - bexp*ln_two1.x;
y = t + three33.x;
base = y - three33.x;
junk2.x = y;
del = (t - base) - eps;
eps = del + del*del*(p3.x*del + p2.x);
i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
j = (junk2.i[LOW_HALF]&511)<<1;
al = coar.x[i]*fine.x[j];
bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
rem=(bet + bet*eps)+al*eps;
res = al + rem;
cor = (al - res) + rem;
if (m>>31) {
ex=junk1.i[LOW_HALF];
if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
if (ex >=-1022) {
binexp.i[HIGH_HALF] = (1023+ex)<<20;
if (res == (res+cor*err_0)) return res*binexp.x;
else return __slowexp(x); /*if error is over bound */
}
ex = -(1022+ex);
binexp.i[HIGH_HALF] = (1023-ex)<<20;
res*=binexp.x;
cor*=binexp.x;
eps=1.0000000001+err_0*binexp.x;
t=1.0+res;
y = ((1.0-t)+res)+cor;
res=t+y;
cor = (t-res)+y;
if (res == (res + eps*cor))
{ binexp.i[HIGH_HALF] = 0x00100000;
return (res-1.0)*binexp.x;
}
else return __slowexp(x); /* if error is over bound */
}
else {
binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
if (res == (res+cor*err_0)) return res*binexp.x*t256.x;
else return __slowexp(x);
}
}
/************************************************************************/
/* Compute e^(x+xx)(Double-Length number) .The routine also receive */
/* bound of error of previous calculation .If after computing exp */
/* error bigger than allows routine return non positive number */
/*else return e^(x + xx) (always positive ) */
/************************************************************************/
double __exp1(double x, double xx, double error) {
double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
mynumber junk1, junk2, binexp = {{0,0}};
#if 0
int4 k;
#endif
int4 i,j,m,n,ex;
junk1.x = x;
m = junk1.i[HIGH_HALF];
n = m&hugeint; /* no sign */
if (n > smallint && n < bigint) {
y = x*log2e.x + three51.x;
bexp = y - three51.x; /* multiply the result by 2**bexp */
junk1.x = y;
eps = bexp*ln_two2.x; /* x = bexp*ln(2) + t - eps */
t = x - bexp*ln_two1.x;
y = t + three33.x;
base = y - three33.x; /* t rounded to a multiple of 2**-18 */
junk2.x = y;
del = (t - base) + (xx-eps); /* x = bexp*ln(2) + base + del */
eps = del + del*del*(p3.x*del + p2.x);
binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+1023)<<20;
i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
j = (junk2.i[LOW_HALF]&511)<<1;
al = coar.x[i]*fine.x[j];
bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
rem=(bet + bet*eps)+al*eps;
res = al + rem;
cor = (al - res) + rem;
if (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
else return -10.0;
}
if (n <= smallint) return 1.0; /* if x->0 e^x=1 */
if (n >= badint) {
if (n > infint) return(zero/zero); /* x is NaN, return invalid */
if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
/* x is finite, cause either overflow or underflow */
if (junk1.i[LOW_HALF] != 0) return (zero/zero); /* x is NaN */
return ((x>0)?inf.x:zero ); /* |x| = inf; return either inf or 0 */
}
y = x*log2e.x + three51.x;
bexp = y - three51.x;
junk1.x = y;
eps = bexp*ln_two2.x;
t = x - bexp*ln_two1.x;
y = t + three33.x;
base = y - three33.x;
junk2.x = y;
del = (t - base) + (xx-eps);
eps = del + del*del*(p3.x*del + p2.x);
i = ((junk2.i[LOW_HALF]>>8)&0xfffffffe)+356;
j = (junk2.i[LOW_HALF]&511)<<1;
al = coar.x[i]*fine.x[j];
bet =(coar.x[i]*fine.x[j+1] + coar.x[i+1]*fine.x[j]) + coar.x[i+1]*fine.x[j+1];
rem=(bet + bet*eps)+al*eps;
res = al + rem;
cor = (al - res) + rem;
if (m>>31) {
ex=junk1.i[LOW_HALF];
if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
if (ex >=-1022) {
binexp.i[HIGH_HALF] = (1023+ex)<<20;
if (res == (res+cor*(1.0+error+err_1))) return res*binexp.x;
else return -10.0;
}
ex = -(1022+ex);
binexp.i[HIGH_HALF] = (1023-ex)<<20;
res*=binexp.x;
cor*=binexp.x;
eps=1.00000000001+(error+err_1)*binexp.x;
t=1.0+res;
y = ((1.0-t)+res)+cor;
res=t+y;
cor = (t-res)+y;
if (res == (res + eps*cor))
{binexp.i[HIGH_HALF] = 0x00100000; return (res-1.0)*binexp.x;}
else return -10.0;
}
else {
binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
if (res == (res+cor*(1.0+error+err_1)))
return res*binexp.x*t256.x;
else return -10.0;
}
}
Reference in New Issue View Git Blame Copy Permalink