lib/libm/shichif.c

   1 /*                                                      shichif.c
   2  *
   3  *      Hyperbolic sine and cosine integrals
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, Chi, Shi;
  10  *
  11  * shichi( x, &Chi, &Shi );
  12  *
  13  *
  14  * DESCRIPTION:
  15  *
  16  * Approximates the integrals
  17  *
  18  *                            x
  19  *                            -
  20  *                           | |   cosh t - 1
  21  *   Chi(x) = eul + ln x +   |    -----------  dt,
  22  *                         | |          t
  23  *                          -
  24  *                          0
  25  *
  26  *               x
  27  *               -
  28  *              | |  sinh t
  29  *   Shi(x) =   |    ------  dt
  30  *            | |       t
  31  *             -
  32  *             0
  33  *
  34  * where eul = 0.57721566490153286061 is Euler's constant.
  35  * The integrals are evaluated by power series for x < 8
  36  * and by Chebyshev expansions for x between 8 and 88.
  37  * For large x, both functions approach exp(x)/2x.
  38  * Arguments greater than 88 in magnitude return MAXNUM.
  39  *
  40  *
  41  * ACCURACY:
  42  *
  43  * Test interval 0 to 88.
  44  *                      Relative error:
  45  * arithmetic   function  # trials      peak         rms
  46  *    IEEE         Shi      20000       3.5e-7      7.0e-8
  47  *        Absolute error, except relative when |Chi| > 1:
  48  *    IEEE         Chi      20000       3.8e-7      7.6e-8
  49  */
  50 \f
  51 /*
  52 Cephes Math Library Release 2.2:  July, 1992
  53 Copyright 1984, 1987, 1992 by Stephen L. Moshier
  54 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  55 */
  56
  57
  58 #include "mconf.h"
  59
  60 /* x exp(-x) shi(x), inverted interval 8 to 18 */
  61 static float S1[] = {
  62 -3.56699611114982536845E-8,
  63  1.44818877384267342057E-7,
  64  7.82018215184051295296E-7,
  65 -5.39919118403805073710E-6,
  66 -3.12458202168959833422E-5,
  67  8.90136741950727517826E-5,
  68  2.02558474743846862168E-3,
  69  2.96064440855633256972E-2,
  70  1.11847751047257036625E0
  71 };
  72
  73 /* x exp(-x) shi(x), inverted interval 18 to 88 */
  74 static float S2[] = {
  75  1.69050228879421288846E-8,
  76  1.25391771228487041649E-7,
  77  1.16229947068677338732E-6,
  78  1.61038260117376323993E-5,
  79  3.49810375601053973070E-4,
  80  1.28478065259647610779E-2,
  81  1.03665722588798326712E0
  82 };
  83
  84
  85 /* x exp(-x) chin(x), inverted interval 8 to 18 */
  86 static float C1[] = {
  87  1.31458150989474594064E-8,
  88 -4.75513930924765465590E-8,
  89 -2.21775018801848880741E-7,
  90  1.94635531373272490962E-6,
  91  4.33505889257316408893E-6,
  92 -6.13387001076494349496E-5,
  93 -3.13085477492997465138E-4,
  94  4.97164789823116062801E-4,
  95  2.64347496031374526641E-2,
  96  1.11446150876699213025E0
  97 };
  98
  99 /* x exp(-x) chin(x), inverted interval 18 to 88 */
 100 static float C2[] = {
 101 -3.00095178028681682282E-9,
 102  7.79387474390914922337E-8,
 103  1.06942765566401507066E-6,
 104  1.59503164802313196374E-5,
 105  3.49592575153777996871E-4,
 106  1.28475387530065247392E-2,
 107  1.03665693917934275131E0
 108 };
 109
 110
 111
 112 /* Sine and cosine integrals */
 113
 114 #define EUL 0.57721566490153286061
 115 extern float MACHEPF, MAXNUMF;
 116
 117 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
 118
 119 #ifdef ANSIC
 120 float logf(float ), expf(float), chbevlf(float, float *, int);
 121 #else
 122 float logf(), expf(), chbevlf();
 123 #endif
 124
 125
 126
 127 #ifdef ANSIC
 128 int shichif( float xx, float *si, float *ci )
 129 #else
 130 int shichif( xx, si, ci )
 131 double xx;
 132 float *si, *ci;
 133 #endif
 134 {
 135 float x, k, z, c, s, a;
 136 short sign;
 137
 138 x = xx;
 139 if( x < 0.0 )
 140         {
 141         sign = -1;
 142         x = -x;
 143         }
 144 else
 145         sign = 0;
 146
 147
 148 if( x == 0.0 )
 149         {
 150         *si = 0.0;
 151         *ci = -MAXNUMF;
 152         return( 0 );
 153         }
 154
 155 if( x >= 8.0 )
 156         goto chb;
 157
 158 z = x * x;
 159
 160 /*      Direct power series expansion   */
 161
 162 a = 1.0;
 163 s = 1.0;
 164 c = 0.0;
 165 k = 2.0;
 166
 167 do
 168         {
 169         a *= z/k;
 170         c += a/k;
 171         k += 1.0;
 172         a /= k;
 173         s += a/k;
 174         k += 1.0;
 175         }
 176 while( fabsf(a/s) > MACHEPF );
 177
 178 s *= x;
 179 goto done;
 180
 181
 182 chb:
 183
 184 if( x < 18.0 )
 185         {
 186         a = (576.0/x - 52.0)/10.0;
 187         k = expf(x) / x;
 188         s = k * chbevlf( a, S1, 9 );
 189         c = k * chbevlf( a, C1, 10 );
 190         goto done;
 191         }
 192
 193 if( x <= 88.0 )
 194         {
 195         a = (6336.0/x - 212.0)/70.0;
 196         k = expf(x) / x;
 197         s = k * chbevlf( a, S2, 7 );
 198         c = k * chbevlf( a, C2, 7 );
 199         goto done;
 200         }
 201 else
 202         {
 203         if( sign )
 204                 *si = -MAXNUMF;
 205         else
 206                 *si = MAXNUMF;
 207         *ci = MAXNUMF;
 208         return(0);
 209         }
 210 done:
 211 if( sign )
 212         s = -s;
 213
 214 *si = s;
 215
 216 *ci = EUL + logf(x) + c;
 217 return(0);
 218 }