1 module digitalnet.implementation;
2 
3 import digitalnet.axiom;
4 import std.traits, std.algorithm, std.array, std.range, std.typecons, std.random, std.exception, std..string;
5 import std.conv : to;
6 debug import std.stdio;
7 
8 ///
9 alias LinearScramble = Tuple!(size_t[], size_t);
10 
11 static assert (isPointSet!(DigitalNet!uint));
12 static assert (hasPrecision!(DigitalNet!uint));
13 static assert (isPointSet!(ShiftedDigitalNet!uint));
14 static assert (hasPrecision!(ShiftedDigitalNet!uint));
15 
16 ///
17 mixin template DigitalNetFunctions(U, Size)
18 {
19 	const U[][] basis;
20 	immutable size_t precision, dimensionF2, dimensionR;
21 	immutable bool bisectable;
22 	this (in U[][] basis, Precision precision = Precision(U.sizeof << 3))
23 	{
24 		this.basis = basis;
25 		this.precision = precision.getPrecision;
26 		this.dimensionF2 = basis.length;
27 		this.dimensionF2.enforce("Zero digital net is not constructible");
28 		this.dimensionR = basis[0].length;
29 		foreach (b; basis)
30 			assert (b.length == dimensionR);
31 		this.front = new U[dimensionR];
32 	}
33 	bool empty;
34 	U[] front;
35 	void popFront()
36 	in
37 	{
38 		assert (!empty);
39 	}
40 	body
41 	{
42 		position += 1;
43 		if (position >> dimensionF2)
44 		{
45 			empty = true;
46 			return;
47 		}
48 		front[] ^= basis[position.bottom_zeros][];
49 	}
50 private:
51 	Size position;
52 	string _toString() const
53 	{
54 		return "%d %d %d %(%(%d %) %)".format(precision, dimensionF2, dimensionR, basis);
55 	}
56 }
57 
58 private U[] scramble(U)(in U[] vector, in size_t[] linearScramble, in size_t precision)
59 body
60 {
61 	static if (is (size_t == ulong))
62 		enum p = 6, q = 63;
63 	static if (is (size_t == uint))
64 		enum p = 5, q = 31;
65 	auto ret = vector.dup;
66 	size_t idx;
67 	foreach (i, ref x; ret)
68 		foreach (j; 1..precision)
69 			foreach (U k; 0..j)
70 			{
71 				x ^= ((linearScramble[idx >> p] >> (idx & q)) & (x >> j) & 1) << k;
72 				idx += 1;
73 			}
74 	return ret;
75 }
76 
77 /// A digital net.
78 struct DigitalNet(U = uint, Size = GreaterInteger!U)
79 	if (isUnsigned!U)
80 {
81 	mixin DigitalNetFunctions!(U, Size);
82 	/// Apply a digital shift.
83 	ShiftedDigitalNet!(U, Size) opBinary(string op)(in U[] shift) const
84 		if (op == "+")
85 	{
86 		return ShiftedDigitalNet!(U, Size)(basis, shift, Precision(precision));
87 	}
88 	DigitalNet!(U, Size) removeShift()
89 	{
90 		return this;
91 	}
92 	ShiftedDigitalNet!(U, Size)[2] bisect() const
93 	{
94 		auto former = ShiftedDigitalNet!(U, Size)(this.basis[1..$], new U[this.dimensionR], Precision(this.precision));
95 		return [former, former + this.basis[0]];
96 	}
97 	public string toString() const
98 	{
99 		return _toString();
100 	}
101 	/// Apply a linear scramble.
102 	DigitalNet!(U, Size) opBinary(string op)(in LinearScramble linearScramble) const
103 		if (op == "*")
104 	{
105 		return DigitalNet!(U, Size)(basis.map!(b => b.scramble(linearScramble[0], precision)).array, Precision(precision));
106 	}
107 }
108 
109 /// A digitally shifted digital net.
110 struct ShiftedDigitalNet(U = uint, Size = GreaterInteger!U)
111 	if (isUnsigned!U)
112 {
113 	mixin DigitalNetFunctions!(U, Size);
114 	const U[] shift;
115 	this (in U[][] basis, Precision precision = Precision(U.sizeof << 3))
116 	{
117 		this.basis = basis;
118 		this.precision = precision.getPrecision;
119 		this.dimensionF2 = basis.length;
120 		this.dimensionF2.enforce("Zero digital net is not constructible");
121 		this.dimensionR = basis[0].length;
122 		foreach (b; basis)
123 			assert (b.length == dimensionR);
124 		this.front = new U[dimensionR];
125 	}
126 	this (in U[][] basis, in U[] shift, Precision precision = Precision(U.sizeof << 3))
127 	{
128 		this (basis, precision);
129 		this.bisectable = bisectThreshold <= this.dimensionF2;
130 		this.shift = shift;
131 		this.front[] = shift[];
132 	}
133 	/// Apply a digital shift.
134 	ShiftedDigitalNet!(U, Size) opBinary(string op)(in U[] shift) const
135 		if (op == "+")
136 	{
137 		auto newShift = new U[dimensionR];
138 		foreach (i, ref s; newShift)
139 			s = this.shift[i] ^ shift[i];
140 		return ShiftedDigitalNet!(U, Size)(basis, newShift, Precision(precision));
141 	}
142 	/// Remove a digital shift.
143 	DigitalNet!(U, Size) removeShift() const
144 	{
145 		return DigitalNet!(U, Size)(basis, Precision(precision));
146 	}
147 	ShiftedDigitalNet!(U, Size)[2] bisect() const
148 	{
149 		auto former = ShiftedDigitalNet!(U, Size)(this.basis[1..$], this.shift, Precision(this.precision));
150 		return [former, former + this.basis[0]];
151 	}
152 	string toString() const
153 	{
154 		return _toString() ~ " %(%d %)".format(shift);
155 	}
156 	/// Apply a linear scramble.
157 	ShiftedDigitalNet!(U, Size) opBinary(string op)(in LinearScramble linearScramble) const
158 		if (op == "*")
159 	{
160 		return ShiftedDigitalNet!(U, Size)(basis.map!(b => b.scramble(linearScramble[0], precision)).array, shift.scramble(linearScramble[0], precision), Precision(precision));
161 	}
162 }
163 
164 private U[] getShift(U)(ref U[][] basis, in size_t dimB, in size_t dimR)
165 {
166 	if (basis.length == dimB)
167 	{
168 		stderr.writeln("run-time warning: add a zero shift");
169 		return new U[dimR];
170 	}
171 	if (basis.length == dimB + 1)
172 	{
173 		auto s = basis[$-1];
174 		basis.length -= 1;
175 		return s;
176 	}
177 	assert (false);
178 }
179 
180 private U[][] noShift(U)(U[][] basis, in size_t dimB)
181 {
182 	if (basis.length == dimB + 1)
183 	{
184 		import std.stdio : stderr;
185 		stderr.writefln("run-time warning: ignoring shift, [%(%d %)]", basis[dimB..$][0]);
186 		return basis[0..$-1];
187 	}
188 	if (basis.length == dimB)
189 		return basis;
190 	assert (false);
191 }
192 
193 private U[][] readDigitalNetParams(U)(const(char)[][] buf, out size_t prec, out size_t dimB, out size_t dimR)
194 {
195 	prec = buf[0].to!size_t;
196 	dimB = buf[1].to!size_t;
197 	dimR = buf[2].to!size_t;
198 	return buf[3..$].map!(to!U).array.chunks(dimR).array;
199 }
200 
201 /// Read a digital net from a string.
202 DigitalNet!U toDigitalNet(U = uint)(const(char)[] x)
203 {
204 	size_t prec, dimB, dimR;
205 	auto basis = x.strip.split(',')[0].split.readDigitalNetParams!U(prec, dimB, dimR);
206 	return DigitalNet!U(noShift(basis, dimB), Precision(prec));
207 }
208 
209 /// Read a digitally shifted digital net from a string.
210 ShiftedDigitalNet!U toShiftedDigitalNet(U = uint)(const(char)[] x)
211 {
212 	size_t prec, dimB, dimR;
213 	auto basis = x.strip.split(',')[0].split.readDigitalNetParams!U(prec, dimB, dimR);
214 	auto shift = getShift(basis, dimB, dimR);
215 	return ShiftedDigitalNet!U(basis, shift, Precision(prec));
216 }
217 
218 ///
219 template GreaterInteger(U)
220 	if (isUnsigned!U)
221 {
222 	import std.bigint;
223 	static if (U.sizeof < ulong.sizeof)
224 		alias GreaterInteger = ulong;
225 	else
226 		alias GreaterInteger = BigInt;
227 }
228 
229 immutable(size_t) getPrecision(Precision n)
230 {
231 	return cast(size_t)n;
232 }
233 immutable(size_t) getDimensionR(DimensionR s)
234 {
235 	return cast(size_t)s;
236 }
237 immutable(size_t) getDimensionF2(DimensionF2 m)
238 {
239 	return cast(size_t)m;
240 }
241 
242 ///
243 alias Precision = Typedef!(size_t, 0, "prec");
244 
245 ///
246 alias DimensionR = Typedef!(size_t, 0, "dimR");
247 
248 ///
249 alias DimensionF2 = Typedef!(size_t, 0, "dimB");
250 
251 /// Construct a digital net by uniform random choice of its basis.
252 DigitalNet!U randomDigitalNet(U = uint)(Precision precision, DimensionR dimensionR, DimensionF2 dimensionF2)
253 	if (isUnsigned!U)
254 {
255 	immutable s = dimensionR.getDimensionR;
256 	auto basis = randomVector!U(precision.getPrecision, s * dimensionF2.getDimensionF2).chunks(s).array;
257 	return DigitalNet!U(basis, precision);
258 }
259 
260 /// Uniformly and randomly pick a digital shift.
261 auto randomDigitalShift(S)(S P)
262 	if (is (S == DigitalNet!(U, Size), U, Size) || is (S == ShiftedDigitalNet!(U, Size), U, Size))
263 {
264 	alias U = ElementType!(typeof (P.front));
265 	return randomVector!U(P.precision, P.dimensionR);
266 }
267 
268 /// ditto
269 U[] randomDigitalShift(U)(Precision precision, DimensionR dimensionR)
270 	if (isUnsigned!U)
271 {
272 	return randomVector!U(precision.getPrecision, dimensionR.getDimensionR);
273 }
274 
275 /// Uniformly and randomly pick a linear scramble.
276 auto randomLinearScramble(S)(S P)
277 	if (is (S == DigitalNet!(U, Size), U, Size) ||
278 		is (S == ShiftedDigitalNet!(U, Size), U, Size))
279 {
280 	immutable size_t
281 		wordSize = size_t.sizeof << 3,
282 		numBits = P.dimensionR * (P.precision * (P.precision - 1) / 2),
283 		backLength = numBits / wordSize + (numBits % wordSize ? 1 : 0);
284 	import std.bitmanip : BitArray;
285 	return LinearScramble(randomVector!size_t(wordSize, backLength), numBits);
286 }
287 
288 /// Shuffle the basis of a digital net.
289 DigitalNet!(U, Size) shuffle(U, Size)(DigitalNet!(U, Size) P)
290 {
291 	auto b = P.basis.map!(a => a.to!(U[])).array;
292 	b.shuffleBasis;
293 	return DigitalNet!(U, Size)(b, Precision(P.precision));
294 }
295 
296 /// Partially shuffle a sequence of vectors.
297 void partialShuffleBasis(U)(U[][] basis)
298 {
299 	import std.random;
300 	auto i = 0.uniform(basis.length);
301 	swap(basis[0], basis[i]);
302 	auto j = 0.uniform(1 << (basis.length - 1));
303 	foreach (k, b; basis[1..$])
304 		if (j >> k & 1)
305 			basis[0][] ^= b[];
306 }
307 
308 /// Shuffle a sequence of vectors.
309 void shuffleBasis(U)(U[][] basis)
310 {
311 	if (basis.length == 1)
312 		return;
313 	basis.partialShuffleBasis;
314 	basis[1..$].shuffleBasis;
315 }
316 
317 /** Obtain another digital net Q from a digital net P where dim P = 1 + dim(P cap Q) = dim Q.
318 */
319 DigitalNet!(U, Size) changeVectorOfBasis(U, Size)(DigitalNet!(U, Size) P)
320 {
321 	auto b = P.basis.map!(a => a.to!(U[])).array;
322 	b.partialShuffleBasis;
323 	b[0] = randomDigitalShift(P);
324 	return DigitalNet!(U, Size)(b, Precision(P.precision));
325 }
326 
327 private:
328 
329 U[] randomVector(U)(size_t precision, size_t length)
330 	if (isUnsigned!U)
331 {
332 	auto ret = new U[length];
333 	randomVector(precision, ret);
334 	return ret;
335 }
336 
337 void randomVector(U)(size_t precision, U[] buffer)
338 {
339 	randomVector(buffer);
340 	if (auto d = (U.sizeof << 3) - precision)
341 		foreach (ref e; buffer)
342 			e >>= d;
343 }
344 
345 void randomVector(U)(U[] buffer)
346 {
347 	foreach (ref e; buffer)
348 		e = uniform!U;
349 }
350 
351 auto bottom_zeros(Size)(Size x)
352 {
353 	assert (x);
354 	size_t ret;
355 	while ((x & 1) == 0)
356 	{
357 		x >>= 1;
358 		ret += 1;
359 	}
360 	return ret;
361 }
362 
363 enum size_t bisectThreshold = 10;