/Python/pymath.c

  1#include "Python.h"
  2
  3#ifdef X87_DOUBLE_ROUNDING
  4/* On x86 platforms using an x87 FPU, this function is called from the
  5   Py_FORCE_DOUBLE macro (defined in pymath.h) to force a floating-point
  6   number out of an 80-bit x87 FPU register and into a 64-bit memory location,
  7   thus rounding from extended precision to double precision. */
  8double _Py_force_double(double x)
  9{
 10	volatile double y;
 11	y = x;
 12	return y;
 13}
 14#endif
 15
 16#ifndef HAVE_HYPOT
 17double hypot(double x, double y)
 18{
 19	double yx;
 20
 21	x = fabs(x);
 22	y = fabs(y);
 23	if (x < y) {
 24		double temp = x;
 25		x = y;
 26		y = temp;
 27	}
 28	if (x == 0.)
 29		return 0.;
 30	else {
 31		yx = y/x;
 32		return x*sqrt(1.+yx*yx);
 33	}
 34}
 35#endif /* HAVE_HYPOT */
 36
 37#ifndef HAVE_COPYSIGN
 38double
 39copysign(double x, double y)
 40{
 41	/* use atan2 to distinguish -0. from 0. */
 42	if (y > 0. || (y == 0. && atan2(y, -1.) > 0.)) {
 43		return fabs(x);
 44	} else {
 45		return -fabs(x);
 46	}
 47}
 48#endif /* HAVE_COPYSIGN */
 49
 50#ifndef HAVE_LOG1P
 51#include <float.h>
 52
 53double
 54log1p(double x)
 55{
 56	/* For x small, we use the following approach.  Let y be the nearest
 57	   float to 1+x, then
 58
 59	     1+x = y * (1 - (y-1-x)/y)
 60
 61	   so log(1+x) = log(y) + log(1-(y-1-x)/y).  Since (y-1-x)/y is tiny,
 62	   the second term is well approximated by (y-1-x)/y.  If abs(x) >=
 63	   DBL_EPSILON/2 or the rounding-mode is some form of round-to-nearest
 64	   then y-1-x will be exactly representable, and is computed exactly
 65	   by (y-1)-x.
 66
 67	   If abs(x) < DBL_EPSILON/2 and the rounding mode is not known to be
 68	   round-to-nearest then this method is slightly dangerous: 1+x could
 69	   be rounded up to 1+DBL_EPSILON instead of down to 1, and in that
 70	   case y-1-x will not be exactly representable any more and the
 71	   result can be off by many ulps.  But this is easily fixed: for a
 72	   floating-point number |x| < DBL_EPSILON/2., the closest
 73	   floating-point number to log(1+x) is exactly x.
 74	*/
 75
 76	double y;
 77	if (fabs(x) < DBL_EPSILON/2.) {
 78		return x;
 79	} else if (-0.5 <= x && x <= 1.) {
 80		/* WARNING: it's possible than an overeager compiler
 81		   will incorrectly optimize the following two lines
 82		   to the equivalent of "return log(1.+x)". If this
 83		   happens, then results from log1p will be inaccurate
 84		   for small x. */
 85		y = 1.+x;
 86		return log(y)-((y-1.)-x)/y;
 87	} else {
 88		/* NaNs and infinities should end up here */
 89		return log(1.+x);
 90	}
 91}
 92#endif /* HAVE_LOG1P */
 93
 94/*
 95 * ====================================================
 96 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 97 *
 98 * Developed at SunPro, a Sun Microsystems, Inc. business.
 99 * Permission to use, copy, modify, and distribute this
100 * software is freely granted, provided that this notice 
101 * is preserved.
102 * ====================================================
103 */
104
105static const double ln2 = 6.93147180559945286227E-01;
106static const double two_pow_m28 = 3.7252902984619141E-09; /* 2**-28 */
107static const double two_pow_p28 = 268435456.0; /* 2**28 */
108static const double zero = 0.0;
109
110/* asinh(x)
111 * Method :
112 *	Based on 
113 *		asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
114 *	we have
115 *	asinh(x) := x  if  1+x*x=1,
116 *		 := sign(x)*(log(x)+ln2)) for large |x|, else
117 *		 := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
118 *		 := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))  
119 */
120
121#ifndef HAVE_ASINH
122double
123asinh(double x)
124{	
125	double w;
126	double absx = fabs(x);
127
128	if (Py_IS_NAN(x) || Py_IS_INFINITY(x)) {
129		return x+x;
130	}
131	if (absx < two_pow_m28) {	/* |x| < 2**-28 */
132		return x;	/* return x inexact except 0 */
133	} 
134	if (absx > two_pow_p28) {	/* |x| > 2**28 */
135		w = log(absx)+ln2;
136	}
137	else if (absx > 2.0) {		/* 2 < |x| < 2**28 */
138		w = log(2.0*absx + 1.0 / (sqrt(x*x + 1.0) + absx));
139	}
140	else {				/* 2**-28 <= |x| < 2= */
141		double t = x*x;
142		w = log1p(absx + t / (1.0 + sqrt(1.0 + t)));
143	}
144	return copysign(w, x);
145	
146}
147#endif /* HAVE_ASINH */
148
149/* acosh(x)
150 * Method :
151 *      Based on
152 *	      acosh(x) = log [ x + sqrt(x*x-1) ]
153 *      we have
154 *	      acosh(x) := log(x)+ln2, if x is large; else
155 *	      acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
156 *	      acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
157 *
158 * Special cases:
159 *      acosh(x) is NaN with signal if x<1.
160 *      acosh(NaN) is NaN without signal.
161 */
162
163#ifndef HAVE_ACOSH
164double
165acosh(double x)
166{
167	if (Py_IS_NAN(x)) {
168		return x+x;
169	}
170	if (x < 1.) {			/* x < 1;  return a signaling NaN */
171		errno = EDOM;
172#ifdef Py_NAN
173		return Py_NAN;
174#else
175		return (x-x)/(x-x);
176#endif
177	}
178	else if (x >= two_pow_p28) {	/* x > 2**28 */
179		if (Py_IS_INFINITY(x)) {
180			return x+x;
181		} else {
182			return log(x)+ln2;	/* acosh(huge)=log(2x) */
183		}
184	}
185	else if (x == 1.) {
186		return 0.0;			/* acosh(1) = 0 */
187	}
188	else if (x > 2.) {			/* 2 < x < 2**28 */
189		double t = x*x;
190		return log(2.0*x - 1.0 / (x + sqrt(t - 1.0)));
191	}
192	else {				/* 1 < x <= 2 */
193		double t = x - 1.0;
194		return log1p(t + sqrt(2.0*t + t*t));
195	}
196}
197#endif /* HAVE_ACOSH */
198
199/* atanh(x)
200 * Method :
201 *    1.Reduced x to positive by atanh(-x) = -atanh(x)
202 *    2.For x>=0.5
203 *		  1	      2x			  x
204 *      atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
205 *		  2	     1 - x		      1 - x
206 *
207 *      For x<0.5
208 *      atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
209 *
210 * Special cases:
211 *      atanh(x) is NaN if |x| >= 1 with signal;
212 *      atanh(NaN) is that NaN with no signal;
213 *
214 */
215
216#ifndef HAVE_ATANH
217double
218atanh(double x)
219{
220	double absx;
221	double t;
222
223	if (Py_IS_NAN(x)) {
224		return x+x;
225	}
226	absx = fabs(x);
227	if (absx >= 1.) {		/* |x| >= 1 */
228		errno = EDOM;
229#ifdef Py_NAN
230		return Py_NAN;
231#else
232		return x/zero;
233#endif
234	}
235	if (absx < two_pow_m28) {	/* |x| < 2**-28 */
236		return x;
237	}
238	if (absx < 0.5) {		/* |x| < 0.5 */
239		t = absx+absx;
240		t = 0.5 * log1p(t + t*absx / (1.0 - absx));
241	} 
242	else {				/* 0.5 <= |x| <= 1.0 */
243		t = 0.5 * log1p((absx + absx) / (1.0 - absx));
244	}
245	return copysign(t, x);
246}
247#endif /* HAVE_ATANH */