lib/libm/sicif.c

   1 /*                                                      sicif.c
   2  *
   3  *      Sine and cosine integrals
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * float x, Ci, Si;
  10  *
  11  * sicif( x, &Si, &Ci );
  12  *
  13  *
  14  * DESCRIPTION:
  15  *
  16  * Evaluates the integrals
  17  *
  18  *                          x
  19  *                          -
  20  *                         |  cos t - 1
  21  *   Ci(x) = eul + ln x +  |  --------- dt,
  22  *                         |      t
  23  *                        -
  24  *                         0
  25  *             x
  26  *             -
  27  *            |  sin t
  28  *   Si(x) =  |  ----- dt
  29  *            |    t
  30  *           -
  31  *            0
  32  *
  33  * where eul = 0.57721566490153286061 is Euler's constant.
  34  * The integrals are approximated by rational functions.
  35  * For x > 8 auxiliary functions f(x) and g(x) are employed
  36  * such that
  37  *
  38  * Ci(x) = f(x) sin(x) - g(x) cos(x)
  39  * Si(x) = pi/2 - f(x) cos(x) - g(x) sin(x)
  40  *
  41  *
  42  * ACCURACY:
  43  *    Test interval = [0,50].
  44  * Absolute error, except relative when > 1:
  45  * arithmetic   function   # trials      peak         rms
  46  *    IEEE        Si        30000       2.1e-7      4.3e-8
  47  *    IEEE        Ci        30000       3.9e-7      2.2e-8
  48  */
  49 \f
  50 /*
  51 Cephes Math Library Release 2.1:  January, 1989
  52 Copyright 1984, 1987, 1989 by Stephen L. Moshier
  53 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  54 */
  55
  56 #include "mconf.h"
  57
  58 static float SN[] = {
  59 -8.39167827910303881427E-11,
  60  4.62591714427012837309E-8,
  61 -9.75759303843632795789E-6,
  62  9.76945438170435310816E-4,
  63 -4.13470316229406538752E-2,
  64  1.00000000000000000302E0,
  65 };
  66 static float SD[] = {
  67   2.03269266195951942049E-12,
  68   1.27997891179943299903E-9,
  69   4.41827842801218905784E-7,
  70   9.96412122043875552487E-5,
  71   1.42085239326149893930E-2,
  72   9.99999999999999996984E-1,
  73 };
  74
  75 static float CN[] = {
  76  2.02524002389102268789E-11,
  77 -1.35249504915790756375E-8,
  78  3.59325051419993077021E-6,
  79 -4.74007206873407909465E-4,
  80  2.89159652607555242092E-2,
  81 -1.00000000000000000080E0,
  82 };
  83 static float CD[] = {
  84   4.07746040061880559506E-12,
  85   3.06780997581887812692E-9,
  86   1.23210355685883423679E-6,
  87   3.17442024775032769882E-4,
  88   5.10028056236446052392E-2,
  89   4.00000000000000000080E0,
  90 };
  91
  92
  93 static float FN4[] = {
  94   4.23612862892216586994E0,
  95   5.45937717161812843388E0,
  96   1.62083287701538329132E0,
  97   1.67006611831323023771E-1,
  98   6.81020132472518137426E-3,
  99   1.08936580650328664411E-4,
 100   5.48900223421373614008E-7,
 101 };
 102 static float FD4[] = {
 103 /*  1.00000000000000000000E0,*/
 104   8.16496634205391016773E0,
 105   7.30828822505564552187E0,
 106   1.86792257950184183883E0,
 107   1.78792052963149907262E-1,
 108   7.01710668322789753610E-3,
 109   1.10034357153915731354E-4,
 110   5.48900252756255700982E-7,
 111 };
 112
 113
 114 static float FN8[] = {
 115   4.55880873470465315206E-1,
 116   7.13715274100146711374E-1,
 117   1.60300158222319456320E-1,
 118   1.16064229408124407915E-2,
 119   3.49556442447859055605E-4,
 120   4.86215430826454749482E-6,
 121   3.20092790091004902806E-8,
 122   9.41779576128512936592E-11,
 123   9.70507110881952024631E-14,
 124 };
 125 static float FD8[] = {
 126 /*  1.00000000000000000000E0,*/
 127   9.17463611873684053703E-1,
 128   1.78685545332074536321E-1,
 129   1.22253594771971293032E-2,
 130   3.58696481881851580297E-4,
 131   4.92435064317881464393E-6,
 132   3.21956939101046018377E-8,
 133   9.43720590350276732376E-11,
 134   9.70507110881952025725E-14,
 135 };
 136
 137 static float GN4[] = {
 138   8.71001698973114191777E-2,
 139   6.11379109952219284151E-1,
 140   3.97180296392337498885E-1,
 141   7.48527737628469092119E-2,
 142   5.38868681462177273157E-3,
 143   1.61999794598934024525E-4,
 144   1.97963874140963632189E-6,
 145   7.82579040744090311069E-9,
 146 };
 147 static float GD4[] = {
 148 /*  1.00000000000000000000E0,*/
 149   1.64402202413355338886E0,
 150   6.66296701268987968381E-1,
 151   9.88771761277688796203E-2,
 152   6.22396345441768420760E-3,
 153   1.73221081474177119497E-4,
 154   2.02659182086343991969E-6,
 155   7.82579218933534490868E-9,
 156 };
 157
 158 static float GN8[] = {
 159   6.97359953443276214934E-1,
 160   3.30410979305632063225E-1,
 161   3.84878767649974295920E-2,
 162   1.71718239052347903558E-3,
 163   3.48941165502279436777E-5,
 164   3.47131167084116673800E-7,
 165   1.70404452782044526189E-9,
 166   3.85945925430276600453E-12,
 167   3.14040098946363334640E-15,
 168 };
 169 static float GD8[] = {
 170 /*  1.00000000000000000000E0,*/
 171   1.68548898811011640017E0,
 172   4.87852258695304967486E-1,
 173   4.67913194259625806320E-2,
 174   1.90284426674399523638E-3,
 175   3.68475504442561108162E-5,
 176   3.57043223443740838771E-7,
 177   1.72693748966316146736E-9,
 178   3.87830166023954706752E-12,
 179   3.14040098946363335242E-15,
 180 };
 181
 182 #define EUL 0.57721566490153286061
 183 extern float MAXNUMF, PIO2F, MACHEPF;
 184
 185
 186
 187 #ifdef ANSIC
 188 float logf(float), sinf(float), cosf(float);
 189 float polevlf(float, float *, int);
 190 float p1evlf(float, float *, int);
 191 #else
 192 float logf(), sinf(), cosf(), polevlf(), p1evlf();
 193 #endif
 194
 195
 196 #ifdef ANSIC
 197 int sicif( float xx, float *si, float *ci )
 198 #else
 199 int sicif( xx, si, ci )
 200 double xx;
 201 float *si, *ci;
 202 #endif
 203 {
 204 float x, z, c, s, f, g;
 205 int sign;
 206
 207 x = xx;
 208 if( x < 0.0 )
 209         {
 210         sign = -1;
 211         x = -x;
 212         }
 213 else
 214         sign = 0;
 215
 216
 217 if( x == 0.0 )
 218         {
 219         *si = 0.0;
 220         *ci = -MAXNUMF;
 221         return( 0 );
 222         }
 223
 224
 225 if( x > 1.0e9 )
 226         {
 227         *si = PIO2F - cosf(x)/x;
 228         *ci = sinf(x)/x;
 229         return( 0 );
 230         }
 231
 232
 233
 234 if( x > 4.0 )
 235         goto asympt;
 236
 237 z = x * x;
 238 s = x * polevlf( z, SN, 5 ) / polevlf( z, SD, 5 );
 239 c = z * polevlf( z, CN, 5 ) / polevlf( z, CD, 5 );
 240
 241 if( sign )
 242         s = -s;
 243 *si = s;
 244 *ci = EUL + logf(x) + c;        /* real part if x < 0 */
 245 return(0);
 246
 247
 248
 249 /* The auxiliary functions are:
 250  *
 251  *
 252  * *si = *si - PIO2;
 253  * c = cos(x);
 254  * s = sin(x);
 255  *
 256  * t = *ci * s - *si * c;
 257  * a = *ci * c + *si * s;
 258  *
 259  * *si = t;
 260  * *ci = -a;
 261  */
 262
 263
 264 asympt:
 265
 266 s = sinf(x);
 267 c = cosf(x);
 268 z = 1.0/(x*x);
 269 if( x < 8.0 )
 270         {
 271         f = polevlf( z, FN4, 6 ) / (x * p1evlf( z, FD4, 7 ));
 272         g = z * polevlf( z, GN4, 7 ) / p1evlf( z, GD4, 7 );
 273         }
 274 else
 275         {
 276         f = polevlf( z, FN8, 8 ) / (x * p1evlf( z, FD8, 8 ));
 277         g = z * polevlf( z, GN8, 8 ) / p1evlf( z, GD8, 9 );
 278         }
 279 *si = PIO2F - f * c - g * s;
 280 if( sign )
 281         *si = -( *si );
 282 *ci = f * s - g * c;
 283
 284 return(0);
 285 }