1 // Written in the D programming language.
2 
3 /**
4 This module provides numerical integration by double-exponential (DE) formula.
5 */
6 module deint.deint;
7 
8 import std.algorithm;
9 import std.array;
10 import std.functional;
11 import std.math;
12 import std.typecons;
13 import std.traits;
14 
15 
16 /**
17 This template checks that a type `T` is a kind of NumInt.
18 */
19 enum bool isNumericalIntegrate(T)
20     = is(typeof((T integ){
21         foreach(xw; lockstep(integ.xs, integ.ws)) {}
22     }));
23 
24 
25 
26 /** This structure stores computing points and weights for numerical integration.
27  * In this type, a numerical integration is computed as
28  * int f(x) dx = sum_i f(xs[i]) ws[i]
29  */
30 struct NumInt(X, W)
31 {
32     /** Constructor
33      * Params:
34      *     xs = computing points
35      *     ws = weights for each computing points
36      */
37     this(immutable(X)[] xs, immutable(W)[] ws)
38     in{
39         assert(xs.length == ws.length);
40     }
41     do {
42         _xs = xs;
43         _ws = ws;
44     }
45 
46 
47   const @property
48   {
49     immutable(X)[] xs() { return _xs; }
50     immutable(W)[] ws() { return _ws; }
51   }
52 
53 
54     /** Execute integration
55     */
56     auto integrate(Fn)(scope Fn func) const
57     {
58         alias T = Unqual!(typeof(func(_xs[0]) * _ws[0]));
59 
60         T sum = 0;
61         foreach(i; 0 .. _xs.length)
62             sum += func(_xs[i]) * _ws[i];
63 
64         return sum;
65     }
66 
67 
68   private:
69     immutable(X)[] _xs;
70     immutable(W)[] _ws;
71 }
72 
73 
74 /** This function create a instance of NumInt!(F, F) for the DE formula.
75  * It is also known as "Tanh-sinh quadrature".
76  * In the DE formula, the integration (1) is converted to (2).
77  * (1) int_{xa}^{xb} f(x) dx
78  * (2) int_{ta}^{tb} f(g(t)) g'(t) dt
79  * 
80  * The type of DE formula is automatically decided from the given interval of the integration.
81  *
82  * Params:
83  *      xa = starting value of original integration.
84  *      xb = end value of original integration.
85  *      isExpDecay = if the integration is formed as int_a^b f(x) exp(-x) dx, this value is Yes. otherwise No.
86  *      trapN = division points of trapezoidal quadrature.
87  *      ta = starting value of integration transformed by DE-formula.
88  *      tb = starting value of integration transformed by DE-formula.
89  *
90  * Reference:
91  * $(HTTP en.wikipedia.org/wiki/Tanh-sinh_quadrature)
92 */
93 NumInt!(F, F) makeDEInt(F)(
94     F xa, F xb,
95     Flag!"isExpDecay" isExpDecay = No.isExpDecay,
96     size_t trapN = 100,
97     F ta = -5, F tb = 5)
98 {
99     NumInt!(F, F) _makeParamsImpl(scope F[2] delegate(F) fn)
100     {
101         immutable(F)[] xs;
102         immutable(F)[] ws;
103         if (!__ctfe)
104         {
105             xs.reserve(trapN);
106             ws.reserve(trapN);
107         }
108         immutable F h = (tb - ta) / (trapN-1);
109         foreach(i; 0 .. trapN) {
110             immutable xWt = fn(i * h + ta);
111             xs ~= xWt[0];
112             ws ~= xWt[1] * h;
113         }
114 
115         return NumInt!(F, F)(xs, ws);
116     }
117 
118 
119     if(xa > xb) {
120         auto invInt = .makeDEInt!F(xb, xa, isExpDecay, trapN, ta, tb);
121         auto newws = invInt.ws.map!"cast(immutable)(-a)".array;
122         return NumInt!(F, F)(invInt.xs, newws);
123     } else if(xa == -F.infinity && xb != F.infinity) {
124         auto tmpInt = .makeDEInt!F(-xb, F.infinity, isExpDecay, trapN, ta, tb);
125         auto newxs = tmpInt.xs.map!"cast(immutable)(-a)".array;
126         return NumInt!(F, F)(newxs, tmpInt.ws);
127     }
128 
129 
130     if(xa == -F.infinity && xb == F.infinity){
131         assert(!isExpDecay);
132 
133         return _makeParamsImpl(delegate(F t){
134             immutable F
135                 sinht = sinh(t),
136                 x = sinh(PI / 2 * sinht),
137                 dx = cosh(PI / 2 * sinht) * PI / 2 * cosh(t);
138 
139             return cast(F[2])[x, dx];
140         });
141     }else if(xb == F.infinity) {
142         return _makeParamsImpl(delegate(F t){
143             if(!isExpDecay){
144                 real x = exp(PI / 2 * sinh(t)),
145                         dx = x * PI / 2 * cosh(t);
146 
147                 return cast(F[2])[x + xa, dx];
148             }else{
149                 real expmt = exp(-t),
150                         x = exp(t - expmt),
151                         dx = (1 + expmt) * x;
152 
153 
154                 return cast(F[2])[x + xa, dx]; 
155             }
156         });
157     }else{
158         if(xa == 0 || xb == 0){
159             return _makeParamsImpl(delegate(F t){
160                 immutable F
161                     cosht = cosh(t),
162                     sinht = sinh(t),
163                     epsinht = exp(PI * sinht),
164                     x = (xa + xb*epsinht)/(1 + epsinht),
165                     cosh2 = cosh(PI / 2 * sinht)^^2,
166                     dx = PI / 2 * cosht / cosh2;
167 
168                 return cast(F[2])[x, dx * (xb - xa)/2];
169             });
170         }else{
171             immutable F diff2 = (xb - xa)/2;
172             immutable F avg2 = (xb + xa)/2;
173 
174             return _makeParamsImpl(delegate(F t){
175                 immutable F
176                     cosht = cosh(t),
177                     sinht = sinh(t),
178                     x = tanh(PI / 2 * sinht) * diff2 + avg2,
179                     cosh2 = cosh(PI / 2 * sinht)^^2,
180                     dx = PI / 2 * cosht / cosh2;
181 
182                 return cast(F[2])[x, dx * diff2];
183             });
184         }
185     }
186 }
187 
188 
189 ///
190 unittest
191 {
192     // integration on [0, 1]
193     auto int01 = makeDEInt!real(0, 1);
194 
195     // int_0^1 x dx = 0.5
196     assert(int01.integrate((real x) => x).approxEqual(0.5));
197 
198     // int_0^1 x^^2 dx = 1/3
199     assert(int01.integrate((real x) => x^^2).approxEqual(1/3.0));
200 
201     // int_0^1 log(x)/sqrt(x) dx = -4
202     assert(int01.integrate((real x) => log(x)/sqrt(x)).approxEqual(-4, 1E-12, 1E-12));
203 
204     // int_{-1}^0 log(-x)/sqrt(-x) dx = -4
205     assert(makeDEInt!real(-1, 0).integrate((real x) => log(-x)/sqrt(-x)).approxEqual(-4, 1E-12, 1E-12));
206 
207     // integration on [-1, 1]
208     auto int11 = makeDEInt!real(-1, 1);
209 
210     // int_{-1}^1 1/sqrt(1-x^2) dx = pi
211     assert(int11.integrate((real x) => 1/(1 + x^^2)).approxEqual(PI/2, 1E-12, 1E-12));
212 
213     // integration on [-inf, inf]
214     auto intII = makeDEInt!real(-real.infinity, real.infinity);
215 
216     // Gaussian integral
217     assert(intII.integrate((real x) => exp(-x^^2)).approxEqual(sqrt(PI)));
218 
219 
220     import std.mathspecial;
221     // integration on [1, inf] and integrand is expected to decay exponentially
222     auto intFI = makeDEInt!real(1, real.infinity, Yes.isExpDecay);
223 
224     // compute incomplete gamma function: int_1^inf x * exp(-x) dx = Gamma(2, 1)
225     assert(intFI.integrate((real x) => x * exp(-x)).approxEqual(gammaIncompleteCompl(2, 1) * gamma(2)));
226 }
227 
228 // Test of \int_a^b exp(x) dx
229 unittest 
230 {
231     // [a, b]
232     real[2][] dataset = [
233         [0.0, 1.0],
234         [-1.0, 1.0],
235         [0.0, -2.0],
236         [2.0, 0.0],
237         [-10.0, 0.0],
238         [0.0, 10.0]
239     ];
240 
241     foreach(data; dataset) {
242         real intDE = makeDEInt!real(data[0], data[1]).integrate((real x) => exp(x));
243         real truevalue = exp(data[1]) - exp(data[0]);
244 
245         assert(approxEqual(intDE, truevalue));
246     }
247 
248 
249     // reverse interval
250     foreach(data; dataset) {
251         //real intDE = intAToB_DE(data[0], data[1], 100, (real x) => exp(x));
252         real intDE = makeDEInt!real(data[1], data[0]).integrate((real x) => exp(x));
253         real truevalue = -(exp(data[1]) - exp(data[0]));
254 
255         assert(approxEqual(intDE, truevalue));
256     }
257 }
258 
259 // Complementary error function: erfc(a) = 2/sqrt(pi) * int_a^inf e^{-x^2} dx
260 unittest 
261 {
262     import std.mathspecial;
263 
264     real[] dataset = [
265         0.0,
266         0.01,
267         0.1,
268         1,
269         10,
270         100,
271     ];
272 
273     foreach(a; dataset) {
274         auto intDE = makeDEInt!real(a, real.infinity).integrate((real x) => exp(-x^^2));
275         intDE *= M_2_SQRTPI;
276 
277         auto truevalue = erfc(a);
278         assert(approxEqual(intDE, truevalue));
279     }
280 }
281 
282 // Complementary error function: erfc(a) = -2/sqrt(pi) * int_(-inf)^(-a) e^{-x^2} dx
283 unittest 
284 {
285     import std.mathspecial;
286 
287     real[] dataset = [
288         0.0,
289         0.01,
290         0.1,
291         1,
292         10,
293         100,
294     ];
295 
296     foreach(a; dataset) {
297         auto intDE = makeDEInt!real(-real.infinity, -a).integrate((real x) => exp(-x^^2));
298         intDE *= M_2_SQRTPI;
299 
300         auto truevalue = erfc(a);
301         assert(approxEqual(intDE, truevalue));
302     }
303 }
304 
305 // int_inf^inf 1/(1+x^2) dx = pi
306 unittest
307 {
308     auto intDE1 = makeDEInt!real(-real.infinity, real.infinity)
309                     .integrate((real x) => 1/(1 + x^^2));
310     assert(approxEqual(intDE1, PI));
311 }
312 
313 // int exp(-x^2) = sqrt(pi)
314 unittest
315 {
316     immutable val1 = makeDEInt!real(0, real.infinity, Yes.isExpDecay, 100)
317         .integrate((real x) => 1.0/(2*sqrt(x)) * exp(-x) );
318 
319     immutable val2 = makeDEInt!real(real.infinity, 0, Yes.isExpDecay, 100)
320         .integrate((real x) => -1.0/(2*sqrt(x)) * exp(-x) );
321 
322     immutable val3 = makeDEInt!real(-real.infinity, 0, Yes.isExpDecay, 100)
323         .integrate((real x) => 1.0/(2*sqrt(-x)) * exp(x) );
324 
325     immutable val4 = makeDEInt!real(0, -real.infinity, Yes.isExpDecay, 100)
326         .integrate((real x) => -1.0/(2*sqrt(-x)) * exp(x) );
327 
328     assert(approxEqual(val1, sqrt(PI)/2));
329     assert(approxEqual(val2, sqrt(PI)/2));
330     assert(approxEqual(val3, sqrt(PI)/2));
331     assert(approxEqual(val4, sqrt(PI)/2));
332 }
333 
334 // int exp(-x^2) * 1i = 1i * sqrt(pi)
335 unittest
336 {
337     import std.complex : Complex;
338 
339     immutable val1 = makeDEInt!real(0, real.infinity, Yes.isExpDecay)
340         .integrate((real x) => 1.0/(2*sqrt(x)) * exp(-x) * Complex!real(0, 1));
341 
342     assert(val1.re.approxEqual(0));
343     assert(val1.im.approxEqual(sqrt(PI)/2));
344 }
345 
346 // unittest for README.md
347 unittest
348 {
349     import std.mathspecial;
350     import deint.utils;
351 
352     // integrate f(x)exp(-x) on (1, inf)
353     // Now, we know that the integrand f(x)exp(-x) decay exponentially.
354     auto intFI = makeDEInt!real(1, real.infinity, Yes.isExpDecay);
355 
356     // incomplete gamma function
357     assert(intFI.integrate((real x) => x * exp(-x)).approxEqual(gammaIncompleteCompl(2, 1) * gamma(2)));
358 
359     // Also, you can use `withWeight` which pre-computes and stores weight exp(-x).
360     // The `withWeight` is useful when integrands have same weights.
361     auto intFIW = intFI.withWeight((real x) => exp(-x));
362 
363     // incomplete gamma function
364     assert(intFIW.integrate((real x) => x).approxEqual(gammaIncompleteCompl(2, 1) * gamma(2)));
365     assert(intFIW.integrate((real x) => x^^2).approxEqual(gammaIncompleteCompl(3, 1) * gamma(3)));
366     assert(intFIW.integrate((real x) => x^^3).approxEqual(gammaIncompleteCompl(4, 1) * gamma(4)));
367 }
368 
369 unittest
370 {
371     // integration on [-inf, inf]
372     auto intII = makeDEInt!real(-real.infinity, real.infinity);
373 
374     // Gaussian integral
375     assert(intII.integrate((real x) => exp(-x^^2)).approxEqual(sqrt(PI)));
376 }
377 
378 unittest
379 {
380     // DEInt!real is a struct which computes x_k and w_k in advance.
381     auto int01 = makeDEInt!real(0, 1);
382 
383     // When f(x) = x, int_0^1 x dx = 0.5
384     assert(int01.integrate((real x) => x).approxEqual(0.5));
385 
386     // `DEInt!real` is reusable.
387     // When f(x) = x^^2, int_0^1 x^^2 dx = 1/3
388     assert(int01.integrate((real x) => x^^2).approxEqual(1/3.0));
389 }
390 
391 
392 
393 /** This function create a instance of NumInt!(F, F) for computing fourtier-type integration by the DE formula.
394  * On the DE formula for the fourier-type integration, the integration (1) is converted to (2).
395  * (1) int_0^{inf} g(x) sin(omega * x) dx
396  * (2) sum_{k=-Nlow}^{Nhigh} g(M*phi(h*k)) sin(M*phi(h*k)/omega) M*phi'(h*k)/omega * h
397  * 
398  * Params:
399  *      isSine = `Yes` if the type of integration is "sine". `No` if the type is "cosine".
400  *      omega = angular frequency
401  *      stepH = step value of each interval
402  *      Nlow = number of negative-valued computing points
403  *      Nhigh = number of positive-valued computing points
404  *
405  * Reference:
406  * $(HTTP www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1040-20.pdf)
407 */
408 NumInt!(F, F) makeDEIntFourier(F)(
409     Flag!"isSine" isSine,
410     F omega,
411     F stepH = 0.1,
412     size_t Nlow = 49, size_t Nhigh = 50)
413 {
414     immutable F beta = 1.0L/4;
415     immutable F h = stepH;
416     immutable F M = PI / h;
417     immutable F alpha = beta / sqrt(1 + M*log(1 + M)/(4*PI));
418 
419     F[2] phi(F t) {
420         F num = t;
421         F den = 1 - exp(-2*t - alpha*(1 - exp(-t)) - beta*(exp(t) - 1));
422         F dden = (den - 1) * (-2 - alpha*exp(-t) - beta * exp(t));
423 
424         F num2 = den - num * dden;
425         F den2 = den^^2;
426 
427         if(abs(t) < 1E-3) {
428             F ddden = (den-1) * (-2 - alpha*exp(-t) - beta * exp(t))^^2
429                     + (den-1) * (   + alpha*exp(-t) - beta * exp(t));
430 
431             num = 1;
432             den = dden;
433 
434             num2 = dden - (dden + num * ddden);
435             den2 = 2*den*dden;
436         }
437 
438         return [num/den, num2/den2];
439     }
440 
441     immutable(F)[] xs;
442     immutable(F)[] ws;
443     if (!__ctfe)
444     {
445         xs.reserve(Nlow + Nhigh);
446         ws.reserve(Nlow + Nhigh);
447     }
448 
449     foreach(long i; -cast(long)Nlow .. cast(long)Nhigh) {
450         F t = h * i;
451         if(isSine) {
452             auto p = phi(t);
453             xs ~= M * p[0] / omega;
454             ws ~= M * p[1] / omega * h;
455         } else {
456             auto p = phi(t - PI/(2*M));
457             xs ~= M * p[0] / omega;
458             ws ~= M * p[1] / omega * h;
459         }
460     }
461 
462     return NumInt!(F, F)(xs, ws);
463 }
464 
465 unittest {
466     auto de = makeDEIntFourier!real(Yes.isSine, 1, 0.2);
467     auto dirichlet = de.integrate((real x) => sin(x)/x);
468 
469     assert(dirichlet.approxEqual(PI/2));
470 }