1 module digitalnet.integration;
2 
3 import digitalnet.axiom;
4 import std.traits, std.math;
5 import std.range : ElementType;
6 import std.typecons : Flag;
7 alias Centering = Flag!"centering";
8 
9 private auto _integral(Centering centering=Centering.yes, S, Function)(S P, Function f)
10 	if (isPointSet!S && hasPrecision!S && hasDimensionR!S &&
11 		isArray!(ParameterTypeTuple!Function) && isFloatingPoint!(ElementType!(ParameterTypeTuple!Function)))
12 {
13 	pragma(msg, "_integral!(", typeid(S), ", ", typeid(Function), ")");
14 	alias F = ElementType!(ParameterTypeTuple!Function);
15 	F result = 0;
16 	foreach (x; P.toReals!(F, centering, S))
17 		result += f(x);
18 	return result * exp2(-cast(ptrdiff_t)P.dimensionF2);
19 }
20 
21 /// Calculate the integral of a function.
22 auto integral(Centering centering=Centering.yes, S, Function)(S P, Function f)
23 	if (isPointSet!S && hasPrecision!S && hasDimensionR!S
24 	 && isArray!(ParameterTypeTuple!Function) && isFloatingPoint!(ElementType!(ParameterTypeTuple!Function))
25 	 && !(isBisectable!S))
26 {
27 	pragma(msg, "integralNonbisectable!(", typeid(S), ", ", typeid(Function), ")");
28 	return _integral!(centering, S, Function)(P, f);
29 }
30 
31 /// ditto
32 auto integral(Centering centering=Centering.yes, S, Function)(S P, Function f)
33 	if (isPointSet!S && hasPrecision!S && hasDimensionR!S
34 	 && isArray!(ParameterTypeTuple!Function) && isFloatingPoint!(ElementType!(ParameterTypeTuple!Function))
35 	 && (isBisectable!S))
36 {
37 	pragma(msg, "integralBisectable!(", typeid(S), ", ", typeid(Function), ")");
38 	if (!P.bisectable)
39 		return _integral!(centering, S, Function)(P, f);
40 	auto Q = P.bisect;
41 	return (integral!(centering, S, Function)(Q[0], f) +
42 	        integral!(centering, S, Function)(Q[1], f)) / 2;
43 }
44 
45 /// Generate points in the hypercube from a digital net.
46 auto toReals(F, Centering centering=Centering.yes, S)(S P)
47 	if (isFloatingPoint!F && isPointSet!S && hasPrecision!S && hasDimensionR!S)
48 {
49 	struct R
50 	{
51 		static if (centering)
52 			enum F h = 0.5;
53 		else
54 			enum F h = 0;
55 		immutable F f;
56 		S P;
57 		alias P this;
58 		F[] front;
59 		void popFront()
60 		{
61 			P.popFront;
62 			if (empty)
63 				return;
64 			foreach (i, ref e; front)
65 				e = (P.front[i] + h) * f;
66 		}
67 		this (S P)
68 		{
69 			f = exp2(-cast(ptrdiff_t)P.precision);
70 			this.P = P;
71 			front.length = P.dimensionR;
72 			foreach (i, ref e; front)
73 				e = (P.front[i] + h) * f;
74 		}
75 	}
76 	return R(P);
77 }