(git:ed6f26b)
Loading...
Searching...
No Matches
grpp_specfunc_dawson.c
Go to the documentation of this file.
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 //
46static long double xChebyshev_Tn_Series(long double x, long double a[],
47 int degree);
48
49// Internally Defined Routines //
50double libgrpp_Dawsons_Integral(double x);
51long double xDawsons_Integral(long double x);
52
53static long double Dawson_Power_Series(long double x);
54static long double Dawson_Chebyshev_Expansion_1_175(long double x);
55static long double Dawson_Chebyshev_Expansion_175_250(long double x);
56static long double Dawson_Chebyshev_Expansion_250_325(long double x);
57static long double Dawson_Chebyshev_Expansion_325_425(long double x);
58static long double Dawson_Chebyshev_Expansion_425_550(long double x);
59static long double Dawson_Chebyshev_Expansion_550_725(long double x);
60static 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////////////////////////////////////////////////////////////////////////////////
83double libgrpp_Dawsons_Integral(double x) {
84 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
109long double xDawsons_Integral(long double x) {
110 long double abs_x = fabsl(x);
111
112 if (abs_x <= 1.0L)
113 return Dawson_Power_Series(x);
114 if (abs_x <= 1.75L)
116 if (abs_x <= 2.50L)
118 if (abs_x <= 3.25L)
120 if (abs_x <= 4.25L)
122 if (abs_x <= 5.50L)
124 if (abs_x <= 7.25L)
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
153static long double Dawson_Power_Series(long double x) {
154 long double two_x2 = -2.0L * x * x;
155 long double sum = 0.0;
156 long double term = 1.0L;
157 long double factorial = 1.0L;
158 long double xn = 1.0L;
159 const long double epsilon = LDBL_EPSILON / 2.0L;
160 int y = 0;
161
162 if (x == 0.0L)
163 return 0.0L;
164 do {
165 sum += term;
166 y += 1;
167 factorial *= (long double)(y + y + 1);
168 xn *= two_x2;
169 term = xn / factorial;
170 } while (fabsl(term) > epsilon * fabsl(sum));
171 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
197static long double Dawson_Chebyshev_Expansion_1_175(long double x) {
198 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 static const int degree = sizeof(c) / sizeof(long double) - 1;
212 static const long double midpoint = (1.75L + 1.0L) / 2.0L;
213 static const long double half_length = (1.75L - 1.0L) / 2.0L;
214 long double daw =
215 xChebyshev_Tn_Series((fabsl(x) - midpoint) / half_length, c, degree);
216 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
242static long double Dawson_Chebyshev_Expansion_175_250(long double x) {
243 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 static const int degree = sizeof(c) / sizeof(long double) - 1;
257 static const long double midpoint = (2.50L + 1.75L) / 2.0L;
258 static const long double half_length = (2.50L - 1.75L) / 2.0;
259 long double daw =
260 xChebyshev_Tn_Series((fabsl(x) - midpoint) / half_length, c, degree);
261 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
287static long double Dawson_Chebyshev_Expansion_250_325(long double x) {
288 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 static const int degree = sizeof(c) / sizeof(long double) - 1;
302 static const long double midpoint = (3.25L + 2.50L) / 2.0L;
303 static const long double half_length = (3.25L - 2.50L) / 2.0L;
304 long double daw =
305 xChebyshev_Tn_Series((fabsl(x) - midpoint) / half_length, c, degree);
306 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
332static long double Dawson_Chebyshev_Expansion_325_425(long double x) {
333 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 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 static const long double midpoint = (4.25L + 3.25L) / 2.0L;
350 static const long double half_length = (4.25L - 3.25L) / 2.0L;
351 long double daw =
352 xChebyshev_Tn_Series((fabsl(x) - midpoint) / half_length, c, degree);
353 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
379static long double Dawson_Chebyshev_Expansion_425_550(long double x) {
380 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 static const int degree = sizeof(c) / sizeof(long double) - 1;
394 static const long double midpoint = (5.50L + 4.25L) / 2.0L;
395 static const long double half_length = (5.50L - 4.25L) / 2.0L;
396 long double daw =
397 xChebyshev_Tn_Series((fabsl(x) - midpoint) / half_length, c, degree);
398 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
424static long double Dawson_Chebyshev_Expansion_550_725(long double x) {
425 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 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 static const long double midpoint = (7.25L + 5.50L) / 2.0L;
451 static const long double half_length = (7.25L - 5.50L) / 2.0L;
452 long double daw =
453 xChebyshev_Tn_Series((fabsl(x) - midpoint) / half_length, c, degree);
454 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////////////////////////////////////////////////////////////////////////////////
479static long double Dawson_Asymptotic_Expansion(long double x) {
480 long double term[Asymptotic_Expansion_Cutoff + 1];
481 long double x2 = x * x;
482 long double two_x = x + x;
483 long double two_x2 = x2 + x2;
484 long double xn = two_x2;
485 long double Sn = 0.0L;
486 long double factorial = 1.0L;
487 int n;
488
489 term[0] = 1.0L;
490 term[1] = 1.0L / xn;
491 for (n = 2; n <= Asymptotic_Expansion_Cutoff; n++) {
492 xn *= two_x2;
493 factorial *= (long double)(n + n - 1);
494 term[n] = factorial / xn;
495 if (term[n] < LDBL_EPSILON / 2.0L)
496 break;
497 }
498
501 for (; n >= 0; n--)
502 Sn += term[n];
503 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
549static long double xChebyshev_Tn_Series(long double x, long double a[],
550 int degree) {
551 long double yp2 = 0.0L;
552 long double yp1 = 0.0L;
553 long double y = 0.0L;
554 long double two_x = x + x;
555 int k;
556
557 // Check that degree >= 0. If not, then return 0. //
558
559 if (degree < 0)
560 return 0.0L;
561
562 // Apply Clenshaw's recursion save the last iteration. //
563
564 for (k = degree; k >= 1; k--, yp2 = yp1, yp1 = y)
565 y = two_x * yp1 - yp2 + a[k];
566
567 // Now apply the last iteration and return the result. //
568
569 return x * yp1 - yp2 + a[0];
570}
static long double Dawson_Chebyshev_Expansion_175_250(long double x)
static long double Dawson_Chebyshev_Expansion_250_325(long double x)
static long double Dawson_Chebyshev_Expansion_550_725(long double x)
#define Asymptotic_Expansion_Cutoff
static long double xChebyshev_Tn_Series(long double x, long double a[], int degree)
static long double Dawson_Chebyshev_Expansion_325_425(long double x)
static long double Dawson_Asymptotic_Expansion(long double x)
long double xDawsons_Integral(long double x)
static long double Dawson_Chebyshev_Expansion_1_175(long double x)
double libgrpp_Dawsons_Integral(double x)
static long double Dawson_Power_Series(long double x)
static long double Dawson_Chebyshev_Expansion_425_550(long double x)