summaryrefslogtreecommitdiff
path: root/libm/float/shichif.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/float/shichif.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/float/shichif.c')
-rw-r--r--libm/float/shichif.c212
1 files changed, 0 insertions, 212 deletions
diff --git a/libm/float/shichif.c b/libm/float/shichif.c
deleted file mode 100644
index ae98021a9..000000000
--- a/libm/float/shichif.c
+++ /dev/null
@@ -1,212 +0,0 @@
-/* shichif.c
- *
- * Hyperbolic sine and cosine integrals
- *
- *
- *
- * SYNOPSIS:
- *
- * float x, Chi, Shi;
- *
- * shichi( x, &Chi, &Shi );
- *
- *
- * DESCRIPTION:
- *
- * Approximates the integrals
- *
- * x
- * -
- * | | cosh t - 1
- * Chi(x) = eul + ln x + | ----------- dt,
- * | | t
- * -
- * 0
- *
- * x
- * -
- * | | sinh t
- * Shi(x) = | ------ dt
- * | | t
- * -
- * 0
- *
- * where eul = 0.57721566490153286061 is Euler's constant.
- * The integrals are evaluated by power series for x < 8
- * and by Chebyshev expansions for x between 8 and 88.
- * For large x, both functions approach exp(x)/2x.
- * Arguments greater than 88 in magnitude return MAXNUM.
- *
- *
- * ACCURACY:
- *
- * Test interval 0 to 88.
- * Relative error:
- * arithmetic function # trials peak rms
- * IEEE Shi 20000 3.5e-7 7.0e-8
- * Absolute error, except relative when |Chi| > 1:
- * IEEE Chi 20000 3.8e-7 7.6e-8
- */
-
-/*
-Cephes Math Library Release 2.2: July, 1992
-Copyright 1984, 1987, 1992 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-*/
-
-
-#include <math.h>
-
-/* x exp(-x) shi(x), inverted interval 8 to 18 */
-static float S1[] = {
--3.56699611114982536845E-8,
- 1.44818877384267342057E-7,
- 7.82018215184051295296E-7,
--5.39919118403805073710E-6,
--3.12458202168959833422E-5,
- 8.90136741950727517826E-5,
- 2.02558474743846862168E-3,
- 2.96064440855633256972E-2,
- 1.11847751047257036625E0
-};
-
-/* x exp(-x) shi(x), inverted interval 18 to 88 */
-static float S2[] = {
- 1.69050228879421288846E-8,
- 1.25391771228487041649E-7,
- 1.16229947068677338732E-6,
- 1.61038260117376323993E-5,
- 3.49810375601053973070E-4,
- 1.28478065259647610779E-2,
- 1.03665722588798326712E0
-};
-
-
-/* x exp(-x) chin(x), inverted interval 8 to 18 */
-static float C1[] = {
- 1.31458150989474594064E-8,
--4.75513930924765465590E-8,
--2.21775018801848880741E-7,
- 1.94635531373272490962E-6,
- 4.33505889257316408893E-6,
--6.13387001076494349496E-5,
--3.13085477492997465138E-4,
- 4.97164789823116062801E-4,
- 2.64347496031374526641E-2,
- 1.11446150876699213025E0
-};
-
-/* x exp(-x) chin(x), inverted interval 18 to 88 */
-static float C2[] = {
--3.00095178028681682282E-9,
- 7.79387474390914922337E-8,
- 1.06942765566401507066E-6,
- 1.59503164802313196374E-5,
- 3.49592575153777996871E-4,
- 1.28475387530065247392E-2,
- 1.03665693917934275131E0
-};
-
-
-
-/* Sine and cosine integrals */
-
-#define EUL 0.57721566490153286061
-extern float MACHEPF, MAXNUMF;
-
-#define fabsf(x) ( (x) < 0 ? -(x) : (x) )
-
-#ifdef ANSIC
-float logf(float ), expf(float), chbevlf(float, float *, int);
-#else
-float logf(), expf(), chbevlf();
-#endif
-
-
-
-int shichif( float xx, float *si, float *ci )
-{
-float x, k, z, c, s, a;
-short sign;
-
-x = xx;
-if( x < 0.0 )
- {
- sign = -1;
- x = -x;
- }
-else
- sign = 0;
-
-
-if( x == 0.0 )
- {
- *si = 0.0;
- *ci = -MAXNUMF;
- return( 0 );
- }
-
-if( x >= 8.0 )
- goto chb;
-
-z = x * x;
-
-/* Direct power series expansion */
-
-a = 1.0;
-s = 1.0;
-c = 0.0;
-k = 2.0;
-
-do
- {
- a *= z/k;
- c += a/k;
- k += 1.0;
- a /= k;
- s += a/k;
- k += 1.0;
- }
-while( fabsf(a/s) > MACHEPF );
-
-s *= x;
-goto done;
-
-
-chb:
-
-if( x < 18.0 )
- {
- a = (576.0/x - 52.0)/10.0;
- k = expf(x) / x;
- s = k * chbevlf( a, S1, 9 );
- c = k * chbevlf( a, C1, 10 );
- goto done;
- }
-
-if( x <= 88.0 )
- {
- a = (6336.0/x - 212.0)/70.0;
- k = expf(x) / x;
- s = k * chbevlf( a, S2, 7 );
- c = k * chbevlf( a, C2, 7 );
- goto done;
- }
-else
- {
- if( sign )
- *si = -MAXNUMF;
- else
- *si = MAXNUMF;
- *ci = MAXNUMF;
- return(0);
- }
-done:
-if( sign )
- s = -s;
-
-*si = s;
-
-*ci = EUL + logf(x) + c;
-return(0);
-}