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 }