2001-03-12 01:04:52 +01:00
|
|
|
/*
|
|
|
|
* IBM Accurate Mathematical Library
|
2002-07-06 08:36:39 +02:00
|
|
|
* written by International Business Machines Corp.
|
|
|
|
* Copyright (C) 2001 Free Software Foundation
|
2001-03-12 01:04:52 +01:00
|
|
|
*
|
|
|
|
* 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
|
2002-08-27 00:40:48 +02:00
|
|
|
* the Free Software Foundation; either version 2.1 of the License, or
|
2001-03-12 01:04:52 +01:00
|
|
|
* (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
|
2002-08-20 23:51:55 +02:00
|
|
|
* GNU Lesser General Public License for more details.
|
2001-03-12 01:04:52 +01:00
|
|
|
*
|
|
|
|
* 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"
|
2001-05-12 16:32:12 +02:00
|
|
|
#include "uexp.h"
|
2001-03-12 01:04:52 +01:00
|
|
|
#include "mydefs.h"
|
|
|
|
#include "MathLib.h"
|
|
|
|
#include "uexp.tbl"
|
2001-05-12 16:32:12 +02:00
|
|
|
#include "math_private.h"
|
|
|
|
|
2001-03-12 01:04:52 +01:00
|
|
|
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;
|
2001-03-12 08:57:09 +01:00
|
|
|
mynumber junk1, junk2, binexp = {{0,0}};
|
|
|
|
#if 0
|
|
|
|
int4 k;
|
|
|
|
#endif
|
|
|
|
int4 i,j,m,n,ex;
|
2001-03-12 01:04:52 +01:00
|
|
|
|
|
|
|
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) {
|
2003-12-28 19:57:35 +01:00
|
|
|
if (n > infint) return(x+x); /* x is NaN */
|
2001-03-12 01:04:52 +01:00
|
|
|
if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
|
|
|
|
/* x is finite, cause either overflow or underflow */
|
2003-12-28 19:57:35 +01:00
|
|
|
if (junk1.i[LOW_HALF] != 0) return (x+x); /* x is NaN */
|
2001-03-12 01:04:52 +01:00
|
|
|
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 */
|
Update.
1998-02-13 17:39 Ulrich Drepper <drepper@cygnus.com>
* elf/Makefile: Don't use --version-script parameter to link ld.so
unconditionally.
1998-01-02 04:19 Geoff Keating <geoffk@ozemail.com.au>
* math/Makefile: Add t_exp.
* math/libm-test.c: Tighten accuracy bounds for exp(), correct
constants.
* math/test-reduce.c: Remove temporarily, it seems to be broken.
* sysdeps/libm-ieee754/e_exp.c: Use accurate table method.
* sysdeps/libm-ieee754/e_expf.c: Use table & double precision for
better accuracy.
* sysdeps/libm-ieee754/s_exp2.c: Use better polynomial; correct
algorithm for very large/very small arguments.
* sysdeps/libm-ieee754/s_exp2f.c: Use slightly better polynomial;
correct algorithm for very large/very small arguments; adjust for
new table.
* sysdeps/libm-ieee754/t_exp.c: New file.
* sysdeps/libm-ieee754/t_exp2f.h: Use table with smaller deltas.
* sysdeps/unix/sysv/linux/powerpc/dl-sysdep.c: Put 'strange test'
back, with comment that explains what breaks when you remove it :-(.
* localedata/xfrm-test.c: Avoid integer overflow.
* stdlib/strfmon.c: char is unsigned, sometimes.
*sysdeps/powerpc
* sysdeps/powerpc/Makefile: Remove quad float support.
* sysdeps/powerpc/q_*.c: Remove, they will become an add-on.
* sysdeps/powerpc/quad_float.h: Likewise.
* sysdeps/powerpc/test-arith.c: Likewise.
* sysdeps/powerpc/test-arithf.c: Likewise.
* sysdeps/generic/s_exp2.c: Remove, we have this implemented now.
* sysdeps/generic/s_exp2f.c: Likewise.
* sysdeps/powerpc/bits/mathinline.h: Use underscores around __asm__,
don't try anything if _SOFT_FLOAT.
1997-12-31 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>
* locale/C-ctype.c (_nl_C_LC_CTYPE_class32): Undo last change.
* locale/programs/ld-ctype.c (CHAR_CLASS32_TRANS): Likewise.
* wctype/wctype.c: Likewise.
* wctype/wctype.h (_ISwxxx): Renamed from _ISxxx, all uses
changed. They are incompatible with the _ISxxx values from
<ctype.h> on little endian machines.
(_ISwbit) [__BYTE_ORDER == __LITTLE_ENDIAN]: Correctly transform
bit number. This fixes the real bug and restores the integrity of
the ctype locale file.
* wctype/wcfuncs.c: Change all _ISxxx to _ISwxxx.
* wctype/wcfuncs_l.c: Likewise.
* wctype/wcextra.c: Likewise.
* wctype/wctype_l.c [__BYTE_ORDER == __LITTLE_ENDIAN]: Use correct
byte swapping.
1998-02-09 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>
* sysdeps/unix/sysv/linux/sparc/sparc32/sysdep.S (errno): Put it
into .bss segment instead of .common, so that aliases on it work.
* sysdeps/unix/sysv/linux/i386/sysdep.S (errno): Add .type and
.size directives, put into .bss segment instead of initializing it
to 4.
1998-02-12 08:00 H.J. Lu <hjl@gnu.org>
* libc.map (gnu_get_libc_release, gnu_get_libc_version): Added.
* version.c (__gnu_get_libc_release, __gnu_get_libc_version): New
functions.
Make names without __ weak aliases.
(__libc_release, __libc_version): Make them static.
* include/gnu/libc-version.h: New file.
* Makefile (headers): Add gnu/libc-version.h.
1998-02-13 Ulrich Drepper <drepper@cygnus.com>
* stdlib/stdlib.h (struct drand48_data): Leave X to user macros
and use x for member name.
Reported by Daniel Lyddy <daniell@cs.berkeley.edu>.
* stdlib/drand48.c: Change according to member name change.
* stdlib/drand48_r.c: Likewise.
* stdlib/lcong48_r.c: Likewise.
* stdlib/lrand48.c: Likewise.
* stdlib/lrand48_r.c: Likewise.
* stdlib/mrand48.c: Likewise.
* stdlib/mrand48_r.c: Likewise.
* stdlib/seed48.c: Likewise.
* stdlib/seed48_r.c: Likewise.
* stdlib/srand48_r.c: Likewise.
1998-02-11 Andreas Jaeger <aj@arthur.rhein-neckar.de>
* nss/test-netdb.c: Add some more test cases.
1998-02-13 11:39 Ulrich Drepper <drepper@cygnus.com>
* libio/iovsscanf.c: Undo last change modifying errno.
1998-02-12 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>
* stdio-common/vfscanf.c: Never try to read another character
after EOF. Don't decrement read_in after EOF, it wasn't
incremented in the first place.
(NEXT_WIDE_CHAR): Set First, not first.
1998-02-06 07:48 H.J. Lu <hjl@gnu.org>
* db/Makefile ($(inst_libdir)/libndbm.a,
$(inst_libdir)/libndbm.so): New targets.
* db2/Makefile: Likewise.
1998-02-12 08:20 H.J. Lu <hjl@gnu.org>
* sysdeps/gnu/errlist.awk (sys_errlist, sys_nerr): Create weak
aliases if HAVE_ELF or PIC or DO_VERSIONING is not defined.
1998-02-12 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>
* sysdeps/generic/_G_config.h: Define _G_wchar_t, for C++
<streambuf.h>.
* sysdeps/unix/sysv/linux/_G_config.h: Likewise.
1998-02-11 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>
* sysdeps/unix/make-syscalls.sh: Fix sed pattern when dealing with
versioned symbols.
1998-02-13 08:14 H.J. Lu <hjl@gnu.org>
* libc.map (_dl_global_scope, _dl_lookup_symbol_skip,
_dl_lookup_versioned_symbol, _dl_lookup_versioned_symbol_skip):
Added for libdl.so.
* elf/rtld.map: New file. Needed to define the GLIBC_2.*
* manual/socket.texi (Host Address Functions): Clarify description
* sysdeps/unix/sysv/linux/alpha/bits/time.h (struct timeval):
1998-02-13 18:54:15 +01:00
|
|
|
}
|
2001-03-12 01:04:52 +01:00
|
|
|
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;
|
Update.
1998-02-13 17:39 Ulrich Drepper <drepper@cygnus.com>
* elf/Makefile: Don't use --version-script parameter to link ld.so
unconditionally.
1998-01-02 04:19 Geoff Keating <geoffk@ozemail.com.au>
* math/Makefile: Add t_exp.
* math/libm-test.c: Tighten accuracy bounds for exp(), correct
constants.
* math/test-reduce.c: Remove temporarily, it seems to be broken.
* sysdeps/libm-ieee754/e_exp.c: Use accurate table method.
* sysdeps/libm-ieee754/e_expf.c: Use table & double precision for
better accuracy.
* sysdeps/libm-ieee754/s_exp2.c: Use better polynomial; correct
algorithm for very large/very small arguments.
* sysdeps/libm-ieee754/s_exp2f.c: Use slightly better polynomial;
correct algorithm for very large/very small arguments; adjust for
new table.
* sysdeps/libm-ieee754/t_exp.c: New file.
* sysdeps/libm-ieee754/t_exp2f.h: Use table with smaller deltas.
* sysdeps/unix/sysv/linux/powerpc/dl-sysdep.c: Put 'strange test'
back, with comment that explains what breaks when you remove it :-(.
* localedata/xfrm-test.c: Avoid integer overflow.
* stdlib/strfmon.c: char is unsigned, sometimes.
*sysdeps/powerpc
* sysdeps/powerpc/Makefile: Remove quad float support.
* sysdeps/powerpc/q_*.c: Remove, they will become an add-on.
* sysdeps/powerpc/quad_float.h: Likewise.
* sysdeps/powerpc/test-arith.c: Likewise.
* sysdeps/powerpc/test-arithf.c: Likewise.
* sysdeps/generic/s_exp2.c: Remove, we have this implemented now.
* sysdeps/generic/s_exp2f.c: Likewise.
* sysdeps/powerpc/bits/mathinline.h: Use underscores around __asm__,
don't try anything if _SOFT_FLOAT.
1997-12-31 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>
* locale/C-ctype.c (_nl_C_LC_CTYPE_class32): Undo last change.
* locale/programs/ld-ctype.c (CHAR_CLASS32_TRANS): Likewise.
* wctype/wctype.c: Likewise.
* wctype/wctype.h (_ISwxxx): Renamed from _ISxxx, all uses
changed. They are incompatible with the _ISxxx values from
<ctype.h> on little endian machines.
(_ISwbit) [__BYTE_ORDER == __LITTLE_ENDIAN]: Correctly transform
bit number. This fixes the real bug and restores the integrity of
the ctype locale file.
* wctype/wcfuncs.c: Change all _ISxxx to _ISwxxx.
* wctype/wcfuncs_l.c: Likewise.
* wctype/wcextra.c: Likewise.
* wctype/wctype_l.c [__BYTE_ORDER == __LITTLE_ENDIAN]: Use correct
byte swapping.
1998-02-09 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>
* sysdeps/unix/sysv/linux/sparc/sparc32/sysdep.S (errno): Put it
into .bss segment instead of .common, so that aliases on it work.
* sysdeps/unix/sysv/linux/i386/sysdep.S (errno): Add .type and
.size directives, put into .bss segment instead of initializing it
to 4.
1998-02-12 08:00 H.J. Lu <hjl@gnu.org>
* libc.map (gnu_get_libc_release, gnu_get_libc_version): Added.
* version.c (__gnu_get_libc_release, __gnu_get_libc_version): New
functions.
Make names without __ weak aliases.
(__libc_release, __libc_version): Make them static.
* include/gnu/libc-version.h: New file.
* Makefile (headers): Add gnu/libc-version.h.
1998-02-13 Ulrich Drepper <drepper@cygnus.com>
* stdlib/stdlib.h (struct drand48_data): Leave X to user macros
and use x for member name.
Reported by Daniel Lyddy <daniell@cs.berkeley.edu>.
* stdlib/drand48.c: Change according to member name change.
* stdlib/drand48_r.c: Likewise.
* stdlib/lcong48_r.c: Likewise.
* stdlib/lrand48.c: Likewise.
* stdlib/lrand48_r.c: Likewise.
* stdlib/mrand48.c: Likewise.
* stdlib/mrand48_r.c: Likewise.
* stdlib/seed48.c: Likewise.
* stdlib/seed48_r.c: Likewise.
* stdlib/srand48_r.c: Likewise.
1998-02-11 Andreas Jaeger <aj@arthur.rhein-neckar.de>
* nss/test-netdb.c: Add some more test cases.
1998-02-13 11:39 Ulrich Drepper <drepper@cygnus.com>
* libio/iovsscanf.c: Undo last change modifying errno.
1998-02-12 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>
* stdio-common/vfscanf.c: Never try to read another character
after EOF. Don't decrement read_in after EOF, it wasn't
incremented in the first place.
(NEXT_WIDE_CHAR): Set First, not first.
1998-02-06 07:48 H.J. Lu <hjl@gnu.org>
* db/Makefile ($(inst_libdir)/libndbm.a,
$(inst_libdir)/libndbm.so): New targets.
* db2/Makefile: Likewise.
1998-02-12 08:20 H.J. Lu <hjl@gnu.org>
* sysdeps/gnu/errlist.awk (sys_errlist, sys_nerr): Create weak
aliases if HAVE_ELF or PIC or DO_VERSIONING is not defined.
1998-02-12 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>
* sysdeps/generic/_G_config.h: Define _G_wchar_t, for C++
<streambuf.h>.
* sysdeps/unix/sysv/linux/_G_config.h: Likewise.
1998-02-11 Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>
* sysdeps/unix/make-syscalls.sh: Fix sed pattern when dealing with
versioned symbols.
1998-02-13 08:14 H.J. Lu <hjl@gnu.org>
* libc.map (_dl_global_scope, _dl_lookup_symbol_skip,
_dl_lookup_versioned_symbol, _dl_lookup_versioned_symbol_skip):
Added for libdl.so.
* elf/rtld.map: New file. Needed to define the GLIBC_2.*
* manual/socket.texi (Host Address Functions): Clarify description
* sysdeps/unix/sysv/linux/alpha/bits/time.h (struct timeval):
1998-02-13 18:54:15 +01:00
|
|
|
}
|
2001-03-12 01:04:52 +01:00
|
|
|
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;
|
2001-03-12 08:57:09 +01:00
|
|
|
mynumber junk1, junk2, binexp = {{0,0}};
|
|
|
|
#if 0
|
|
|
|
int4 k;
|
|
|
|
#endif
|
|
|
|
int4 i,j,m,n,ex;
|
2001-03-12 01:04:52 +01:00
|
|
|
|
|
|
|
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;
|
|
|
|
}
|
1996-03-05 22:41:30 +01:00
|
|
|
}
|