1 /*
2 * Rys quadrature
3
4 Code is edited based on
5 * PyQuante quantum chemistry program suite http://pyquante.sourceforge.net
6 * BDF program package
7
8 Coefficients of Chebyshev polynomials are generated based on Toru Shiozaki's
9 integral project libslater
10 */
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <string.h>
15 #include <float.h>
16 #include <math.h>
17 #include "config.h"
18 #include "rys_roots.h"
19 #include "roots_for_x0.dat"
20
21 #ifdef HAVE_QUADMATH_H
22 #include <quadmath.h>
23 #endif
24 #define PIE4 0.78539816339744827900
25 #define THRESHOLD_ZERO (DBL_EPSILON * 8)
26
27 static void rys_root1(double x, double *roots, double *weights);
28 static void rys_root2(double x, double *roots, double *weights);
29 static void rys_root3(double x, double *roots, double *weights);
30 static void rys_root4(double x, double *roots, double *weights);
31 static void rys_root5(double x, double *roots, double *weights);
32 typedef int QuadratureFunction(int n, double x, double lower, double *roots, double *weights);
33 #ifndef HAVE_QUADMATH_H
34 #define CINTqrys_schmidt CINTlrys_schmidt
35 #define CINTqrys_laguerre CINTlrys_laguerre
36 #define CINTqrys_jacobi CINTlrys_jacobi
37 #endif
38
segment_solve(int n,double x,double lower,double * u,double * w,double breakpoint,QuadratureFunction fn1,QuadratureFunction fn2)39 static int segment_solve(int n, double x, double lower, double *u, double *w,
40 double breakpoint, QuadratureFunction fn1, QuadratureFunction fn2)
41 {
42 int error;
43 if (x <= breakpoint) {
44 error = fn1(n, x, lower, u, w);
45 } else {
46 error = fn2(n, x, lower, u, w);
47 }
48 if (error) {
49 error = CINTqrys_schmidt(n, x, lower, u, w);
50 }
51 if (error) {
52 fprintf(stderr, "libcint rys_roots failed. nroots=%d\n", n);
53 #ifndef KEEP_GOING
54 exit(1);
55 #endif
56 }
57 return error;
58 }
59
CINTrys_roots(int nroots,double x,double * u,double * w)60 void CINTrys_roots(int nroots, double x, double *u, double *w)
61 {
62 if (x < THRESHOLD_ZERO) {
63 int off = nroots * (nroots - 1) / 2;
64 int i;
65 for (i = 0; i < nroots; i++) {
66 u[i] = ROOTS_FOR_X0[off + i];
67 w[i] = WEIGHTS_FOR_X0[off + i];
68 }
69 return;
70 }
71
72 switch (nroots) {
73 case 1:
74 rys_root1(x, u, w);
75 break;
76 case 2:
77 rys_root2(x, u, w);
78 break;
79 case 3:
80 rys_root3(x, u, w);
81 break;
82 case 4:
83 rys_root4(x, u, w);
84 break;
85 case 5:
86 rys_root5(x, u, w);
87 break;
88 case 6: case 7:
89 segment_solve(nroots, x, 0., u, w, 11, CINTrys_jacobi, CINTrys_schmidt);
90 break;
91 case 8:
92 segment_solve(nroots, x, 0., u, w, 11, CINTrys_jacobi, CINTlrys_schmidt);
93 break;
94 case 9:
95 segment_solve(nroots, x, 0., u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre);
96 break;
97 case 10: case 11:
98 segment_solve(nroots, x, 0., u, w, 18, CINTlrys_jacobi, CINTlrys_laguerre);
99 break;
100 case 12:
101 segment_solve(nroots, x, 0., u, w, 22, CINTlrys_jacobi, CINTlrys_laguerre);
102 break;
103 default:
104 segment_solve(nroots, x, 0., u, w, 50, CINTqrys_jacobi, CINTqrys_laguerre);
105 break;
106 }
107 }
108
109 #ifdef WITH_RANGE_COULOMB
110 /*
111 * lower is the lower bound of the sr integral
112 */
segment_solve1(int n,double x,double lower,double * u,double * w,double lower_bp1,double lower_bp2,double breakpoint,QuadratureFunction fn1,QuadratureFunction fn2,QuadratureFunction fn3)113 static int segment_solve1(int n, double x, double lower, double *u, double *w,
114 double lower_bp1, double lower_bp2, double breakpoint,
115 QuadratureFunction fn1, QuadratureFunction fn2, QuadratureFunction fn3)
116 {
117 int error;
118 if (lower < lower_bp1) {
119 if (x <= breakpoint) {
120 error = fn1(n, x, lower, u, w);
121 } else {
122 error = fn2(n, x, lower, u, w);
123 }
124 } else if (lower < lower_bp2) {
125 error = fn3(n, x, lower, u, w);
126 } else {
127 fprintf(stderr, "libcint SR-rys_roots does not support nroots=%d x=%g lower=%g\n",
128 n, x, lower);
129 #ifndef KEEP_GOING
130 exit(1);
131 #endif
132 return 1;
133 }
134 if (error) {
135 error = CINTqrys_schmidt(n, x, lower, u, w);
136 if (error) {
137 fprintf(stderr, "libcint SR-rys_roots failed. nroots=%d\n", n);
138 #ifndef KEEP_GOING
139 exit(1);
140 #endif
141 }
142 }
143 return error;
144 }
145
CINTsr_rys_roots(int nroots,double x,double lower,double * u,double * w)146 void CINTsr_rys_roots(int nroots, double x, double lower, double *u, double *w)
147 {
148 switch (nroots) {
149 case 1:
150 CINTrys_schmidt(nroots, x, lower, u, w);
151 break;
152 case 2:
153 if (lower < 0.99) {
154 CINTrys_schmidt(nroots, x, lower, u, w);
155 } else {
156 CINTqrys_jacobi(nroots, x, lower, u, w);
157 }
158 break;
159 case 3:
160 if (lower < 0.93) {
161 CINTrys_schmidt(nroots, x, lower, u, w);
162 } else if (lower < 0.97) {
163 segment_solve(nroots, x, lower, u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre);
164 } else {
165 CINTqrys_jacobi(nroots, x, lower, u, w);
166 }
167 break;
168 case 4:
169 if (lower < 0.85) {
170 CINTrys_schmidt(nroots, x, lower, u, w);
171 } else if (lower < 0.9) {
172 segment_solve(nroots, x, lower, u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre);
173 } else {
174 CINTqrys_jacobi(nroots, x, lower, u, w);
175 }
176 break;
177 case 5:
178 if (lower < 0.45) {
179 CINTrys_schmidt(nroots, x, lower, u, w);
180 } else if (lower < 0.8) {
181 segment_solve(nroots, x, lower, u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre);
182 } else {
183 CINTqrys_jacobi(nroots, x, lower, u, w);
184 }
185 break;
186 case 6:
187 if (lower < 0.35) {
188 CINTrys_schmidt(nroots, x, lower, u, w);
189 } else if (lower < 0.8) {
190 segment_solve(nroots, x, lower, u, w, 10, CINTlrys_jacobi, CINTlrys_laguerre);
191 } else {
192 CINTqrys_jacobi(nroots, x, lower, u, w);
193 }
194 break;
195 case 7:
196 segment_solve1(nroots, x, lower, u, w, 0.55, 1., 60, CINTlrys_jacobi, CINTlrys_laguerre, CINTqrys_jacobi);
197 break;
198 case 8: case 9: case 10:
199 //CINTqrys_jacobi(nroots, x, lower, u, w);
200 segment_solve1(nroots, x, lower, u, w, 0.15, 1., 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi);
201 break;
202 case 11: case 12:
203 segment_solve1(nroots, x, lower, u, w, 0.15, 1., 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi);
204 break;
205 case 13: case 14:
206 segment_solve1(nroots, x, lower, u, w, 0.25, 1., 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi);
207 break;
208 case 15: case 16:
209 segment_solve1(nroots, x, lower, u, w, 0.25, 0.75, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi);
210 break;
211 case 17:
212 segment_solve1(nroots, x, lower, u, w, 0.25, 0.65, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi);
213 break;
214 case 18:
215 segment_solve1(nroots, x, lower, u, w, 0.15, 0.65, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi);
216 break;
217 case 19:
218 segment_solve1(nroots, x, lower, u, w, 0.15, 0.55, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi);
219 break;
220 case 20: case 21:
221 segment_solve1(nroots, x, lower, u, w, 0.25, 0.45, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi);
222 break;
223 case 22: case 23: case 24:
224 segment_solve1(nroots, x, lower, u, w, 0.25, 0.35, 60, CINTqrys_jacobi, CINTqrys_laguerre, CINTqrys_jacobi);
225 break;
226 default:
227 fprintf(stderr, "libcint SR-rys_roots does not support nroots=%d\n", nroots);
228 #ifndef KEEP_GOING
229 exit(1);
230 #endif
231 }
232 }
233 #endif
234
rys_root1(double X,double * roots,double * weights)235 static void rys_root1(double X, double *roots, double *weights)
236 {
237 double Y, F1;
238
239 if (X > 33.) {
240 weights[0] = sqrt(PIE4/X);
241 roots[0] = 0.5E+00/(X-0.5E+00);
242 return;
243 } else if (X < 3.e-7) {
244 weights[0] = 1.0E+00 -X/3.0E+00;
245 roots[0] = 0.5E+00 -X/5.0E+00;
246 return;
247 }
248
249 double E = exp(-X);
250 if (X > 15.) {
251 Y = 1./X;
252 F1 = ((( 1.9623264149430E-01*Y-4.9695241464490E-01)*Y -
253 6.0156581186481E-05)* E + sqrt(PIE4/X) - E)*Y;
254 F1 *= .5;
255 } else if (X > 10.) {
256 Y = 1./X;
257 F1 = ((((-1.8784686463512E-01*Y+2.2991849164985E-01)*Y -
258 4.9893752514047E-01)*Y-2.1916512131607E-05)* E
259 + sqrt(PIE4/X) - E)*Y;
260 F1 *= .5;
261 } else if (X > 5.) {
262 Y = 1./X;
263 F1 = ((((((( 4.6897511375022E-01*Y-6.9955602298985E-01)*Y +
264 5.3689283271887E-01)*Y-3.2883030418398E-01)*Y +
265 2.4645596956002E-01)*Y-4.9984072848436E-01)*Y -
266 3.1501078774085E-06)* E + sqrt(PIE4/X) - E)*Y;
267 F1 *= .5;
268 } else if (X > 3.){
269 Y = X-4.0E+00;
270 F1 = ((((((((((-2.62453564772299E-11*Y+3.24031041623823E-10 )*Y-
271 3.614965656163E-09)*Y+3.760256799971E-08)*Y-
272 3.553558319675E-07)*Y+3.022556449731E-06)*Y-
273 2.290098979647E-05)*Y+1.526537461148E-04)*Y-
274 8.81947375894379E-04 )*Y+4.33207949514611E-03 )*Y-
275 1.75257821619926E-02 )*Y+5.28406320615584E-02;
276 } else if (X > 1.) {
277 Y = X-2.0E+00;
278 F1 = ((((((((((-1.61702782425558E-10*Y+1.96215250865776E-09 )*Y-
279 2.14234468198419E-08 )*Y+2.17216556336318E-07 )*Y-
280 1.98850171329371E-06 )*Y+1.62429321438911E-05 )*Y-
281 1.16740298039895E-04 )*Y+7.24888732052332E-04 )*Y-
282 3.79490003707156E-03 )*Y+1.61723488664661E-02 )*Y-
283 5.29428148329736E-02 )*Y+1.15702180856167E-01;
284 } else {
285 F1 = ((((((((-8.36313918003957E-08*X+1.21222603512827E-06 )*X-
286 1.15662609053481E-05 )*X+9.25197374512647E-05 )*X-
287 6.40994113129432E-04 )*X+3.78787044215009E-03 )*X-
288 1.85185172458485E-02 )*X+7.14285713298222E-02 )*X-
289 1.99999999997023E-01 )*X+3.33333333333318E-01;
290 }
291
292 double WW1 = 2. * X * F1 + E;
293 weights[0] = WW1;
294 roots[0] = F1 / (WW1 - F1);
295 }
296
rys_root2(double X,double * roots,double * weights)297 static void rys_root2(double X, double *roots, double *weights)
298 {
299
300 double R12, R22, W22;
301 double RT1, RT2, WW1, WW2;
302 double F1, E, Y;
303
304 R12 = 2.75255128608411E-01;
305 R22 = 2.72474487139158E+00;
306 W22 = 9.17517095361369E-02;
307
308 if (X < 3.e-7){
309 RT1 = 1.30693606237085E-01 -2.90430236082028E-02 *X;
310 RT2 = 2.86930639376291E+00 -6.37623643058102E-01 *X;
311 WW1 = 6.52145154862545E-01 -1.22713621927067E-01 *X;
312 WW2 = 3.47854845137453E-01 -2.10619711404725E-01 *X;
313 } else if (X < 1.) {
314 F1 = ((((((((-8.36313918003957E-08*X+1.21222603512827E-06 )*X-
315 1.15662609053481E-05 )*X+9.25197374512647E-05 )*X-
316 6.40994113129432E-04 )*X+3.78787044215009E-03 )*X-
317 1.85185172458485E-02 )*X+7.14285713298222E-02 )*X-
318 1.99999999997023E-01 )*X+3.33333333333318E-01;
319 WW1 = (X+X)*F1+exp(-X);
320 RT1 = (((((((-2.35234358048491E-09*X+2.49173650389842E-08)*X-
321 4.558315364581E-08)*X-2.447252174587E-06)*X+
322 4.743292959463E-05)*X-5.33184749432408E-04 )*X+
323 4.44654947116579E-03 )*X-2.90430236084697E-02 )*X+
324 1.30693606237085E-01;
325 RT2 = (((((((-2.47404902329170E-08*X+2.36809910635906E-07)*X+
326 1.835367736310E-06)*X-2.066168802076E-05)*X-
327 1.345693393936E-04)*X-5.88154362858038E-05 )*X+
328 5.32735082098139E-02 )*X-6.37623643056745E-01 )*X+
329 2.86930639376289E+00;
330 WW2 = ((F1-WW1)*RT1+F1)*(1.0E+00+RT2)/(RT2-RT1);
331 WW1 = WW1-WW2;
332 } else if (X < 3.) {
333 Y = X-2.0E+00;
334 F1 = ((((((((((-1.61702782425558E-10*Y+1.96215250865776E-09 )*Y-
335 2.14234468198419E-08 )*Y+2.17216556336318E-07 )*Y-
336 1.98850171329371E-06 )*Y+1.62429321438911E-05 )*Y-
337 1.16740298039895E-04 )*Y+7.24888732052332E-04 )*Y-
338 3.79490003707156E-03 )*Y+1.61723488664661E-02 )*Y-
339 5.29428148329736E-02 )*Y+1.15702180856167E-01;
340 WW1 = (X+X)*F1+exp(-X);
341 RT1 = (((((((((-6.36859636616415E-12*Y+8.47417064776270E-11)*Y-
342 5.152207846962E-10)*Y-3.846389873308E-10)*Y+
343 8.472253388380E-08)*Y-1.85306035634293E-06 )*Y+
344 2.47191693238413E-05 )*Y-2.49018321709815E-04 )*Y+
345 2.19173220020161E-03 )*Y-1.63329339286794E-02 )*Y+
346 8.68085688285261E-02;
347 RT2 = ((((((((( 1.45331350488343E-10*Y+2.07111465297976E-09)*Y-
348 1.878920917404E-08)*Y-1.725838516261E-07)*Y+
349 2.247389642339E-06)*Y+9.76783813082564E-06 )*Y-
350 1.93160765581969E-04 )*Y-1.58064140671893E-03 )*Y+
351 4.85928174507904E-02 )*Y-4.30761584997596E-01 )*Y+
352 1.80400974537950E+00;
353 WW2 = ((F1-WW1)*RT1+F1)*(1.0E+00+RT2)/(RT2-RT1);
354 WW1 = WW1-WW2;
355 } else if (X < 5.){
356 Y = X-4.0E+00;
357 F1 = ((((((((((-2.62453564772299E-11*Y+3.24031041623823E-10 )*Y-
358 3.614965656163E-09)*Y+3.760256799971E-08)*Y-
359 3.553558319675E-07)*Y+3.022556449731E-06)*Y-
360 2.290098979647E-05)*Y+1.526537461148E-04)*Y-
361 8.81947375894379E-04 )*Y+4.33207949514611E-03 )*Y-
362 1.75257821619926E-02 )*Y+5.28406320615584E-02;
363 WW1 = (X+X)*F1+exp(-X);
364 RT1 = ((((((((-4.11560117487296E-12*Y+7.10910223886747E-11)*Y-
365 1.73508862390291E-09 )*Y+5.93066856324744E-08 )*Y-
366 9.76085576741771E-07 )*Y+1.08484384385679E-05 )*Y-
367 1.12608004981982E-04 )*Y+1.16210907653515E-03 )*Y-
368 9.89572595720351E-03 )*Y+6.12589701086408E-02;
369 RT2 = (((((((((-1.80555625241001E-10*Y+5.44072475994123E-10)*Y+
370 1.603498045240E-08)*Y-1.497986283037E-07)*Y-
371 7.017002532106E-07)*Y+1.85882653064034E-05 )*Y-
372 2.04685420150802E-05 )*Y-2.49327728643089E-03 )*Y+
373 3.56550690684281E-02 )*Y-2.60417417692375E-01 )*Y+
374 1.12155283108289E+00;
375 WW2 = ((F1-WW1)*RT1+F1)*(1.0E+00+RT2)/(RT2-RT1);
376 WW1 = WW1-WW2;
377 } else if (X < 10) {
378 E = exp(-X);
379 WW1 = (((((( 4.6897511375022E-01/X-6.9955602298985E-01)/X +
380 5.3689283271887E-01)/X-3.2883030418398E-01)/X +
381 2.4645596956002E-01)/X-4.9984072848436E-01)/X -
382 3.1501078774085E-06)*E + sqrt(PIE4/X);
383 F1 = (WW1-E)/(X+X);
384 Y = X-7.5E+00;
385 RT1 = (((((((((((((-1.43632730148572E-16*Y+2.38198922570405E-16)*
386 Y+1.358319618800E-14)*Y-7.064522786879E-14)*Y-
387 7.719300212748E-13)*Y+7.802544789997E-12)*Y+
388 6.628721099436E-11)*Y-1.775564159743E-09)*Y+
389 1.713828823990E-08)*Y-1.497500187053E-07)*Y+
390 2.283485114279E-06)*Y-3.76953869614706E-05 )*Y+
391 4.74791204651451E-04 )*Y-4.60448960876139E-03 )*Y+
392 3.72458587837249E-02;
393 RT2 = (((((((((((( 2.48791622798900E-14*Y-1.36113510175724E-13)*Y-
394 2.224334349799E-12)*Y+4.190559455515E-11)*Y-
395 2.222722579924E-10)*Y-2.624183464275E-09)*Y+
396 6.128153450169E-08)*Y-4.383376014528E-07)*Y-
397 2.49952200232910E-06 )*Y+1.03236647888320E-04 )*Y-
398 1.44614664924989E-03 )*Y+1.35094294917224E-02 )*Y-
399 9.53478510453887E-02 )*Y+5.44765245686790E-01;
400 WW2 = ((F1-WW1)*RT1+F1)*(1.0E+00+RT2)/(RT2-RT1);
401 WW1 = WW1-WW2;
402 } else if (X < 15) {
403 E = exp(-X);
404 WW1 = (((-1.8784686463512E-01/X+2.2991849164985E-01)/X -
405 4.9893752514047E-01)/X-2.1916512131607E-05)*E
406 + sqrt(PIE4/X);
407 F1 = (WW1-E)/(X+X);
408 RT1 = ((((-1.01041157064226E-05*X+1.19483054115173E-03)*X -
409 6.73760231824074E-02)*X+1.25705571069895E+00)*X +
410 (((-8.57609422987199E+03/X+5.91005939591842E+03)/X -
411 1.70807677109425E+03)/X+2.64536689959503E+02)/X -
412 2.38570496490846E+01)*E + R12/(X-R12);
413 RT2 = ((( 3.39024225137123E-04*X-9.34976436343509E-02)*X -
414 4.22216483306320E+00)*X +
415 (((-2.08457050986847E+03/X -
416 1.04999071905664E+03)/X+3.39891508992661E+02)/X -
417 1.56184800325063E+02)/X+8.00839033297501E+00)*E + R22/(X-R22);
418 WW2 = ((F1-WW1)*RT1+F1)*(1.0E+00+RT2)/(RT2-RT1);
419 WW1 = WW1-WW2;
420 } else if (X < 33) {
421 E = exp(-X);
422 WW1 = (( 1.9623264149430E-01/X-4.9695241464490E-01)/X -
423 6.0156581186481E-05)*E + sqrt(PIE4/X);
424 F1 = (WW1-E)/(X+X);
425 RT1 = ((((-1.14906395546354E-06*X+1.76003409708332E-04)*X -
426 1.71984023644904E-02)*X-1.37292644149838E-01)*X +
427 (-4.75742064274859E+01/X+9.21005186542857E+00)/X -
428 2.31080873898939E-02)*E + R12/(X-R12);
429 RT2 = ((( 3.64921633404158E-04*X-9.71850973831558E-02)*X -
430 4.02886174850252E+00)*X +
431 (-1.35831002139173E+02/X -
432 8.66891724287962E+01)/X+2.98011277766958E+00)*E + R22/(X-R22);
433 WW2 = ((F1-WW1)*RT1+F1)*(1.0E+00+RT2)/(RT2-RT1);
434 WW1 = WW1-WW2;
435 } else if (X < 40) {
436 WW1 = sqrt(PIE4/X);
437 E = exp(-X);
438 RT1 = (-8.78947307498880E-01*X+1.09243702330261E+01)*E + R12/(X-R12);
439 RT2 = (-9.28903924275977E+00*X+8.10642367843811E+01)*E + R22/(X-R22);
440 WW2 = ( 4.46857389308400E+00*X-7.79250653461045E+01)*E + W22*WW1;
441 WW1 = WW1-WW2;
442 } else {
443 WW1 = sqrt(PIE4/X);
444 RT1 = R12/(X-R12);
445 RT2 = R22/(X-R22);
446 WW2 = W22*WW1;
447 WW1 = WW1-WW2;
448 }
449 roots[0] = RT1;
450 roots[1] = RT2;
451 weights[0] = WW1;
452 weights[1] = WW2;
453 }
454
rys_root3(double X,double * roots,double * weights)455 static void rys_root3(double X, double *roots, double *weights)
456 {
457
458 double R13, R23, W23, R33, W33;
459 double RT1, RT2, RT3, WW1, WW2, WW3;
460 double F1, F2, E, T1, T2, T3, A1, A2, Y;
461
462 R13 = 1.90163509193487E-01;
463 R23 = 1.78449274854325E+00;
464 W23 = 1.77231492083829E-01;
465 R33 = 5.52534374226326E+00;
466 W33 = 5.11156880411248E-03;
467
468 if (X < 3.e-7){
469 RT1 = 6.03769246832797E-02 -9.28875764357368E-03 *X;
470 RT2 = 7.76823355931043E-01 -1.19511285527878E-01 *X;
471 RT3 = 6.66279971938567E+00 -1.02504611068957E+00 *X;
472 WW1 = 4.67913934572691E-01 -5.64876917232519E-02 *X;
473 WW2 = 3.60761573048137E-01 -1.49077186455208E-01 *X;
474 WW3 = 1.71324492379169E-01 -1.27768455150979E-01 *X;
475 } else if (X < 1.) {
476 RT1 = ((((((-5.10186691538870E-10*X+2.40134415703450E-08)*X-
477 5.01081057744427E-07 )*X+7.58291285499256E-06 )*X-
478 9.55085533670919E-05 )*X+1.02893039315878E-03 )*X-
479 9.28875764374337E-03 )*X+6.03769246832810E-02;
480 RT2 = ((((((-1.29646524960555E-08*X+7.74602292865683E-08)*X+
481 1.56022811158727E-06 )*X-1.58051990661661E-05 )*X-
482 3.30447806384059E-04 )*X+9.74266885190267E-03 )*X-
483 1.19511285526388E-01 )*X+7.76823355931033E-01;
484 RT3 = ((((((-9.28536484109606E-09*X-3.02786290067014E-07)*X-
485 2.50734477064200E-06 )*X-7.32728109752881E-06 )*X+
486 2.44217481700129E-04 )*X+4.94758452357327E-02 )*X-
487 1.02504611065774E+00 )*X+6.66279971938553E+00;
488 F2 = ((((((((-7.60911486098850E-08*X+1.09552870123182E-06 )*X-
489 1.03463270693454E-05 )*X+8.16324851790106E-05 )*X-
490 5.55526624875562E-04 )*X+3.20512054753924E-03 )*X-
491 1.51515139838540E-02 )*X+5.55555554649585E-02 )*X-
492 1.42857142854412E-01 )*X+1.99999999999986E-01;
493 E = exp(-X);
494 F1 = ((X+X)*F2+E)/3.0E+00;
495 WW1 = (X+X)*F1+E;
496 T1 = RT1/(RT1+1.0E+00);
497 T2 = RT2/(RT2+1.0E+00);
498 T3 = RT3/(RT3+1.0E+00);
499 A2 = F2-T1*F1;
500 A1 = F1-T1*WW1;
501 WW3 = (A2-T2*A1)/((T3-T2)*(T3-T1));
502 WW2 = (T3*A1-A2)/((T3-T2)*(T2-T1));
503 WW1 = WW1-WW2-WW3;
504 } else if (X < 3.) {
505 Y = X-2.0E+00;
506 RT1 = (((((((( 1.44687969563318E-12*Y+4.85300143926755E-12)*Y-
507 6.55098264095516E-10 )*Y+1.56592951656828E-08 )*Y-
508 2.60122498274734E-07 )*Y+3.86118485517386E-06 )*Y-
509 5.13430986707889E-05 )*Y+6.03194524398109E-04 )*Y-
510 6.11219349825090E-03 )*Y+4.52578254679079E-02;
511 RT2 = ((((((( 6.95964248788138E-10*Y-5.35281831445517E-09)*Y-
512 6.745205954533E-08)*Y+1.502366784525E-06)*Y+
513 9.923326947376E-07)*Y-3.89147469249594E-04 )*Y+
514 7.51549330892401E-03 )*Y-8.48778120363400E-02 )*Y+
515 5.73928229597613E-01;
516 RT3 = ((((((((-2.81496588401439E-10*Y+3.61058041895031E-09)*Y+
517 4.53631789436255E-08 )*Y-1.40971837780847E-07 )*Y-
518 6.05865557561067E-06 )*Y-5.15964042227127E-05 )*Y+
519 3.34761560498171E-05 )*Y+5.04871005319119E-02 )*Y-
520 8.24708946991557E-01 )*Y+4.81234667357205E+00;
521 F2 = ((((((((((-1.48044231072140E-10*Y+1.78157031325097E-09 )*Y-
522 1.92514145088973E-08 )*Y+1.92804632038796E-07 )*Y-
523 1.73806555021045E-06 )*Y+1.39195169625425E-05 )*Y-
524 9.74574633246452E-05 )*Y+5.83701488646511E-04 )*Y-
525 2.89955494844975E-03 )*Y+1.13847001113810E-02 )*Y-
526 3.23446977320647E-02 )*Y+5.29428148329709E-02;
527 E = exp(-X);
528 F1 = ((X+X)*F2+E)/3.0E+00;
529 WW1 = (X+X)*F1+E;
530 T1 = RT1/(RT1+1.0E+00);
531 T2 = RT2/(RT2+1.0E+00);
532 T3 = RT3/(RT3+1.0E+00);
533 A2 = F2-T1*F1;
534 A1 = F1-T1*WW1;
535 WW3 = (A2-T2*A1)/((T3-T2)*(T3-T1));
536 WW2 = (T3*A1-A2)/((T3-T2)*(T2-T1));
537 WW1 = WW1-WW2-WW3;
538 } else if (X < 5.){
539 Y = X-4.0E+00;
540 RT1 = ((((((( 1.44265709189601E-11*Y-4.66622033006074E-10)*Y+
541 7.649155832025E-09)*Y-1.229940017368E-07)*Y+
542 2.026002142457E-06)*Y-2.87048671521677E-05 )*Y+
543 3.70326938096287E-04 )*Y-4.21006346373634E-03 )*Y+
544 3.50898470729044E-02;
545 RT2 = ((((((((-2.65526039155651E-11*Y+1.97549041402552E-10)*Y+
546 2.15971131403034E-09 )*Y-7.95045680685193E-08 )*Y+
547 5.15021914287057E-07 )*Y+1.11788717230514E-05 )*Y-
548 3.33739312603632E-04 )*Y+5.30601428208358E-03 )*Y-
549 5.93483267268959E-02 )*Y+4.31180523260239E-01;
550 RT3 = ((((((((-3.92833750584041E-10*Y-4.16423229782280E-09)*Y+
551 4.42413039572867E-08 )*Y+6.40574545989551E-07 )*Y-
552 3.05512456576552E-06 )*Y-1.05296443527943E-04 )*Y-
553 6.14120969315617E-04 )*Y+4.89665802767005E-02 )*Y-
554 6.24498381002855E-01 )*Y+3.36412312243724E+00;
555 F2 = ((((((((((-2.36788772599074E-11*Y+2.89147476459092E-10 )*Y-
556 3.18111322308846E-09 )*Y+3.25336816562485E-08 )*Y-
557 3.00873821471489E-07 )*Y+2.48749160874431E-06 )*Y-
558 1.81353179793672E-05 )*Y+1.14504948737066E-04 )*Y-
559 6.10614987696677E-04 )*Y+2.64584212770942E-03 )*Y-
560 8.66415899015349E-03 )*Y+1.75257821619922E-02;
561 E = exp(-X);
562 F1 = ((X+X)*F2+E)/3.0E+00;
563 WW1 = (X+X)*F1+E;
564 T1 = RT1/(RT1+1.0E+00);
565 T2 = RT2/(RT2+1.0E+00);
566 T3 = RT3/(RT3+1.0E+00);
567 A2 = F2-T1*F1;
568 A1 = F1-T1*WW1;
569 WW3 = (A2-T2*A1)/((T3-T2)*(T3-T1));
570 WW2 = (T3*A1-A2)/((T3-T2)*(T2-T1));
571 WW1 = WW1-WW2-WW3;
572 } else if (X < 10) {
573 E = exp(-X);
574 WW1 = (((((( 4.6897511375022E-01/X-6.9955602298985E-01)/X +
575 5.3689283271887E-01)/X-3.2883030418398E-01)/X +
576 2.4645596956002E-01)/X-4.9984072848436E-01)/X -
577 3.1501078774085E-06)*E + sqrt(PIE4/X);
578 F1 = (WW1-E)/(X+X);
579 F2 = (F1+F1+F1-E)/(X+X);
580 Y = X-7.5E+00;
581 RT1 = ((((((((((( 5.74429401360115E-16*Y+7.11884203790984E-16)*Y-
582 6.736701449826E-14)*Y-6.264613873998E-13)*Y+
583 1.315418927040E-11)*Y-4.23879635610964E-11 )*Y+
584 1.39032379769474E-09 )*Y-4.65449552856856E-08 )*Y+
585 7.34609900170759E-07 )*Y-1.08656008854077E-05 )*Y+
586 1.77930381549953E-04 )*Y-2.39864911618015E-03 )*Y+
587 2.39112249488821E-02;
588 RT2 = ((((((((((( 1.13464096209120E-14*Y+6.99375313934242E-15)*Y-
589 8.595618132088E-13)*Y-5.293620408757E-12)*Y-
590 2.492175211635E-11)*Y+2.73681574882729E-09 )*Y-
591 1.06656985608482E-08 )*Y-4.40252529648056E-07 )*Y+
592 9.68100917793911E-06 )*Y-1.68211091755327E-04 )*Y+
593 2.69443611274173E-03 )*Y-3.23845035189063E-02 )*Y+
594 2.75969447451882E-01;
595 RT3 = (((((((((((( 6.66339416996191E-15*Y+1.84955640200794E-13)*Y-
596 1.985141104444E-12)*Y-2.309293727603E-11)*Y+
597 3.917984522103E-10)*Y+1.663165279876E-09)*Y-
598 6.205591993923E-08)*Y+8.769581622041E-09)*Y+
599 8.97224398620038E-06 )*Y-3.14232666170796E-05 )*Y-
600 1.83917335649633E-03 )*Y+3.51246831672571E-02 )*Y-
601 3.22335051270860E-01 )*Y+1.73582831755430E+00;
602 T1 = RT1/(RT1+1.0E+00);
603 T2 = RT2/(RT2+1.0E+00);
604 T3 = RT3/(RT3+1.0E+00);
605 A2 = F2-T1*F1;
606 A1 = F1-T1*WW1;
607 WW3 = (A2-T2*A1)/((T3-T2)*(T3-T1));
608 WW2 = (T3*A1-A2)/((T3-T2)*(T2-T1));
609 WW1 = WW1-WW2-WW3;
610 } else if (X < 15) {
611 E = exp(-X);
612 WW1 = (((-1.8784686463512E-01/X+2.2991849164985E-01)/X -
613 4.9893752514047E-01)/X-2.1916512131607E-05)*E
614 + sqrt(PIE4/X);
615 F1 = (WW1-E)/(X+X);
616 F2 = (F1+F1+F1-E)/(X+X);
617 Y = X-12.5E+00;
618 RT1 = ((((((((((( 4.42133001283090E-16*Y-2.77189767070441E-15)*Y-
619 4.084026087887E-14)*Y+5.379885121517E-13)*Y+
620 1.882093066702E-12)*Y-8.67286219861085E-11 )*Y+
621 7.11372337079797E-10 )*Y-3.55578027040563E-09 )*Y+
622 1.29454702851936E-07 )*Y-4.14222202791434E-06 )*Y+
623 8.04427643593792E-05 )*Y-1.18587782909876E-03 )*Y+
624 1.53435577063174E-02;
625 RT2 = ((((((((((( 6.85146742119357E-15*Y-1.08257654410279E-14)*Y-
626 8.579165965128E-13)*Y+6.642452485783E-12)*Y+
627 4.798806828724E-11)*Y-1.13413908163831E-09 )*Y+
628 7.08558457182751E-09 )*Y-5.59678576054633E-08 )*Y+
629 2.51020389884249E-06 )*Y-6.63678914608681E-05 )*Y+
630 1.11888323089714E-03 )*Y-1.45361636398178E-02 )*Y+
631 1.65077877454402E-01;
632 RT3 = (((((((((((( 3.20622388697743E-15*Y-2.73458804864628E-14)*Y-
633 3.157134329361E-13)*Y+8.654129268056E-12)*Y-
634 5.625235879301E-11)*Y-7.718080513708E-10)*Y+
635 2.064664199164E-08)*Y-1.567725007761E-07)*Y-
636 1.57938204115055E-06 )*Y+6.27436306915967E-05 )*Y-
637 1.01308723606946E-03 )*Y+1.13901881430697E-02 )*Y-
638 1.01449652899450E-01 )*Y+7.77203937334739E-01;
639 T1 = RT1/(RT1+1.0E+00);
640 T2 = RT2/(RT2+1.0E+00);
641 T3 = RT3/(RT3+1.0E+00);
642 A2 = F2-T1*F1;
643 A1 = F1-T1*WW1;
644 WW3 = (A2-T2*A1)/((T3-T2)*(T3-T1));
645 WW2 = (T3*A1-A2)/((T3-T2)*(T2-T1));
646 WW1 = WW1-WW2-WW3;
647 } else if (X < 33) {
648 E = exp(-X);
649 WW1 = (( 1.9623264149430E-01/X-4.9695241464490E-01)/X -
650 6.0156581186481E-05)*E + sqrt(PIE4/X);
651 F1 = (WW1-E)/(X+X);
652 F2 = (F1+F1+F1-E)/(X+X);
653 if (X < 20) {
654 RT1 = ((((((-2.43270989903742E-06*X+3.57901398988359E-04)*X -
655 2.34112415981143E-02)*X+7.81425144913975E-01)*X -
656 1.73209218219175E+01)*X+2.43517435690398E+02)*X +
657 (-1.97611541576986E+04/X+9.82441363463929E+03)/X -
658 2.07970687843258E+03)*E + R13/(X-R13);
659 RT2 = (((((-2.62627010965435E-04*X+3.49187925428138E-02)*X -
660 3.09337618731880E+00)*X+1.07037141010778E+02)*X -
661 2.36659637247087E+03)*X +
662 ((-2.91669113681020E+06/X +
663 1.41129505262758E+06)/X-2.91532335433779E+05)/X +
664 3.35202872835409E+04)*E + R23/(X-R23);
665 RT3 = ((((( 9.31856404738601E-05*X-2.87029400759565E-02)*X -
666 7.83503697918455E-01)*X-1.84338896480695E+01)*X +
667 4.04996712650414E+02)*X +
668 (-1.89829509315154E+05/X +
669 5.11498390849158E+04)/X-6.88145821789955E+03)*E
670 + R33/(X-R33);
671 } else {
672 RT1 = ((((-4.97561537069643E-04*X-5.00929599665316E-02)*X +
673 1.31099142238996E+00)*X-1.88336409225481E+01)*X -
674 6.60344754467191E+02 /X+1.64931462413877E+02)*E
675 + R13/(X-R13);
676 RT2 = ((((-4.48218898474906E-03*X-5.17373211334924E-01)*X +
677 1.13691058739678E+01)*X-1.65426392885291E+02)*X -
678 6.30909125686731E+03 /X+1.52231757709236E+03)*E
679 + R23/(X-R23);
680 RT3 = ((((-1.38368602394293E-02*X-1.77293428863008E+00)*X +
681 1.73639054044562E+01)*X-3.57615122086961E+02)*X -
682 1.45734701095912E+04 /X+2.69831813951849E+03)*E
683 + R33/(X-R33);
684 }
685 T1 = RT1/(RT1+1.0E+00);
686 T2 = RT2/(RT2+1.0E+00);
687 T3 = RT3/(RT3+1.0E+00);
688 A2 = F2-T1*F1;
689 A1 = F1-T1*WW1;
690 WW3 = (A2-T2*A1)/((T3-T2)*(T3-T1));
691 WW2 = (T3*A1-A2)/((T3-T2)*(T2-T1));
692 WW1 = WW1-WW2-WW3;
693 } else if (X < 47) {
694 WW1 = sqrt(PIE4/X);
695 E = exp(-X);
696 RT1 = ((-7.39058467995275E+00*X+3.21318352526305E+02)*X -
697 3.99433696473658E+03)*E + R13/(X-R13);
698 RT2 = ((-7.38726243906513E+01*X+3.13569966333873E+03)*X -
699 3.86862867311321E+04)*E + R23/(X-R23);
700 RT3 = ((-2.63750565461336E+02*X+1.04412168692352E+04)*X -
701 1.28094577915394E+05)*E + R33/(X-R33);
702 WW3 = ((( 1.52258947224714E-01*X-8.30661900042651E+00)*X +
703 1.92977367967984E+02)*X-1.67787926005344E+03)*E
704 + W33*WW1;
705 WW2 = (( 6.15072615497811E+01*X-2.91980647450269E+03)*X +
706 3.80794303087338E+04)*E + W23*WW1;
707 WW1 = WW1-WW2-WW3;
708 } else {
709 WW1 = sqrt(PIE4/X);
710 RT1 = R13/(X-R13);
711 RT2 = R23/(X-R23);
712 RT3 = R33/(X-R33);
713 WW2 = W23*WW1;
714 WW3 = W33*WW1;
715 WW1 = WW1-WW2-WW3;
716 }
717 roots[0] = RT1;
718 roots[1] = RT2;
719 roots[2] = RT3;
720 weights[0] = WW1;
721 weights[1] = WW2;
722 weights[2] = WW3;
723 }
724
rys_root4(double X,double * roots,double * weights)725 static void rys_root4(double X, double *roots, double *weights)
726 {
727 double R14, R24, W24, R34, W34, R44, W44;
728 double RT1, RT2, RT3, RT4, WW1, WW2, WW3, WW4;
729 double Y, E;
730
731 R14 = 1.45303521503316E-01;
732 R24 = 1.33909728812636E+00;
733 W24 = 2.34479815323517E-01;
734 R34 = 3.92696350135829E+00;
735 W34 = 1.92704402415764E-02;
736 R44 = 8.58863568901199E+00;
737 W44 = 2.25229076750736E-04;
738
739 if (X <= 3.0E-7) {
740 RT1 = 3.48198973061471E-02 -4.09645850660395E-03 *X;
741 RT2 = 3.81567185080042E-01 -4.48902570656719E-02 *X;
742 RT3 = 1.73730726945891E+00 -2.04389090547327E-01 *X;
743 RT4 = 1.18463056481549E+01 -1.39368301742312E+00 *X;
744 WW1 = 3.62683783378362E-01 -3.13844305713928E-02 *X;
745 WW2 = 3.13706645877886E-01 -8.98046242557724E-02 *X;
746 WW3 = 2.22381034453372E-01 -1.29314370958973E-01 *X;
747 WW4 = 1.01228536290376E-01 -8.28299075414321E-02 *X;
748 } else if (X <= 1.0) {
749 RT1 = ((((((-1.95309614628539E-10*X+5.19765728707592E-09)*X-
750 1.01756452250573E-07 )*X+1.72365935872131E-06 )*X-
751 2.61203523522184E-05 )*X+3.52921308769880E-04 )*X-
752 4.09645850658433E-03 )*X+3.48198973061469E-02;
753 RT2 = (((((-1.89554881382342E-08*X+3.07583114342365E-07)*X+
754 1.270981734393E-06)*X-1.417298563884E-04)*X+
755 3.226979163176E-03)*X-4.48902570678178E-02 )*X+
756 3.81567185080039E-01;
757 RT3 = (((((( 1.77280535300416E-09*X+3.36524958870615E-08)*X-
758 2.58341529013893E-07 )*X-1.13644895662320E-05 )*X-
759 7.91549618884063E-05 )*X+1.03825827346828E-02 )*X-
760 2.04389090525137E-01 )*X+1.73730726945889E+00;
761 RT4 = (((((-5.61188882415248E-08*X-2.49480733072460E-07)*X+
762 3.428685057114E-06)*X+1.679007454539E-04)*X+
763 4.722855585715E-02)*X-1.39368301737828E+00 )*X+
764 1.18463056481543E+01;
765 WW1 = ((((((-1.14649303201279E-08*X+1.88015570196787E-07)*X-
766 2.33305875372323E-06 )*X+2.68880044371597E-05 )*X-
767 2.94268428977387E-04 )*X+3.06548909776613E-03 )*X-
768 3.13844305680096E-02 )*X+3.62683783378335E-01;
769 WW2 = ((((((((-4.11720483772634E-09*X+6.54963481852134E-08)*X-
770 7.20045285129626E-07 )*X+6.93779646721723E-06 )*X-
771 6.05367572016373E-05 )*X+4.74241566251899E-04 )*X-
772 3.26956188125316E-03 )*X+1.91883866626681E-02 )*X-
773 8.98046242565811E-02 )*X+3.13706645877886E-01;
774 WW3 = ((((((((-3.41688436990215E-08*X+5.07238960340773E-07)*X-
775 5.01675628408220E-06 )*X+4.20363420922845E-05 )*X-
776 3.08040221166823E-04 )*X+1.94431864731239E-03 )*X-
777 1.02477820460278E-02 )*X+4.28670143840073E-02 )*X-
778 1.29314370962569E-01 )*X+2.22381034453369E-01;
779 WW4 = ((((((((( 4.99660550769508E-09*X-7.94585963310120E-08)*X+
780 8.359072409485E-07)*X-7.422369210610E-06)*X+
781 5.763374308160E-05)*X-3.86645606718233E-04 )*X+
782 2.18417516259781E-03 )*X-9.99791027771119E-03 )*X+
783 3.48791097377370E-02 )*X-8.28299075413889E-02 )*X+
784 1.01228536290376E-01;
785 } else if (X <= 5) {
786 Y = X-3.0E+00;
787 RT1 = (((((((((-1.48570633747284E-15*Y-1.33273068108777E-13)*Y+
788 4.068543696670E-12)*Y-9.163164161821E-11)*Y+
789 2.046819017845E-09)*Y-4.03076426299031E-08 )*Y+
790 7.29407420660149E-07 )*Y-1.23118059980833E-05 )*Y+
791 1.88796581246938E-04 )*Y-2.53262912046853E-03 )*Y+
792 2.51198234505021E-02;
793 RT2 = ((((((((( 1.35830583483312E-13*Y-2.29772605964836E-12)*Y-
794 3.821500128045E-12)*Y+6.844424214735E-10)*Y-
795 1.048063352259E-08)*Y+1.50083186233363E-08 )*Y+
796 3.48848942324454E-06 )*Y-1.08694174399193E-04 )*Y+
797 2.08048885251999E-03 )*Y-2.91205805373793E-02 )*Y+
798 2.72276489515713E-01;
799 RT3 = ((((((((( 5.02799392850289E-13*Y+1.07461812944084E-11)*Y-
800 1.482277886411E-10)*Y-2.153585661215E-09)*Y+
801 3.654087802817E-08)*Y+5.15929575830120E-07 )*Y-
802 9.52388379435709E-06 )*Y-2.16552440036426E-04 )*Y+
803 9.03551469568320E-03 )*Y-1.45505469175613E-01 )*Y+
804 1.21449092319186E+00;
805 RT4 = (((((((((-1.08510370291979E-12*Y+6.41492397277798E-11)*Y+
806 7.542387436125E-10)*Y-2.213111836647E-09)*Y-
807 1.448228963549E-07)*Y-1.95670833237101E-06 )*Y-
808 1.07481314670844E-05 )*Y+1.49335941252765E-04 )*Y+
809 4.87791531990593E-02 )*Y-1.10559909038653E+00 )*Y+
810 8.09502028611780E+00;
811 WW1 = ((((((((((-4.65801912689961E-14*Y+7.58669507106800E-13)*Y-
812 1.186387548048E-11)*Y+1.862334710665E-10)*Y-
813 2.799399389539E-09)*Y+4.148972684255E-08)*Y-
814 5.933568079600E-07)*Y+8.168349266115E-06)*Y-
815 1.08989176177409E-04 )*Y+1.41357961729531E-03 )*Y-
816 1.87588361833659E-02 )*Y+2.89898651436026E-01;
817 WW2 = ((((((((((((-1.46345073267549E-14*Y+2.25644205432182E-13)*Y-
818 3.116258693847E-12)*Y+4.321908756610E-11)*Y-
819 5.673270062669E-10)*Y+7.006295962960E-09)*Y-
820 8.120186517000E-08)*Y+8.775294645770E-07)*Y-
821 8.77829235749024E-06 )*Y+8.04372147732379E-05 )*Y-
822 6.64149238804153E-04 )*Y+4.81181506827225E-03 )*Y-
823 2.88982669486183E-02 )*Y+1.56247249979288E-01;
824 WW3 = ((((((((((((( 9.06812118895365E-15*Y-1.40541322766087E-13)*
825 Y+1.919270015269E-12)*Y-2.605135739010E-11)*Y+
826 3.299685839012E-10)*Y-3.86354139348735E-09 )*Y+
827 4.16265847927498E-08 )*Y-4.09462835471470E-07 )*Y+
828 3.64018881086111E-06 )*Y-2.88665153269386E-05 )*Y+
829 2.00515819789028E-04 )*Y-1.18791896897934E-03 )*Y+
830 5.75223633388589E-03 )*Y-2.09400418772687E-02 )*Y+
831 4.85368861938873E-02;
832 WW4 = ((((((((((((((-9.74835552342257E-16*Y+1.57857099317175E-14)*
833 Y-2.249993780112E-13)*Y+3.173422008953E-12)*Y-
834 4.161159459680E-11)*Y+5.021343560166E-10)*Y-
835 5.545047534808E-09)*Y+5.554146993491E-08)*Y-
836 4.99048696190133E-07 )*Y+3.96650392371311E-06 )*Y-
837 2.73816413291214E-05 )*Y+1.60106988333186E-04 )*Y-
838 7.64560567879592E-04 )*Y+2.81330044426892E-03 )*Y-
839 7.16227030134947E-03 )*Y+9.66077262223353E-03;
840 } else if (X <= 10.0) {
841 Y = X-7.5E+00;
842 RT1 = ((((((((( 4.64217329776215E-15*Y-6.27892383644164E-15)*Y+
843 3.462236347446E-13)*Y-2.927229355350E-11)*Y+
844 5.090355371676E-10)*Y-9.97272656345253E-09 )*Y+
845 2.37835295639281E-07 )*Y-4.60301761310921E-06 )*Y+
846 8.42824204233222E-05 )*Y-1.37983082233081E-03 )*Y+
847 1.66630865869375E-02;
848 RT2 = ((((((((( 2.93981127919047E-14*Y+8.47635639065744E-13)*Y-
849 1.446314544774E-11)*Y-6.149155555753E-12)*Y+
850 8.484275604612E-10)*Y-6.10898827887652E-08 )*Y+
851 2.39156093611106E-06 )*Y-5.35837089462592E-05 )*Y+
852 1.00967602595557E-03 )*Y-1.57769317127372E-02 )*Y+
853 1.74853819464285E-01;
854 RT3 = (((((((((( 2.93523563363000E-14*Y-6.40041776667020E-14)*Y-
855 2.695740446312E-12)*Y+1.027082960169E-10)*Y-
856 5.822038656780E-10)*Y-3.159991002539E-08)*Y+
857 4.327249251331E-07)*Y+4.856768455119E-06)*Y-
858 2.54617989427762E-04 )*Y+5.54843378106589E-03 )*Y-
859 7.95013029486684E-02 )*Y+7.20206142703162E-01;
860 RT4 = (((((((((((-1.62212382394553E-14*Y+7.68943641360593E-13)*Y+
861 5.764015756615E-12)*Y-1.380635298784E-10)*Y-
862 1.476849808675E-09)*Y+1.84347052385605E-08 )*Y+
863 3.34382940759405E-07 )*Y-1.39428366421645E-06 )*Y-
864 7.50249313713996E-05 )*Y-6.26495899187507E-04 )*Y+
865 4.69716410901162E-02 )*Y-6.66871297428209E-01 )*Y+
866 4.11207530217806E+00;
867 WW1 = ((((((((((-1.65995045235997E-15*Y+6.91838935879598E-14)*Y-
868 9.131223418888E-13)*Y+1.403341829454E-11)*Y-
869 3.672235069444E-10)*Y+6.366962546990E-09)*Y-
870 1.039220021671E-07)*Y+1.959098751715E-06)*Y-
871 3.33474893152939E-05 )*Y+5.72164211151013E-04 )*Y-
872 1.05583210553392E-02 )*Y+2.26696066029591E-01;
873 WW2 = ((((((((((((-3.57248951192047E-16*Y+6.25708409149331E-15)*Y-
874 9.657033089714E-14)*Y+1.507864898748E-12)*Y-
875 2.332522256110E-11)*Y+3.428545616603E-10)*Y-
876 4.698730937661E-09)*Y+6.219977635130E-08)*Y-
877 7.83008889613661E-07 )*Y+9.08621687041567E-06 )*Y-
878 9.86368311253873E-05 )*Y+9.69632496710088E-04 )*Y-
879 8.14594214284187E-03 )*Y+8.50218447733457E-02;
880 WW3 = ((((((((((((( 1.64742458534277E-16*Y-2.68512265928410E-15)*
881 Y+3.788890667676E-14)*Y-5.508918529823E-13)*Y+
882 7.555896810069E-12)*Y-9.69039768312637E-11 )*Y+
883 1.16034263529672E-09 )*Y-1.28771698573873E-08 )*Y+
884 1.31949431805798E-07 )*Y-1.23673915616005E-06 )*Y+
885 1.04189803544936E-05 )*Y-7.79566003744742E-05 )*Y+
886 5.03162624754434E-04 )*Y-2.55138844587555E-03 )*Y+
887 1.13250730954014E-02;
888 WW4 = ((((((((((((((-1.55714130075679E-17*Y+2.57193722698891E-16)*
889 Y-3.626606654097E-15)*Y+5.234734676175E-14)*Y-
890 7.067105402134E-13)*Y+8.793512664890E-12)*Y-
891 1.006088923498E-10)*Y+1.050565098393E-09)*Y-
892 9.91517881772662E-09 )*Y+8.35835975882941E-08 )*Y-
893 6.19785782240693E-07 )*Y+3.95841149373135E-06 )*Y-
894 2.11366761402403E-05 )*Y+9.00474771229507E-05 )*Y-
895 2.78777909813289E-04 )*Y+5.26543779837487E-04;
896 } else if (X <= 15) {
897 Y = X-12.5E+00;
898 RT1 = ((((((((((( 4.94869622744119E-17*Y+8.03568805739160E-16)*Y-
899 5.599125915431E-15)*Y-1.378685560217E-13)*Y+
900 7.006511663249E-13)*Y+1.30391406991118E-11 )*Y+
901 8.06987313467541E-11 )*Y-5.20644072732933E-09 )*Y+
902 7.72794187755457E-08 )*Y-1.61512612564194E-06 )*Y+
903 4.15083811185831E-05 )*Y-7.87855975560199E-04 )*Y+
904 1.14189319050009E-02;
905 RT2 = ((((((((((( 4.89224285522336E-16*Y+1.06390248099712E-14)*Y-
906 5.446260182933E-14)*Y-1.613630106295E-12)*Y+
907 3.910179118937E-12)*Y+1.90712434258806E-10 )*Y+
908 8.78470199094761E-10 )*Y-5.97332993206797E-08 )*Y+
909 9.25750831481589E-07 )*Y-2.02362185197088E-05 )*Y+
910 4.92341968336776E-04 )*Y-8.68438439874703E-03 )*Y+
911 1.15825965127958E-01;
912 RT3 = (((((((((( 6.12419396208408E-14*Y+1.12328861406073E-13)*Y-
913 9.051094103059E-12)*Y-4.781797525341E-11)*Y+
914 1.660828868694E-09)*Y+4.499058798868E-10)*Y-
915 2.519549641933E-07)*Y+4.977444040180E-06)*Y-
916 1.25858350034589E-04 )*Y+2.70279176970044E-03 )*Y-
917 3.99327850801083E-02 )*Y+4.33467200855434E-01;
918 RT4 = ((((((((((( 4.63414725924048E-14*Y-4.72757262693062E-14)*Y-
919 1.001926833832E-11)*Y+6.074107718414E-11)*Y+
920 1.576976911942E-09)*Y-2.01186401974027E-08 )*Y-
921 1.84530195217118E-07 )*Y+5.02333087806827E-06 )*Y+
922 9.66961790843006E-06 )*Y-1.58522208889528E-03 )*Y+
923 2.80539673938339E-02 )*Y-2.78953904330072E-01 )*Y+
924 1.82835655238235E+00;
925 WW4 = ((((((((((((( 2.90401781000996E-18*Y-4.63389683098251E-17)*
926 Y+6.274018198326E-16)*Y-8.936002188168E-15)*Y+
927 1.194719074934E-13)*Y-1.45501321259466E-12 )*Y+
928 1.64090830181013E-11 )*Y-1.71987745310181E-10 )*Y+
929 1.63738403295718E-09 )*Y-1.39237504892842E-08 )*Y+
930 1.06527318142151E-07 )*Y-7.27634957230524E-07 )*Y+
931 4.12159381310339E-06 )*Y-1.74648169719173E-05 )*Y+
932 8.50290130067818E-05;
933 WW3 = ((((((((((((-4.19569145459480E-17*Y+5.94344180261644E-16)*Y-
934 1.148797566469E-14)*Y+1.881303962576E-13)*Y-
935 2.413554618391E-12)*Y+3.372127423047E-11)*Y-
936 4.933988617784E-10)*Y+6.116545396281E-09)*Y-
937 6.69965691739299E-08 )*Y+7.52380085447161E-07 )*Y-
938 8.08708393262321E-06 )*Y+6.88603417296672E-05 )*Y-
939 4.67067112993427E-04 )*Y+5.42313365864597E-03;
940 WW2 = ((((((((((-6.22272689880615E-15*Y+1.04126809657554E-13)*Y-
941 6.842418230913E-13)*Y+1.576841731919E-11)*Y-
942 4.203948834175E-10)*Y+6.287255934781E-09)*Y-
943 8.307159819228E-08)*Y+1.356478091922E-06)*Y-
944 2.08065576105639E-05 )*Y+2.52396730332340E-04 )*Y-
945 2.94484050194539E-03 )*Y+6.01396183129168E-02;
946 WW1 = (((-1.8784686463512E-01/X+2.2991849164985E-01)/X -
947 4.9893752514047E-01)/X-2.1916512131607E-05)*exp(-X) +
948 sqrt(PIE4/X)-WW4-WW3-WW2;
949 } else if (X <= 20) {
950 WW1 = sqrt(PIE4/X);
951 Y = X-17.5E+00;
952 RT1 = ((((((((((( 4.36701759531398E-17*Y-1.12860600219889E-16)*Y-
953 6.149849164164E-15)*Y+5.820231579541E-14)*Y+
954 4.396602872143E-13)*Y-1.24330365320172E-11 )*Y+
955 6.71083474044549E-11 )*Y+2.43865205376067E-10 )*Y+
956 1.67559587099969E-08 )*Y-9.32738632357572E-07 )*Y+
957 2.39030487004977E-05 )*Y-4.68648206591515E-04 )*Y+
958 8.34977776583956E-03;
959 RT2 = ((((((((((( 4.98913142288158E-16*Y-2.60732537093612E-16)*Y-
960 7.775156445127E-14)*Y+5.766105220086E-13)*Y+
961 6.432696729600E-12)*Y-1.39571683725792E-10 )*Y+
962 5.95451479522191E-10 )*Y+2.42471442836205E-09 )*Y+
963 2.47485710143120E-07 )*Y-1.14710398652091E-05 )*Y+
964 2.71252453754519E-04 )*Y-4.96812745851408E-03 )*Y+
965 8.26020602026780E-02;
966 RT3 = ((((((((((( 1.91498302509009E-15*Y+1.48840394311115E-14)*Y-
967 4.316925145767E-13)*Y+1.186495793471E-12)*Y+
968 4.615806713055E-11)*Y-5.54336148667141E-10 )*Y+
969 3.48789978951367E-10 )*Y-2.79188977451042E-09 )*Y+
970 2.09563208958551E-06 )*Y-6.76512715080324E-05 )*Y+
971 1.32129867629062E-03 )*Y-2.05062147771513E-02 )*Y+
972 2.88068671894324E-01;
973 RT4 = (((((((((((-5.43697691672942E-15*Y-1.12483395714468E-13)*Y+
974 2.826607936174E-12)*Y-1.266734493280E-11)*Y-
975 4.258722866437E-10)*Y+9.45486578503261E-09 )*Y-
976 5.86635622821309E-08 )*Y-1.28835028104639E-06 )*Y+
977 4.41413815691885E-05 )*Y-7.61738385590776E-04 )*Y+
978 9.66090902985550E-03 )*Y-1.01410568057649E-01 )*Y+
979 9.54714798156712E-01;
980 WW4 = ((((((((((((-7.56882223582704E-19*Y+7.53541779268175E-18)*Y-
981 1.157318032236E-16)*Y+2.411195002314E-15)*Y-
982 3.601794386996E-14)*Y+4.082150659615E-13)*Y-
983 4.289542980767E-12)*Y+5.086829642731E-11)*Y-
984 6.35435561050807E-10 )*Y+6.82309323251123E-09 )*Y-
985 5.63374555753167E-08 )*Y+3.57005361100431E-07 )*Y-
986 2.40050045173721E-06 )*Y+4.94171300536397E-05;
987 WW3 = (((((((((((-5.54451040921657E-17*Y+2.68748367250999E-16)*Y+
988 1.349020069254E-14)*Y-2.507452792892E-13)*Y+
989 1.944339743818E-12)*Y-1.29816917658823E-11 )*Y+
990 3.49977768819641E-10 )*Y-8.67270669346398E-09 )*Y+
991 1.31381116840118E-07 )*Y-1.36790720600822E-06 )*Y+
992 1.19210697673160E-05 )*Y-1.42181943986587E-04 )*Y+
993 4.12615396191829E-03;
994 WW2 = (((((((((((-1.86506057729700E-16*Y+1.16661114435809E-15)*Y+
995 2.563712856363E-14)*Y-4.498350984631E-13)*Y+
996 1.765194089338E-12)*Y+9.04483676345625E-12 )*Y+
997 4.98930345609785E-10 )*Y-2.11964170928181E-08 )*Y+
998 3.98295476005614E-07 )*Y-5.49390160829409E-06 )*Y+
999 7.74065155353262E-05 )*Y-1.48201933009105E-03 )*Y+
1000 4.97836392625268E-02;
1001 WW1 = (( 1.9623264149430E-01/X-4.9695241464490E-01)/X -
1002 6.0156581186481E-05)*exp(-X)+WW1-WW2-WW3-WW4;
1003 } else if (X <= 35) {
1004 WW1 = sqrt(PIE4/X);
1005 E = exp(-X);
1006 RT1 = ((((((-4.45711399441838E-05*X+1.27267770241379E-03)*X -
1007 2.36954961381262E-01)*X+1.54330657903756E+01)*X -
1008 5.22799159267808E+02)*X+1.05951216669313E+04)*X +
1009 (-2.51177235556236E+06/X+8.72975373557709E+05)/X -
1010 1.29194382386499E+05)*E + R14/(X-R14);
1011 RT2 = (((((-7.85617372254488E-02*X+6.35653573484868E+00)*X -
1012 3.38296938763990E+02)*X+1.25120495802096E+04)*X -
1013 3.16847570511637E+05)*X +
1014 ((-1.02427466127427E+09/X +
1015 3.70104713293016E+08)/X-5.87119005093822E+07)/X +
1016 5.38614211391604E+06)*E + R24/(X-R24);
1017 RT3 = (((((-2.37900485051067E-01*X+1.84122184400896E+01)*X -
1018 1.00200731304146E+03)*X+3.75151841595736E+04)*X -
1019 9.50626663390130E+05)*X +
1020 ((-2.88139014651985E+09/X +
1021 1.06625915044526E+09)/X-1.72465289687396E+08)/X +
1022 1.60419390230055E+07)*E + R34/(X-R34);
1023 RT4 = ((((((-6.00691586407385E-04*X-3.64479545338439E-01)*X +
1024 1.57496131755179E+01)*X-6.54944248734901E+02)*X +
1025 1.70830039597097E+04)*X-2.90517939780207E+05)*X +
1026 (3.49059698304732E+07/X-1.64944522586065E+07)/X +
1027 2.96817940164703E+06)*E + R44/(X-R44);
1028 if (X <= 25)
1029 WW4 = ((((((( 2.33766206773151E-07*X-
1030 3.81542906607063E-05)*X +3.51416601267000E-03)*X-
1031 1.66538571864728E-01)*X +4.80006136831847E+00)*X-
1032 8.73165934223603E+01)*X +9.77683627474638E+02)*X +
1033 1.66000945117640E+04/X -6.14479071209961E+03)*E + W44*WW1;
1034 else
1035 WW4 = (((((( 5.74245945342286E-06*X-
1036 7.58735928102351E-05)*X +2.35072857922892E-04)*X-
1037 3.78812134013125E-03)*X +3.09871652785805E-01)*X-
1038 7.11108633061306E+00)*X +5.55297573149528E+01)*E + W44*WW1;
1039 WW3 = (((((( 2.36392855180768E-04*X-9.16785337967013E-03)*X +
1040 4.62186525041313E-01)*X-1.96943786006540E+01)*X +
1041 4.99169195295559E+02)*X-6.21419845845090E+03)*X +
1042 ((+5.21445053212414E+07/X-1.34113464389309E+07)/X +
1043 1.13673298305631E+06)/X-2.81501182042707E+03)*E + W34*WW1;
1044 WW2 = (((((( 7.29841848989391E-04*X-3.53899555749875E-02)*X +
1045 2.07797425718513E+00)*X-1.00464709786287E+02)*X +
1046 3.15206108877819E+03)*X-6.27054715090012E+04)*X +
1047 (+1.54721246264919E+07/X-5.26074391316381E+06)/X +
1048 7.67135400969617E+05)*E + W24*WW1;
1049 WW1 = (( 1.9623264149430E-01/X-4.9695241464490E-01)/X -
1050 6.0156581186481E-05)*E + WW1-WW2-WW3-WW4;
1051 } else if (X <= 53) {
1052 WW1 = sqrt(PIE4/X);
1053 E = exp(-X)*pow(X,4);
1054 RT4 = ((-2.19135070169653E-03*X-1.19108256987623E-01)*X -
1055 7.50238795695573E-01)*E + R44/(X-R44);
1056 RT3 = ((-9.65842534508637E-04*X-4.49822013469279E-02)*X +
1057 6.08784033347757E-01)*E + R34/(X-R34);
1058 RT2 = ((-3.62569791162153E-04*X-9.09231717268466E-03)*X +
1059 1.84336760556262E-01)*E + R24/(X-R24);
1060 RT1 = ((-4.07557525914600E-05*X-6.88846864931685E-04)*X +
1061 1.74725309199384E-02)*E + R14/(X-R14);
1062 WW4 = (( 5.76631982000990E-06*X-7.89187283804890E-05)*X +
1063 3.28297971853126E-04)*E + W44*WW1;
1064 WW3 = (( 2.08294969857230E-04*X-3.77489954837361E-03)*X +
1065 2.09857151617436E-02)*E + W34*WW1;
1066 WW2 = (( 6.16374517326469E-04*X-1.26711744680092E-02)*X +
1067 8.14504890732155E-02)*E + W24*WW1;
1068 WW1 = WW1-WW2-WW3-WW4;
1069 } else {
1070 WW1 = sqrt(PIE4/X);
1071 RT1 = R14/(X-R14);
1072 RT2 = R24/(X-R24);
1073 RT3 = R34/(X-R34);
1074 RT4 = R44/(X-R44);
1075 WW4 = W44*WW1;
1076 WW3 = W34*WW1;
1077 WW2 = W24*WW1;
1078 WW1 = WW1-WW2-WW3-WW4;
1079 }
1080 roots[0] = RT1;
1081 roots[1] = RT2;
1082 roots[2] = RT3;
1083 roots[3] = RT4;
1084 weights[0] = WW1;
1085 weights[1] = WW2;
1086 weights[2] = WW3;
1087 weights[3] = WW4;
1088 }
1089
rys_root5(double X,double * roots,double * weights)1090 static void rys_root5(double X, double *roots, double *weights)
1091 {
1092 double R15,R25,W25,R35,W35,R45,W45,R55,W55;
1093 double RT1, RT2, RT3, RT4, RT5, WW1, WW2, WW3, WW4, WW5;
1094 double Y, E, XXX;
1095
1096 R15 = 1.17581320211778E-01;
1097 R25 = 1.07456201243690E+00;
1098 W25 = 2.70967405960535E-01;
1099 R35 = 3.08593744371754E+00;
1100 W35 = 3.82231610015404E-02;
1101 R45 = 6.41472973366203E+00;
1102 W45 = 1.51614186862443E-03;
1103 R55 = 1.18071894899717E+01;
1104 W55 = 8.62130526143657E-06;
1105
1106 if (X < 3.e-7){
1107 RT1 = 2.26659266316985E-02 -2.15865967920897E-03 *X;
1108 RT2 = 2.31271692140903E-01 -2.20258754389745E-02 *X;
1109 RT3 = 8.57346024118836E-01 -8.16520023025515E-02 *X;
1110 RT4 = 2.97353038120346E+00 -2.83193369647137E-01 *X;
1111 RT5 = 1.84151859759051E+01 -1.75382723579439E+00 *X;
1112 WW1 = 2.95524224714752E-01 -1.96867576909777E-02 *X;
1113 WW2 = 2.69266719309995E-01 -5.61737590184721E-02 *X;
1114 WW3 = 2.19086362515981E-01 -9.71152726793658E-02 *X;
1115 WW4 = 1.49451349150580E-01 -1.02979262193565E-01 *X;
1116 WW5 = 6.66713443086877E-02 -5.73782817488315E-02 *X;
1117 } else if (X < 1.0){
1118 RT1 = ((((((-4.46679165328413E-11*X+1.21879111988031E-09)*X-
1119 2.62975022612104E-08 )*X+5.15106194905897E-07 )*X-
1120 9.27933625824749E-06 )*X+1.51794097682482E-04 )*X-
1121 2.15865967920301E-03 )*X+2.26659266316985E-02;
1122 RT2 = (((((( 1.93117331714174E-10*X-4.57267589660699E-09)*X+
1123 2.48339908218932E-08 )*X+1.50716729438474E-06 )*X-
1124 6.07268757707381E-05 )*X+1.37506939145643E-03 )*X-
1125 2.20258754419939E-02 )*X+2.31271692140905E-01;
1126 RT3 = ((((( 4.84989776180094E-09*X+1.31538893944284E-07)*X-
1127 2.766753852879E-06)*X-7.651163510626E-05)*X+
1128 4.033058545972E-03)*X-8.16520022916145E-02 )*X+
1129 8.57346024118779E-01;
1130 RT4 = ((((-2.48581772214623E-07*X-4.34482635782585E-06)*X-
1131 7.46018257987630E-07 )*X+1.01210776517279E-02 )*X-
1132 2.83193369640005E-01 )*X+2.97353038120345E+00;
1133 RT5 = (((((-8.92432153868554E-09*X+1.77288899268988E-08)*X+
1134 3.040754680666E-06)*X+1.058229325071E-04)*X+
1135 4.596379534985E-02)*X-1.75382723579114E+00 )*X+
1136 1.84151859759049E+01;
1137 WW1 = ((((((-2.03822632771791E-09*X+3.89110229133810E-08)*X-
1138 5.84914787904823E-07 )*X+8.30316168666696E-06 )*X-
1139 1.13218402310546E-04 )*X+1.49128888586790E-03 )*X-
1140 1.96867576904816E-02 )*X+2.95524224714749E-01;
1141 WW2 = ((((((( 8.62848118397570E-09*X-1.38975551148989E-07)*X+
1142 1.602894068228E-06)*X-1.646364300836E-05)*X+
1143 1.538445806778E-04)*X-1.28848868034502E-03 )*X+
1144 9.38866933338584E-03 )*X-5.61737590178812E-02 )*X+
1145 2.69266719309991E-01;
1146 WW3 = ((((((((-9.41953204205665E-09*X+1.47452251067755E-07)*X-
1147 1.57456991199322E-06 )*X+1.45098401798393E-05 )*X-
1148 1.18858834181513E-04 )*X+8.53697675984210E-04 )*X-
1149 5.22877807397165E-03 )*X+2.60854524809786E-02 )*X-
1150 9.71152726809059E-02 )*X+2.19086362515979E-01;
1151 WW4 = ((((((((-3.84961617022042E-08*X+5.66595396544470E-07)*X-
1152 5.52351805403748E-06 )*X+4.53160377546073E-05 )*X-
1153 3.22542784865557E-04 )*X+1.95682017370967E-03 )*X-
1154 9.77232537679229E-03 )*X+3.79455945268632E-02 )*X-
1155 1.02979262192227E-01 )*X+1.49451349150573E-01;
1156 WW5 = ((((((((( 4.09594812521430E-09*X-6.47097874264417E-08)*X+
1157 6.743541482689E-07)*X-5.917993920224E-06)*X+
1158 4.531969237381E-05)*X-2.99102856679638E-04 )*X+
1159 1.65695765202643E-03 )*X-7.40671222520653E-03 )*X+
1160 2.50889946832192E-02 )*X-5.73782817487958E-02 )*X+
1161 6.66713443086877E-02;
1162 } else if (X < 5.0) {
1163 Y = X-3.0E+00;
1164 RT1 = ((((((((-2.58163897135138E-14*Y+8.14127461488273E-13)*Y-
1165 2.11414838976129E-11 )*Y+5.09822003260014E-10 )*Y-
1166 1.16002134438663E-08 )*Y+2.46810694414540E-07 )*Y-
1167 4.92556826124502E-06 )*Y+9.02580687971053E-05 )*Y-
1168 1.45190025120726E-03 )*Y+1.73416786387475E-02;
1169 RT2 = ((((((((( 1.04525287289788E-14*Y+5.44611782010773E-14)*Y-
1170 4.831059411392E-12)*Y+1.136643908832E-10)*Y-
1171 1.104373076913E-09)*Y-2.35346740649916E-08 )*Y+
1172 1.43772622028764E-06 )*Y-4.23405023015273E-05 )*Y+
1173 9.12034574793379E-04 )*Y-1.52479441718739E-02 )*Y+
1174 1.76055265928744E-01;
1175 RT3 = (((((((((-6.89693150857911E-14*Y+5.92064260918861E-13)*Y+
1176 1.847170956043E-11)*Y-3.390752744265E-10)*Y-
1177 2.995532064116E-09)*Y+1.57456141058535E-07 )*Y-
1178 3.95859409711346E-07 )*Y-9.58924580919747E-05 )*Y+
1179 3.23551502557785E-03 )*Y-5.97587007636479E-02 )*Y+
1180 6.46432853383057E-01;
1181 RT4 = ((((((((-3.61293809667763E-12*Y-2.70803518291085E-11)*Y+
1182 8.83758848468769E-10 )*Y+1.59166632851267E-08 )*Y-
1183 1.32581997983422E-07 )*Y-7.60223407443995E-06 )*Y-
1184 7.41019244900952E-05 )*Y+9.81432631743423E-03 )*Y-
1185 2.23055570487771E-01 )*Y+2.21460798080643E+00;
1186 RT5 = ((((((((( 7.12332088345321E-13*Y+3.16578501501894E-12)*Y-
1187 8.776668218053E-11)*Y-2.342817613343E-09)*Y-
1188 3.496962018025E-08)*Y-3.03172870136802E-07 )*Y+
1189 1.50511293969805E-06 )*Y+1.37704919387696E-04 )*Y+
1190 4.70723869619745E-02 )*Y-1.47486623003693E+00 )*Y+
1191 1.35704792175847E+01;
1192 WW1 = ((((((((( 1.04348658616398E-13*Y-1.94147461891055E-12)*Y+
1193 3.485512360993E-11)*Y-6.277497362235E-10)*Y+
1194 1.100758247388E-08)*Y-1.88329804969573E-07 )*Y+
1195 3.12338120839468E-06 )*Y-5.04404167403568E-05 )*Y+
1196 8.00338056610995E-04 )*Y-1.30892406559521E-02 )*Y+
1197 2.47383140241103E-01;
1198 WW2 = ((((((((((( 3.23496149760478E-14*Y-5.24314473469311E-13)*Y+
1199 7.743219385056E-12)*Y-1.146022750992E-10)*Y+
1200 1.615238462197E-09)*Y-2.15479017572233E-08 )*Y+
1201 2.70933462557631E-07 )*Y-3.18750295288531E-06 )*Y+
1202 3.47425221210099E-05 )*Y-3.45558237388223E-04 )*Y+
1203 3.05779768191621E-03 )*Y-2.29118251223003E-02 )*Y+
1204 1.59834227924213E-01;
1205 WW3 = ((((((((((((-3.42790561802876E-14*Y+5.26475736681542E-13)*Y-
1206 7.184330797139E-12)*Y+9.763932908544E-11)*Y-
1207 1.244014559219E-09)*Y+1.472744068942E-08)*Y-
1208 1.611749975234E-07)*Y+1.616487851917E-06)*Y-
1209 1.46852359124154E-05 )*Y+1.18900349101069E-04 )*Y-
1210 8.37562373221756E-04 )*Y+4.93752683045845E-03 )*Y-
1211 2.25514728915673E-02 )*Y+6.95211812453929E-02;
1212 WW4 = ((((((((((((( 1.04072340345039E-14*Y-1.60808044529211E-13)*
1213 Y+2.183534866798E-12)*Y-2.939403008391E-11)*Y+
1214 3.679254029085E-10)*Y-4.23775673047899E-09 )*Y+
1215 4.46559231067006E-08 )*Y-4.26488836563267E-07 )*Y+
1216 3.64721335274973E-06 )*Y-2.74868382777722E-05 )*Y+
1217 1.78586118867488E-04 )*Y-9.68428981886534E-04 )*Y+
1218 4.16002324339929E-03 )*Y-1.28290192663141E-02 )*Y+
1219 2.22353727685016E-02;
1220 WW5 = ((((((((((((((-8.16770412525963E-16*Y+1.31376515047977E-14)*
1221 Y-1.856950818865E-13)*Y+2.596836515749E-12)*Y-
1222 3.372639523006E-11)*Y+4.025371849467E-10)*Y-
1223 4.389453269417E-09)*Y+4.332753856271E-08)*Y-
1224 3.82673275931962E-07 )*Y+2.98006900751543E-06 )*Y-
1225 2.00718990300052E-05 )*Y+1.13876001386361E-04 )*Y-
1226 5.23627942443563E-04 )*Y+1.83524565118203E-03 )*Y-
1227 4.37785737450783E-03 )*Y+5.36963805223095E-03;
1228 } else if (X < 10.0) {
1229 Y = X-7.5E+00;
1230 RT1 = ((((((((-1.13825201010775E-14*Y+1.89737681670375E-13)*Y-
1231 4.81561201185876E-12 )*Y+1.56666512163407E-10 )*Y-
1232 3.73782213255083E-09 )*Y+9.15858355075147E-08 )*Y-
1233 2.13775073585629E-06 )*Y+4.56547356365536E-05 )*Y-
1234 8.68003909323740E-04 )*Y+1.22703754069176E-02;
1235 RT2 = (((((((((-3.67160504428358E-15*Y+1.27876280158297E-14)*Y-
1236 1.296476623788E-12)*Y+1.477175434354E-11)*Y+
1237 5.464102147892E-10)*Y-2.42538340602723E-08 )*Y+
1238 8.20460740637617E-07 )*Y-2.20379304598661E-05 )*Y+
1239 4.90295372978785E-04 )*Y-9.14294111576119E-03 )*Y+
1240 1.22590403403690E-01;
1241 RT3 = ((((((((( 1.39017367502123E-14*Y-6.96391385426890E-13)*Y+
1242 1.176946020731E-12)*Y+1.725627235645E-10)*Y-
1243 3.686383856300E-09)*Y+2.87495324207095E-08 )*Y+
1244 1.71307311000282E-06 )*Y-7.94273603184629E-05 )*Y+
1245 2.00938064965897E-03 )*Y-3.63329491677178E-02 )*Y+
1246 4.34393683888443E-01;
1247 RT4 = ((((((((((-1.27815158195209E-14*Y+1.99910415869821E-14)*Y+
1248 3.753542914426E-12)*Y-2.708018219579E-11)*Y-
1249 1.190574776587E-09)*Y+1.106696436509E-08)*Y+
1250 3.954955671326E-07)*Y-4.398596059588E-06)*Y-
1251 2.01087998907735E-04 )*Y+7.89092425542937E-03 )*Y-
1252 1.42056749162695E-01 )*Y+1.39964149420683E+00;
1253 RT5 = ((((((((((-1.19442341030461E-13*Y-2.34074833275956E-12)*Y+
1254 6.861649627426E-12)*Y+6.082671496226E-10)*Y+
1255 5.381160105420E-09)*Y-6.253297138700E-08)*Y-
1256 2.135966835050E-06)*Y-2.373394341886E-05)*Y+
1257 2.88711171412814E-06 )*Y+4.85221195290753E-02 )*Y-
1258 1.04346091985269E+00 )*Y+7.89901551676692E+00;
1259 WW1 = ((((((((( 7.95526040108997E-15*Y-2.48593096128045E-13)*Y+
1260 4.761246208720E-12)*Y-9.535763686605E-11)*Y+
1261 2.225273630974E-09)*Y-4.49796778054865E-08 )*Y+
1262 9.17812870287386E-07 )*Y-1.86764236490502E-05 )*Y+
1263 3.76807779068053E-04 )*Y-8.10456360143408E-03 )*Y+
1264 2.01097936411496E-01;
1265 WW2 = ((((((((((( 1.25678686624734E-15*Y-2.34266248891173E-14)*Y+
1266 3.973252415832E-13)*Y-6.830539401049E-12)*Y+
1267 1.140771033372E-10)*Y-1.82546185762009E-09 )*Y+
1268 2.77209637550134E-08 )*Y-4.01726946190383E-07 )*Y+
1269 5.48227244014763E-06 )*Y-6.95676245982121E-05 )*Y+
1270 8.05193921815776E-04 )*Y-8.15528438784469E-03 )*Y+
1271 9.71769901268114E-02;
1272 WW3 = ((((((((((((-8.20929494859896E-16*Y+1.37356038393016E-14)*Y-
1273 2.022863065220E-13)*Y+3.058055403795E-12)*Y-
1274 4.387890955243E-11)*Y+5.923946274445E-10)*Y-
1275 7.503659964159E-09)*Y+8.851599803902E-08)*Y-
1276 9.65561998415038E-07 )*Y+9.60884622778092E-06 )*Y-
1277 8.56551787594404E-05 )*Y+6.66057194311179E-04 )*Y-
1278 4.17753183902198E-03 )*Y+2.25443826852447E-02;
1279 WW4 = ((((((((((((((-1.08764612488790E-17*Y+1.85299909689937E-16)*
1280 Y-2.730195628655E-15)*Y+4.127368817265E-14)*Y-
1281 5.881379088074E-13)*Y+7.805245193391E-12)*Y-
1282 9.632707991704E-11)*Y+1.099047050624E-09)*Y-
1283 1.15042731790748E-08 )*Y+1.09415155268932E-07 )*Y-
1284 9.33687124875935E-07 )*Y+7.02338477986218E-06 )*Y-
1285 4.53759748787756E-05 )*Y+2.41722511389146E-04 )*Y-
1286 9.75935943447037E-04 )*Y+2.57520532789644E-03;
1287 WW5 = ((((((((((((((( 7.28996979748849E-19*Y-1.26518146195173E-17)
1288 *Y+1.886145834486E-16)*Y-2.876728287383E-15)*Y+
1289 4.114588668138E-14)*Y-5.44436631413933E-13 )*Y+
1290 6.64976446790959E-12 )*Y-7.44560069974940E-11 )*Y+
1291 7.57553198166848E-10 )*Y-6.92956101109829E-09 )*Y+
1292 5.62222859033624E-08 )*Y-3.97500114084351E-07 )*Y+
1293 2.39039126138140E-06 )*Y-1.18023950002105E-05 )*Y+
1294 4.52254031046244E-05 )*Y-1.21113782150370E-04 )*Y+
1295 1.75013126731224E-04;
1296 } else if (X < 15.0) {
1297 Y = X-12.5E+00;
1298 RT1 = ((((((((((-4.16387977337393E-17*Y+7.20872997373860E-16)*Y+
1299 1.395993802064E-14)*Y+3.660484641252E-14)*Y-
1300 4.154857548139E-12)*Y+2.301379846544E-11)*Y-
1301 1.033307012866E-09)*Y+3.997777641049E-08)*Y-
1302 9.35118186333939E-07 )*Y+2.38589932752937E-05 )*Y-
1303 5.35185183652937E-04 )*Y+8.85218988709735E-03;
1304 RT2 = ((((((((((-4.56279214732217E-16*Y+6.24941647247927E-15)*Y+
1305 1.737896339191E-13)*Y+8.964205979517E-14)*Y-
1306 3.538906780633E-11)*Y+9.561341254948E-11)*Y-
1307 9.772831891310E-09)*Y+4.240340194620E-07)*Y-
1308 1.02384302866534E-05 )*Y+2.57987709704822E-04 )*Y-
1309 5.54735977651677E-03 )*Y+8.68245143991948E-02;
1310 RT3 = ((((((((((-2.52879337929239E-15*Y+2.13925810087833E-14)*Y+
1311 7.884307667104E-13)*Y-9.023398159510E-13)*Y-
1312 5.814101544957E-11)*Y-1.333480437968E-09)*Y-
1313 2.217064940373E-08)*Y+1.643290788086E-06)*Y-
1314 4.39602147345028E-05 )*Y+1.08648982748911E-03 )*Y-
1315 2.13014521653498E-02 )*Y+2.94150684465425E-01;
1316 RT4 = ((((((((((-6.42391438038888E-15*Y+5.37848223438815E-15)*Y+
1317 8.960828117859E-13)*Y+5.214153461337E-11)*Y-
1318 1.106601744067E-10)*Y-2.007890743962E-08)*Y+
1319 1.543764346501E-07)*Y+4.520749076914E-06)*Y-
1320 1.88893338587047E-04 )*Y+4.73264487389288E-03 )*Y-
1321 7.91197893350253E-02 )*Y+8.60057928514554E-01;
1322 RT5 = (((((((((((-2.24366166957225E-14*Y+4.87224967526081E-14)*Y+
1323 5.587369053655E-12)*Y-3.045253104617E-12)*Y-
1324 1.223983883080E-09)*Y-2.05603889396319E-09 )*Y+
1325 2.58604071603561E-07 )*Y+1.34240904266268E-06 )*Y-
1326 5.72877569731162E-05 )*Y-9.56275105032191E-04 )*Y+
1327 4.23367010370921E-02 )*Y-5.76800927133412E-01 )*Y+
1328 3.87328263873381E+00;
1329 WW1 = ((((((((( 8.98007931950169E-15*Y+7.25673623859497E-14)*Y+
1330 5.851494250405E-14)*Y-4.234204823846E-11)*Y+
1331 3.911507312679E-10)*Y-9.65094802088511E-09 )*Y+
1332 3.42197444235714E-07 )*Y-7.51821178144509E-06 )*Y+
1333 1.94218051498662E-04 )*Y-5.38533819142287E-03 )*Y+
1334 1.68122596736809E-01;
1335 WW2 = ((((((((((-1.05490525395105E-15*Y+1.96855386549388E-14)*Y-
1336 5.500330153548E-13)*Y+1.003849567976E-11)*Y-
1337 1.720997242621E-10)*Y+3.533277061402E-09)*Y-
1338 6.389171736029E-08)*Y+1.046236652393E-06)*Y-
1339 1.73148206795827E-05 )*Y+2.57820531617185E-04 )*Y-
1340 3.46188265338350E-03 )*Y+7.03302497508176E-02;
1341 WW3 = ((((((((((( 3.60020423754545E-16*Y-6.24245825017148E-15)*Y+
1342 9.945311467434E-14)*Y-1.749051512721E-12)*Y+
1343 2.768503957853E-11)*Y-4.08688551136506E-10 )*Y+
1344 6.04189063303610E-09 )*Y-8.23540111024147E-08 )*Y+
1345 1.01503783870262E-06 )*Y-1.20490761741576E-05 )*Y+
1346 1.26928442448148E-04 )*Y-1.05539461930597E-03 )*Y+
1347 1.15543698537013E-02;
1348 WW4 = ((((((((((((( 2.51163533058925E-18*Y-4.31723745510697E-17)*
1349 Y+6.557620865832E-16)*Y-1.016528519495E-14)*Y+
1350 1.491302084832E-13)*Y-2.06638666222265E-12 )*Y+
1351 2.67958697789258E-11 )*Y-3.23322654638336E-10 )*Y+
1352 3.63722952167779E-09 )*Y-3.75484943783021E-08 )*Y+
1353 3.49164261987184E-07 )*Y-2.92658670674908E-06 )*Y+
1354 2.12937256719543E-05 )*Y-1.19434130620929E-04 )*Y+
1355 6.45524336158384E-04;
1356 WW5 = ((((((((((((((-1.29043630202811E-19*Y+2.16234952241296E-18)*
1357 Y-3.107631557965E-17)*Y+4.570804313173E-16)*Y-
1358 6.301348858104E-15)*Y+8.031304476153E-14)*Y-
1359 9.446196472547E-13)*Y+1.018245804339E-11)*Y-
1360 9.96995451348129E-11 )*Y+8.77489010276305E-10 )*Y-
1361 6.84655877575364E-09 )*Y+4.64460857084983E-08 )*Y-
1362 2.66924538268397E-07 )*Y+1.24621276265907E-06 )*Y-
1363 4.30868944351523E-06 )*Y+9.94307982432868E-06;
1364 } else if (X < 20.0){
1365 Y = X-17.5E+00;
1366 RT1 = (((((((((( 1.91875764545740E-16*Y+7.8357401095707E-16)*Y-
1367 3.260875931644E-14)*Y-1.186752035569E-13)*Y+
1368 4.275180095653E-12)*Y+3.357056136731E-11)*Y-
1369 1.123776903884E-09)*Y+1.231203269887E-08)*Y-
1370 3.99851421361031E-07 )*Y+1.45418822817771E-05 )*Y-
1371 3.49912254976317E-04 )*Y+6.67768703938812E-03;
1372 RT2 = (((((((((( 2.02778478673555E-15*Y+1.01640716785099E-14)*Y-
1373 3.385363492036E-13)*Y-1.615655871159E-12)*Y+
1374 4.527419140333E-11)*Y+3.853670706486E-10)*Y-
1375 1.184607130107E-08)*Y+1.347873288827E-07)*Y-
1376 4.47788241748377E-06 )*Y+1.54942754358273E-04 )*Y-
1377 3.55524254280266E-03 )*Y+6.44912219301603E-02;
1378 RT3 = (((((((((( 7.79850771456444E-15*Y+6.00464406395001E-14)*Y-
1379 1.249779730869E-12)*Y-1.020720636353E-11)*Y+
1380 1.814709816693E-10)*Y+1.766397336977E-09)*Y-
1381 4.603559449010E-08)*Y+5.863956443581E-07)*Y-
1382 2.03797212506691E-05 )*Y+6.31405161185185E-04 )*Y-
1383 1.30102750145071E-02 )*Y+2.10244289044705E-01;
1384 RT4 = (((((((((((-2.92397030777912E-15*Y+1.94152129078465E-14)*Y+
1385 4.859447665850E-13)*Y-3.217227223463E-12)*Y-
1386 7.484522135512E-11)*Y+7.19101516047753E-10 )*Y+
1387 6.88409355245582E-09 )*Y-1.44374545515769E-07 )*Y+
1388 2.74941013315834E-06 )*Y-1.02790452049013E-04 )*Y+
1389 2.59924221372643E-03 )*Y-4.35712368303551E-02 )*Y+
1390 5.62170709585029E-01;
1391 RT5 = ((((((((((( 1.17976126840060E-14*Y+1.24156229350669E-13)*Y-
1392 3.892741622280E-12)*Y-7.755793199043E-12)*Y+
1393 9.492190032313E-10)*Y-4.98680128123353E-09 )*Y-
1394 1.81502268782664E-07 )*Y+2.69463269394888E-06 )*Y+
1395 2.50032154421640E-05 )*Y-1.33684303917681E-03 )*Y+
1396 2.29121951862538E-02 )*Y-2.45653725061323E-01 )*Y+
1397 1.89999883453047E+00;
1398 WW1 = (((((((((( 1.74841995087592E-15*Y-6.95671892641256E-16)*Y-
1399 3.000659497257E-13)*Y+2.021279817961E-13)*Y+
1400 3.853596935400E-11)*Y+1.461418533652E-10)*Y-
1401 1.014517563435E-08)*Y+1.132736008979E-07)*Y-
1402 2.86605475073259E-06 )*Y+1.21958354908768E-04 )*Y-
1403 3.86293751153466E-03 )*Y+1.45298342081522E-01;
1404 WW2 = ((((((((((-1.11199320525573E-15*Y+1.85007587796671E-15)*Y+
1405 1.220613939709E-13)*Y+1.275068098526E-12)*Y-
1406 5.341838883262E-11)*Y+6.161037256669E-10)*Y-
1407 1.009147879750E-08)*Y+2.907862965346E-07)*Y-
1408 6.12300038720919E-06 )*Y+1.00104454489518E-04 )*Y-
1409 1.80677298502757E-03 )*Y+5.78009914536630E-02;
1410 WW3 = ((((((((((-9.49816486853687E-16*Y+6.67922080354234E-15)*Y+
1411 2.606163540537E-15)*Y+1.983799950150E-12)*Y-
1412 5.400548574357E-11)*Y+6.638043374114E-10)*Y-
1413 8.799518866802E-09)*Y+1.791418482685E-07)*Y-
1414 2.96075397351101E-06 )*Y+3.38028206156144E-05 )*Y-
1415 3.58426847857878E-04 )*Y+8.39213709428516E-03;
1416 WW4 = ((((((((((( 1.33829971060180E-17*Y-3.44841877844140E-16)*Y+
1417 4.745009557656E-15)*Y-6.033814209875E-14)*Y+
1418 1.049256040808E-12)*Y-1.70859789556117E-11 )*Y+
1419 2.15219425727959E-10 )*Y-2.52746574206884E-09 )*Y+
1420 3.27761714422960E-08 )*Y-3.90387662925193E-07 )*Y+
1421 3.46340204593870E-06 )*Y-2.43236345136782E-05 )*Y+
1422 3.54846978585226E-04;
1423 WW5 = ((((((((((((( 2.69412277020887E-20*Y-4.24837886165685E-19)*
1424 Y+6.030500065438E-18)*Y-9.069722758289E-17)*Y+
1425 1.246599177672E-15)*Y-1.56872999797549E-14 )*Y+
1426 1.87305099552692E-13 )*Y-2.09498886675861E-12 )*Y+
1427 2.11630022068394E-11 )*Y-1.92566242323525E-10 )*Y+
1428 1.62012436344069E-09 )*Y-1.23621614171556E-08 )*Y+
1429 7.72165684563049E-08 )*Y-3.59858901591047E-07 )*Y+
1430 2.43682618601000E-06;
1431 } else if (X < 25.0) {
1432 Y = X-22.5E+00;
1433 RT1 = (((((((((-1.13927848238726E-15*Y+7.39404133595713E-15)*Y+
1434 1.445982921243E-13)*Y-2.676703245252E-12)*Y+
1435 5.823521627177E-12)*Y+2.17264723874381E-10 )*Y+
1436 3.56242145897468E-09 )*Y-3.03763737404491E-07 )*Y+
1437 9.46859114120901E-06 )*Y-2.30896753853196E-04 )*Y+
1438 5.24663913001114E-03;
1439 RT2 = (((((((((( 2.89872355524581E-16*Y-1.22296292045864E-14)*Y+
1440 6.184065097200E-14)*Y+1.649846591230E-12)*Y-
1441 2.729713905266E-11)*Y+3.709913790650E-11)*Y+
1442 2.216486288382E-09)*Y+4.616160236414E-08)*Y-
1443 3.32380270861364E-06 )*Y+9.84635072633776E-05 )*Y-
1444 2.30092118015697E-03 )*Y+5.00845183695073E-02;
1445 RT3 = (((((((((( 1.97068646590923E-15*Y-4.89419270626800E-14)*Y+
1446 1.136466605916E-13)*Y+7.546203883874E-12)*Y-
1447 9.635646767455E-11)*Y-8.295965491209E-11)*Y+
1448 7.534109114453E-09)*Y+2.699970652707E-07)*Y-
1449 1.42982334217081E-05 )*Y+3.78290946669264E-04 )*Y-
1450 8.03133015084373E-03 )*Y+1.58689469640791E-01;
1451 RT4 = (((((((((( 1.33642069941389E-14*Y-1.55850612605745E-13)*Y-
1452 7.522712577474E-13)*Y+3.209520801187E-11)*Y-
1453 2.075594313618E-10)*Y-2.070575894402E-09)*Y+
1454 7.323046997451E-09)*Y+1.851491550417E-06)*Y-
1455 6.37524802411383E-05 )*Y+1.36795464918785E-03 )*Y-
1456 2.42051126993146E-02 )*Y+3.97847167557815E-01;
1457 RT5 = ((((((((((-6.07053986130526E-14*Y+1.04447493138843E-12)*Y-
1458 4.286617818951E-13)*Y-2.632066100073E-10)*Y+
1459 4.804518986559E-09)*Y-1.835675889421E-08)*Y-
1460 1.068175391334E-06)*Y+3.292234974141E-05)*Y-
1461 5.94805357558251E-04 )*Y+8.29382168612791E-03 )*Y-
1462 9.93122509049447E-02 )*Y+1.09857804755042E+00;
1463 WW1 = (((((((((-9.10338640266542E-15*Y+1.00438927627833E-13)*Y+
1464 7.817349237071E-13)*Y-2.547619474232E-11)*Y+
1465 1.479321506529E-10)*Y+1.52314028857627E-09 )*Y+
1466 9.20072040917242E-09 )*Y-2.19427111221848E-06 )*Y+
1467 8.65797782880311E-05 )*Y-2.82718629312875E-03 )*Y+
1468 1.28718310443295E-01;
1469 WW2 = ((((((((( 5.52380927618760E-15*Y-6.43424400204124E-14)*Y-
1470 2.358734508092E-13)*Y+8.261326648131E-12)*Y+
1471 9.229645304956E-11)*Y-5.68108973828949E-09 )*Y+
1472 1.22477891136278E-07 )*Y-2.11919643127927E-06 )*Y+
1473 4.23605032368922E-05 )*Y-1.14423444576221E-03 )*Y+
1474 5.06607252890186E-02;
1475 WW3 = ((((((((( 3.99457454087556E-15*Y-5.11826702824182E-14)*Y-
1476 4.157593182747E-14)*Y+4.214670817758E-12)*Y+
1477 6.705582751532E-11)*Y-3.36086411698418E-09 )*Y+
1478 6.07453633298986E-08 )*Y-7.40736211041247E-07 )*Y+
1479 8.84176371665149E-06 )*Y-1.72559275066834E-04 )*Y+
1480 7.16639814253567E-03;
1481 WW4 = (((((((((((-2.14649508112234E-18*Y-2.45525846412281E-18)*Y+
1482 6.126212599772E-16)*Y-8.526651626939E-15)*Y+
1483 4.826636065733E-14)*Y-3.39554163649740E-13 )*Y+
1484 1.67070784862985E-11 )*Y-4.42671979311163E-10 )*Y+
1485 6.77368055908400E-09 )*Y-7.03520999708859E-08 )*Y+
1486 6.04993294708874E-07 )*Y-7.80555094280483E-06 )*Y+
1487 2.85954806605017E-04;
1488 WW5 = ((((((((((((-5.63938733073804E-21*Y+6.92182516324628E-20)*Y-
1489 1.586937691507E-18)*Y+3.357639744582E-17)*Y-
1490 4.810285046442E-16)*Y+5.386312669975E-15)*Y-
1491 6.117895297439E-14)*Y+8.441808227634E-13)*Y-
1492 1.18527596836592E-11 )*Y+1.36296870441445E-10 )*Y-
1493 1.17842611094141E-09 )*Y+7.80430641995926E-09 )*Y-
1494 5.97767417400540E-08 )*Y+1.65186146094969E-06;
1495 } else if (X < 40) {
1496 WW1 = sqrt(PIE4/X);
1497 E = exp(-X);
1498 RT1 = ((((((((-1.73363958895356E-06*X+1.19921331441483E-04)*X -
1499 1.59437614121125E-02)*X+1.13467897349442E+00)*X -
1500 4.47216460864586E+01)*X+1.06251216612604E+03)*X -
1501 1.52073917378512E+04)*X+1.20662887111273E+05)*X -
1502 4.07186366852475E+05)*E + R15/(X-R15);
1503 RT2 = ((((((((-1.60102542621710E-05*X+1.10331262112395E-03)*X -
1504 1.50043662589017E-01)*X+1.05563640866077E+01)*X -
1505 4.10468817024806E+02)*X+9.62604416506819E+03)*X -
1506 1.35888069838270E+05)*X+1.06107577038340E+06)*X -
1507 3.51190792816119E+06)*E + R25/(X-R25);
1508 RT3 = ((((((((-4.48880032128422E-05*X+2.69025112122177E-03)*X -
1509 4.01048115525954E-01)*X+2.78360021977405E+01)*X -
1510 1.04891729356965E+03)*X+2.36985942687423E+04)*X -
1511 3.19504627257548E+05)*X+2.34879693563358E+06)*X -
1512 7.16341568174085E+06)*E + R35/(X-R35);
1513 RT4 = ((((((((-6.38526371092582E-05*X-2.29263585792626E-03)*X -
1514 7.65735935499627E-02)*X+9.12692349152792E+00)*X -
1515 2.32077034386717E+02)*X+2.81839578728845E+02)*X +
1516 9.59529683876419E+04)*X-1.77638956809518E+06)*X +
1517 1.02489759645410E+07)*E + R45/(X-R45);
1518 RT5 = ((((((((-3.59049364231569E-05*X-2.25963977930044E-02)*X +
1519 1.12594870794668E+00)*X-4.56752462103909E+01)*X +
1520 1.05804526830637E+03)*X-1.16003199605875E+04)*X -
1521 4.07297627297272E+04)*X+2.22215528319857E+06)*X -
1522 1.61196455032613E+07)*E + R55/(X-R55);
1523 WW5 = (((((((((-4.61100906133970E-10*X+1.43069932644286E-07)*X -
1524 1.63960915431080E-05)*X+1.15791154612838E-03)*X -
1525 5.30573476742071E-02)*X+1.61156533367153E+00)*X -
1526 3.23248143316007E+01)*X+4.12007318109157E+02)*X -
1527 3.02260070158372E+03)*X+9.71575094154768E+03)*E + W55*WW1;
1528 WW4 = (((((((((-2.40799435809950E-08*X+8.12621667601546E-06)*X -
1529 9.04491430884113E-04)*X+6.37686375770059E-02)*X -
1530 2.96135703135647E+00)*X+9.15142356996330E+01)*X -
1531 1.86971865249111E+03)*X+2.42945528916947E+04)*X -
1532 1.81852473229081E+05)*X+5.96854758661427E+05)*E + W45*WW1;
1533 WW3 = (((((((( 1.83574464457207E-05*X-1.54837969489927E-03)*X +
1534 1.18520453711586E-01)*X-6.69649981309161E+00)*X +
1535 2.44789386487321E+02)*X-5.68832664556359E+03)*X +
1536 8.14507604229357E+04)*X-6.55181056671474E+05)*X +
1537 2.26410896607237E+06)*E + W35*WW1;
1538 WW2 = (((((((( 2.77778345870650E-05*X-2.22835017655890E-03)*X +
1539 1.61077633475573E-01)*X-8.96743743396132E+00)*X +
1540 3.28062687293374E+02)*X-7.65722701219557E+03)*X +
1541 1.10255055017664E+05)*X-8.92528122219324E+05)*X +
1542 3.10638627744347E+06)*E + W25*WW1;
1543 WW1 = WW1-0.01962E+00*E-WW2-WW3-WW4-WW5;
1544 } else if (X < 59.0) {
1545 WW1 = sqrt(PIE4/X);
1546 XXX = X * X * X;
1547 E = XXX*exp(-X);
1548 RT1 = (((-2.43758528330205E-02*X+2.07301567989771E+00)*X -
1549 6.45964225381113E+01)*X+7.14160088655470E+02)*E + R15/(X-R15);
1550 RT2 = (((-2.28861955413636E-01*X+1.93190784733691E+01)*X -
1551 5.99774730340912E+02)*X+6.61844165304871E+03)*E + R25/(X-R25);
1552 RT3 = (((-6.95053039285586E-01*X+5.76874090316016E+01)*X -
1553 1.77704143225520E+03)*X+1.95366082947811E+04)*E + R35/(X-R35);
1554 RT4 = (((-1.58072809087018E+00*X+1.27050801091948E+02)*X -
1555 3.86687350914280E+03)*X+4.23024828121420E+04)*E + R45/(X-R45);
1556 RT5 = (((-3.33963830405396E+00*X+2.51830424600204E+02)*X -
1557 7.57728527654961E+03)*X+8.21966816595690E+04)*E + R55/(X-R55);
1558 E = XXX*E;
1559 WW5 = (( 1.35482430510942E-08*X-3.27722199212781E-07)*X +
1560 2.41522703684296E-06)*E + W55*WW1;
1561 WW4 = (( 1.23464092261605E-06*X-3.55224564275590E-05)*X +
1562 3.03274662192286E-04)*E + W45*WW1;
1563 WW3 = (( 1.34547929260279E-05*X-4.19389884772726E-04)*X +
1564 3.87706687610809E-03)*E + W35*WW1;
1565 WW2 = (( 2.09539509123135E-05*X-6.87646614786982E-04)*X +
1566 6.68743788585688E-03)*E + W25*WW1;
1567 WW1 = WW1-WW2-WW3-WW4-WW5;
1568 } else {
1569 WW1 = sqrt(PIE4/X);
1570 RT1 = R15/(X-R15);
1571 RT2 = R25/(X-R25);
1572 RT3 = R35/(X-R35);
1573 RT4 = R45/(X-R45);
1574 RT5 = R55/(X-R55);
1575 WW2 = W25*WW1;
1576 WW3 = W35*WW1;
1577 WW4 = W45*WW1;
1578 WW5 = W55*WW1;
1579 WW1 = WW1-WW2-WW3-WW4-WW5;
1580 }
1581 roots[0] = RT1;
1582 roots[1] = RT2;
1583 roots[2] = RT3;
1584 roots[3] = RT4;
1585 roots[4] = RT5;
1586 weights[0] = WW1;
1587 weights[1] = WW2;
1588 weights[2] = WW3;
1589 weights[3] = WW4;
1590 weights[4] = WW5;
1591 }
1592
1593 #define POLYNOMIAL_VALUE1(p, a, order, x) \
1594 p = a[order]; \
1595 for (i = 1; i <= order; i++) { \
1596 p = p * x + a[order-i]; \
1597 }
1598
R_dnode(double * a,double * roots,int order)1599 static int R_dnode(double *a, double *roots, int order)
1600 {
1601 const double accrt = 1e-15;
1602 double x0, x1, xi, x1init, p0, p1, pi, p1init;
1603 int i, m, n;
1604
1605 x1init = 0;
1606 p1init = a[0];
1607 for (m = 0; m < order; ++m) {
1608 x0 = x1init;
1609 p0 = p1init;
1610 x1init = roots[m];
1611 POLYNOMIAL_VALUE1(p1init, a, order, x1init);
1612
1613 // When all coefficients a are 0, short-circuit the rest code to
1614 // ensure the roots from the lower order polynomials are preserved
1615 if (p1init == 0) {
1616 // roots[m] = x1init;
1617 continue;
1618 }
1619 if (p0 * p1init > 0) {
1620 fprintf(stderr, "ROOT NUMBER %d WAS NOT FOUND FOR POLYNOMIAL OF ORDER %d\n", m, order);
1621 return 1;
1622 }
1623 if (x0 <= x1init) {
1624 x1 = x1init;
1625 p1 = p1init;
1626 } else {
1627 x1 = x0;
1628 p1 = p0;
1629 x0 = x1init;
1630 p0 = p1init;
1631 }
1632 // interpolate/extrapolate between [x0,x1]
1633 if (p1 == 0) {
1634 roots[m] = x1;
1635 continue;
1636 } else if (p0 == 0) {
1637 roots[m] = x0;
1638 continue;
1639 } else {
1640 xi = x0 + (x0 - x1) / (p1 - p0) * p0;
1641 }
1642 n = 0;
1643 while (x1 > accrt+x0 || x0 > x1+accrt) {
1644 n++;
1645 if (n > 200) {
1646 fprintf(stderr, "libcint::rys_roots NO CONV. IN R_dnode\n");
1647 return 1;
1648 }
1649 POLYNOMIAL_VALUE1(pi, a, order, xi);
1650 if (pi == 0) {
1651 break;
1652 } else if (p0 * pi <= 0) {
1653 x1 = xi;
1654 p1 = pi;
1655 xi = x0 * .25 + xi * .75;
1656 } else {
1657 x0 = xi;
1658 p0 = pi;
1659 xi = xi * .75 + x1 * .25;
1660 }
1661 POLYNOMIAL_VALUE1(pi, a, order, xi);
1662 if (pi == 0) {
1663 break;
1664 } else if (p0 * pi <= 0) {
1665 x1 = xi;
1666 p1 = pi;
1667 } else {
1668 x0 = xi;
1669 p0 = pi;
1670 }
1671
1672 xi = x0 + (x0 - x1) / (p1 - p0) * p0;
1673 }
1674 roots[m] = xi;
1675 }
1676 return 0;
1677 }
1678
1679 #define SET_ZERO(a, n, start) \
1680 for (k = start; k < n; ++k) { \
1681 for (i = 0; i < n; ++i) { \
1682 a[i + k * n] = 0; \
1683 } \
1684 } \
1685
R_dsmit(double * cs,double * fmt_ints,int n)1686 static int R_dsmit(double *cs, double *fmt_ints, int n)
1687 {
1688 int i, j, k;
1689 double fac, dot, tmp;
1690 double v[MXRYSROOTS];
1691
1692 fac = -fmt_ints[1] / fmt_ints[0];
1693 tmp = fmt_ints[2] + fac * fmt_ints[1];
1694 if (tmp <= 0) {
1695 fprintf(stderr, "libcint::rys_roots negative value in sqrt for roots %d (j=1)\n", n-1);
1696 SET_ZERO(cs, n, 1);
1697 return 1;
1698 }
1699 tmp = 1 / sqrt(tmp);
1700 cs[0+0*n] = 1 / sqrt(fmt_ints[0]);
1701 cs[0+1*n] = fac * tmp;
1702 cs[1+1*n] = tmp;
1703
1704 for (j = 2; j < n; ++j) {
1705 for (k = 0; k < j; ++k) {
1706 v[k] = 0;
1707 }
1708 fac = fmt_ints[j + j];
1709 for (k = 0; k < j; ++k) {
1710 dot = 0;
1711 for (i = 0; i <= k; ++i) {
1712 dot += cs[i + k * n] * fmt_ints[i+j];
1713 }
1714 for (i = 0; i <= k; ++i) {
1715 v[i] -= dot * cs[i + k * n];
1716 }
1717 fac -= dot * dot;
1718 }
1719
1720 if (fac <= 0) {
1721 fprintf(stderr, "libcint::rys_roots negative value in sqrt for roots %d (j=%d)\n", n-1, j);
1722 // set rest coefficients to 0
1723 SET_ZERO(cs, n, j);
1724 return j;
1725 }
1726 fac = 1 / sqrt(fac);
1727 cs[j + j * n] = fac;
1728 for (k = 0; k < j; ++k) {
1729 cs[k + j * n] = fac * v[k];
1730 }
1731 }
1732 return 0;
1733 }
1734
1735 /*
1736 * Using RDK algorithm (based on Schmidt orthogonalization to search rys roots
1737 * and weights
1738 */
_rdk_rys_roots(int nroots,double * fmt_ints,double * roots,double * weights)1739 static int _rdk_rys_roots(int nroots, double *fmt_ints,
1740 double *roots, double *weights)
1741 {
1742 int i, k, j, order;
1743 int nroots1 = nroots + 1;
1744 double rt[MXRYSROOTS + MXRYSROOTS * MXRYSROOTS];
1745 double *cs = rt + nroots1;
1746 double *a;
1747 double root, poly, dum;
1748
1749 // to avoid numerical instability for very small fmt integrals
1750 if (fmt_ints[0] == 0) {
1751 for (k = 0; k < nroots; ++k) {
1752 roots[k] = 0;
1753 weights[k] = 0;
1754 }
1755 return 0;
1756 }
1757 if (nroots == 1) {
1758 roots[0] = fmt_ints[1] / (fmt_ints[0] - fmt_ints[1]);
1759 weights[0] = fmt_ints[0];
1760 return 0;
1761 }
1762
1763 int error = R_dsmit(cs, fmt_ints, nroots1);
1764 #ifndef KEEP_GOING
1765 if (error) {
1766 exit(error);
1767 }
1768 #else
1769 if (error == 1) {
1770 return 1;
1771 }
1772 #endif
1773
1774 dum = sqrt(cs[2*nroots1+1] * cs[2*nroots1+1] - 4 * cs[2*nroots1+0] * cs[2*nroots1+2]);
1775 rt[0] = .5 * (-cs[2*nroots1+1] - dum) / cs[2*nroots1+2];
1776 rt[1] = .5 * (-cs[2*nroots1+1] + dum) / cs[2*nroots1+2];
1777 for (i = 2; i < nroots; i++) {
1778 rt[i] = 1;
1779 }
1780
1781 for (k = 2; k < nroots; ++k) {
1782 order = k + 1;
1783 a = cs + order * nroots1;
1784 error = R_dnode(a, rt, order);
1785 if (error) {
1786 #ifndef KEEP_GOING
1787 exit(error);
1788 #else
1789 return error;
1790 #endif
1791 }
1792 }
1793
1794 for (k = 0; k < nroots; ++k) {
1795 root = rt[k];
1796 // When singularity was caught in R_dsmit, they are typically
1797 // caused by high order Rys polynomials. We assume the contributions
1798 // from these high order Rys polynomials are negligible. Only the
1799 // roots obtained from lower polynomials are used.
1800 if (root == 1) {
1801 roots[k] = 0;
1802 weights[k] = 0;
1803 continue;
1804 }
1805
1806 dum = 1 / fmt_ints[0];
1807 for (j = 1; j < nroots; ++j) {
1808 order = j;
1809 a = cs + j * nroots1;
1810 // poly = poly_value1(cs[:order+1,j], order, root);
1811 POLYNOMIAL_VALUE1(poly, a, order, root);
1812 dum += poly * poly;
1813 }
1814 roots[k] = root / (1 - root);
1815 weights[k] = 1 / dum;
1816 }
1817 return 0;
1818 }
1819
CINTrys_schmidt(int nroots,double x,double lower,double * roots,double * weights)1820 int CINTrys_schmidt(int nroots, double x, double lower, double *roots, double *weights)
1821 {
1822 double fmt_ints[MXRYSROOTS*2];
1823 if (lower == 0) {
1824 gamma_inc_like(fmt_ints, x, nroots*2);
1825 } else {
1826 fmt_erfc_like(fmt_ints, x, lower, nroots*2);
1827 }
1828 return _rdk_rys_roots(nroots, fmt_ints, roots, weights);
1829 }
1830
1831 /*
1832 ******************************************************
1833 * 80 bit double
1834 */
1835 #ifdef HAVE_SQRTL
1836 #define SQRTL sqrtl
1837 #else
c99_sqrtl(long double x)1838 static long double c99_sqrtl(long double x)
1839 {
1840 long double z = sqrt(x);
1841 // ref. Babylonian method
1842 // http://en.wikipedia.org/wiki/Methods_of_computing_square_roots
1843 // An extra update should be enough due to the quadratic convergence
1844 return (z*z + x)/(z * 2);
1845 }
1846 #define SQRTL c99_sqrtl
1847 #endif
1848
1849 #ifdef HAVE_EXPL
1850 #define EXPL expl
1851 #else
1852 // Does it need to swith to 128 bit expl algorithm?
1853 // ref https://github.com/JuliaLang/openlibm/ld128/e_expl.c
c99_expl(long double x)1854 static long double c99_expl(long double x)
1855 {
1856 return exp(x);
1857 }
1858 #define EXPL c99_expl
1859 #endif
1860
R_lsmit(long double * cs,long double * fmt_ints,int n)1861 static int R_lsmit(long double *cs, long double *fmt_ints, int n)
1862 {
1863 int i, j, k;
1864 long double fac, dot, tmp;
1865 long double v[MXRYSROOTS];
1866
1867 fac = -fmt_ints[1] / fmt_ints[0];
1868 tmp = fmt_ints[2] + fac * fmt_ints[1];
1869 if (tmp <= 0) {
1870 fprintf(stderr, "libcint::rys_roots negative value in sqrtl for roots %d (j=1)\n", n-1);
1871 SET_ZERO(cs, n, 1);
1872 return 1;
1873 }
1874 tmp = 1 / SQRTL(tmp);
1875 cs[0+0*n] = 1 / SQRTL(fmt_ints[0]);
1876 cs[0+1*n] = fac * tmp;
1877 cs[1+1*n] = tmp;
1878
1879 for (j = 2; j < n; ++j) {
1880 for (k = 0; k < j; ++k) {
1881 v[k] = 0;
1882 }
1883 fac = fmt_ints[j + j];
1884 for (k = 0; k < j; ++k) {
1885 dot = 0;
1886 for (i = 0; i <= k; ++i) {
1887 dot += cs[i + k * n] * fmt_ints[i+j];
1888 }
1889 for (i = 0; i <= k; ++i) {
1890 v[i] -= dot * cs[i + k * n];
1891 }
1892 fac -= dot * dot;
1893 }
1894
1895 if (fac <= 0) {
1896 fprintf(stderr, "libcint::rys_roots negative value in sqrtl for roots %d (j=%d)\n", n-1, j);
1897 // set rest coefficients to 0
1898 SET_ZERO(cs, n, j);
1899 return j;
1900 }
1901 fac = 1 / SQRTL(fac);
1902 cs[j + j * n] = fac;
1903 for (k = 0; k < j; ++k) {
1904 cs[k + j * n] = fac * v[k];
1905 }
1906 }
1907 return 0;
1908 }
1909
CINTlrys_schmidt(int nroots,double x,double lower,double * roots,double * weights)1910 int CINTlrys_schmidt(int nroots, double x, double lower, double *roots, double *weights)
1911 {
1912 int i, k, j, order, error;
1913 int nroots1 = nroots + 1;
1914 long double fmt_ints[MXRYSROOTS * 2 + MXRYSROOTS * MXRYSROOTS];
1915 long double *qcs = fmt_ints + nroots1 * 2;
1916 double rt[MXRYSROOTS + MXRYSROOTS * MXRYSROOTS];
1917 double *cs = rt + nroots;
1918 double *a;
1919 double root, poly, dum, dum0;
1920
1921 if (lower == 0) {
1922 lgamma_inc_like(fmt_ints, x, nroots*2);
1923 } else {
1924 fmt_lerfc_like(fmt_ints, x, lower, nroots*2);
1925 }
1926
1927 if (fmt_ints[0] == 0) {
1928 for (k = 0; k < nroots; ++k) {
1929 roots[k] = 0;
1930 weights[k] = 0;
1931 }
1932 return 0;
1933 }
1934
1935 if (nroots == 1) {
1936 rt[0] = fmt_ints[1] / fmt_ints[0];
1937 } else {
1938 error = R_lsmit(qcs, fmt_ints, nroots1);
1939 #ifndef KEEP_GOING
1940 if (error) {
1941 exit(error);
1942 }
1943 #else
1944 if (error == 1) {
1945 return 1;
1946 }
1947 #endif
1948 for (k = 1; k < nroots1; k++) {
1949 for (i = 0; i <= k; i++) {
1950 cs[k * nroots1 + i] = qcs[k * nroots1 + i];
1951 }
1952 }
1953 dum = sqrt(cs[2*nroots1+1] * cs[2*nroots1+1] - 4 * cs[2*nroots1+0] * cs[2*nroots1+2]);
1954 rt[0] = .5 * (-cs[2*nroots1+1] - dum) / cs[2*nroots1+2];
1955 rt[1] = .5 * (-cs[2*nroots1+1] + dum) / cs[2*nroots1+2];
1956 for (i = 2; i < nroots; i++) {
1957 rt[i] = 1;
1958 }
1959 }
1960
1961 for (k = 2; k < nroots; ++k) {
1962 order = k + 1;
1963 a = cs + order * nroots1;
1964 error = R_dnode(a, rt, order);
1965 if (error) {
1966 #ifndef KEEP_GOING
1967 exit(error);
1968 #else
1969 return error;
1970 #endif
1971 }
1972 }
1973
1974 dum0 = 1 / fmt_ints[0];
1975 for (k = 0; k < nroots; ++k) {
1976 root = rt[k];
1977 if (root == 1) {
1978 roots[k] = 0;
1979 weights[k] = 0;
1980 continue;
1981 }
1982
1983 dum = dum0;
1984 for (j = 1; j < nroots; ++j) {
1985 order = j;
1986 a = cs + j * nroots1;
1987 // poly = poly_value1(cs[:order+1,j], order, root);
1988 POLYNOMIAL_VALUE1(poly, a, order, root);
1989 dum += poly * poly;
1990 }
1991 roots[k] = root / (1 - root);
1992 weights[k] = 1 / dum;
1993 }
1994 return 0;
1995 }
1996
1997 #ifdef HAVE_QUADMATH_H
R_qsmit(__float128 * cs,__float128 * fmt_ints,int n)1998 static int R_qsmit(__float128 *cs, __float128 *fmt_ints, int n)
1999 {
2000 int i, j, k;
2001 __float128 fac, dot, tmp;
2002 __float128 v[MXRYSROOTS];
2003
2004 fac = -fmt_ints[1] / fmt_ints[0];
2005 tmp = fmt_ints[2] + fac * fmt_ints[1];
2006 if (tmp <= 0) {
2007 fprintf(stderr, "libcint::rys_roots negative value in sqrtq for roots %d (j=1)\n", n-1);
2008 SET_ZERO(cs, n, 1);
2009 return 1;
2010 }
2011 tmp = 1 / sqrtq(tmp);
2012 cs[0+0*n] = 1 / sqrtq(fmt_ints[0]);
2013 cs[0+1*n] = fac * tmp;
2014 cs[1+1*n] = tmp;
2015
2016 for (j = 2; j < n; ++j) {
2017 for (k = 0; k < j; ++k) {
2018 v[k] = 0;
2019 }
2020 fac = fmt_ints[j + j];
2021 for (k = 0; k < j; ++k) {
2022 dot = 0;
2023 for (i = 0; i <= k; ++i) {
2024 dot += cs[i + k * n] * fmt_ints[i+j];
2025 }
2026 for (i = 0; i <= k; ++i) {
2027 v[i] -= dot * cs[i + k * n];
2028 }
2029 fac -= dot * dot;
2030 }
2031
2032 if (fac <= 0) {
2033 fprintf(stderr, "libcint::rys_roots negative value in sqrtq for roots %d (j=%d)\n", n-1, j);
2034 // set rest coefficients to 0
2035 SET_ZERO(cs, n, j);
2036 return j;
2037 }
2038 fac = 1.q / sqrtq(fac);
2039 cs[j + j * n] = fac;
2040 for (k = 0; k < j; ++k) {
2041 cs[k + j * n] = fac * v[k];
2042 }
2043 }
2044 return 0;
2045 }
2046
CINTqrys_schmidt(int nroots,double x,double lower,double * roots,double * weights)2047 int CINTqrys_schmidt(int nroots, double x, double lower, double *roots, double *weights)
2048 {
2049 int i, k, j, order, error;
2050 int nroots1 = nroots + 1;
2051 __float128 fmt_ints[MXRYSROOTS * 2 + MXRYSROOTS * MXRYSROOTS];
2052 __float128 *qcs = fmt_ints + nroots1 * 2;
2053 double rt[MXRYSROOTS + MXRYSROOTS * MXRYSROOTS];
2054 double *cs = rt + nroots;
2055 double *a;
2056 double root, poly, dum, dum0;
2057
2058 if (lower == 0) {
2059 qgamma_inc_like(fmt_ints, x, nroots*2);
2060 } else {
2061 fmt_qerfc_like(fmt_ints, x, lower, nroots*2);
2062 }
2063
2064 if (fmt_ints[0] == 0) {
2065 for (k = 0; k < nroots; ++k) {
2066 roots[k] = 0;
2067 weights[k] = 0;
2068 }
2069 return 0;
2070 }
2071
2072 if (nroots == 1) {
2073 rt[0] = fmt_ints[1] / fmt_ints[0];
2074 } else {
2075 error = R_qsmit(qcs, fmt_ints, nroots1);
2076 #ifndef KEEP_GOING
2077 if (error) {
2078 exit(error);
2079 }
2080 #else
2081 if (error == 1) {
2082 return 1;
2083 }
2084 #endif
2085 for (k = 1; k < nroots1; k++) {
2086 for (i = 0; i <= k; i++) {
2087 cs[k * nroots1 + i] = qcs[k * nroots1 + i];
2088 }
2089 }
2090 dum = sqrt(cs[2*nroots1+1] * cs[2*nroots1+1] - 4 * cs[2*nroots1+0] * cs[2*nroots1+2]);
2091 rt[0] = .5 * (-cs[2*nroots1+1] - dum) / cs[2*nroots1+2];
2092 rt[1] = .5 * (-cs[2*nroots1+1] + dum) / cs[2*nroots1+2];
2093 for (i = 2; i < nroots; i++) {
2094 rt[i] = 1;
2095 }
2096 }
2097
2098 for (k = 2; k < nroots; ++k) {
2099 order = k + 1;
2100 a = cs + order * nroots1;
2101 error = R_dnode(a, rt, order);
2102 if (error) {
2103 #ifndef KEEP_GOING
2104 exit(error);
2105 #else
2106 return error;
2107 #endif
2108 }
2109 }
2110
2111 dum0 = 1 / fmt_ints[0];
2112 for (k = 0; k < nroots; ++k) {
2113 root = rt[k];
2114 if (root == 1) {
2115 roots[k] = 0;
2116 weights[k] = 0;
2117 continue;
2118 }
2119
2120 dum = dum0;
2121 for (j = 1; j < nroots; ++j) {
2122 order = j;
2123 a = cs + j * nroots1;
2124 // poly = poly_value1(cs[:order+1,j], order, root);
2125 POLYNOMIAL_VALUE1(poly, a, order, root);
2126 dum += poly * poly;
2127 }
2128 roots[k] = root / (1 - root);
2129 weights[k] = 1 / dum;
2130 }
2131 return 0;
2132 }
2133
2134 #endif
2135