exp.c

操作系统SunOS 4.1.3版本的源码

#ifndef lint
static	char sccsid[] = "@(#)exp.c 1.1 92/07/30 SMI";
#endif

/*
 * Copyright (c) 1988 by Sun Microsystems, Inc.
 */

/* exp(x)
 * Hybrid alogrithm of Peter Tang's Table driven method (for 
 * large arguments) and rational approximation (for small arguments).
 * Written by K.C. Ng, November 1988.
 * Method (Table driven):
 *	1. Argument Reduction: given the input x, find r and integer k 
 *	   and j such that
 *	             x = (32k+j)*ln2 + r,  |r| <= (1/64)*ln2 .  
 *
 *	2. exp(x) = 2^k * (2^(j/32) + 2^(j/32)*expm1(r)) 
 *	   Note:
 *	   a. expm1(r) is computed by 
 *		(i)  If divide is fast:
 *		     expm1(r) = (2r)/(2-R), R = r - r^2*(t1 + t2*r^2)
 *		(ii) Otherwise (divide is more than 7 multiplications time)
 *		     expm1(r) = r + t1*r^2 + t2*r^3 + ... + t5*r^6
 *	   b. 2^(j/32) is represented as S[j]+S_trail[j] where
 *	      S[j] = 2^(j/32) rounded and S_trail[j] = 2^(j/32) - S[j].
 *	      We also define S2[j]=2*S[j] for computing a(i)'s expm1: 
 *			expm1(r) = (S2[j]*r)/(2.0-R)
 *			
 * Special cases:
 *	exp(INF) is INF, exp(NaN) is NaN;
 *	exp(-INF)=  0;
 *	for finite argument, only exp(0)=1 is exact.
 *
 * Accuracy:
 *	according to an error analysis, the error is always less than
 *	an ulp (unit in the last place).
 *
 * Misc. info.
 *	For IEEE double 
 *		if x >  7.09782712893383973096e+02 then exp(x) overflow
 *		if x < -7.45133219101941108420e+02 then exp(x) underflow
 *
 * Constants:
 * The hexadecimal values are the intended ones for the following constants.
 * The decimal values may be used, provided that the compiler will convert
 * from decimal to binary accurately enough to produce the hexadecimal values
 * shown.
 */

#include <math.h>

extern double SVID_libm_err();

static double
threshold1	= 7.09782712893383973096e+02,	/* 0x40862E42, 0xFEFA39EF */
threshold2	= 7.45133219101941108420e+02,	/* 0x40874910, 0xD52D3051 */
huge		= 1.0e30,
one		= 1.0,
ln2_onehalf     = 1.03972077083991796412e+00,   /* 0x3FF0A2B2, 0x3F3BAB73 */
ln2		= 6.93147180559945286227e-01,	/* 0x3fe62e42, 0xfefa39ef */
ln2_2		= 3.46573590279972643113e-01,	/* 0x3fd62e42, 0xfefa39ef */
ln2hi           = 6.93147180369123816490e-01,   /* 0x3fe62e42, 0xfee00000 */
ln2lo           = 1.90821492927058770002e-10,   /* 0x3dea39ef, 0x35793c76 */
ln2_64		= 1.08304246962491450973e-02,   /* 0x3f862e42, 0xfefa39ef */
ln2_32hi	= 2.16608493865351192653e-02,	/* 0x3f962e42, 0xfee00000 */
ln2_32lo	= 5.96317165397058656257e-12,	/* 0x3d9a39ef, 0x35793c76 */
invln2		= 1.44269504088896338700e+00,	/* 0x3ff71547, 0x652b82fe */
invln2_32       = 4.61662413084468283841e+01,   /* 0x40471547, 0x652b82fe */
twom18		= 3.81469726562500000000e-06,	/* 0x3ed00000, 0x00000000 */
twom28		= 3.72529029846191406250e-09,	/* 0x3e300000, 0x00000000 */
twom54		= 5.55111512312578270212e-17;	/* 0x3c900000, 0x00000000 */

/* rational or polynomial approximation on [0,(ln2)/2] */
#ifdef SLOWDIV
static double
p1 = 5.00000000000000000000e-1,
p2 = 1.66666666666666784982e-1,
p3 = 4.16666666666197204123e-2,
p4 = 8.33333333331901796590e-3,
p5 = 1.38888889206249588980e-3,
p6 = 1.98412698941833241823e-4,
p7 = 2.48015161960545453685e-5,
p8 = 2.75572352833991096674e-6,
p9 = 2.76222881768990971035e-7,
p10= 2.51123092975208043742e-8;
#else
static double
p1     =  1.6666666666666601904E-1    , /*Hex  2^-3    *  1.555555555553E */
p2     = -2.7777777777015593384E-3    , /*Hex  2^-9    * -1.6C16C16BEBD93 */
p3     =  6.6137563214379343612E-5    , /*Hex  2^-14   *  1.1566AAF25DE2C */
p4     = -1.6533902205465251539E-6    , /*Hex  2^-20   * -1.BBD41C5D26BF1 */
p5     =  4.1381367970572384604E-8    ; /*Hex  2^-25   *  1.6376972BEA4D0 */
#endif

/* rational or polynomial approximation on [0,(ln2)/64] */
#ifdef SLOWDIV
static double
t1    =   5.0000000000000000000e-1    , /* 3fe0000000000000 */ 
t2    =   1.6666666666526086527e-1    , /* 3fc5555555548f7c */ 
t3    =   4.1666666666226079285e-2    , /* 3fa5555555545d4e */ 
t4    =   8.3333679843421958056e-3    , /* 3f811115b7aa905e */ 
t5    =   1.3888949086377719040e-3    ; /* 3f56c1728d739765 */ 
#else
static double		
t1    =   1.6666666666627465031E-1    , /* 3FC5555555551E29 */
t2    =  -2.7777669763261076298E-3    ; /* BF66C1664A3720A8 */
#endif

static double S[] = {
 1.00000000000000000000e+00	, /* 3FF0000000000000 */
 1.02189714865411662714e+00	, /* 3FF059B0D3158574 */
 1.04427378242741375480e+00	, /* 3FF0B5586CF9890F */
 1.06714040067682369717e+00	, /* 3FF11301D0125B51 */
 1.09050773266525768967e+00	, /* 3FF172B83C7D517B */
 1.11438674259589243221e+00	, /* 3FF1D4873168B9AA */
 1.13878863475669156458e+00	, /* 3FF2387A6E756238 */
 1.16372485877757747552e+00	, /* 3FF29E9DF51FDEE1 */
 1.18920711500272102690e+00	, /* 3FF306FE0A31B715 */
 1.21524735998046895524e+00	, /* 3FF371A7373AA9CB */
 1.24185781207348400201e+00	, /* 3FF3DEA64C123422 */
 1.26905095719173321989e+00	, /* 3FF44E086061892D */
 1.29683955465100964055e+00	, /* 3FF4BFDAD5362A27 */
 1.32523664315974132322e+00	, /* 3FF5342B569D4F82 */
 1.35425554693689265129e+00	, /* 3FF5AB07DD485429 */
 1.38390988196383202258e+00	, /* 3FF6247EB03A5585 */
 1.41421356237309514547e+00	, /* 3FF6A09E667F3BCD */
 1.44518080697704665027e+00	, /* 3FF71F75E8EC5F74 */
 1.47682614593949934623e+00	, /* 3FF7A11473EB0187 */
 1.50916442759342284141e+00	, /* 3FF82589994CCE13 */
 1.54221082540794074411e+00	, /* 3FF8ACE5422AA0DB */
 1.57598084510788649659e+00	, /* 3FF93737B0CDC5E5 */
 1.61049033194925428347e+00	, /* 3FF9C49182A3F090 */
 1.64575547815396494578e+00	, /* 3FFA5503B23E255D */
 1.68179283050742900407e+00	, /* 3FFAE89F995AD3AD */
 1.71861929812247793414e+00	, /* 3FFB7F76F2FB5E47 */
 1.75625216037329945351e+00	, /* 3FFC199BDD85529C */
 1.79470907500310716820e+00	, /* 3FFCB720DCEF9069 */
 1.83400808640934243066e+00	, /* 3FFD5818DCFBA487 */
 1.87416763411029996256e+00	, /* 3FFDFC97337B9B5F */
 1.91520656139714740007e+00	, /* 3FFEA4AFA2A490DA */
 1.95714412417540017941e+00	, /* 3FFF50765B6E4540 */
};

static double S_trail[] = {
 0.00000000000000000000e+00	, 
 5.10922502897344389359e-17	, /* 3C8D73E2A475B465 */
 8.55188970553796365958e-17	, /* 3C98A62E4ADC610A */
-7.89985396684158212226e-17	, /* BC96C51039449B3A */
-3.04678207981247114697e-17	, /* BC819041B9D78A76 */
 1.04102784568455709549e-16	, /* 3C9E016E00A2643C */
 8.91281267602540777782e-17	, /* 3C99B07EB6C70573 */
 3.82920483692409349872e-17	, /* 3C8612E8AFAD1255 */
 3.98201523146564611098e-17	, /* 3C86F46AD23182E4 */
-7.71263069268148813091e-17	, /* BC963AEABF42EAE2 */
 4.65802759183693679123e-17	, /* 3C8ADA0911F09EBC */
 2.66793213134218609523e-18	, /* 3C489B7A04EF80D0 */
 2.53825027948883149593e-17	, /* 3C7D4397AFEC42E2 */
-2.85873121003886075697e-17	, /* BC807ABE1DB13CAC */
 7.70094837980298946162e-17	, /* 3C96324C054647AD */
-6.77051165879478628716e-17	, /* BC9383C17E40B497 */
-9.66729331345291345105e-17	, /* BC9BDD3413B26456 */
-3.02375813499398731940e-17	, /* BC816E4786887A99 */
-3.48399455689279579579e-17	, /* BC841577EE04992F */
-1.01645532775429503911e-16	, /* BC9D4C1DD41532D8 */
 7.94983480969762085616e-17	, /* 3C96E9F156864B27 */
-1.01369164712783039808e-17	, /* BC675FC781B57EBC */
 2.47071925697978878522e-17	, /* 3C7C7C46B071F2BE */
-1.01256799136747726038e-16	, /* BC9D2F6EDB8D41E1 */
 8.19901002058149652013e-17	, /* 3C97A1CD345DCC81 */
-1.85138041826311098821e-17	, /* BC75584F7E54AC3B */
 2.96014069544887330703e-17	, /* 3C811065895048DD */
 1.82274584279120867698e-17	, /* 3C7503CBD1E949DB */
 3.28310722424562658722e-17	, /* 3C82ED02D75B3706 */
-6.12276341300414256164e-17	, /* BC91A5CD4F184B5C */
-1.06199460561959626376e-16	, /* BC9E9C23179C2893 */
 8.96076779103666776760e-17	, /* 3C99D3E12DD8A18B */
};

static double S2[] = {
 2.00000000000000000000e+00	, /* 4000000000000000 */
 2.04379429730823325428e+00	, /* 400059B0D3158574 */
 2.08854756485482750961e+00	, /* 4000B5586CF9890F */
 2.13428080135364739434e+00	, /* 40011301D0125B51 */
 2.18101546533051537935e+00	, /* 400172B83C7D517B */
 2.22877348519178486441e+00	, /* 4001D4873168B9AA */
 2.27757726951338312915e+00	, /* 4002387A6E756238 */
 2.32744971755515495104e+00	, /* 40029E9DF51FDEE1 */
 2.37841423000544205379e+00	, /* 400306FE0A31B715 */
 2.43049471996093791049e+00	, /* 400371A7373AA9CB */
 2.48371562414696800403e+00	, /* 4003DEA64C123422 */
 2.53810191438346643977e+00	, /* 40044E086061892D */
 2.59367910930201928110e+00	, /* 4004BFDAD5362A27 */
 2.65047328631948264643e+00	, /* 4005342B569D4F82 */
 2.70851109387378530258e+00	, /* 4005AB07DD485429 */
 2.76781976392766404516e+00	, /* 4006247EB03A5585 */
 2.82842712474619029095e+00	, /* 4006A09E667F3BCD */
 2.89036161395409330055e+00	, /* 40071F75E8EC5F74 */
 2.95365229187899869245e+00	, /* 4007A11473EB0187 */
 3.01832885518684568282e+00	, /* 40082589994CCE13 */
 3.08442165081588148823e+00	, /* 4008ACE5422AA0DB */
 3.15196169021577299318e+00	, /* 40093737B0CDC5E5 */
 3.22098066389850856694e+00	, /* 4009C49182A3F090 */
 3.29151095630792989155e+00	, /* 400A5503B23E255D */
 3.36358566101485800814e+00	, /* 400AE89F995AD3AD */
 3.43723859624495586829e+00	, /* 400B7F76F2FB5E47 */
 3.51250432074659890702e+00	, /* 400C199BDD85529C */
 3.58941815000621433640e+00	, /* 400CB720DCEF9069 */
 3.66801617281868486131e+00	, /* 400D5818DCFBA487 */
 3.74833526822059992512e+00	, /* 400DFC97337B9B5F */
 3.83041312279429480014e+00	, /* 400EA4AFA2A490DA */
 3.91428824835080035882e+00	, /* 400F50765B6E4540 */
};

double exp(x)
double x;
{
	double y,z,c,t,hi,lo;
	int k,xsb,j,m;

	y = fabs(x);

    /* filter out non-finite arugment */
	if(!finite(x)) {
	    if(x!=x) return x+x; else return (signbit(x)==0)? x:0.0;
	}

	if(y > ln2_64) {
	    xsb = signbit(x);
	    if (y<=ln2_onehalf) {
		k=0;
	        if(y>ln2_2) 
		    if(xsb==0)
			{hi = x - ln2hi; lo =  ln2lo; x = hi - lo; k =  1;}
		    else
			{hi = x + ln2hi; lo = -ln2lo; x = hi - lo; k = -1;}
		t  = x*x;
#ifdef SLOWDIV
		c  = (-t)*(p1+x*(p2+x*(p3+x*(p4+x*(p5+
			x*(p6+x*(p7+x*(p8+x*(p9+x*p10)))))))));
		if(k==0) return one-(c-x); 
		else {
		    y = one-((lo+c)-hi);
#else
		c  = x - t*(p1+t*(p2+t*(p3+t*(p4+t*p5))));
		if(k==0) return one-((x*c)/(c-2.0)-x); 
		else {
		    y = one-((lo-(x*c)/(2.0-c))-hi);
#endif
		    if(k==1) return y+y; else return 0.5*y;
		}
	    } else if(y>((xsb==0)?threshold1:threshold2)) {
		if(xsb==0) return SVID_libm_err(x,x,6); /* overflow */
		else 	   return SVID_libm_err(x,x,7); /* underflow */
	    } else {
		k  = invln2_32*x+((xsb==0)?0.5:-0.5);
		j  = k&0x1f; m = k>>5;
		t  = k;
		z  = (x - t*ln2_32hi) - t*ln2_32lo;
	    }
	} else if(y < twom18) {
	    dummy(x+huge);	/* raise inexact flags */
	    if(y < twom28)
		return one+x;
	    else
		return one+x*(one+0.5*x);
	} else {
	    z = x;
	    m = 0;
	    j = 0;
	}

    /* z is now in primary range */
#ifdef SLOWDIV
	t = (-z)*(1.0+z*(t1+z*(t2+z*(t3+z*(t4+z*t5)))));
	y = S[j]-(S[j]*t-S_trail[j]);
#else
	t = z*z;
	y = S[j]-((S2[j]*z)/((z-t*(t1+t*t2))-2.0)-S_trail[j]);
#endif
	if(m==0) 	return y; 
	else if (m < -1021)	
		return twom54*scalbn(y,m+54); 	/* subnormal output */
	else
		return scalbn(y,m);		/* guaranteed normal output */
}

static dummy(x)
double x;
{
	return 0;
}