summaryrefslogtreecommitdiff
path: root/libm/float/shichif.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/float/shichif.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/float/shichif.c')
-rw-r--r--libm/float/shichif.c212
1 files changed, 212 insertions, 0 deletions
diff --git a/libm/float/shichif.c b/libm/float/shichif.c
new file mode 100644
index 000000000..ae98021a9
--- /dev/null
+++ b/libm/float/shichif.c
@@ -0,0 +1,212 @@
+/* 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);
+}