1 module chunker.polynomials;
2 
3 /// Represents a polynomial from `F_2[X]`.
4 struct Pol
5 {
6 	alias Base = ulong;
7 	Base value;
8 
9 	/// Returns `this+y`.
10 	public Pol opBinary(string op)(Pol y) const
11 	if (op == "+")
12 	{
13 		return Pol(this.value ^ y.value);
14 	}
15 
16 	unittest
17 	{
18 		struct Test
19 		{
20 			Pol x, y;
21 			Pol sum;
22 		}
23 		Test[] tests =
24 		[
25 			{Pol(23), Pol(16), Pol(23 ^ 16)},
26 			{Pol(0x9a7e30d1e855e0a0), Pol(0x670102a1f4bcd414), Pol(0xfd7f32701ce934b4)},
27 			{Pol(0x9a7e30d1e855e0a0), Pol(0x9a7e30d1e855e0a0), Pol(0)},
28 		];
29 
30 		foreach (i, test; tests)
31 		{
32 			assert(test.sum == test.x + test.y,
33 				format!"test %d failed: sum != x+y"(i));
34 
35 			assert(test.sum == test.y + test.x,
36 				format!"test %d failed: sum != y+x"(i));
37 		}
38 	}
39 
40 
41 	/// Returns true if the multiplication would overflow `Pol.Base`.
42 	/// Code by Rob Pike, see
43 	/// https://groups.google.com/d/msg/golang-nuts/h5oSN5t3Au4/KaNQREhZh0QJ
44 	private static bool mulOverflows(Pol a, Pol b)
45 	{
46 		if (a.value <= 1 || b.value <= 1)
47 			return false;
48 		auto c = mulImpl(a, b);
49 		auto d = c / b;
50 		if (d != a)
51 			return true;
52 
53 		return false;
54 	}
55 
56 	private static Pol mulImpl(Pol x, Pol y)
57 	{
58 		if (x.value == 0 || y.value == 0)
59 			return Pol(0);
60 
61 		Pol res;
62 		foreach (i; 0 .. y.deg + 1)
63 			if ((y.value & (1L << uint(i))) > 0)
64 				res = res + Pol(x.value << uint(i));
65 
66 		return res;
67 	}
68 
69 	/// Returns `this*y`.
70 	/// When an overflow occurs, throws an exception.
71 	public Pol opBinary(string op)(Pol y) const
72 	if (op == "*")
73 	{
74 		if (mulOverflows(this, y))
75 			throw new Exception("multiplication would overflow ulong");
76 
77 		return mulImpl(this, y);
78 	}
79 
80 	version(unittest) static private Pol parseBin(string s)
81 	{
82 		import std.conv : to;
83 		return Pol(s.to!Base(2));
84 	}
85 
86 	unittest
87 	{
88 		struct Test
89 		{
90 			Pol x, y;
91 			Pol res;
92 		}
93 		Test[] tests =
94 		[
95 			{Pol(1), Pol(2), Pol(2)},
96 			{
97 				parseBin("1101"),
98 				parseBin("10"),
99 				parseBin("11010"),
100 			},
101 			{
102 				parseBin("1101"),
103 				parseBin("11"),
104 				parseBin("10111"),
105 			},
106 			{
107 				Pol(0x40000000),
108 				Pol(0x40000000),
109 				Pol(0x1000000000000000),
110 			},
111 			{
112 				parseBin("1010"),
113 				parseBin("100100"),
114 				parseBin("101101000"),
115 			},
116 			{
117 				parseBin("100"),
118 				parseBin("11"),
119 				parseBin("1100"),
120 			},
121 			{
122 				parseBin("11"),
123 				parseBin("110101"),
124 				parseBin("1011111"),
125 			},
126 			{
127 				parseBin("10011"),
128 				parseBin("110101"),
129 				parseBin("1100001111"),
130 			},
131 		];
132 
133 		foreach (i, test; tests)
134 		{
135 			auto m = test.x * test.y;
136 			assert(test.res == m,
137 				format!"TestPolMul failed for test %d: %s * %s: want %s, got %s"
138 				(i, test.x, test.y, test.res, m));
139 
140 			m = test.y * test.x;
141 			assert(test.res == m,
142 				format!"TestPolMul failed for %d: %s * %s: want %s, got %s"
143 				(i, test.x, test.y, test.res, m));
144 		}
145 	}
146 
147 	unittest
148 	{
149 		auto x = Pol(1L << 63);
150 		try
151 		{
152 			x * Pol(2);
153 			assert(false, "overflow test did not panic");
154 		}
155 		catch (Exception e)
156 		{
157 			// try to recover overflow error
158 			if (e.msg == "multiplication would overflow ulong")
159 				return;
160 
161 			stderr.writefln!"invalid error raised: %s"(e);
162 			// re-raise error if not overflow
163 			throw e;
164 		}
165 	}
166 
167 
168 	/// Returns the degree of the polynomial.
169 	/// If `this` is zero, -1 is returned.
170 	@property public int deg() const
171 	{
172 		Base x = this.value;
173 
174 		// the degree of 0 is -1
175 		if (x == 0)
176 			return -1;
177 
178 		// see https://graphics.stanford.edu/~seander/bithacks.html#IntegerLog
179 
180 		auto r = 0;
181 		if ((ulong(x) & 0xffffffff00000000) > 0)
182 		{
183 			x >>= 32;
184 			r |= 32;
185 		}
186 
187 		if ((ulong(x) & 0xffff0000) > 0)
188 		{
189 			x >>= 16;
190 			r |= 16;
191 		}
192 
193 		if ((ulong(x) & 0xff00) > 0)
194 		{
195 			x >>= 8;
196 			r |= 8;
197 		}
198 
199 		if ((ulong(x) & 0xf0) > 0)
200 		{
201 			x >>= 4;
202 			r |= 4;
203 		}
204 
205 		if ((ulong(x) & 0xc) > 0)
206 		{
207 			x >>= 2;
208 			r |= 2;
209 		}
210 
211 		if ((ulong(x) & 0x2) > 0)
212 		{
213 			r |= 1;
214 		}
215 
216 		return r;
217 	}
218 
219 	unittest
220 	{
221 		Pol x;
222 		assert(x.deg == -1,
223 			format!"deg(0) is not -1: %s"(x.deg));
224 
225 		x = Pol(1);
226 		assert(x.deg == 0,
227 			format!"deg(1) is not 0: %s"(x.deg));
228 
229 		foreach (i; 0 .. 64)
230 		{
231 			x = Pol(1L << uint(i));
232 			assert(x.deg == i,
233 				format!"deg(1<<%d) is not %d: %s"(i, i, x.deg));
234 		}
235 	}
236 
237 	version(benchmarkPolynomials) private static void benchmarkPolDeg()
238 	{
239 		auto f = Pol(0x3af4b284899);
240 		auto d = f.deg;
241 		assert(d == 41,
242 			format!"BenchmarkPolDeg: Wrong degree %d returned, expected %d"
243 			(d, 41));
244 
245 		Benchmark.benchmark({
246 			f.deg;
247 		});
248 	}
249 
250 
251 	/// Returns the coefficients in hex.
252 	public string toString() const
253 	{
254 		import std.conv : to;
255 		return "0x" ~ to!string(this.value, 16);
256 	}
257 
258 	/// Returns the string representation of the polynomial x.
259 	public string expand() const
260 	{
261 		import std.format : format;
262 
263 		if (value == 0)
264 			return "0";
265 
266 		auto s = "";
267 		for (auto i = deg; i > 1; i--)
268 			if ((value & (1L<<uint(i))) > 0)
269 				s ~= format!"+x^%d"(i);
270 
271 		if ((value & 2) > 0)
272 			s ~= "+x";
273 
274 		if ((value & 1) > 0)
275 			s ~= "+1";
276 
277 		return s[1 .. $];
278 	}
279 
280 	unittest
281 	{
282 		auto pol = Pol(0x3DA3358B4DC173);
283 		auto s = pol.expand();
284 		assert(s == "x^53+x^52+x^51+x^50+x^48+x^47+x^45+x^41+x^40+x^37+x^36+x^34+x^32+x^31+x^27+x^25+x^24+x^22+x^19+x^18+x^16+x^15+x^14+x^8+x^6+x^5+x^4+x+1",
285 			"wrong result");
286 	}
287 
288 
289 	/// Returns quotient and remainder from division `[x / d, x % d]`,
290 	/// see https://en.wikipedia.org/wiki/Division_algorithm
291 	public static Pol[2] divMod(Pol x, Pol d)
292 	{
293 		if (x.value == 0)
294 			return [Pol(0), Pol(0)];
295 
296 		if (d.value == 0)
297 			assert(false, "division by zero");
298 
299 		auto d_deg = d.deg;
300 		auto diff = x.deg - d_deg;
301 		if (diff < 0)
302 			return [Pol(0), x];
303 
304 		Pol q;
305 		while (diff >= 0)
306 		{
307 			auto m = d.value << uint(diff);
308 			q.value |= (1L << uint(diff));
309 			x = x + Pol(m);
310 
311 			diff = x.deg - d_deg;
312 		}
313 
314 		return [q, x];
315 	}
316 
317 	version(benchmarkPolynomials) private static void benchmarkPolDivMod()
318 	{
319 		auto f = Pol(0x2482734cacca49);
320 		auto g = Pol(0x3af4b284899);
321 
322 		Benchmark.benchmark({
323 			divMod(g, f);
324 		});
325 	}
326 
327 
328 	/// Returns the integer division result `this / d`.
329 	public Pol opBinary(string op)(Pol d) const
330 	if (op == "/")
331 	{
332 		return divMod(this, d)[0];
333 	}
334 
335 	unittest
336 	{
337 		struct Test
338 		{
339 			private Pol x, y;
340 			private Pol res;
341 		}
342 
343 		Test[] tests =
344 		[
345 			{Pol(10), Pol(50), Pol(0)},
346 			{Pol(0), Pol(1), Pol(0)},
347 			{
348 				parseBin("101101000"), // 0x168
349 				parseBin("1010"),      // 0xa
350 				parseBin("100100"),    // 0x24
351 			},
352 			{Pol(2), Pol(2), Pol(1)},
353 			{
354 				Pol(0x8000000000000000),
355 				Pol(0x8000000000000000),
356 				Pol(1),
357 			},
358 			{
359 				parseBin("1100"),
360 				parseBin("100"),
361 				parseBin("11"),
362 			},
363 			{
364 				parseBin("1100001111"),
365 				parseBin("10011"),
366 				parseBin("110101"),
367 			},
368 		];
369 
370 		foreach (i, test; tests)
371 		{
372 			auto m = test.x / test.y;
373 			assert(test.res == m,
374 				format!"TestPolDiv failed for test %d: %s * %s: want %s, got %s"
375 				(i, test.x, test.y, test.res, m));
376 		}
377 	}
378 
379 	version(benchmarkPolynomials) private static void benchmarkPolDiv()
380 	{
381 		auto f = Pol(0x2482734cacca49);
382 		auto g = Pol(0x3af4b284899);
383 
384 		Benchmark.benchmark({
385 			g / f;
386 		});
387 	}
388 
389 
390 	/// Returns the remainder of `this / d`.
391 	public Pol opBinary(string op)(Pol d) const
392 	if (op == "%")
393 	{
394 		return divMod(this, d)[1];
395 	}
396 
397 	unittest
398 	{
399 		struct Test
400 		{
401 			private Pol x, y;
402 			private Pol res;
403 		}
404 		Test[] tests =
405 		[
406 			{Pol(10), Pol(50), Pol(10)},
407 			{Pol(0), Pol(1), Pol(0)},
408 			{
409 				parseBin("101101001"),
410 				parseBin("1010"),
411 				parseBin("1"),
412 			},
413 			{Pol(2), Pol(2), Pol(0)},
414 			{
415 				Pol(0x8000000000000000),
416 				Pol(0x8000000000000000),
417 				Pol(0),
418 			},
419 			{
420 				parseBin("1100"),
421 				parseBin("100"),
422 				parseBin("0"),
423 			},
424 			{
425 				parseBin("1100001111"),
426 				parseBin("10011"),
427 				parseBin("0"),
428 			},
429 		];
430 
431 		foreach (i, test; tests)
432 		{
433 			auto res = test.x % test.y;
434 			assert(test.res == res,
435 				format!"test %d failed: want %s, got %s"(i, test.res, res));
436 		}
437 	}
438 
439 	version(benchmarkPolynomials) private static void benchmarkPolMod()
440 	{
441 		auto f = Pol(0x2482734cacca49);
442 		auto g = Pol(0x3af4b284899);
443 
444 		Benchmark.benchmark({
445 			g % f;
446 		});
447 	}
448 
449 
450 	/// I really dislike having a function that does not terminate, so specify a
451 	/// really large upper bound for finding a new irreducible polynomial, and
452 	/// return an error when no irreducible polynomial has been found within
453 	/// randPolMaxTries.
454 	private enum randPolMaxTries = 1e6;
455 
456 	/// Returns a new random irreducible polynomial of degree 53 using
457 	/// the default `rndGen` as source.  It is equivalent to calling
458 	/// `derive(generate!(() => uniform!ubyte))`.
459 	public static Pol getRandom()
460 	{
461 		import std.random : uniform;
462 		import std.range : generate;
463 		return derive(generate!(() => uniform!ubyte));
464 	}
465 
466 	unittest
467 	{
468 		getRandom();
469 	}
470 
471 	version(benchmarkPolynomials) private static void benchmarkPolGetRandom()
472 	{
473 		Benchmark.benchmark({
474 			getRandom();
475 		});
476 	}
477 
478 
479 	/// Returns an irreducible polynomial of degree 53
480 	/// (largest prime number below 64-8) by reading bytes from source.
481 	/// There are (2^53-2/53) irreducible polynomials of degree 53 in
482 	/// `F_2[X]`, c.f. Michael O. Rabin (1981): "Fingerprinting by Random
483 	/// Polynomials", page 4. If no polynomial could be found in one
484 	/// million tries, an error is returned.
485 	public static Pol derive(R)(R source)
486 	if (is (typeof(source.front) == ubyte))
487 	{
488 		foreach (i; 0 .. randPolMaxTries)
489 		{
490 			Pol f;
491 
492 			// choose polynomial at (pseudo)random
493 			f.value = 0;
494 			foreach (b; 0 .. Base.sizeof)
495 			{
496 				f.value = (f.value << 8) | source.front;
497 				source.popFront();
498 			}
499 
500 			// mask away bits above bit 53
501 			f.value &= Base((1L << 54) - 1);
502 
503 			// set highest and lowest bit so that the degree is 53 and the
504 			// polynomial is not trivially reducible
505 			f.value |= (1L << 53) | 1;
506 
507 			// test if f is irreducible
508 			if (f.irreducible)
509 				return f;
510 		}
511 
512 		// If this is reached, we haven't found an irreducible polynomial in
513 		// randPolMaxTries. This error is very unlikely to occur.
514 		throw new Exception("unable to find new random irreducible polynomial");
515 	}
516 
517 	/// Computes the Greatest Common Divisor of `x` and `f`.
518 	public static Pol gcd(Pol x, Pol f)
519 	{
520 		if (f.value == 0)
521 			return x;
522 
523 		if (x.value == 0)
524 			return f;
525 
526 		if (x.deg < f.deg)
527 		{
528 			import std.algorithm.mutation : swap;
529 			swap(x, f);
530 		}
531 
532 		return gcd(f, x % f);
533 	}
534 
535 	unittest
536 	{
537 		struct Test
538 		{
539 			private Pol f1;
540 			private Pol f2;
541 			private Pol gcd;
542 		}
543 		Test[] tests =
544 		[
545 			{Pol(10), Pol(50), Pol(2)},
546 			{Pol(0), Pol(1), Pol(1)},
547 			{
548 				parseBin("101101001"),
549 				parseBin("1010"),
550 				parseBin("1"),
551 			},
552 			{Pol(2), Pol(2), Pol(2)},
553 			{
554 				parseBin("1010"),
555 				parseBin("11"),
556 				parseBin("11"),
557 			},
558 			{
559 				Pol(0x8000000000000000),
560 				Pol(0x8000000000000000),
561 				Pol(0x8000000000000000),
562 			},
563 			{
564 				parseBin("1100"),
565 				parseBin("101"),
566 				parseBin("11"),
567 			},
568 			{
569 				parseBin("1100001111"),
570 				parseBin("10011"),
571 				parseBin("10011"),
572 			},
573 			{
574 				Pol(0x3DA3358B4DC173),
575 				Pol(0x3DA3358B4DC173),
576 				Pol(0x3DA3358B4DC173),
577 			},
578 			{
579 				Pol(0x3DA3358B4DC173),
580 				Pol(0x230d2259defd),
581 				Pol(1),
582 			},
583 			{
584 				Pol(0x230d2259defd),
585 				Pol(0x51b492b3eff2),
586 				parseBin("10011"),
587 			},
588 		];
589 
590 		foreach (i, test; tests)
591 		{
592 			auto r = gcd(test.f1, test.f2);
593 			assert(test.gcd == r,
594 				format!"GCD test %d (%+s) failed: got %s, wanted %s"
595 				(i, test, r, test.gcd));
596 
597 			r = gcd(test.f2, test.f1);
598 			assert(test.gcd == r,
599 				format!"GCD test %d (%+s) failed: got %s, wanted %s"
600 				(i, test, r, test.gcd));
601 		}
602 	}
603 
604 
605 	/// Returns `true` iff `x` is irreducible over `F_2`.
606 	/// This function uses Ben Or's reducibility test.
607 	//
608 	/// For details see "Tests and Constructions of Irreducible
609 	/// Polynomials over Finite Fields".
610 	@property public bool irreducible() const
611 	{
612 		foreach (i; 1 .. this.deg/2 + 1)
613 			if (gcd(this, qp(uint(i), this)).value != 1)
614 				return false;
615 
616 		return true;
617 	}
618 
619 	private
620 	{
621 		struct IrredTest
622 		{
623 			private Pol f;
624 			private bool irred;
625 		}
626 		static immutable IrredTest[] irredTests =
627 		[
628 			{Pol(0x38f1e565e288df), false},
629 			{Pol(0x3DA3358B4DC173), true},
630 			{Pol(0x30a8295b9d5c91), false},
631 			{Pol(0x255f4350b962cb), false},
632 			{Pol(0x267f776110a235), false},
633 			{Pol(0x2f4dae10d41227), false},
634 			{Pol(0x2482734cacca49), true},
635 			{Pol(0x312daf4b284899), false},
636 			{Pol(0x29dfb6553d01d1), false},
637 			{Pol(0x3548245eb26257), false},
638 			{Pol(0x3199e7ef4211b3), false},
639 			{Pol(0x362f39017dae8b), false},
640 			{Pol(0x200d57aa6fdacb), false},
641 			{Pol(0x35e0a4efa1d275), false},
642 			{Pol(0x2ced55b026577f), false},
643 			{Pol(0x260b012010893d), false},
644 			{Pol(0x2df29cbcd59e9d), false},
645 			{Pol(0x3f2ac7488bd429), false},
646 			{Pol(0x3e5cb1711669fb), false},
647 			{Pol(0x226d8de57a9959), false},
648 			{Pol(0x3c8de80aaf5835), false},
649 			{Pol(0x2026a59efb219b), false},
650 			{Pol(0x39dfa4d13fb231), false},
651 			{Pol(0x3143d0464b3299), false},
652 		];
653 	}
654 
655 	unittest
656 	{
657 		foreach (_, test; irredTests)
658 			assert(test.f.irreducible == test.irred,
659 				format!"Irreducibility test for Polynomial %s failed: got %s, wanted %s"
660 				(test.f, test.f.irreducible, test.irred));
661 	}
662 
663 	version(benchmarkPolynomials) private static void benchmarkPolIrreducible()
664 	{
665 		// find first irreducible polynomial
666 		Pol pol;
667 		foreach (_, test; irredTests)
668 			if (test.irred)
669 			{
670 				pol = test.f;
671 				break;
672 			}
673 
674 		Benchmark.benchmark({
675 			if (!pol.irreducible)
676 				assert(false, format!"Irreducibility test for Polynomial %s failed"(pol));
677 		});
678 	}
679 
680 
681 	/// Computes `this*f mod g`.
682 	public Pol mulMod(Pol f, Pol g)
683 	{
684 		if (value == 0 || f.value == 0)
685 			return Pol(0);
686 
687 		Pol res;
688 		foreach (i; 0 .. f.deg + 1)
689 			if ((f.value & (1L << uint(i))) > 0)
690 			{
691 				auto a = this;
692 				foreach (j; 0 .. i)
693 					a = (a * Pol(2)) % g;
694 				res = (res + a) % g;
695 			}
696 
697 		return res;
698 	}
699 
700 	unittest
701 	{
702 		struct Test
703 		{
704 			private Pol f1;
705 			private Pol f2;
706 			private Pol g;
707 			private Pol mod;
708 		}
709 		Test[] tests =
710 		[
711 			{
712 				Pol(0x1230),
713 				Pol(0x230),
714 				Pol(0x55),
715 				Pol(0x22),
716 			},
717 			{
718 				Pol(0x0eae8c07dbbb3026),
719 				Pol(0xd5d6db9de04771de),
720 				Pol(0xdd2bda3b77c9),
721 				Pol(0x425ae8595b7a),
722 			},
723 		];
724 
725 		foreach (i, test; tests)
726 		{
727 			auto mod = test.f1.mulMod(test.f2, test.g);
728 			assert(mod == test.mod,
729 				format!"mulMod test %d (%+s) failed: got %s, wanted %s"
730 				(i, test, mod, test.mod));
731 		}
732 	}
733 
734 
735 	/// Computes the polynomial `(x^(2^p)-x) mod g`.
736 	/// This is needed for the reducibility test.
737 	static private Pol qp(uint p, Pol g)
738 	{
739 		auto num = (1L << p);
740 		auto i = 1;
741 
742 		// start with x
743 		auto res = Pol(2);
744 
745 		while (i < num)
746 		{
747 			// repeatedly square res
748 			res = res.mulMod(res, g);
749 			i *= 2;
750 		}
751 
752 		// add x
753 		return (res + Pol(2)) % g;
754 	}
755 }
756 
757 version (benchmarkPolynomials)
758 {
759 	import chunker.internal.benchmark;
760 	mixin BenchmarkThisModule;
761 
762 	static foreach (name; __traits(allMembers, Pol))
763 		static if (name.length > 9 && name[0..9] == "benchmark")
764 			mixin(`private void ` ~ name ~ `() { Pol.` ~ name ~ `(); }`);
765 	version = test;
766 }
767 version(unittest) version = test;
768 version(test) import std.format : format;
769 version(test) import std.stdio : stderr;
770