SiPixelRecHits/src/sicif.h

0001 /*                          sicif.c
0002  *
0003  *  Sine and cosine integrals
0004  *
0005  *
0006  *
0007  * SYNOPSIS:
0008  *
0009  * float x, Ci, Si;
0010  *
0011  * sicif( x, &Si, &Ci );
0012  *
0013  *
0014  * DESCRIPTION:
0015  *
0016  * Evaluates the integrals
0017  *
0018  *                          x
0019  *                          -
0020  *                         |  cos t - 1
0021  *   Ci(x) = eul + ln x +  |  --------- dt,
0022  *                         |      t
0023  *                        -
0024  *                         0
0025  *             x
0026  *             -
0027  *            |  sin t
0028  *   Si(x) =  |  ----- dt
0029  *            |    t
0030  *           -
0031  *            0
0032  *
0033  * where eul = 0.57721566490153286061 is Euler's constant.
0034  * The integrals are approximated by rational functions.
0035  * For x > 8 auxiliary functions f(x) and g(x) are employed
0036  * such that
0037  *
0038  * Ci(x) = f(x) sin(x) - g(x) cos(x)
0039  * Si(x) = pi/2 - f(x) cos(x) - g(x) sin(x)
0040  *
0041  *
0042  * ACCURACY:
0043  *    Test interval = [0,50].
0044  * Absolute error, except relative when > 1:
0045  * arithmetic   function   # trials      peak         rms
0046  *    IEEE        Si        30000       2.1e-7      4.3e-8
0047  *    IEEE        Ci        30000       3.9e-7      2.2e-8
0048  */
0049
0050 /*
0051  Cephes Math Library Release 2.1:  January, 1989
0052  Copyright 1984, 1987, 1989 by Stephen L. Moshier
0053  Direct inquiries to 30 Frost Street, Cambridge, MA 02140
0054  */
0055
0056 #include "vdt/vdtMath.h"
0057
0058 static const float SN[] = {
0059     -8.39167827910303881427E-11,
0060     4.62591714427012837309E-8,
0061     -9.75759303843632795789E-6,
0062     9.76945438170435310816E-4,
0063     -4.13470316229406538752E-2,
0064     1.00000000000000000302E0,
0065 };
0066 static const float SD[] = {
0067     2.03269266195951942049E-12,
0068     1.27997891179943299903E-9,
0069     4.41827842801218905784E-7,
0070     9.96412122043875552487E-5,
0071     1.42085239326149893930E-2,
0072     9.99999999999999996984E-1,
0073 };
0074
0075 static const float CN[] = {
0076     2.02524002389102268789E-11,
0077     -1.35249504915790756375E-8,
0078     3.59325051419993077021E-6,
0079     -4.74007206873407909465E-4,
0080     2.89159652607555242092E-2,
0081     -1.00000000000000000080E0,
0082 };
0083 static const float CD[] = {
0084     4.07746040061880559506E-12,
0085     3.06780997581887812692E-9,
0086     1.23210355685883423679E-6,
0087     3.17442024775032769882E-4,
0088     5.10028056236446052392E-2,
0089     4.00000000000000000080E0,
0090 };
0091
0092 static const float FN4[] = {
0093     4.23612862892216586994E0,
0094     5.45937717161812843388E0,
0095     1.62083287701538329132E0,
0096     1.67006611831323023771E-1,
0097     6.81020132472518137426E-3,
0098     1.08936580650328664411E-4,
0099     5.48900223421373614008E-7,
0100 };
0101 static const float FD4[] = {
0102     /*  1.00000000000000000000E0,*/
0103     8.16496634205391016773E0,
0104     7.30828822505564552187E0,
0105     1.86792257950184183883E0,
0106     1.78792052963149907262E-1,
0107     7.01710668322789753610E-3,
0108     1.10034357153915731354E-4,
0109     5.48900252756255700982E-7,
0110 };
0111
0112 static const float FN8[] = {
0113     4.55880873470465315206E-1,
0114     7.13715274100146711374E-1,
0115     1.60300158222319456320E-1,
0116     1.16064229408124407915E-2,
0117     3.49556442447859055605E-4,
0118     4.86215430826454749482E-6,
0119     3.20092790091004902806E-8,
0120     9.41779576128512936592E-11,
0121     9.70507110881952024631E-14,
0122 };
0123 static const float FD8[] = {
0124     /*  1.00000000000000000000E0,*/
0125     9.17463611873684053703E-1,
0126     1.78685545332074536321E-1,
0127     1.22253594771971293032E-2,
0128     3.58696481881851580297E-4,
0129     4.92435064317881464393E-6,
0130     3.21956939101046018377E-8,
0131     9.43720590350276732376E-11,
0132     9.70507110881952025725E-14,
0133 };
0134
0135 static const float GN4[] = {
0136     8.71001698973114191777E-2,
0137     6.11379109952219284151E-1,
0138     3.97180296392337498885E-1,
0139     7.48527737628469092119E-2,
0140     5.38868681462177273157E-3,
0141     1.61999794598934024525E-4,
0142     1.97963874140963632189E-6,
0143     7.82579040744090311069E-9,
0144 };
0145 static const float GD4[] = {
0146     /*  1.00000000000000000000E0,*/
0147     1.64402202413355338886E0,
0148     6.66296701268987968381E-1,
0149     9.88771761277688796203E-2,
0150     6.22396345441768420760E-3,
0151     1.73221081474177119497E-4,
0152     2.02659182086343991969E-6,
0153     7.82579218933534490868E-9,
0154 };
0155
0156 static const float GN8[] = {
0157     6.97359953443276214934E-1,
0158     3.30410979305632063225E-1,
0159     3.84878767649974295920E-2,
0160     1.71718239052347903558E-3,
0161     3.48941165502279436777E-5,
0162     3.47131167084116673800E-7,
0163     1.70404452782044526189E-9,
0164     3.85945925430276600453E-12,
0165     3.14040098946363334640E-15,
0166 };
0167 static const float GD8[] = {
0168     /*  1.00000000000000000000E0,*/
0169     1.68548898811011640017E0,
0170     4.87852258695304967486E-1,
0171     4.67913194259625806320E-2,
0172     1.90284426674399523638E-3,
0173     3.68475504442561108162E-5,
0174     3.57043223443740838771E-7,
0175     1.72693748966316146736E-9,
0176     3.87830166023954706752E-12,
0177     3.14040098946363335242E-15,
0178 };
0179
0180 inline float polevlf(float xx, const float *coef, int N) {
0181   float ans, x;
0182   const float *p;
0183   int i;
0184
0185   x = xx;
0186   p = coef;
0187   ans = *p++;
0188
0189   i = N;
0190   do
0191     ans = ans * x + *p++;
0192   while (--i);
0193
0194   return (ans);
0195 }
0196
0197 /*                          p1evl() */
0198 /*                                          N
0199  * Evaluate polynomial when coefficient of x  is 1.0.
0200  * Otherwise same as polevl.
0201  */
0202 inline float p1evlf(float xx, const float *coef, int N) {
0203   float ans, x;
0204   const float *p;
0205   int i;
0206
0207   x = xx;
0208   p = coef;
0209   ans = x + *p++;
0210   i = N - 1;
0211
0212   do
0213     ans = ans * x + *p++;
0214   while (--i);
0215
0216   return (ans);
0217 }
0218
0219 inline int sicif(float xx, float &si, float &ci) {
0220   const float MAXNUMF = 1.7014117331926442990585209174225846272e38;
0221   const float PIO2F = 1.5707963267948966192;
0222   // const float MACHEPF = 5.9604644775390625E-8;
0223   const float EUL = 0.57721566490153286061;
0224
0225   float x, z, c, s, f, g;
0226   int sign;
0227
0228   x = xx;
0229   if (x < 0.0f) {
0230     sign = -1;
0231     x = -x;
0232   } else
0233     sign = 0;
0234
0235   if (x == 0.0f) {
0236     si = 0.0;
0237     ci = -MAXNUMF;
0238     return (0);
0239   }
0240
0241   if (x > 1.0e9f) {
0242     float su, cu;
0243     vdt::fast_sincosf(x, su, cu);
0244     si = PIO2F - cu / x;
0245     ci = su / x;
0246     return (0);
0247   }
0248
0249   if (x > 4.0f)
0250     goto asympt;
0251
0252   z = x * x;
0253   s = x * polevlf(z, SN, 5) / polevlf(z, SD, 5);
0254   c = z * polevlf(z, CN, 5) / polevlf(z, CD, 5);
0255
0256   if (sign)
0257     s = -s;
0258   si = s;
0259   ci = EUL + vdt::fast_logf(x) + c; /* real part if x < 0 */
0260   return (0);
0261
0262   /* The auxiliary functions are:
0263     *
0264     *
0265     * *si = *si - PIO2;
0266     * c = cos(x);
0267     * s = sin(x);
0268     *
0269     * t = *ci * s - *si * c;
0270     * a = *ci * c + *si * s;
0271     *
0272     * *si = t;
0273     * *ci = -a;
0274     */
0275
0276 asympt:
0277   vdt::fast_sincosf(x, s, c);
0278   z = 1.0f / (x * x);
0279   if (x < 8.0f) {
0280     f = polevlf(z, FN4, 6) / (x * p1evlf(z, FD4, 7));
0281     g = z * polevlf(z, GN4, 7) / p1evlf(z, GD4, 7);
0282   } else {
0283     f = polevlf(z, FN8, 8) / (x * p1evlf(z, FD8, 8));
0284     g = z * polevlf(z, GN8, 8) / p1evlf(z, GD8, 9);
0285   }
0286   si = PIO2F - f * c - g * s;
0287   if (sign)
0288     si = -(si);
0289   ci = f * s - g * c;
0290
0291   return (0);
0292 }