1 /++
2 Note:
3 The module doesn't provide full arithmetic API for now.
4 +/
5 module mir.bignum.integer;
6
7 import mir.bitop;
8 import mir.serde: serdeProxy, serdeScoped;
9 import mir.utility;
10 import std.traits;
11
12 /++
13 Stack-allocated big signed integer.
14 Params:
15 maxSize64 = count of 64bit words in coefficient
16 +/
17 @serdeScoped @serdeProxy!(const(char)[])
18 struct BigInt(size_t maxSize64)
19 if (maxSize64 && maxSize64 <= ushort.max)
20 {
21 import mir.bignum.low_level_view;
22 import mir.bignum.fixed;
23
24 ///
25 bool sign;
26 ///
27 uint length;
28 ///
29 size_t[ulong.sizeof / size_t.sizeof * maxSize64] data = void;
30
31 ///
32 this(size_t size)(UInt!size fixedInt)
33 {
34 this(fixedInt.data);
35 }
36
37 ///
38 this(size_t N)(size_t[N] data)
39 if (N <= this.data.length)
40 {
41 sign = false;
42 version(LittleEndian)
43 this.data[0 .. N] = data;
44 else
45 this.data[$ - N .. $] = data;
46 length = data.length;
47 normalize;
48 }
49
50 ///
51 this(ulong data)
52 {
53 sign = false;
54 static if (size_t.sizeof == ulong.sizeof)
55 {
56 length = 1;
57 view.leastSignificantFirst[0] = data;
58 }
59 else
60 {
61 length = 2;
62 auto d = view.leastSignificantFirst;
63 d[0] = cast(uint) data;
64 d[1] = cast(uint) (data >> 32);
65 }
66 normalize;
67 }
68
69 ///
70 this()(scope const(char)[] str) @safe pure @nogc
71 // if (isSomeChar!C)
72 {
73 if (fromStringImpl(str))
74 return;
75 static if (__traits(compiles, () @nogc { throw new Exception("Can't parse BigInt."); }))
76 {
77 import mir.exception: MirException;
78 throw new MirException("Can't parse BigInt!" ~ maxSize64.stringof ~ " from string `", str , "`.");
79 }
80 else
81 {
82 static immutable exception = new Exception("Can't parse BigInt!" ~ maxSize64.stringof ~ ".");
83 throw exception;
84 }
85 }
86
87 ///
88 ref opAssign(ulong data) return
89 {
90 static if (size_t.sizeof == ulong.sizeof)
91 {
92 length = 1;
93 view.leastSignificantFirst[0] = data;
94 }
95 else
96 {
97 length = 2;
98 auto d = view.leastSignificantFirst;
99 d[0] = cast(uint) data;
100 d[1] = cast(uint) (data >> 32);
101 }
102 normalize;
103 return this;
104 }
105
106 static if (maxSize64 == 3)
107 ///
108 version(mir_bignum_test) @safe pure @nogc unittest
109 {
110 import mir.math.constant: PI;
111 BigInt!4 integer = "-34010447314490204552169750449563978034784726557588085989975288830070948234680"; // constructor
112 assert(integer.sign);
113 integer.sign = false;
114 assert(integer == BigInt!4.fromHexString("4b313b23aa560e1b0985f89cbe6df5460860e39a64ba92b4abdd3ee77e4e05b8"));
115 }
116
117 ///
118 ref opAssign(size_t rhsMaxSize64)(auto ref scope const BigInt!rhsMaxSize64 rhs) return
119 if (rhsMaxSize64 < maxSize64)
120 {
121 this.sign = rhs.sign;
122 this.length = rhs.length;
123 version(LittleEndian)
124 {
125 data[0 .. length] = rhs.data[0 .. length];
126 }
127 else
128 {
129 data[$ - length .. $] = rhs.data[$ - length .. $];
130 }
131 return this;
132 }
133
134 /++
135 Returns: false in case of overflow or incorrect string.
136 Precondition: non-empty coefficients.
137 +/
138 bool fromStringImpl(C)(scope const(C)[] str)
139 scope @trusted pure @nogc nothrow
140 if (isSomeChar!C)
141 {
142 auto work = BigIntView!size_t(data[]);
143 if (work.fromStringImpl(str))
144 {
145 length = cast(uint) work.coefficients.length;
146 sign = work.sign;
147 return true;
148 }
149 return false;
150 }
151
152 ///
153 BigInt copy() @property
154 {
155 BigInt ret;
156 ret.sign = sign;
157 ret.length = length;
158 ret.data = data;
159 return ret;
160 }
161
162 ///
163 bool opEquals()(auto ref const BigInt rhs)
164 const @safe pure nothrow @nogc
165 {
166 return view == rhs.view;
167 }
168
169 ///
170 bool opEquals()(size_t rhs, bool rhsSign = false)
171 const @safe pure nothrow @nogc
172 {
173 return rhs == 0 && length == 0 || length == 1 && sign == rhsSign && view.unsigned.leastSignificant == rhs;
174 }
175
176 ///
177 bool opEquals()(sizediff_t rhs)
178 const @safe pure nothrow @nogc
179 {
180 auto sign = rhs < 0;
181 return opEquals(sign ? ulong(-rhs) : ulong(rhs), sign);
182 }
183
184 /++
185 +/
186 auto opCmp()(auto ref const BigInt rhs)
187 const @safe pure nothrow @nogc
188 {
189 return view.opCmp(rhs.view);
190 }
191
192 ///
193 BigIntView!size_t view()() @property
194 {
195 version (LittleEndian)
196 return typeof(return)(data[0 .. length], sign);
197 else
198 return typeof(return)(data[$ - length .. $], sign);
199 }
200
201 ///
202 BigIntView!(const size_t) view()() const @property
203 {
204 version (LittleEndian)
205 return typeof(return)(data[0 .. length], sign);
206 else
207 return typeof(return)(data[$ - length .. $], sign);
208 }
209
210 ///
211 void normalize()()
212 {
213 auto norm = view.normalized;
214 this.length = cast(uint) norm.unsigned.coefficients.length;
215 this.sign = norm.sign;
216 }
217
218 /++
219 +/
220 void putCoefficient(size_t value)
221 {
222 assert(length < data.length);
223 version (LittleEndian)
224 data[length++] = value;
225 else
226 data[$ - ++length] = value;
227 }
228
229 /++
230 Performs `size_t overflow = (big += overflow) *= scalar` operatrion.
231 Params:
232 rhs = unsigned value to multiply by
233 overflow = initial overflow value
234 Returns:
235 unsigned overflow value
236 +/
237 size_t opOpAssign(string op : "*")(size_t rhs, size_t overflow = 0u)
238 @safe pure nothrow @nogc
239 {
240 if (length == 0)
241 goto L;
242 overflow = view.unsigned.opOpAssign!op(rhs, overflow);
243 if (overflow && length < data.length)
244 {
245 L:
246 putCoefficient(overflow);
247 overflow = 0;
248 }
249 return overflow;
250 }
251
252 /++
253 Performs `uint remainder = (overflow$big) /= scalar` operatrion, where `$` denotes big-endian concatenation.
254 Precondition: `overflow < rhs`
255 Params:
256 rhs = unsigned value to devide by
257 overflow = initial unsigned overflow
258 Returns:
259 unsigned remainder value (evaluated overflow)
260 +/
261 uint opOpAssign(string op : "/")(uint rhs, uint overflow = 0)
262 @safe pure nothrow @nogc
263 {
264 assert(overflow < rhs);
265 if (length)
266 return view.unsigned.opOpAssign!op(rhs, overflow);
267 return overflow;
268 }
269
270 /++
271 Performs `size_t overflow = (big += overflow) *= fixed` operatrion.
272 Params:
273 rhs = unsigned value to multiply by
274 overflow = initial overflow value
275 Returns:
276 unsigned overflow value
277 +/
278 UInt!size opOpAssign(string op : "*", size_t size)(UInt!size rhs, UInt!size overflow = 0)
279 @safe pure nothrow @nogc
280 {
281 if (length == 0)
282 goto L;
283 overflow = view.unsigned.opOpAssign!op(rhs, overflow);
284 if (overflow && length < data.length)
285 {
286 L:
287 static if (size <= 64)
288 {
289 auto o = cast(ulong)overflow;
290 static if (size_t.sizeof == ulong.sizeof)
291 {
292 putCoefficient(o);
293 overflow = UInt!size.init;
294 }
295 else
296 {
297 putCoefficient(cast(uint)o);
298 o >>= 32;
299 if (length < data.length)
300 {
301 putCoefficient(cast(uint)o);
302 o = 0;
303 }
304 overflow = UInt!size(o);
305 }
306 }
307 else
308 {
309 do
310 {
311 putCoefficient(cast(size_t)overflow);
312 overflow >>= size_t.sizeof * 8;
313 }
314 while(overflow && length < data.length);
315 }
316 }
317 return overflow;
318 }
319
320 /++
321 Performs `size_t overflow = big *= fixed` operatrion.
322 Params:
323 rhs = unsigned value to multiply by
324 Returns:
325 overflow
326 +/
327 bool opOpAssign(string op, size_t rhsMaxSize64)(ref const BigInt!rhsMaxSize64 rhs)
328 @safe pure nothrow @nogc
329 if (op == "+" || op == "-")
330 {
331 return opOpAssign!op(rhs.view);
332 }
333
334 /// ditto
335 bool opOpAssign(string op)(BigIntView!(const size_t) rhs)
336 @safe pure nothrow @nogc
337 if (op == "+" || op == "-")
338 {
339 sizediff_t diff = length - rhs.coefficients.length;
340 if (diff < 0)
341 {
342 auto oldLength = length;
343 length = cast(int)rhs.coefficients.length;
344 view.unsigned.leastSignificantFirst[oldLength .. $] = 0;
345 }
346 else
347 if (rhs.coefficients.length == 0)
348 return false;
349 auto thisView = view;
350 auto overflow = thisView.opOpAssign!op(rhs);
351 this.sign = thisView.sign;
352 if (overflow)
353 {
354 if (length < data.length)
355 {
356 putCoefficient(overflow);
357 overflow = false;
358 }
359 }
360 else
361 {
362 normalize;
363 }
364 return overflow;
365 }
366
367 /++
368 +/
369 static BigInt fromHexString(bool allowUnderscores = false)(scope const(char)[] str)
370 @trusted pure
371 {
372 BigInt ret;
373 if (ret.fromHexStringImpl!(char, allowUnderscores)(str))
374 return ret;
375 version(D_Exceptions)
376 throw hexStringException;
377 else
378 assert(0, hexStringErrorMsg);
379 }
380
381 /++
382 +/
383 bool fromHexStringImpl(C, bool allowUnderscores = false)(scope const(C)[] str)
384 @safe pure @nogc nothrow
385 if (isSomeChar!C)
386 {
387 auto work = BigIntView!size_t(data);
388 auto ret = work.fromHexStringImpl!(C, allowUnderscores)(str);
389 if (ret)
390 {
391 length = cast(uint)work.unsigned.coefficients.length;
392 sign = work.sign;
393 }
394 return ret;
395 }
396
397 ///
398 bool mulPow5(size_t degree)
399 {
400 // assert(approxCanMulPow5(degree));
401 if (length == 0)
402 return false;
403 enum n = MaxWordPow5!size_t;
404 enum wordInit = size_t(5) ^^ n;
405 size_t word = wordInit;
406 bool of;
407 while(degree)
408 {
409 if (degree >= n)
410 {
411 degree -= n;
412 }
413 else
414 {
415 word = 1;
416 do word *= 5;
417 while(--degree);
418 }
419 if (auto overflow = view *= word)
420 {
421 of = length >= data.length;
422 if (!of)
423 putCoefficient(overflow);
424 }
425 }
426 return of;
427 }
428
429 ///
430 ref BigInt opOpAssign(string op)(size_t shift)
431 @safe pure nothrow @nogc return
432 if (op == "<<" || op == ">>")
433 {
434 auto index = shift / (size_t.sizeof * 8);
435 auto bs = shift % (size_t.sizeof * 8);
436 auto ss = size_t.sizeof * 8 - bs;
437 static if (op == ">>")
438 {
439 if (index >= length)
440 {
441 length = 0;
442 return this;
443 }
444 auto d = view.leastSignificantFirst;
445 if (bs)
446 {
447 foreach (j; 0 .. d.length - (index + 1))
448 {
449 d[j] = (d[j + index] >>> bs) | (d[j + (index + 1)] << ss);
450 }
451 }
452 else
453 {
454 foreach (j; 0 .. d.length - (index + 1))
455 {
456 d[j] = d[j + index];
457 }
458 }
459 auto most = d[$ - (index + 1)] = d.back >>> bs;
460 length -= index + (most == 0);
461 }
462 else
463 {
464 if (index >= data.length || length == 0)
465 {
466 length = 0;
467 return this;
468 }
469
470 if (bs)
471 {
472 auto most = view.unsigned.mostSignificant >> ss;
473 length += index;
474 if (length < data.length)
475 {
476 if (most)
477 {
478 length++;
479 view.unsigned.mostSignificant = most;
480 length--;
481 }
482 }
483 else
484 {
485 length = data.length;
486 }
487
488 auto d = view.leastSignificantFirst;
489 foreach_reverse (j; index + 1 .. length)
490 {
491 d[j] = (d[j - index] << bs) | (d[j - (index + 1)] >> ss);
492 }
493 d[index] = d.front << bs;
494 if (length < data.length)
495 length += most != 0;
496 }
497 else
498 {
499 length = cast(uint) min(length + index, cast(uint)data.length);
500 auto d = view.leastSignificantFirst;
501 foreach_reverse (j; index .. length)
502 {
503 d[j] = d[j - index];
504 }
505 }
506 view.leastSignificantFirst[0 .. index] = 0;
507 }
508 return this;
509 }
510
511 ///
512 T opCast(T, bool wordNormalized = false, bool nonZero = false)() const
513 if (isFloatingPoint!T && isMutable!T)
514 {
515 return view.opCast!(T, wordNormalized, nonZero);
516 }
517
518 ///
519 T opCast(T, bool nonZero = false)() const
520 if (is(T == long) || is(T == int))
521 {
522 return this.view.opCast!(T, nonZero);
523 }
524
525 /++
526 Returns: overflow flag
527 +/
528 bool copyFrom(W, WordEndian endian)(BigIntView!(const W, endian) view)
529 {
530 static if (W.sizeof > size_t.sizeof && endian == TargetEndian)
531 {
532 return this.copyFrom(cast(BigIntView!(const size_t))view);
533 }
534 else
535 {
536 this.sign = view.sign;
537 auto lhs = BigUIntView!W(cast(W[])data);
538 auto rhs = view;
539 auto overflow = lhs.coefficients.length < rhs.coefficients.length;
540 auto n = overflow ? lhs.coefficients.length : rhs.coefficients.length;
541 lhs.leastSignificantFirst[0 .. n] = rhs.leastSignificantFirst[0 .. n];
542 this.length = cast(uint)(n / (size_t.sizeof / W.sizeof));
543 if (auto tail = n % (size_t.sizeof / W.sizeof))
544 {
545 this.length++;
546 auto shift = ((size_t.sizeof / W.sizeof) - tail) * (W.sizeof * 8);
547 auto value = this.view.unsigned.mostSignificant;
548 value <<= shift;
549 value >>= shift;
550 this.view.unsigned.mostSignificant = value;
551 }
552 return overflow;
553 }
554 }
555
556 /// ditto
557 bool copyFrom(W, WordEndian endian)(BigUIntView!(const W, endian) view)
558 {
559 return this.copyFrom(BigIntView!(const W, endian)(view));
560 }
561
562 ///
563 immutable(C)[] toString(C = char)() const @safe pure nothrow
564 if(isSomeChar!C && isMutable!C)
565 {
566 C[ceilLog10Exp2(data.length * (size_t.sizeof * 8)) + 1] buffer = void;
567 BigInt copy = this;
568 auto len = copy.view.unsigned.toStringImpl(buffer);
569 if (sign)
570 buffer[$ - ++len] = '-';
571 return buffer[$ - len .. $].idup;
572 }
573
574 static if (maxSize64 == 3)
575 ///
576 version(mir_bignum_test) @safe pure unittest
577 {
578 auto str = "-34010447314490204552169750449563978034784726557588085989975288830070948234680";
579 auto integer = BigInt!4(str);
580 assert(integer.toString == str);
581
582 integer = BigInt!4.init;
583 assert(integer.toString == "0");
584 }
585
586 ///
587 void toString(C = char, W)(scope ref W w) const
588 if(isSomeChar!C && isMutable!C)
589 {
590 C[ceilLog10Exp2(data.length * (size_t.sizeof * 8)) + 1] buffer = void;
591 BigInt copy = this;
592 auto len = copy.view.unsigned.toStringImpl(buffer);
593 if (sign)
594 buffer[$ - ++len] = '-';
595 w.put(buffer[$ - len .. $]);
596 }
597
598 static if (maxSize64 == 3)
599 /// Check @nogc toString impl
600 version(mir_bignum_test) @safe pure @nogc unittest
601 {
602 import mir.format: stringBuf;
603 auto str = "-34010447314490204552169750449563978034784726557588085989975288830070948234680";
604 auto integer = BigInt!4(str);
605 stringBuf buffer;
606 buffer << integer;
607 assert(buffer.data == str, buffer.data);
608 }
609 }
610
611
612 ///
613 version(mir_bignum_test)
614 unittest
615 {
616 import mir.bignum.fixed;
617 import mir.bignum.low_level_view;
618
619 auto a = BigInt!4.fromHexString("4b313b23aa560e1b0985f89cbe6df5460860e39a64ba92b4abdd3ee77e4e05b8");
620 auto b = BigInt!4.fromHexString("c39b18a9f06fd8e962d99935cea0707f79a222050aaeaaaed17feb7aa76999d7");
621 auto c = BigInt!4.fromHexString("7869dd864619cace5953a09910327b3971413e6aa5f417fa25a2ac93291b941f");
622 c.sign = true;
623 assert(a != b);
624 assert(a < b);
625 a -= b;
626 assert(a.sign);
627 assert(a == c);
628 a -= a;
629 assert(!a.sign);
630 assert(!a.length);
631
632 auto d = BigInt!4.fromHexString("0de1a911c6dc8f90a7169a148e65d22cf34f6a8254ae26362b064f26ac44218a");
633 assert((b *= 0x7869dd86) == 0x5c019770);
634 assert(b == d);
635
636 d = BigInt!4.fromHexString("856eeb23e68cc73f2a517448862cdc97e83f9dfa23768296724bf00fda7df32a");
637 auto o = b *= UInt!128.fromHexString("f79a222050aaeaaa417fa25a2ac93291");
638 assert(o == UInt!128.fromHexString("d6d15b99499b73e68c3331eb0f7bf16"));
639 assert(b == d);
640
641 d = BigInt!4.fromHexString("d"); // initial value
642 d.mulPow5(60);
643 c = BigInt!4.fromHexString("81704fcef32d3bd8117effd5c4389285b05d");
644 assert(d == c);
645
646 d >>= 80;
647 c = BigInt!4.fromHexString("81704fcef32d3bd8");
648 assert(d == c);
649
650 c = BigInt!4.fromHexString("c39b18a9f06fd8e962d99935cea0707f79a222050aaeaaaed17feb7aa76999d7");
651 d = BigInt!4.fromHexString("9935cea0707f79a222050aaeaaaed17feb7aa76999d700000000000000000000");
652 c <<= 80;
653 assert(d == c);
654 c >>= 80;
655 c <<= 84;
656 d <<= 4;
657 assert(d == c);
658 assert(c != b);
659 b.sign = true;
660 assert(!c.copyFrom(b.view));
661 assert(c == b);
662 b >>= 18;
663 auto bView = cast(BigIntView!ushort)b.view;
664 assert(!c.copyFrom(bView.topLeastSignificantPart(bView.unsigned.coefficients.length - 1)));
665 assert(c == b);
666 }
667
668 version(mir_bignum_test)
669 @safe pure @nogc unittest
670 {
671 BigInt!4 i = "-0";
672 assert(i.view.coefficients.length == 0);
673 assert(cast(long) i == 0);
674 }