Line data Source code
1 : /*----------------------------------------------------------------------------*/
2 : /* CP2K: A general program to perform molecular dynamics simulations */
3 : /* Copyright 2000-2025 CP2K developers group <https://cp2k.org> */
4 : /* */
5 : /* SPDX-License-Identifier: MIT */
6 : /*----------------------------------------------------------------------------*/
7 :
8 : /*
9 : * libgrpp - a library for the evaluation of integrals over
10 : * generalized relativistic pseudopotentials.
11 : *
12 : * Copyright (C) 2021-2023 Alexander Oleynichenko
13 : */
14 :
15 : #include <float.h> // required for LDBL_EPSILON
16 : #include <math.h> // required for fabsl()
17 : #ifndef M_PI
18 : #define M_PI 3.1415926535897932384626433
19 : #endif
20 :
21 : #include "grpp_specfunc.h"
22 :
23 : /*
24 : * This code is taken from the website:
25 : * http://www.mymathlib.com/functions/dawsons_integral.html
26 : */
27 :
28 : ////////////////////////////////////////////////////////////////////////////////
29 : // File: specfunc_dawson.c //
30 : // Routine(s): //
31 : // Dawsons_Integral //
32 : // xDawsons_Integral //
33 : ////////////////////////////////////////////////////////////////////////////////
34 : ////////////////////////////////////////////////////////////////////////////////
35 : // Description: //
36 : // Dawson's integral, Daw(x), is the integral from 0 to x of the //
37 : // integrand: //
38 : // exp(-x^2) exp(t^2) dt. //
39 : // I.e. //
40 : // Daw(x) = exp(-x^2) * I[0,x] (exp(t^2)) dt, //
41 : // where I[0,x] indicates the integral from 0 to x. //
42 : ////////////////////////////////////////////////////////////////////////////////
43 : #define Asymptotic_Expansion_Cutoff 50
44 :
45 : // Externally Defined Routines //
46 : static long double xChebyshev_Tn_Series(long double x, long double a[],
47 : int degree);
48 :
49 : // Internally Defined Routines //
50 : double libgrpp_Dawsons_Integral(double x);
51 : long double xDawsons_Integral(long double x);
52 :
53 : static long double Dawson_Power_Series(long double x);
54 : static long double Dawson_Chebyshev_Expansion_1_175(long double x);
55 : static long double Dawson_Chebyshev_Expansion_175_250(long double x);
56 : static long double Dawson_Chebyshev_Expansion_250_325(long double x);
57 : static long double Dawson_Chebyshev_Expansion_325_425(long double x);
58 : static long double Dawson_Chebyshev_Expansion_425_550(long double x);
59 : static long double Dawson_Chebyshev_Expansion_550_725(long double x);
60 : static long double Dawson_Asymptotic_Expansion(long double x);
61 :
62 : ////////////////////////////////////////////////////////////////////////////////
63 : // double libgrpp_Dawsons_Integral( double x ) //
64 : // //
65 : // Description: //
66 : // Dawson's integral, Daw(x), is the integral with integrand //
67 : // exp(-x^2) exp(t^2) dt //
68 : // where the integral extends from 0 to x. //
69 : // //
70 : // Arguments: //
71 : // double x The argument of Dawson's integral Daw(). //
72 : // //
73 : // Return Value: //
74 : // The value of Dawson's integral, Daw(), evaluated at x. //
75 : // //
76 : // Example: //
77 : // double y, x; //
78 : // //
79 : // ( code to initialize x ) //
80 : // //
81 : // y = Dawsons_Integral( x ); //
82 : ////////////////////////////////////////////////////////////////////////////////
83 0 : double libgrpp_Dawsons_Integral(double x) {
84 0 : return (double)xDawsons_Integral((long double)x);
85 : }
86 :
87 : ////////////////////////////////////////////////////////////////////////////////
88 : // long double xDawsons_Integral( long double x ) //
89 : // //
90 : // Description: //
91 : // Dawson's integral, Daw(x), is the integral with integrand //
92 : // exp(-x^2) exp(t^2) dt //
93 : // where the integral extends from 0 to x. //
94 : // //
95 : // Arguments: //
96 : // long double x The argument of Dawson's integral Daw(). //
97 : // //
98 : // Return Value: //
99 : // The value of Dawson's integral, Daw(), evaluated at x. //
100 : // //
101 : // Example: //
102 : // long double y, x; //
103 : // //
104 : // ( code to initialize x ) //
105 : // //
106 : // y = xDawsons_Integral( x ); //
107 : ////////////////////////////////////////////////////////////////////////////////
108 :
109 0 : long double xDawsons_Integral(long double x) {
110 0 : long double abs_x = fabsl(x);
111 :
112 0 : if (abs_x <= 1.0L)
113 0 : return Dawson_Power_Series(x);
114 0 : if (abs_x <= 1.75L)
115 0 : return Dawson_Chebyshev_Expansion_1_175(x);
116 0 : if (abs_x <= 2.50L)
117 0 : return Dawson_Chebyshev_Expansion_175_250(x);
118 0 : if (abs_x <= 3.25L)
119 0 : return Dawson_Chebyshev_Expansion_250_325(x);
120 0 : if (abs_x <= 4.25L)
121 0 : return Dawson_Chebyshev_Expansion_325_425(x);
122 0 : if (abs_x <= 5.50L)
123 0 : return Dawson_Chebyshev_Expansion_425_550(x);
124 0 : if (abs_x <= 7.25L)
125 0 : return Dawson_Chebyshev_Expansion_550_725(x);
126 0 : return Dawson_Asymptotic_Expansion(x);
127 : }
128 :
129 : ////////////////////////////////////////////////////////////////////////////////
130 : // static long double Dawson_Power_Series( long double x ) //
131 : // //
132 : // Description: //
133 : // Evaluate Dawsons integral for -1 <= x <= 1 using the power series //
134 : // expansion: //
135 : // Daw(x) = x Sum [(-2x^2)^j / (2j + 1)!!, //
136 : // where the sum extends over j = 0, ... . //
137 : // //
138 : // Arguments: //
139 : // long double x The argument of Dawson's integral Daw(), where //
140 : // where |x| <= 1. //
141 : // //
142 : // Return Value: //
143 : // The value of Dawson's integral, Daw(), evaluated at x where |x| <= 1. //
144 : // //
145 : // Example: //
146 : // long double y, x; //
147 : // //
148 : // ( code to initialize x ) //
149 : // //
150 : // y = Dawson_Power_Series(x); //
151 : ////////////////////////////////////////////////////////////////////////////////
152 :
153 0 : static long double Dawson_Power_Series(long double x) {
154 0 : long double two_x2 = -2.0L * x * x;
155 0 : long double sum = 0.0;
156 0 : long double term = 1.0L;
157 0 : long double factorial = 1.0L;
158 0 : long double xn = 1.0L;
159 0 : const long double epsilon = LDBL_EPSILON / 2.0L;
160 0 : int y = 0;
161 :
162 0 : if (x == 0.0L)
163 : return 0.0L;
164 0 : do {
165 0 : sum += term;
166 0 : y += 1;
167 0 : factorial *= (long double)(y + y + 1);
168 0 : xn *= two_x2;
169 0 : term = xn / factorial;
170 0 : } while (fabsl(term) > epsilon * fabsl(sum));
171 0 : return x * sum;
172 : }
173 :
174 : ////////////////////////////////////////////////////////////////////////////////
175 : // static long double Dawson_Chebyshev_Expansion_1_175( long double x ) //
176 : // //
177 : // Description: //
178 : // Evaluate Dawsons integral for 1 <= |x| <= 1.75 using the Chebyshev //
179 : // expansion of Daw(x) for x in the interval [1,1.75]. //
180 : // //
181 : // Arguments: //
182 : // long double x The argument of Dawson's integral Daw(), where //
183 : // where 1 <= |x| <= 1.75. //
184 : // //
185 : // Return Value: //
186 : // The value of Dawson's integral, Daw(), evaluated at x where //
187 : // 1 <= x <= 1.75. //
188 : // //
189 : // Example: //
190 : // long double y, x; //
191 : // //
192 : // ( code to initialize x ) //
193 : // //
194 : // y = Dawson_Chebyshev_Expansion_1_175(x); //
195 : ////////////////////////////////////////////////////////////////////////////////
196 :
197 0 : static long double Dawson_Chebyshev_Expansion_1_175(long double x) {
198 0 : static long double c[] = {
199 : +4.563960711239483142081e-1L, -9.268566100670767619861e-2L,
200 : -7.334392170021420220239e-3L, +3.379523740404396755124e-3L,
201 : -3.085898448678595090813e-4L, -1.519846724619319512311e-5L,
202 : +4.903955822454009397182e-6L, -2.106910538629224721838e-7L,
203 : -2.930676220603996192089e-8L, +3.326790071774057337673e-9L,
204 : +3.335593457695769191326e-11L, -2.279104036721012221982e-11L,
205 : +7.877561633156348806091e-13L, +9.173158167107974472228e-14L,
206 : -7.341175636102869400671e-15L, -1.763370444125849029511e-16L,
207 : +3.792946298506435014290e-17L, -4.251969162435936250171e-19L,
208 : -1.358295820818448686821e-19L, +5.268740962820224108235e-21L,
209 : +3.414939674304748094484e-22L};
210 :
211 0 : static const int degree = sizeof(c) / sizeof(long double) - 1;
212 0 : static const long double midpoint = (1.75L + 1.0L) / 2.0L;
213 0 : static const long double half_length = (1.75L - 1.0L) / 2.0L;
214 0 : long double daw =
215 0 : xChebyshev_Tn_Series((fabsl(x) - midpoint) / half_length, c, degree);
216 0 : return (x > 0.0L) ? daw : -daw;
217 : }
218 :
219 : ////////////////////////////////////////////////////////////////////////////////
220 : // static long double Dawson_Chebyshev_Expansion_175_250( long double x ) //
221 : // //
222 : // Description: //
223 : // Evaluate Dawsons integral for 1.75 <= |x| <= 2.50 using the Chebyshev //
224 : // expansion of Daw(x) for x in the interval [1.75,2.50]. //
225 : // //
226 : // Arguments: //
227 : // long double x The argument of Dawson's integral Daw(), where //
228 : // where 1.75 <= |x| <= 2.50. //
229 : // //
230 : // Return Value: //
231 : // The value of Dawson's integral, Daw(), evaluated at x where //
232 : // 1.75 <= x <= 2.50. //
233 : // //
234 : // Example: //
235 : // long double y, x; //
236 : // //
237 : // ( code to initialize x ) //
238 : // //
239 : // y = Dawson_Chebyshev_Expansion_175_250(x); //
240 : ////////////////////////////////////////////////////////////////////////////////
241 :
242 0 : static long double Dawson_Chebyshev_Expansion_175_250(long double x) {
243 0 : static long double c[] = {
244 : +2.843711194548592808550e-1L, -6.791774139166808940530e-2L,
245 : +6.955211059379384327814e-3L, -2.726366582146839486784e-4L,
246 : -6.516682485087925163874e-5L, +1.404387911504935155228e-5L,
247 : -1.103288540946056915318e-6L, -1.422154597293404846081e-8L,
248 : +1.102714664312839585330e-8L, -8.659211557383544255053e-10L,
249 : -8.048589443963965285748e-12L, +6.092061709996351761426e-12L,
250 : -3.580977611213519234324e-13L, -1.085173558590137965737e-14L,
251 : +2.411707924175380740802e-15L, -7.760751294610276598631e-17L,
252 : -6.701490147030045891595e-18L, +6.350145841254563572100e-19L,
253 : -2.034625734538917052251e-21L, -2.260543651146274653910e-21L,
254 : +9.782419961387425633151e-23L};
255 :
256 0 : static const int degree = sizeof(c) / sizeof(long double) - 1;
257 0 : static const long double midpoint = (2.50L + 1.75L) / 2.0L;
258 0 : static const long double half_length = (2.50L - 1.75L) / 2.0;
259 0 : long double daw =
260 0 : xChebyshev_Tn_Series((fabsl(x) - midpoint) / half_length, c, degree);
261 0 : return (x > 0.0L) ? daw : -daw;
262 : }
263 :
264 : ////////////////////////////////////////////////////////////////////////////////
265 : // static long double Dawson_Chebyshev_Expansion_250_325( long double x ) //
266 : // //
267 : // Description: //
268 : // Evaluate Dawsons integral for 2.5 <= |x| <= 3.25 using the Chebyshev //
269 : // expansion of Daw(x) for x in the interval [2.50,3.25]. //
270 : // //
271 : // Arguments: //
272 : // long double x The argument of Dawson's integral Daw(), where //
273 : // where 2.50 <= |x| <= 3.25. //
274 : // //
275 : // Return Value: //
276 : // The value of Dawson's integral, Daw(), evaluated at x where //
277 : // 2.50 <= x <= 3.25. //
278 : // //
279 : // Example: //
280 : // long double y, x; //
281 : // //
282 : // ( code to initialize x ) //
283 : // //
284 : // y = Dawson_Chebyshev_Expansion_250_325(x); //
285 : ////////////////////////////////////////////////////////////////////////////////
286 :
287 0 : static long double Dawson_Chebyshev_Expansion_250_325(long double x) {
288 0 : static long double c[] = {
289 : +1.901351274204578126827e-1L, -3.000575522193632460118e-2L,
290 : +2.672138524890489432579e-3L, -2.498237548675235150519e-4L,
291 : +2.013483163459701593271e-5L, -8.454663603108548182962e-7L,
292 : -8.036589636334016432368e-8L, +2.055498509671357933537e-8L,
293 : -2.052151324060186596995e-9L, +8.584315967075483822464e-11L,
294 : +5.062689357469596748991e-12L, -1.038671167196342609090e-12L,
295 : +6.367962851860231236238e-14L, +3.084688422647419767229e-16L,
296 : -3.417946142546575188490e-16L, +2.311567730100119302160e-17L,
297 : -6.170132546983726244716e-20L, -9.133176920944950460847e-20L,
298 : +5.712092431423316128728e-21L, +1.269641078369737220790e-23L,
299 : -2.072659711527711312699e-23L};
300 :
301 0 : static const int degree = sizeof(c) / sizeof(long double) - 1;
302 0 : static const long double midpoint = (3.25L + 2.50L) / 2.0L;
303 0 : static const long double half_length = (3.25L - 2.50L) / 2.0L;
304 0 : long double daw =
305 0 : xChebyshev_Tn_Series((fabsl(x) - midpoint) / half_length, c, degree);
306 0 : return (x > 0.0L) ? daw : -daw;
307 : }
308 :
309 : ////////////////////////////////////////////////////////////////////////////////
310 : // static long double Dawson_Chebyshev_Expansion_325_425( long double x ) //
311 : // //
312 : // Description: //
313 : // Evaluate Dawsons integral for 3.25 <= |x| <= 4.25 using the Chebyshev //
314 : // expansion of Daw(x) for x in the interval [3.25,4.75]. //
315 : // //
316 : // Arguments: //
317 : // long double x The argument of Dawson's integral Daw(), where //
318 : // where 3.25 <= |x| <= 4.25. //
319 : // //
320 : // Return Value: //
321 : // The value of Dawson's integral, Daw(), evaluated at x where //
322 : // 3.25 <= x <= 4.25. //
323 : // //
324 : // Example: //
325 : // long double y, x; //
326 : // //
327 : // ( code to initialize x ) //
328 : // //
329 : // y = Dawson_Chebyshev_Expansion_325_425(x); //
330 : ////////////////////////////////////////////////////////////////////////////////
331 :
332 0 : static long double Dawson_Chebyshev_Expansion_325_425(long double x) {
333 0 : static long double c[] = {
334 : +1.402884974484995678749e-1L, -2.053975371995937033959e-2L,
335 : +1.595388628922920119352e-3L, -1.336894584910985998203e-4L,
336 : +1.224903774178156286300e-5L, -1.206856028658387948773e-6L,
337 : +1.187997233269528945503e-7L, -1.012936061496824448259e-8L,
338 : +5.244408240062370605664e-10L, +2.901444759022254846562e-11L,
339 : -1.168987502493903926906e-11L, +1.640096995420504465839e-12L,
340 : -1.339190668554209618318e-13L, +3.643815972666851044790e-15L,
341 : +6.922486581126169160232e-16L, -1.158761251467106749752e-16L,
342 : +8.164320395639210093180e-18L, -5.397918405779863087588e-20L,
343 : -5.052069908100339242896e-20L, +5.322512674746973445361e-21L,
344 : -1.869294542789169825747e-22L};
345 :
346 0 : static const int degree = sizeof(c) / sizeof(long double) - 1;
347 : // static const long double lower_bound = 3.25L;
348 : // static const long double upper_bound = 4.25L;
349 0 : static const long double midpoint = (4.25L + 3.25L) / 2.0L;
350 0 : static const long double half_length = (4.25L - 3.25L) / 2.0L;
351 0 : long double daw =
352 0 : xChebyshev_Tn_Series((fabsl(x) - midpoint) / half_length, c, degree);
353 0 : return (x > 0.0L) ? daw : -daw;
354 : }
355 :
356 : ////////////////////////////////////////////////////////////////////////////////
357 : // static long double Dawson_Chebyshev_Expansion_425_550( long double x ) //
358 : // //
359 : // Description: //
360 : // Evaluate Dawsons integral for 4.25 <= |x| <= 5.50 using the Chebyshev //
361 : // expansion of Daw(x) for x in the interval [4.25,5.50]. //
362 : // //
363 : // Arguments: //
364 : // long double x The argument of Dawson's integral Daw(), where //
365 : // where 4.25 <= |x| <= 5.50. //
366 : // //
367 : // Return Value: //
368 : // The value of Dawson's integral, Daw(), evaluated at x where //
369 : // 4.25 <= x <= 5.50. //
370 : // //
371 : // Example: //
372 : // long double y, x; //
373 : // //
374 : // ( code to initialize x ) //
375 : // //
376 : // y = Dawson_Chebyshev_Expansion_425_550(x); //
377 : ////////////////////////////////////////////////////////////////////////////////
378 :
379 0 : static long double Dawson_Chebyshev_Expansion_425_550(long double x) {
380 0 : static long double c[] = {
381 : +1.058610209741581514157e-1L, -1.429297757627935191694e-2L,
382 : +9.911301703835545472874e-4L, -7.079903107876049846509e-5L,
383 : +5.229587914675267516134e-6L, -4.016071345964089296212e-7L,
384 : +3.231734714422926453741e-8L, -2.752870944370338482109e-9L,
385 : +2.503059741885009530630e-10L, -2.418699000594890423278e-11L,
386 : +2.410158905786160001792e-12L, -2.327254341132174000949e-13L,
387 : +1.958284411563056492727e-14L, -1.099893145048991004460e-15L,
388 : -2.959085292526991317697e-17L, +1.966366179276295203082e-17L,
389 : -3.314408783993662492621e-18L, +3.635520318133814622089e-19L,
390 : -2.550826919215104648800e-20L, +3.830090587178262542288e-22L,
391 : +1.836693763159216122739e-22L};
392 :
393 0 : static const int degree = sizeof(c) / sizeof(long double) - 1;
394 0 : static const long double midpoint = (5.50L + 4.25L) / 2.0L;
395 0 : static const long double half_length = (5.50L - 4.25L) / 2.0L;
396 0 : long double daw =
397 0 : xChebyshev_Tn_Series((fabsl(x) - midpoint) / half_length, c, degree);
398 0 : return (x > 0.0L) ? daw : -daw;
399 : }
400 :
401 : ////////////////////////////////////////////////////////////////////////////////
402 : // static long double Dawson_Chebyshev_Expansion_550_725( long double x ) //
403 : // //
404 : // Description: //
405 : // Evaluate Dawsons integral for 5.50 <= |x| <= 7.25 using the Chebyshev //
406 : // expansion of Daw(x) for x in the interval [5.50,7.25]. //
407 : // //
408 : // Arguments: //
409 : // long double x The argument of Dawson's integral Daw(), where //
410 : // where 5.50 <= |x| <= 7.25. //
411 : // //
412 : // Return Value: //
413 : // The value of Dawson's integral, Daw(), evaluated at x where //
414 : // 5.50 <= x <= 7.25. //
415 : // //
416 : // Example: //
417 : // long double y, x; //
418 : // //
419 : // ( code to initialize x ) //
420 : // //
421 : // y = Dawson_Chebyshev_Expansion_550_725(x); //
422 : ////////////////////////////////////////////////////////////////////////////////
423 :
424 0 : static long double Dawson_Chebyshev_Expansion_550_725(long double x) {
425 0 : static long double c[] = {+8.024637207807814739314e-2L,
426 : -1.136614891549306029413e-2L,
427 : +8.164249750628661856014e-4L,
428 : -5.951964778701328943018e-5L,
429 : +4.407349502747483429390e-6L,
430 : -3.317746826184531133862e-7L,
431 : +2.541483569880571680365e-8L,
432 : -1.983391157250772649001e-9L,
433 : +1.579050614491277335581e-10L,
434 : -1.284592098551537518322e-11L,
435 : +1.070070857004674207604e-12L,
436 : -9.151832297362522251950e-14L,
437 : +8.065447314948125338081e-15L,
438 : -7.360105847607056315915e-16L,
439 : +6.995966000187407197283e-17L,
440 : -6.964349343411584120055e-18L,
441 : +7.268789359189778223225e-19L,
442 : -7.885125241947769024019e-20L,
443 : +8.689022564130615225208e-21L,
444 : -9.353211304381231554634e-22L +
445 : 9.218280404899298404756e-23L};
446 :
447 0 : static const int degree = sizeof(c) / sizeof(long double) - 1;
448 : // static const long double lower_bound = 5.50L;
449 : // static const long double upper_bound = 7.25L;
450 0 : static const long double midpoint = (7.25L + 5.50L) / 2.0L;
451 0 : static const long double half_length = (7.25L - 5.50L) / 2.0L;
452 0 : long double daw =
453 0 : xChebyshev_Tn_Series((fabsl(x) - midpoint) / half_length, c, degree);
454 0 : return (x > 0.0L) ? daw : -daw;
455 : }
456 :
457 : ////////////////////////////////////////////////////////////////////////////////
458 : // static long double Dawson_Asymptotic_Expansion( long double x ) //
459 : // //
460 : // Description: //
461 : // For a large magnitude of the argument x, Dawson's integral can be //
462 : // expressed as the asymptotic series //
463 : // Daw(x) ~ (1/2x) [ 1 + 1 / (2x^2) + ... + (2j - 1)!! / (2x^2)^j + ... ] //
464 : // //
465 : // Arguments: //
466 : // long double x The argument of Dawson's integral Daw(), where //
467 : // |x| > 7. //
468 : // //
469 : // Return Value: //
470 : // The value of Dawson's integral, Daw(), evaluated at x where |x| > 7. //
471 : // //
472 : // Example: //
473 : // long double y, x; //
474 : // //
475 : // ( code to initialize x ) //
476 : // //
477 : // y = Dawson_Asymptotic_Expansion( x ); //
478 : ////////////////////////////////////////////////////////////////////////////////
479 0 : static long double Dawson_Asymptotic_Expansion(long double x) {
480 0 : long double term[Asymptotic_Expansion_Cutoff + 1];
481 0 : long double x2 = x * x;
482 0 : long double two_x = x + x;
483 0 : long double two_x2 = x2 + x2;
484 0 : long double xn = two_x2;
485 0 : long double Sn = 0.0L;
486 0 : long double factorial = 1.0L;
487 0 : int n;
488 :
489 0 : term[0] = 1.0L;
490 0 : term[1] = 1.0L / xn;
491 0 : for (n = 2; n <= Asymptotic_Expansion_Cutoff; n++) {
492 0 : xn *= two_x2;
493 0 : factorial *= (long double)(n + n - 1);
494 0 : term[n] = factorial / xn;
495 0 : if (term[n] < LDBL_EPSILON / 2.0L)
496 : break;
497 : }
498 :
499 0 : if (n > Asymptotic_Expansion_Cutoff)
500 : n = Asymptotic_Expansion_Cutoff;
501 0 : for (; n >= 0; n--)
502 0 : Sn += term[n];
503 0 : return Sn / two_x;
504 : }
505 :
506 : ////////////////////////////////////////////////////////////////////////////////
507 : // File: xchebyshev_Tn_series.c //
508 : // Routine(s): //
509 : // xChebyshev_Tn_Series //
510 : ////////////////////////////////////////////////////////////////////////////////
511 : ////////////////////////////////////////////////////////////////////////////////
512 : // long double xChebyshev_Tn_Series(long double x, long double a[],int degree)//
513 : // //
514 : // Description: //
515 : // This routine uses Clenshaw's recursion algorithm to evaluate a given //
516 : // polynomial p(x) expressed as a linear combination of Chebyshev //
517 : // polynomials of the first kind, Tn, at a point x, //
518 : // p(x) = a[0] + a[1]*T[1](x) + a[2]*T[2](x) + ... + a[deg]*T[deg](x). //
519 : // //
520 : // Clenshaw's recursion formula applied to Chebyshev polynomials of the //
521 : // first kind is: //
522 : // Set y[degree + 2] = 0, y[degree + 1] = 0, then for k = degree, ..., 1 //
523 : // set y[k] = 2 * x * y[k+1] - y[k+2] + a[k]. Finally //
524 : // set y[0] = x * y[1] - y[2] + a[0]. Then p(x) = y[0]. //
525 : // //
526 : // Arguments: //
527 : // long double x //
528 : // The point at which to evaluate the polynomial. //
529 : // long double a[] //
530 : // The coefficients of the expansion in terms of Chebyshev polynomials,//
531 : // i.e. a[k] is the coefficient of T[k](x). Note that in the calling //
532 : // routine a must be defined double a[N] where N >= degree + 1. //
533 : // int degree //
534 : // The degree of the polynomial p(x). //
535 : // //
536 : // Return Value: //
537 : // The value of the polynomial at x. //
538 : // If degree is negative, then 0.0 is returned. //
539 : // //
540 : // Example: //
541 : // long double x, a[N], p; //
542 : // int deg = N - 1; //
543 : // //
544 : // ( code to initialize x, and a[i] i = 0, ... , a[deg] ) //
545 : // //
546 : // p = xChebyshev_Tn_Series(x, a, deg); //
547 : ////////////////////////////////////////////////////////////////////////////////
548 :
549 0 : static long double xChebyshev_Tn_Series(long double x, long double a[],
550 : int degree) {
551 0 : long double yp2 = 0.0L;
552 0 : long double yp1 = 0.0L;
553 0 : long double y = 0.0L;
554 0 : long double two_x = x + x;
555 0 : int k;
556 :
557 : // Check that degree >= 0. If not, then return 0. //
558 :
559 0 : if (degree < 0)
560 : return 0.0L;
561 :
562 : // Apply Clenshaw's recursion save the last iteration. //
563 :
564 0 : for (k = degree; k >= 1; k--, yp2 = yp1, yp1 = y)
565 0 : y = two_x * yp1 - yp2 + a[k];
566 :
567 : // Now apply the last iteration and return the result. //
568 :
569 0 : return x * yp1 - yp2 + a[0];
570 : }
|
