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;