Private version of the cosine and sine integral.
202 const float zero = 0.;
203 const float q2[7] = { .10340013040487,3.319092135933,
204 20.449478501379,41.280784189142,32.426421069514,10.041164382905,
206 const float p3[6] = { -2.3909964453136,-147.98219500504,
207 -254.3763397689,-119.55761038372,-19.630408535939,-.9999999999036
209 const float q3[6] = { 177.60070940351,530.68509610812,
210 462.23027156148,156.81843364539,21.630408494238,1. };
211 const float p4[8] = { -8.6693733995107,-549.14226552109,
212 -4210.0161535707,-249301.39345865,-119623.66934925,
213 -22174462.775885,3892804.213112,-391546073.8091 };
214 const float q4[8] = { 34.171875,-1607.0892658722,35730.029805851,
215 -483547.43616216,4285596.2461175,-24903337.574054,89192576.757561,
217 const float a1[8] = { -2.1808638152072,-21.901023385488,
218 9.3081638566217,25.076281129356,-33.184253199722,60.121799083008,
219 -43.253113287813,1.0044310922808 };
220 const float b1[8] = { 0.,3.9370770185272,300.89264837292,
221 -6.2504116167188,1003.6743951673,14.325673812194,2736.2411988933,
223 const float a2[8] = { -3.4833465360285,-18.65454548834,
224 -8.2856199414064,-32.34673303054,17.960168876925,1.7565631546961,
225 -1.9502232128966,.99999429607471 };
226 const float b2[8] = { 0.,69.500065588743,57.283719383732,
227 25.777638423844,760.76114800773,28.951672792514,-3.4394226689987,
229 const float a3[6] = { -27.780928934438,-10.10479081576,
230 -9.1483008216736,-5.0223317461851,-3.0000077799358,
232 const float one = 1.;
233 const float b3[6] = { 0.,122.39993926823,2.7276100778779,
234 -7.1897518395045,-2.9990118065262,1.999999942826 };
235 const float two = 2.;
236 const float three = 3.;
237 const float x0 = .37250741078137;
238 const float xl[6] = { -24.,-12.,-6.,0.,1.,4. };
239 const float p1[5] = { 4.293125234321,39.894153870321,
240 292.52518866921,425.69682638592,-434.98143832952 };
241 const float q1[5] = { 1.,18.899288395003,150.95038744251,
242 568.05252718987,753.58564359843 };
243 const float p2[7] = { .43096783946939,6.9052252278444,
244 23.019255939133,24.378408879132,9.0416155694633,.99997957705159,
248 float v,
y, ap, bp, aq,
dp, bq, dq;
252 for (
int i__ = 2; i__ <= 5; ++i__) {
253 ap = a3[i__ - 1] - x + b3[i__ - 1] / ap;
255 y = vdt::fast_expf(-x) / x * (one - (a3[5] + b3[5] / ap) / x);
256 }
else if (
x <= xl[1]) {
258 for (
int i__ = 2; i__ <= 7; ++i__) {
259 ap = a2[i__ - 1] - x + b2[i__ - 1] / ap;
261 y = vdt::fast_expf(-x) / x * (a2[7] + b2[7] / ap);
262 }
else if (x <= xl[2]) {
264 for (
int i__ = 2; i__ <= 7; ++i__) {
265 ap = a1[i__ - 1] - x + b1[i__ - 1] / ap;
267 y = vdt::fast_expf(-x) / x * (a1[7] + b1[7] / ap);
268 }
else if (x < xl[3]) {
269 v = -two * (x / three + one);
272 for (
int i__ = 2; i__ <= 8; ++i__) {
275 dp = p4[i__ - 1] - ap + v * bp;
279 for (
int i__ = 2; i__ <= 8; ++i__) {
282 dq = q4[i__ - 1] - aq + v * bq;
284 y = -vdt::fast_logf(-x / x0) + (x + x0) * (dp - ap) / (dq - aq);
285 }
else if (x == xl[3]) {
287 }
else if (x < xl[4]) {
290 for (
int i__ = 2; i__ <= 5; ++i__) {
291 ap = p1[i__ - 1] + x * ap;
292 aq = q1[i__ - 1] + x * aq;
294 y = -vdt::fast_logf(x) + ap / aq;
295 }
else if (x <= xl[5]) {
299 for (
int i__ = 2; i__ <= 7; ++i__) {
300 ap = p2[i__ - 1] + y * ap;
301 aq = q2[i__ - 1] + y * aq;
303 y = vdt::fast_expf(-x) * ap / aq;
308 for (
int i__ = 2; i__ <= 6; ++i__) {
309 ap = p3[i__ - 1] + y * ap;
310 aq = q3[i__ - 1] + y * aq;
312 y = vdt::fast_expf(-x) * y * (one + y * ap / aq);