wafom.original source code

1 module wafom.original;
2 
3 import digitalnet.implementation;
4 import std.algorithm, std.math, std.functional;
5 public import std.math : lg = log2;
6 
7 real WAFOM(U, Size)(DigitalNet!(U, Size) P, real c = 1, real exponent = 1)
8 {
9 	return (P + new U[P.dimensionR]).WAFOMimpl(getWAFOMMultiplicand(P.precision, c, exponent));
10 }
11 
12 real RMSWAFOM(U, Size)(DigitalNet!(U, Size) P, real c = 1)
13 {
14 	return (P + new U[P.dimensionR]).WAFOMimpl(getWAFOMMultiplicand(P.precision, c, 2));
15 }
16 
17 private class WAFOMMultiplicand
18 {
19 	real[2][] t;
20 	this (size_t n, real c, real exponent)
21 	{
22 		t = new real[2][n];
23 		foreach (j, ref u; t)
24 		{
25 			immutable eps = exp2(exponent * (c - 1 - (n - j)));
26 			u[0] = 1 + eps;
27 			u[1] = 1 - eps;
28 		}
29 	}
30 }
31 
32 auto _getWAFOMMultiplicand(size_t n, real c, real exponent)
33 {
34 	return new WAFOMMultiplicand(n, c, exponent);
35 }
36 
37 alias getWAFOMMultiplicand = memoize!_getWAFOMMultiplicand;
38 
39 private real WAFOMimpl(U, Size)(ShiftedDigitalNet!(U, Size) P, WAFOMMultiplicand wafomMultiplicand)
40 {
41 	if (P.bisectable)
42 	{
43 		auto Q = P.bisect;
44 		return (Q[0].WAFOMimpl(wafomMultiplicand) + Q[1].WAFOMimpl(wafomMultiplicand)) / 2;
45 	}
46 	real ret = 0;
47 	foreach (X; P)
48 	{
49 		real current = 1;
50 		foreach (x; X)
51 			foreach (j, u; wafomMultiplicand.t)
52 				current *= u[x >> j & 1];
53 		ret += current - 1;
54 	}
55 	return ret * exp2(-cast(ptrdiff_t)P.dimensionF2);
56 }