Unofficial Dub Package Documentation 
1 /** High-level wrapper for GNU Multiple Precision (MP) library.
2 
3     TODO _ccc Make _ccc two => uint _initCCC and uint _opCCC, mutatingCallCount
4     return their sum. Use these counters to make more precise asserts on how
5     assignments and operators are lowered to calls to C __gmpz-calls
6 
7     TODO disallow construction and assignment from floating point? Check with
8     other GMP interfaces and std.bigint.
9 
10     TODO Should `in` operator be used for anything good?
11 
12     TODO Support operators
13           <<      >>      >>>
14           &= |= ^=
15 */
16 module gmp;
17 
18 debug import core.stdc.stdio : printf;
19 import std.typecons : Unsigned, Unqual;
20 import std.meta : AliasSeq;
21 import std.traits : isInstanceOf;
22 import std.algorithm.mutation : move, moveEmplace;
23 
24 /// Call unittests taking long to execute.
25 enum unittestLong = false;
26 
27 version(unittest)
28 {
29     version = ccc;              // do C mutation call count
30 }
31 
32 /** Is `true` iff `T` is a GNU MP arithmetic type (`long`, `ulong` or `double`). */
33 enum isGMPArithmetic(T) = is(T == long) && is(T == ulong) && is(T == double);
34 
35 /// Is `true` if type `T` can be evaluated to a `MpZ` value.
36 enum isMpZExpr(T) = (is(Unqual!(typeof(T.eval())) == MpZ)); // which returns an `MpZ`
37 /// Is `true` if type `T` is lazy (not yet evaluated) `MpZ`-expression.
38 enum isLazyMpZExpr(T) = (!is(Unqual!T == MpZ) &&            // exclude direct value
39                          isMpZExpr!T);
40 
41 version(unittest)
42 {
43     static assert(!isMpZExpr!int);
44 }
45 
46 // TODO use these imports instead of the ones below
47 // import deimos.gmp.gmp;
48 // import deimos.gmp.integer;
49 
50 /** Arbitrary (multi) precision signed integer (Z).
51     Wrapper for GNU MP (GMP)'s type `mpz_t` and functions `mpz_.*`.
52  */
53 struct MpZ
54 {
55     /// Default conversion base.
56     private enum defaultBase = 10;
57 
58     pure nothrow pragma(inline, true):
59 
60     /// Convert to `string` in base `base`.
61     string toString(uint base = defaultBase, bool upperCaseDigits = false) const @trusted
62     {
63         assert((base >= -2 && base <= -36) ||
64                (base >= 2 && base <= 62));
65 
66         immutable size = sizeInBase(base);
67 
68         string str = new string(size + 1); // one extra for minus sign
69         __gmpz_get_str(cast(char*)str.ptr, base, _ptr); // fill it
70 
71         if (upperCaseDigits)
72         {
73             foreach (ref e; cast(char[])str)
74             {
75                 import std.ascii : isLower, toUpper;
76                 if (e.isLower)
77                 {
78                     e = e.toUpper;
79                 }
80             }
81         }
82 
83         return str[0] == '-' ? str : str.ptr[0 .. size];
84     }
85 
86     // TODO toRCString wrapped in UniqueRange
87 
88     // Returns: A unique hash of the `MpZ` value suitable for use in a hash table.
89     size_t toHash() const
90     {
91         import core.internal.hash : hashOf;
92         typeof(return) hash = _limbCount;
93         foreach (immutable i; 0 .. _limbCount)
94         {
95             hash ^= _limbs[i].hashOf;
96         }
97         return hash;
98     }
99 
100     @nogc:
101 
102     /// No default construction.
103     @disable this();
104 
105     /// Construct empty (undefined) from explicit `null`.
106     this(typeof(null)) @trusted
107     {
108         initialize();             // TODO is there a faster way?
109         assert(this == MpZ.init); // if this is same as default
110     }
111 
112     /// Construct from expression `expr`.
113     this(Expr)(Expr expr)
114         if (isLazyMpZExpr!Expr)
115     {
116         // TODO ok to just assume zero-initialized contents at `_z` before...
117         this = expr;            // ...calling this.opAssign ?x
118     }
119 
120     /** Construct from `value`. */
121     this(T)(T value) @trusted
122         if (isArithmetic!T)
123     {
124         version(ccc) ++_ccc;
125         static      if (isUnsigned!T)
126             __gmpz_init_set_ui(_ptr, value);
127         else static if (isFloating!T)
128             __gmpz_init_set_d(_ptr, value);
129         else                    // isSigned integral
130             __gmpz_init_set_si(_ptr, value);
131     }
132 
133     /** Construct from `value` in base `base`.
134         If `base` is 0 it's guessed from contents of `value`.
135         */
136     pragma(inline, false)
137     this(in string value, uint base = 0) @trusted // TODO Use Optional/Nullable when value is nan, or inf
138     {
139         assert(base == 0 || (base >= 2 && base <= 62));
140         char* stringz = _allocStringzCopyOf(value);
141         immutable int status = __gmpz_init_set_str(_ptr, stringz, base); version(ccc) ++_ccc;
142         qualifiedFree(stringz);
143         assert(status == 0, "Parameter `value` does not contain an integer");
144     }
145 
146     /// Returns: the Mersenne prime, M(p) = 2 ^^ p - 1
147     static MpZ mersennePrime(Integral)(Integral p)
148         if (isIntegral!Integral)
149     {
150         return typeof(this).pow(2UL, p) - 1;
151     }
152 
153     /// Use copy construction.
154     enum useCopy = false;       // disable copy construction for now
155 
156     static if (useCopy)
157     {
158         /// Construct copy of `value`.
159         this(ref const MpZ value) @trusted
160         {
161             __gmpz_init_set(_ptr, value._ptr); version(ccc) ++_ccc;
162         }
163         /// Construct copy of `value`.
164         this(MpZ value) @trusted
165         {
166             moveEmplace(value, this); // fast
167         }
168     }
169     else
170     {
171         /// Disable copy construction.
172         @disable this(this);
173     }
174 
175     /** Initialize internal struct. */
176     private void initialize() @trusted // cannot be called `init` as that will override builtin type property
177     {
178         __gmpz_init(_ptr); version(ccc) ++_ccc;
179     }
180 
181     /// Swap content of `this` with `rhs`.
182     void swap(ref MpZ rhs) @safe
183     {
184         import std.algorithm.mutation : swap;
185         swap(this, rhs); // faster than __gmpz_swap(_ptr, rhs._ptr); version(ccc) ++_ccc;
186     }
187 
188     /// Returns: (duplicate) copy of `this`.
189     MpZ dup() const @trusted
190     {
191         typeof(return) y = void;
192         __gmpz_init_set(y._ptr, _ptr); ++y._ccc;
193         return y;
194     }
195 
196     /// Assign from `rhs`.
197     ref MpZ opAssign()(auto ref const MpZ rhs) return @trusted // TODO scope
198     {
199         __gmpz_set(_ptr, rhs._ptr); version(ccc) ++_ccc;
200         return this;
201     }
202     /// ditto
203     ref MpZ opAssign(Expr)(auto ref Expr rhs) return @trusted // TODO scope
204         if (isLazyMpZExpr!Expr)
205     {
206         static      if (isInstanceOf!(MpzAddExpr, Expr))
207         {
208             __gmpz_add(this._ptr,
209                        rhs.e1.eval()._ptr,
210                        rhs.e2.eval()._ptr); version(ccc) ++_ccc;
211         }
212         else static if (isInstanceOf!(MpzSubExpr, Expr))
213         {
214             __gmpz_sub(this._ptr,
215                        rhs.e1.eval()._ptr,
216                        rhs.e2.eval()._ptr); version(ccc) ++_ccc;
217         }
218         else static if (isInstanceOf!(MpzMulExpr, Expr))
219         {
220             __gmpz_mul(this._ptr,
221                        rhs.e1.eval()._ptr,
222                        rhs.e2.eval()._ptr); version(ccc) ++_ccc;
223         }
224         else static if (isInstanceOf!(MpzDivExpr, Expr))
225         {
226             __gmpz_tdiv_q(this._ptr,
227                           rhs.e1.eval()._ptr,
228                           rhs.e2.eval()._ptr); version(ccc) ++_ccc;
229         }
230         else static if (isInstanceOf!(MpzModExpr, Expr))
231         {
232             __gmpz_tdiv_r(this._ptr,
233                           rhs.e1.eval()._ptr,
234                           rhs.e2.eval()._ptr); version(ccc) ++_ccc;
235         }
236         else static if (isInstanceOf!(MpzNegExpr, Expr))
237         {
238             __gmpz_neg(this._ptr,
239                        rhs.e1.eval()._ptr); version(ccc) ++_ccc;
240         }
241         else
242         {
243             this = rhs.eval();     // evaluate and move
244         }
245         return this;
246     }
247     /// ditto
248     ref MpZ opAssign(T)(T rhs) return @trusted // TODO scope
249         if (isArithmetic!T)
250     {
251         static if      (isUnsigned!T)
252         {
253             __gmpz_set_ui(_ptr, rhs);
254         }
255         else static if (isFloating!T)
256         {
257             __gmpz_set_d(_ptr, rhs);
258         }
259         else static if (isSigned!T)
260         {
261             pragma(msg, T);
262             __gmpz_set_si(_ptr, rhs);
263         }
264         else
265         {
266             static assert(false);
267         }
268         version(ccc) ++_ccc;
269         return this;
270     }
271 
272     /** Assign `this` from `string` `rhs` interpreted in base `base`.
273         If `base` is 0 it's guessed from contents of `value`.
274     */
275     ref MpZ fromString(in string rhs, uint base = 0) return @trusted // TODO scope
276     {
277         assert(base == 0 || (base >= 2 && base <= 62));
278         char* stringz = _allocStringzCopyOf(rhs);
279         immutable int status = __gmpz_set_str(_ptr, stringz, base); version(ccc) ++_ccc;
280         qualifiedFree(stringz);
281         assert(status == 0, "Parameter `rhs` does not contain an integer");
282         return this;
283     }
284 
285     /// Destruct `this`.
286     ~this() @trusted
287     {
288         if (_ptr)
289         {
290             __gmpz_clear(_ptr); version(ccc) ++_ccc;
291         }
292     }
293 
294     /// Returns: `true` iff `this` equals `rhs`.
295     bool opEquals()(auto ref const MpZ rhs) const @trusted // TODO scope
296     {
297         if (_ptr == rhs._ptr)   // fast equality
298         {
299             return true;        // fast bailout
300         }
301         return __gmpz_cmp(_ptr, rhs._ptr) == 0;
302     }
303     /// ditto
304     bool opEquals(Rhs)(Rhs rhs) const @trusted // TODO scope
305         if (isArithmetic!Rhs)
306     {
307         if (rhs == 0)
308         {
309             return isZero;      // optimization
310         }
311         static      if (isUnsigned!Rhs)
312         {
313             return __gmpz_cmp_ui(_ptr, cast(ulong)rhs) == 0;
314         }
315         else static if (isFloating!Rhs)
316         {
317             return __gmpz_cmp_d(_ptr, cast(double)rhs) == 0; // TODO correct to do this cast here?
318         }
319         else                    // isSigned integral
320         {
321             return __gmpz_cmp_si(_ptr, cast(long)rhs) == 0;
322         }
323     }
324 
325     /// Compare `this` to `rhs`.
326     int opCmp()(auto ref const MpZ rhs) const @trusted // TODO scope
327     {
328         if (rhs == 0)
329         {
330             return sgn();       // optimization
331         }
332         return __gmpz_cmp(_ptr, rhs._ptr);
333     }
334 
335     /// ditto
336     int opCmp(T)(T rhs) const @trusted // TODO scope
337         if (isArithmetic!T)
338     {
339         if (rhs == 0)
340         {
341             return sgn();       // optimization
342         }
343         static      if (isUnsigned!T)
344         {
345             return __gmpz_cmp_ui(_ptr, rhs);
346         }
347         else static if (isFloating!T)
348         {
349             return __gmpz_cmp_d(_ptr, rhs);
350         }
351         else                    // isSigned integral
352         {
353             return __gmpz_cmp_si(_ptr, rhs);
354         }
355     }
356 
357     /// Cast to `bool`.
358     bool opCast(T : bool)() const // TODO scope
359     {
360         return !isZero;
361     }
362 
363     /// Cast to arithmetic type `T`.
364     T opCast(T)() const @trusted // TODO scope
365         if (isArithmetic!T)
366     {
367         static      if (isUnsigned!T)
368         {
369             return cast(T)__gmpz_get_ui(_ptr);
370         }
371         else static if (isFloating!T)
372         {
373             return cast(T)__gmpz_get_d(_ptr);
374         }
375         else                    // isSigned integral
376         {
377             return cast(T)__gmpz_get_si(_ptr);
378         }
379     }
380 
381     /** Returns: The value of this as a `long`, or +/- `long.max` if outside
382         the representable range.
383     */
384     long toLong() const @trusted
385     {
386         // TODO can probably be optimized
387         if (this <= long.min)
388         {
389             return long.min;
390         }
391         else if (this >= long.max)
392         {
393             return long.max;
394         }
395         else
396         {
397             return cast(long)__gmpz_get_si(_ptr);
398         }
399     }
400 
401     /** Returns: The value of this as a `int`, or +/- `int.max` if outside
402         the representable range.
403     */
404     int toInt() const @trusted
405     {
406         // TODO can probably be optimized
407         if (this <= int.min)
408         {
409             return int.min;
410         }
411         else if (this >= int.max)
412         {
413             return int.max;
414         }
415         else
416         {
417             return cast(int)__gmpz_get_si(_ptr);
418         }
419     }
420 
421     /** Returns: `this` `s` `rhs`. */
422     MpZ opBinary(string s)(auto ref const MpZ rhs) const @trusted // direct value
423         if ((s == "+" || s == "-" ||
424              s == "*" || s == "/" || s == "%" ||
425              s == "&" || s == "|" || s == "^"))
426     {
427         static if (!__traits(isRef, rhs)) // r-value `rhs`
428         {
429             MpZ* mut_rhs = (cast(MpZ*)(&rhs)); // @trusted because `MpZ` has no aliased indirections
430             static      if (s == "+")
431             {
432                 __gmpz_add(mut_rhs._ptr, _ptr, rhs._ptr); version(ccc) ++mut_rhs._ccc;
433             }
434             else static if (s == "-")
435             {
436                 __gmpz_sub(mut_rhs._ptr, _ptr, rhs._ptr); version(ccc) ++mut_rhs._ccc;
437             }
438             else static if (s == "*")
439             {
440                 __gmpz_mul(mut_rhs._ptr, _ptr, rhs._ptr); version(ccc) ++mut_rhs._ccc;
441             }
442             else static if (s == "/")
443             {
444                 assert(rhs != 0, "Divison by zero");
445                 __gmpz_tdiv_q(mut_rhs._ptr, _ptr, rhs._ptr); version(ccc) ++mut_rhs._ccc;
446             }
447             else static if (s == "%")
448             {
449                 assert(rhs != 0, "Divison by zero");
450                 __gmpz_tdiv_r(mut_rhs._ptr, _ptr, rhs._ptr); version(ccc) ++mut_rhs._ccc;
451             }
452             else static if (s == "&")
453             {
454                 __gmpz_and(mut_rhs._ptr, _ptr, rhs._ptr); version(ccc) ++mut_rhs._ccc;
455             }
456             else static if (s == "|")
457             {
458                 __gmpz_ior(mut_rhs._ptr, _ptr, rhs._ptr); version(ccc) ++mut_rhs._ccc;
459             }
460             else static if (s == "^")
461             {
462                 __gmpz_xor(mut_rhs._ptr, _ptr, rhs._ptr); version(ccc) ++mut_rhs._ccc;
463             }
464             else
465             {
466                 static assert(false);
467             }
468             return move(*mut_rhs); // TODO shouldn't have to call `move` here
469         }
470         else
471         {
472             typeof(return) y = null;
473             static      if (s == "+")
474             {
475                 __gmpz_add(y._ptr, _ptr, rhs._ptr); version(ccc) ++y._ccc;
476             }
477             else static if (s == "-")
478             {
479                 __gmpz_sub(y._ptr, _ptr, rhs._ptr); version(ccc) ++y._ccc;
480             }
481             else static if (s == "*")
482             {
483                 __gmpz_mul(y._ptr, _ptr, rhs._ptr); version(ccc) ++y._ccc;
484             }
485             else static if (s == "/")
486             {
487                 assert(rhs != 0, "Divison by zero");
488                 __gmpz_tdiv_q(y._ptr, _ptr, rhs._ptr); version(ccc) ++y._ccc;
489             }
490             else static if (s == "%")
491             {
492                 assert(rhs != 0, "Divison by zero");
493                 __gmpz_tdiv_r(y._ptr, _ptr, rhs._ptr); version(ccc) ++y._ccc;
494             }
495             else static if (s == "&")
496             {
497                 __gmpz_and(y._ptr, _ptr, rhs._ptr); version(ccc) ++y._ccc;
498             }
499             else static if (s == "|")
500             {
501                 __gmpz_ior(y._ptr, _ptr, rhs._ptr); version(ccc) ++y._ccc;
502             }
503             else static if (s == "^")
504             {
505                 __gmpz_xor(y._ptr, _ptr, rhs._ptr); version(ccc) ++y._ccc;
506             }
507             else
508             {
509                 static assert(false);
510             }
511             return y;
512         }
513     }
514 
515     /// ditto
516     MpZ opBinary(string s, Rhs)(auto ref const Rhs rhs) const
517         if (isLazyMpZExpr!Rhs && // lazy value
518             (s == "+" || s == "-" || s == "*" || s == "/" || s == "%"))
519     {
520         static assert(false, "TODO");
521     }
522 
523     MpZ opBinary(string s, Rhs)(Rhs rhs) const @trusted
524         if ((s == "+" || s == "-" || s == "*" || s == "/" || s == "^^") &&
525             isUnsigned!Rhs)
526     {
527         typeof(return) y = null;
528         version(ccc) ++y._ccc;
529         static      if (s == "+")
530         {
531             __gmpz_add_ui(y._ptr, _ptr, rhs);
532         }
533         else static if (s == "-")
534         {
535             __gmpz_sub_ui(y._ptr, _ptr, rhs);
536         }
537         else static if (s == "*")
538         {
539             __gmpz_mul_ui(y._ptr, _ptr, rhs);
540         }
541         else static if (s == "/")
542         {
543             assert(rhs != 0, "Divison by zero");
544             __gmpz_tdiv_q_ui(y._ptr, _ptr, rhs);
545         }
546         else static if (s == "^^")
547         {
548             __gmpz_pow_ui(y._ptr, _ptr, rhs);
549         }
550         else
551         {
552             static assert(false);
553         }
554         return y;
555     }
556 
557     MpZ opBinary(string s, Rhs)(Rhs rhs) const @trusted
558         if ((s == "+" || s == "-" || s == "*" || s == "/" || s == "^^") &&
559             isSigned!Rhs)
560     {
561         typeof(return) y = null;
562         version(ccc) ++y._ccc;
563         static      if (s == "+")
564         {
565             if (rhs < 0)        // TODO handle `rhs == rhs.min`
566             {
567                 immutable ulong pos_rhs = -rhs; // make it positive
568                 __gmpz_sub_ui(y._ptr, _ptr, pos_rhs);
569             }
570             else
571             {
572                 __gmpz_add_ui(y._ptr, _ptr, rhs);
573             }
574         }
575         else static if (s == "-")
576         {
577             if (rhs < 0)        // TODO handle `rhs == rhs.min`
578             {
579                 __gmpz_add_ui(y._ptr, _ptr, -rhs); // x - (-y) == x + y
580             }
581             else
582             {
583                 __gmpz_sub_ui(y._ptr, _ptr, rhs); // rhs is positive
584             }
585         }
586         else static if (s == "*")
587         {
588             __gmpz_mul_si(y._ptr, _ptr, rhs);
589         }
590         else static if (s == "/")
591         {
592             assert(rhs != 0, "Divison by zero");
593             if (rhs < 0)        // TODO handle `rhs == rhs.min`
594             {
595                 immutable ulong pos_rhs = -rhs; // make it positive
596                 __gmpz_tdiv_q_ui(y._ptr, _ptr, pos_rhs);
597                 y.negate();     // negate result
598             }
599             else
600             {
601                 __gmpz_tdiv_q_ui(y._ptr, _ptr, rhs);
602             }
603         }
604         else static if (s == "^^")
605         {
606             assert(rhs >= 0, "TODO Negative power exponent needs MpQ return");
607             __gmpz_pow_ui(y._ptr, _ptr, rhs);
608         }
609         else
610         {
611             static assert(false);
612         }
613         return y;
614     }
615 
616     /// Remainer propagates modulus type.
617     Unqual!Rhs opBinary(string s, Rhs)(Rhs rhs) const @trusted
618         if ((s == "%") &&
619             isIntegral!Rhs)
620     {
621         assert(rhs != 0, "Divison by zero");
622         MpZ y = null;
623         version(ccc) ++y._ccc;
624         static if (isSigned!Rhs)
625         {
626             if (rhs < 0)        // TODO handle `rhs == rhs.min`
627             {
628                 immutable ulong pos_rhs = -rhs; // make it positive
629                 return -cast(typeof(return))__gmpz_tdiv_r_ui(y._ptr, _ptr, pos_rhs);
630             }
631             else
632             {
633                 return cast(typeof(return))__gmpz_tdiv_r_ui(y._ptr, _ptr, rhs);
634             }
635         }
636         else static if (isUnsigned!Rhs)
637         {
638             return cast(typeof(return))__gmpz_tdiv_r_ui(y._ptr, _ptr, rhs);
639         }
640         else
641         {
642             static assert(false);
643         }
644     }
645 
646     /// Returns: an unsigned type `lhs` divided by `this`.
647     MpZ opBinaryRight(string s, Lhs)(Lhs lhs) const @trusted
648         if ((s == "+" || s == "-" || s == "*" || s == "%") &&
649             isUnsigned!Lhs)
650     {
651         typeof(return) y = null;
652         version(ccc) ++y._ccc;
653         static      if (s == "+")
654         {
655             __gmpz_add_ui(y._ptr, _ptr, lhs); // commutative
656         }
657         else static if (s == "-")
658         {
659             __gmpz_ui_sub(y._ptr, lhs, _ptr);
660         }
661         else static if (s == "*")
662         {
663             __gmpz_mul_ui(y._ptr, _ptr, lhs); // commutative
664         }
665         else static if (s == "%")
666         {
667             assert(this != 0, "Divison by zero");
668             __gmpz_tdiv_r(y._ptr, MpZ(lhs)._ptr, _ptr); // convert `lhs` to MpZ
669         }
670         else
671         {
672             static assert(false);
673         }
674         return y;
675     }
676 
677     /// Returns: a signed type `lhs` divided by `this`.
678     MpZ opBinaryRight(string s, Lhs)(Lhs lhs) const @trusted
679         if ((s == "+" || s == "-" || s == "*" || s == "%") &&
680             isSigned!Lhs)
681     {
682         static if (s == "+" || s == "*")
683         {
684             return opBinary!s(lhs); // commutative reuse
685         }
686         else static if (s == "-")
687         {
688             typeof(return) y = null;
689             version(ccc) ++y._ccc;
690             if (lhs < 0)        // TODO handle `lhs == lhs.min`
691             {
692                 immutable ulong pos_rhs = -lhs; // make it positive
693                 __gmpz_add_ui(y._ptr, _ptr, pos_rhs);
694             }
695             else
696             {
697                 __gmpz_sub_ui(y._ptr, _ptr, lhs);
698             }
699             y.negate();
700             return y;
701         }
702         else static if (s == "%")
703         {
704             typeof(return) y = null;
705             version(ccc) ++y._ccc;
706             assert(this != 0, "Divison by zero");
707             __gmpz_tdiv_r(y._ptr, MpZ(lhs)._ptr, _ptr); // convert `lhs` to MpZ
708             return y;
709         }
710         else
711         {
712             static assert(false);
713         }
714     }
715 
716     /// Dividend propagates quotient type to signed.
717     Unqual!Lhs opBinaryRight(string s, Lhs)(Lhs lhs) const @trusted
718         if ((s == "/") &&
719             isIntegral!Lhs)
720     {
721         MpZ y = null; // TODO avoid if !__traits(isRef, this)
722         version(ccc) ++y._ccc;
723         assert(this != 0, "Divison by zero");
724         __gmpz_tdiv_q(y._ptr, MpZ(lhs)._ptr, _ptr);
725         static      if (isSigned!Lhs)
726         {
727             return cast(typeof(return))y;
728         }
729         else static if (isUnsigned!Lhs)
730         {
731             import std.typecons : Signed;
732             return cast(Signed!(typeof(return)))y;
733         }
734         else
735         {
736             static assert(false);
737         }
738     }
739 
740     /// Exponentation.
741     MpZ opBinaryRight(string s, Base)(Base base) const
742         if ((s == "^^") &&
743             isIntegral!Base)
744     {
745         static assert(false, "Convert `this MpZ` exponent to `ulong` and calculate power via static method `pow()`");
746         // MpZ exp = null;
747         // __gmpz_pow();
748         // return exp;
749     }
750 
751     ref MpZ opOpAssign(string s)(auto ref const MpZ rhs) return @trusted // TODO scope
752         if ((s == "+" || s == "-" ||
753              s == "*" || s == "/" || s == "%" ||
754              s == "&" || s == "|" || s == "^"))
755     {
756         static      if (s == "+")
757         {
758             __gmpz_add(_ptr, _ptr, rhs._ptr); version(ccc) ++_ccc;
759         }
760         else static if (s == "-")
761         {
762             __gmpz_sub(_ptr, _ptr, rhs._ptr); version(ccc) ++_ccc;
763         }
764         else static if (s == "*")
765         {
766             if (rhs == -1)
767             {
768                 negate();
769             }
770             else
771             {
772                 __gmpz_mul(_ptr, _ptr, rhs._ptr); version(ccc) ++_ccc;
773             }
774         }
775         else static if (s == "/")
776         {
777             assert(rhs != 0, "Divison by zero");
778             if (rhs == -1)
779             {
780                 negate();
781             }
782             else
783             {
784                 __gmpz_tdiv_q(_ptr, _ptr, rhs._ptr); version(ccc) ++_ccc;
785             }
786         }
787         else static if (s == "%")
788         {
789             assert(rhs != 0, "Divison by zero");
790             __gmpz_tdiv_r(_ptr, _ptr, rhs._ptr); version(ccc) ++_ccc;
791         }
792         else static if (s == "&")
793         {
794             __gmpz_and(_ptr, _ptr, rhs._ptr); version(ccc) ++_ccc;
795         }
796         else static if (s == "|")
797         {
798             __gmpz_ior(_ptr, _ptr, rhs._ptr); version(ccc) ++_ccc;
799         }
800         else static if (s == "^")
801         {
802             __gmpz_xor(_ptr, _ptr, rhs._ptr); version(ccc) ++_ccc;
803         }
804         else
805         {
806             static assert(false);
807         }
808         return this;
809     }
810 
811     ref MpZ opOpAssign(string s, Rhs)(Rhs rhs) return @trusted // TODO scope
812         if ((s == "+" || s == "-" || s == "*" || s == "/" || s == "%" || s == "^^") &&
813             isUnsigned!Rhs)
814     {
815         static      if (s == "+")
816         {
817             __gmpz_add_ui(_ptr, _ptr, rhs); version(ccc) ++_ccc;
818         }
819         else static if (s == "-")
820         {
821             __gmpz_sub_ui(_ptr, _ptr, rhs); version(ccc) ++_ccc;
822         }
823         else static if (s == "*")
824         {
825             __gmpz_mul_ui(_ptr, _ptr, rhs); version(ccc) ++_ccc;
826         }
827         else static if (s == "/")
828         {
829             assert(rhs != 0, "Divison by zero");
830             __gmpz_tdiv_q_ui(_ptr, _ptr, rhs); version(ccc) ++_ccc;
831         }
832         else static if (s == "%")
833         {
834             assert(rhs != 0, "Divison by zero");
835             __gmpz_tdiv_r_ui(_ptr, _ptr, rhs); version(ccc) ++_ccc;
836         }
837         else static if (s == "^^")
838         {
839             __gmpz_pow_ui(_ptr, _ptr, rhs); version(ccc) ++_ccc;
840         }
841         else
842         {
843             static assert(false);
844         }
845         return this;
846     }
847 
848     ref MpZ opOpAssign(string s, Rhs)(Rhs rhs) return @trusted // TODO scope
849         if ((s == "+" || s == "-" || s == "*" || s == "/" || s == "%" || s == "^^") &&
850             isSigned!Rhs)
851     {
852         static      if (s == "+")
853         {
854             if (rhs < 0)        // TODO handle `rhs == rhs.min`
855             {
856                 assert(rhs != rhs.min);
857                 immutable ulong pos_rhs = -rhs; // make it positive
858                 __gmpz_sub_ui(_ptr, _ptr, pos_rhs); version(ccc) ++_ccc;
859             }
860             else
861             {
862                 __gmpz_add_ui(_ptr, _ptr, rhs); version(ccc) ++_ccc;
863             }
864         }
865         else static if (s == "-")
866         {
867             if (rhs < 0)        // TODO handle `rhs == rhs.min`
868             {
869                 assert(rhs != rhs.min);
870                 immutable ulong pos_rhs = -rhs; // make it positive
871                 __gmpz_add_ui(_ptr, _ptr, pos_rhs); version(ccc) ++_ccc;
872             }
873             else
874             {
875                 __gmpz_sub_ui(_ptr, _ptr, rhs); version(ccc) ++_ccc;
876             }
877         }
878         else static if (s == "*")
879         {
880             if (rhs == -1)
881             {
882                 negate();       // optimization
883             }
884             else
885             {
886                 __gmpz_mul_si(_ptr, _ptr, rhs); version(ccc) ++_ccc;
887             }
888         }
889         else static if (s == "/")
890         {
891             assert(rhs != 0, "Divison by zero");
892             if (rhs < 0)        // TODO handle `rhs == rhs.min`
893             {
894                 immutable ulong pos_rhs = -rhs; // make it positive
895                 __gmpz_tdiv_q_ui(_ptr, _ptr, pos_rhs); version(ccc) ++_ccc;
896                 negate();
897             }
898             else
899             {
900                 __gmpz_tdiv_q_ui(_ptr, _ptr, rhs); version(ccc) ++_ccc;
901             }
902         }
903         else static if (s == "%")
904         {
905             assert(rhs != 0, "Divison by zero");
906             __gmpz_tdiv_r_ui(_ptr, _ptr, rhs); version(ccc) ++_ccc;
907         }
908         else static if (s == "^^")
909         {
910             assert(rhs >= 0, "Negative power exponent");
911             __gmpz_pow_ui(_ptr, _ptr, rhs); version(ccc) ++_ccc;
912         }
913         else
914         {
915             static assert(false);
916         }
917         return this;
918     }
919 
920     /// Returns: `this`.
921     ref inout(MpZ) opUnary(string s)() inout
922         if (s == "+")
923     {
924         return this;
925     }
926 
927     /// Returns: negation of `this`.
928     MpZ opUnary(string s)() const
929         if (s == "-")
930     {
931         typeof(return) y = this.dup;
932         y.negate();
933         return y;
934     }
935 
936     /// Returns: negation of `this`.
937     MpZ unaryMinus() const @safe
938     {
939         pragma(msg, "memberFun:", __traits(isRef, this) ? "ref" : "non-ref", " this");
940         typeof(return) y = this.dup;
941         y.negate();
942         return y;
943     }
944 
945     /** Negate `this` in-place.
946         Returns: `void` to make it obvious `this` is mutated.
947      */
948     void negate() @safe
949     {
950         _z._mp_size = -_z._mp_size; // fast C macro `mpz_neg` in gmp.h
951     }
952 
953     /** Make `this` the absolute value of itself in-place.
954         Returns: `void` to make it obvious `this` is mutated.
955      */
956     void absolute() @trusted
957     {
958         __gmpz_abs(_ptr, _ptr); version(ccc) ++_ccc;
959     }
960 
961     /// Increase `this` by one.
962     ref MpZ opUnary(string s)() return @trusted // TODO scope
963         if (s == "++")
964     {
965         __gmpz_add_ui(_ptr, _ptr, 1); version(ccc) ++_ccc;
966         return this;
967     }
968 
969     /// Decrease `this` by one.
970     ref MpZ opUnary(string s)() return @trusted // TODO scope
971         if (s == "--")
972     {
973         __gmpz_sub_ui(_ptr, _ptr, 1); version(ccc) ++_ccc;
974         return this;
975     }
976 
977     /// Returns: `base` raised to the power of `exp`.
978     static typeof(this) pow(Base, Exp)(Base base, Exp exp) @trusted
979         if (isIntegral!Base &&
980             isIntegral!Exp)
981     {
982         static if (isSigned!Base)
983         {
984             immutable bool negate = base < 0;
985             immutable ubase = cast(ulong)(negate ? -base : base);
986         }
987         else
988         {
989             immutable bool negate = false;
990             immutable ubase = base;
991         }
992 
993         static if (isSigned!Exp)
994         {
995             assert(exp >= 0, "Negative power exponent"); // TODO return mpq?
996             immutable uexp = cast(ulong)exp;
997         }
998         else
999         {
1000             immutable uexp = exp;
1001         }
1002 
1003         typeof(return) y = null;
1004         __gmpz_ui_pow_ui(y._ptr, ubase, uexp); version(ccc) ++y._ccc;
1005         if (negate && exp & 1)  // if negative odd exponent
1006         {
1007             y.negate();
1008         }
1009 
1010         return y;
1011     }
1012 
1013     /// Returns: number of digits in base `base`.
1014     size_t sizeInBase(uint base) const @trusted
1015     {
1016         return __gmpz_sizeinbase(_ptr, base);
1017     }
1018 
1019     /// Returns: `true` iff `this` fits in a `T`.
1020     bool fitsIn(T)() const @trusted
1021         if (isIntegral!T)
1022     {
1023         static      if (is(T == ulong))  { return __gmpz_fits_ulong_p(_ptr) != 0; }
1024         else static if (is(T ==  long))  { return __gmpz_fits_slong_p(_ptr) != 0; }
1025         else static if (is(T ==  uint))  { return __gmpz_fits_uint_p(_ptr) != 0; }
1026         else static if (is(T ==   int))  { return __gmpz_fits_sint_p(_ptr) != 0; }
1027         else static if (is(T == ushort)) { return __gmpz_fits_ushort_p(_ptr) != 0; }
1028         else static if (is(T ==  short)) { return __gmpz_fits_sshort_p(_ptr) != 0; }
1029         else { static assert(false, "Unsupported type " ~ T.stringof); }
1030     }
1031 
1032     /// Check if `this` is zero.
1033     @property bool isZero() const @safe
1034     {
1035         return _z._mp_size == 0; // fast
1036     }
1037 
1038     /// Check if `this` is odd. TODO use as specialcase in: this & 1
1039     @property bool isOdd() const @safe
1040     {
1041         return (_z._mp_size != 0) & cast(int)(_z._mp_d[0]); // fast C macro `mpz_odd_p` in gmp.h
1042     }
1043 
1044     /// Check if `this` is odd.
1045     @property bool isEven() const @safe
1046     {
1047         return !isOdd;            // fast C macro `mpz_even_p` in gmp.h
1048     }
1049 
1050     /// Check if `this` is negative. TODO use as specialcase in opCmp
1051     @property bool isNegative() const @safe
1052     {
1053         return _z._mp_size < 0; // fast
1054     }
1055 
1056     /// Check if `this` is positive.
1057     @property bool isPositive() const @safe
1058     {
1059         return !isNegative;     // fast
1060     }
1061 
1062     /** Returns: sign as either
1063         - -1 (`this` < 0),
1064         -  0 (`this` == 0), or
1065         - +1 (`this` > 0).
1066      */
1067     @property int sgn() const @safe
1068     {
1069         return _z._mp_size < 0 ? -1 : _z._mp_size > 0; // fast C macro `mpz_sgn` in gmp.h
1070     }
1071 
1072     /// Number of significant `uint`s used for storing `this`.
1073     @property size_t uintLength() const
1074     {
1075         assert(false, "TODO use mpz_size");
1076     }
1077 
1078     /// Number of significant `ulong`s used for storing `this`.
1079     @property size_t uintLong() const
1080     {
1081         assert(false, "TODO use mpz_size");
1082     }
1083 
1084 private:
1085 
1086     /// Returns: evaluation of `this` expression which in this is a no-op.
1087     ref inout(MpZ) eval() @safe inout return { return this; }
1088 
1089     /// Type of limb in internal representation.
1090     alias Limb = __mp_limb_t;   // GNU MP alias
1091 
1092     /** Returns: limbs. */
1093     inout(Limb)[] _limbs() inout return @system // TODO scope
1094     {
1095         // import std.math : abs;
1096         return _z._mp_d[0 .. _limbCount];
1097     }
1098 
1099     /// Get number of limbs in internal representation.
1100     @property uint _limbCount() const @safe
1101     {
1102         return _integralAbs(_z._mp_size);
1103     }
1104 
1105     /// @nogc-variant of `toStringz` with heap allocation of null-terminated C-string `stringz`.
1106     char* _allocStringzCopyOf(const string value) @trusted @nogc
1107     {
1108         char* stringz = cast(char*)qualifiedMalloc(value.length + 1); // maximum this many characters
1109         size_t i = 0;
1110         foreach (immutable char ch; value[])
1111         {
1112             if (ch != '_')
1113             {
1114                 stringz[i] = ch;
1115                 i += 1;
1116             }
1117         }
1118         stringz[i] = '\0'; // set C null terminator
1119         return stringz;
1120     }
1121 
1122     /// Returns: pointer to internal GMP C struct.
1123     inout(__mpz_struct)* _ptr() inout return @system // TODO scope
1124     {
1125         return &_z;
1126     }
1127 
1128     __mpz_struct _z;            // internal libgmp C struct
1129 
1130     // qualified C memory managment
1131     static @safe
1132     {
1133         pragma(mangle, "malloc") void* qualifiedMalloc(size_t size);
1134         pragma(mangle, "free") void qualifiedFree(void* ptr);
1135     }
1136 
1137     /** Number of calls made to `__gmpz`--functions that construct or changes
1138         this value. Used to verify correct lowering and evaluation of template
1139         expressions.
1140 
1141         For instance the `x` in `x = y + z` should be assigned only once inside
1142         a call to `mpz_add`.
1143      */
1144     @property size_t mutatingCallCount() const @safe { return _ccc; }
1145 
1146     version(ccc) size_t _ccc;  // C mutation call count. number of calls to C GMP function calls that mutate this object
1147 }
1148 
1149 version(unittest) static assert(isMpZExpr!MpZ);
1150 
1151 pure nothrow pragma(inline, true):
1152 
1153 /** Instantiator for `MpZ`. */
1154 MpZ mpz(Args...)(Args args) @safe @nogc
1155 {
1156     return typeof(return)(args);
1157 }
1158 
1159 /// Swap contents of `x` with contents of `y`.
1160 void swap()(ref MpZ x,
1161             ref MpZ y) @nogc
1162 {
1163     import std.algorithm.mutation : swap;
1164     swap(x, y); // x.swap(y);
1165 }
1166 
1167 /// Returns: `x` as a `string` in decimal base.
1168 string toDecimalString()(auto ref const MpZ x) // for `std.bigint.BigInt` compatibility
1169 {
1170     return x.toString(10);
1171 }
1172 
1173 /// Returns: `x` as a uppercased `string` in hexadecimal base without any base prefix (0x).
1174 string toHex()(auto ref const MpZ x) // for `std.bigint.BigInt` compatibility
1175 {
1176     return x.toString(16, true);
1177 }
1178 
1179 /// Returns: the absolute value of `x` converted to the corresponding unsigned type.
1180 Unsigned!T absUnsign(T)(auto ref const MpZ x) // for `std.bigint.BigInt` compatibility
1181     if (isIntegral!T)
1182 {
1183     return _integralAbs(cast(T)x);
1184 }
1185 
1186 /** Returns: absolute value of `x`.
1187     Written as a free function instead of `MpZ`-member because `__traits(isRef, this)` cannot be used.
1188  */
1189 MpZ abs()(auto ref const MpZ x) @trusted @nogc
1190 {
1191     static if (__traits(isRef, x)) // l-value `x`
1192     {
1193         typeof(return) y = null; // must use temporary
1194         __gmpz_abs(y._ptr, x._ptr); version(ccc) ++y._ccc;
1195         return y;
1196     }
1197     else                        // r-value `x`
1198     {
1199         MpZ* mut_x = (cast(MpZ*)(&x)); // @trusted because `MpZ` has no aliased indirections
1200         mut_x.absolute();
1201         return move(*mut_x);    // TODO shouldn't have to call `move` here
1202     }
1203 }
1204 
1205 /// Returns: comparison of the absolute values of `x` and `y`.
1206 int cmpabs()(auto ref const MpZ x,
1207              auto ref const MpZ y) @trusted @nogc
1208 {
1209     return __gmpz_cmpabs(x._ptr, y._ptr);
1210 }
1211 /// ditto
1212 int cmpabs()(auto ref const MpZ x, double y) @trusted @nogc
1213 {
1214     return __gmpz_cmpabs_d(x._ptr, y);
1215 }
1216 /// ditto
1217 int cmpabs()(auto ref const MpZ x, ulong y) @trusted @nogc
1218 {
1219     return __gmpz_cmpabs_ui(x._ptr, y);
1220 }
1221 
1222 /** Returns: next prime greater than `x`.
1223     Written as a free function instead of `MpZ`-member because `__traits(isRef, this)` cannot be used.
1224 */
1225 MpZ nextPrime()(auto ref const MpZ x) @trusted @nogc
1226 {
1227     static if (__traits(isRef, x)) // l-value `x`
1228     {
1229         typeof(return) y = null; // must use temporary
1230         __gmpz_nextprime(y._ptr, x._ptr); version(ccc) ++y._ccc;
1231         return y;
1232     }
1233     else                        // r-value `x`
1234     {
1235         MpZ* mut_x = (cast(MpZ*)(&x)); // @trusted because `MpZ` has no aliased indirections
1236         __gmpz_nextprime(mut_x._ptr, x._ptr); version(ccc) ++mut_x._ccc;
1237         return move(*mut_x);    // TODO shouldn't have to call `move` here
1238     }
1239 }
1240 
1241 /// Returns: greatest common divisor (gcd) of `x` and `y`.
1242 MpZ gcd()(auto ref const MpZ x,
1243           auto ref const MpZ y) @trusted @nogc
1244 {
1245     static      if (!__traits(isRef, x)) // r-value `x`
1246     {
1247         MpZ* mut_x = (cast(MpZ*)(&x)); // @trusted because `MpZ` has no aliased indirections
1248         __gmpz_gcd(mut_x._ptr, x._ptr, y._ptr); version(ccc) ++mut_x._ccc;
1249         return move(*mut_x);    // TODO shouldn't have to call `move` here
1250     }
1251     else static if (!__traits(isRef, y)) // r-value `y`
1252     {
1253         MpZ* mut_y = (cast(MpZ*)(&y)); // @trusted because `MpZ` has no aliased indirections
1254         __gmpz_gcd(mut_y._ptr, x._ptr, y._ptr); version(ccc) ++mut_y._ccc;
1255         return move(*mut_y);    // TODO shouldn't have to call `move` here
1256     }
1257     else                        // l-value `x` and `y`
1258     {
1259         MpZ z = null;
1260         __gmpz_gcd(z._ptr, x._ptr, y._ptr); version(ccc) ++z._ccc;
1261         return z;
1262     }
1263 }
1264 /// ditto
1265 MpZ gcd()(auto ref const MpZ x,
1266           ulong y) @trusted @nogc
1267 {
1268     static if (__traits(isRef, x)) // l-value `x`
1269     {
1270         MpZ z = null;
1271         const z_ui = __gmpz_gcd_ui(z._ptr, x._ptr, y); version(ccc) ++z._ccc;
1272         return z;
1273     }
1274     else
1275     {
1276         MpZ* mut_x = (cast(MpZ*)(&x)); // @trusted because `MpZ` has no aliased indirections
1277         const z_ui = __gmpz_gcd_ui(mut_x._ptr, x._ptr, y); version(ccc) ++mut_x._ccc;
1278         return move(*mut_x);    // TODO shouldn't have to call `move` here
1279     }
1280 }
1281 
1282 /// Returns: least common multiple (lcm) of `x` and `y`.
1283 MpZ lcm()(auto ref const MpZ x,
1284           auto ref const MpZ y) @trusted @nogc
1285 {
1286     static      if (!__traits(isRef, x)) // r-value `x`
1287     {
1288         MpZ* mut_x = (cast(MpZ*)(&x)); // @trusted because `MpZ` has no aliased indirections
1289         __gmpz_lcm(mut_x._ptr, x._ptr, y._ptr); version(ccc) ++mut_x._ccc;
1290         return move(*mut_x);    // TODO shouldn't have to call `move` here
1291     }
1292     else static if (!__traits(isRef, y)) // r-value `y`
1293     {
1294         MpZ* mut_y = (cast(MpZ*)(&y)); // @trusted because `MpZ` has no aliased indirections
1295         __gmpz_lcm(mut_y._ptr, x._ptr, y._ptr); version(ccc) ++mut_y._ccc;
1296         return move(*mut_y);    // TODO shouldn't have to call `move` here
1297     }
1298     else                        // l-value `x` and `y`
1299     {
1300         MpZ z = null;
1301         __gmpz_lcm(z._ptr, x._ptr, y._ptr); version(ccc) ++z._ccc;
1302         return z;
1303     }
1304 }
1305 /// ditto
1306 MpZ lcm()(auto ref const MpZ x,
1307           ulong y) @trusted @nogc
1308 {
1309     static if (__traits(isRef, x)) // l-value `x`
1310     {
1311         MpZ z = null;
1312         __gmpz_lcm_ui(z._ptr, x._ptr, y);
1313         return z;
1314     }
1315     else
1316     {
1317         MpZ* mut_x = (cast(MpZ*)(&x)); // @trusted because `MpZ` has no aliased indirections
1318         __gmpz_lcm_ui(mut_x._ptr, x._ptr, y); version(ccc) ++mut_x._ccc;
1319         return move(*mut_x);    // TODO shouldn't have to call `move` here
1320     }
1321 }
1322 
1323 /** Returns: `base` ^^ `exp` (modulo `mod`).
1324     Parameter `exp` must be positive.
1325 */
1326 MpZ powm()(auto ref const MpZ base,
1327            auto ref const MpZ exp,
1328            auto ref const MpZ mod) @trusted
1329 {
1330     typeof(return) y = 0; // result, TODO reuse `exp` or `mod` if any is an r-value
1331     assert(exp >= 0, "Negative exponent");
1332     __gmpz_powm(y._ptr, base._ptr, exp._ptr, mod._ptr); version(ccc) ++y._ccc;
1333     return y;
1334 }
1335 /// ditto
1336 MpZ powm()(auto ref const MpZ base,
1337            ulong exp,
1338            auto ref const MpZ mod) @trusted
1339 {
1340     typeof(return) y = 0;       // result, TODO reuse `exp` or `mod` if any is an r-value
1341     __gmpz_powm_ui(y._ptr, base._ptr, exp, mod._ptr); version(ccc) ++y._ccc;
1342     return y;
1343 }
1344 
1345 ///
1346 @safe @nogc unittest
1347 {
1348     const x = 42.Z;
1349     assert(x.unaryMinus() == -42);   // l-value `this`
1350     assert(42.Z.unaryMinus() == -42); // r-value `this`
1351 }
1352 
1353 /// convert to string
1354 @safe unittest
1355 {
1356     assert(mpz(    42).toString ==   `42`);
1357     assert(mpz(   -42).toString ==  `-42`);
1358     assert(mpz(`-101`).toString == `-101`);
1359 
1360     assert(mpz(-42).toDecimalString == `-42`);
1361     assert(mpz( 42).toDecimalString ==  `42`);
1362 
1363     assert(mpz( 0).toHex == `0`);
1364     assert(mpz( 1).toHex == `1`);
1365     assert(mpz( 9).toHex == `9`);
1366     assert(mpz(10).toHex == `A`);
1367     assert(mpz(14).toHex == `E`);
1368     assert(mpz(15).toHex == `F`);
1369     assert(mpz(16).toHex == `10`);
1370 
1371     assert(mpz(-42).absUnsign!ulong == 42);
1372     assert(mpz( 42).absUnsign!ulong == 42);
1373 }
1374 
1375 /// opBinary with r-value right-hand-side
1376 @safe @nogc unittest
1377 {
1378     Z a = 42;
1379     {
1380         Z b = a + 1.Z;              // r-value `rhs`
1381         assert(b.mutatingCallCount == 2);
1382         assert(b == 43);
1383     }
1384 
1385     {
1386         Z b = a - 1.Z;              // r-value `rhs`
1387         assert(b.mutatingCallCount == 2);
1388         assert(b == 41);
1389     }
1390 
1391     {
1392         Z b = a * 2.Z;              // r-value `rhs`
1393         assert(b.mutatingCallCount == 2);
1394         assert(b == 84);
1395     }
1396 
1397     {
1398         Z b = a / 2.Z;              // r-value `rhs`
1399         assert(b.mutatingCallCount == 2);
1400         assert(b == 21);
1401     }
1402 
1403     {
1404         Z b = a % 10.Z;              // r-value `rhs`
1405         assert(b.mutatingCallCount == 2);
1406         assert(b == 2);
1407     }
1408 }
1409 
1410 ///
1411 @safe @nogc unittest
1412 {
1413     const _ = (cast(uint)42).Z;
1414     const a = 42.Z;
1415     const b = 43UL.Z;
1416     const c = 43.0.Z;
1417     const z = 0.Z;
1418 
1419     // `opOpAssign` with `Unsigned`
1420     auto w = 42.Z;
1421     assert(w == 42);
1422 
1423     w += 100UL;
1424     assert(w == 142);
1425 
1426     w -= 100.Z;
1427     assert(w == 42);
1428 
1429     w += 100UL;
1430     assert(w == 142);
1431 
1432     w -= 100UL;
1433     assert(w == 42);
1434 
1435     w *= 100UL;
1436     assert(w == 4200);
1437 
1438     w /= 100UL;
1439     assert(w == 42);
1440 
1441     w %= 10UL;
1442     assert(w == 2);
1443 
1444     w ^^= 6UL;
1445     assert(w == 64);
1446 
1447     w = 42;
1448     assert(w == 42);
1449 
1450     w += 100;
1451     assert(w == 142);
1452 
1453     w -= 100;
1454     assert(w == 42);
1455 
1456     w *= 100;
1457     assert(w == 4200);
1458 
1459     w /= 100;
1460     assert(w == 42);
1461 
1462     w *= 100;
1463     assert(w == 4200);
1464 
1465     w /= -100;
1466     assert(w == -42);
1467 
1468     w *= -1;
1469     assert(w == 42);
1470 
1471     w %= 10;
1472     assert(w == 2);
1473 
1474     w = 2;
1475     w ^^= 6;
1476     assert(w == 64);
1477 
1478     w = 32.0;
1479     assert(w == 32);
1480 
1481     w = 42UL;
1482     assert(w == 42);
1483 
1484     w /= 2.Z;
1485     assert(w == 21);
1486 
1487     w /= -2.Z;
1488     assert(w == -10);
1489 
1490     w *= -2.Z;
1491     assert(w == 20);
1492 
1493     w %= 3.Z;
1494     assert(w == 2);
1495 
1496     w *= -1.Z;
1497     assert(w == -2);
1498 
1499     w /= -1.Z;
1500     assert(w == 2);
1501 
1502     // equality
1503     assert(z == 0);
1504     assert(z == cast(uint)0);
1505     assert(z == 0L);
1506     assert(z == 0UL);
1507     assert(z == 0.0f);
1508     assert(z == 0.0);
1509 
1510     // eval cast
1511 
1512     assert(a);
1513     assert(cast(ulong)a == a);
1514     assert(cast(ulong)a == 42);
1515     assert(cast(long)a == a);
1516     assert(cast(long)a == 42);
1517     assert(cast(double)a == 42.0);
1518 
1519     // binary
1520 
1521     assert(`0b11`.Z == 3);
1522     assert(`0B11`.Z == 3);
1523 
1524     // octal
1525 
1526     assert(`07`.Z == 7);
1527     assert(`010`.Z == 8);
1528 
1529     // hexadecimal
1530 
1531     assert(`0x10`.Z == 16);
1532     assert(`0X10`.Z == 16);
1533 
1534     // decimal
1535 
1536     assert(`101`.Z == 101);
1537     assert(`101`.Z == 101);
1538 
1539     immutable ic = 101UL.Z;
1540 
1541     assert(a == a.dup);
1542     assert(ic == ic.dup);
1543 
1544     // equality
1545 
1546     assert(a == a);
1547     assert(a == 42.Z);
1548     assert(a == 42.0);
1549     assert(a == 42);
1550     assert(a == cast(uint)42);
1551     assert(a == 42UL);
1552     assert(_ == 42);
1553     assert(c == 43.0);
1554 
1555     // non-equality
1556 
1557     assert(a != b);
1558 
1559     // less than
1560 
1561     assert(a < b);
1562     assert(a < 43.Z);
1563     assert(a < 43);
1564     assert(a < cast(uint)43);
1565     assert(a < 43UL);
1566     assert(a < 43.0);
1567 
1568     assert(-1.Z < 0.Z);
1569     assert(-1.Z < 0L);
1570     assert(-1.Z < 0UL);
1571     assert(-1.Z < 0.0);
1572 
1573     // greater than
1574 
1575     assert(b > a);
1576     assert(b > 42.Z);
1577     assert(b > 42);
1578     assert(b > cast(uint)42);
1579     assert(b > 42UL);
1580     assert(b > 42.0);
1581 
1582     assert(+1.Z > 0.Z);
1583     assert(+1.Z > 0L);
1584     assert(+1.Z > 0UL);
1585     assert(+1.Z > 0.0);
1586 
1587     // absolute value
1588 
1589     assert(abs(a) == a);        // free function
1590     assert(a.abs == a);         // UFCS
1591     assert(abs(-42.Z) == 42);
1592     assert(abs(-a) == a);
1593 
1594     // absolute value comparison
1595 
1596     assert(cmpabs(-43.Z,  44.Z) == -1);
1597     assert(cmpabs(-43.Z, -44.Z) == -1);
1598     assert(cmpabs(-44.Z, -43.Z) == +1);
1599     assert(cmpabs(-43.Z, -43.Z) ==  0);
1600 
1601     assert(cmpabs(-43.Z,  44.0) == -1);
1602     assert(cmpabs(-43.Z, -44.0) == -1);
1603     assert(cmpabs(-44.Z, -43.0) == +1);
1604     assert(cmpabs(-43.Z, -43.0) ==  0);
1605 
1606     assert(cmpabs(-43.Z, 44) == -1);
1607     assert(cmpabs( 43.Z, 44) == -1);
1608     assert(cmpabs(-44.Z, 43) == +1);
1609     assert(cmpabs( 44.Z, 43) == +1);
1610     assert(cmpabs(-43.Z, 43) ==  0);
1611 
1612     Z _43 = 43;
1613     Z _4 = 4;
1614     Z _24 = 24;
1615 
1616     // next prime
1617     assert(nextPrime(_4) == 5);
1618     assert(nextPrime(24.Z) == 29);
1619 
1620     assert(nextPrime(_24) == 29);
1621     assert(nextPrime(24.Z) == 29);
1622     assert(24.Z.nextPrime() == 29);
1623 
1624     assert(nextPrime(_43) == 47);
1625     assert(nextPrime(43.Z) == 47);
1626     assert(43.Z.nextPrime() == 47);
1627 
1628     // greatest common divisor
1629 
1630     assert(gcd(43.Z,  44.Z) == 1);
1631     assert(gcd(4.Z, 24.Z) == 4);
1632     assert(gcd(6.Z, 24.Z) == 6);
1633     assert(gcd(10.Z, 100.Z) == 10);
1634 
1635     assert(gcd(43.Z,  44) == 1);
1636     assert(gcd(4.Z, 24) == 4);
1637     assert(gcd(6.Z, 24) == 6);
1638     assert(gcd(10.Z, 100) == 10);
1639 
1640     assert(gcd(_43,  44) == 1);
1641     assert(gcd(_4, 24) == 4);
1642     assert(gcd(_4, _24) == 4);
1643     assert(gcd(_4, 24.Z) == 4);
1644 
1645     // least common multiple
1646 
1647     assert(lcm(43.Z,  44.Z) == 1892);
1648     assert(lcm(4.Z, 24.Z) == 24);
1649     assert(lcm(6.Z, 24.Z) == 24);
1650     assert(lcm(10.Z, 100.Z) == 100);
1651 
1652     assert(lcm(43.Z,  44) == 1892);
1653     assert(lcm(4.Z, 24) == 24);
1654     assert(lcm(6.Z, 24) == 24);
1655     assert(lcm(10.Z, 100) == 100);
1656 
1657     assert(lcm(_43,  44) == 1892);
1658     assert(lcm(_4, 24) == 24);
1659     assert(lcm(_4, _24) == 24);
1660     assert(lcm(_4, 24.Z) == 24);
1661 
1662     // negated value
1663 
1664     assert(-a == -42);
1665     assert(-(-a) == a);
1666 
1667     auto n = 42.Z;
1668     n.negate();
1669     assert(n == -42);
1670     n.negate();
1671     assert(n == 42);
1672     n.negate();
1673     assert(n == -42);
1674     n.absolute();
1675     assert(n == 42);
1676     n.absolute();
1677     assert(n == 42);
1678 
1679     // addition
1680 
1681     assert(a + b == b + a);     // commutative
1682     assert(a + 43.Z == b + a);
1683     assert(a - 43.Z == -(43.Z - a));
1684     assert(a + 0 == a);
1685     assert(a + 1 != a);
1686     assert(0 + a == a);
1687     assert(1 + a != a);
1688     assert(a + 0UL == a);
1689     assert(a + 1UL != a);
1690     assert(a + b == 42 + 43);
1691     assert(1 + a == 43);
1692     assert(a + (-1) == 41);
1693     assert(1UL + a == 43);
1694     assert(a + 1 == 1 + a);       // commutative
1695     assert(a + (-1) == (-1) + a); // commutative
1696     assert(1UL + a == a + 1UL);   // commutative
1697 
1698     // subtraction
1699 
1700     assert(a - 2 == 40);
1701     assert(2 - a == -40);
1702     assert(-2 - a == -44);
1703     assert(a - 2 == -(2 - a));     // commutative
1704     assert(a - (-2) == 44);
1705     assert(44UL - 42.Z == 2);
1706 
1707     // multiplication
1708 
1709     assert(a * 1UL == a);
1710     assert(a * 1 == a);
1711     assert(1 * a == a);
1712     assert(1UL * a == a);
1713     assert((-1) * a == -a);
1714     assert(a * 2 != a);
1715     assert(a * -2 == -(2*a));
1716     assert(a * b == b * a);
1717     assert(a * b == 42UL * 43UL);
1718 
1719     // division
1720 
1721     assert(27.Z / 3.Z == 9);
1722     assert(27.Z /   3  == 9);
1723 
1724     assert(27.Z / 10.Z  == 2);
1725     assert(27.Z /   10   == 2);
1726     assert(27.Z /   10UL == 2);
1727 
1728     assert(27.Z / -3   == -9);
1729     assert(27.Z /  3UL ==  9);
1730 
1731     assert(27.Z / -10   == -2);
1732 
1733     assert(28   /  3.Z ==  9);
1734     assert(28UL /  3.Z ==  9);
1735 
1736     assert(28   / -3.Z == -9);
1737     assert(28UL / -3.Z == -9);
1738 
1739     // modulo/remainder
1740 
1741     assert(27.Z % 3.Z == 0);
1742     assert(27.Z % 10.Z == 7);
1743 
1744     assert(27.Z % 3 == 0);
1745     assert(-27.Z % 3 == 0);
1746 
1747     assert(27.Z % 10 == 7);
1748     assert(27.Z % 10 == 7);
1749 
1750     assert(28   % 3.Z == 1);
1751     assert(28UL % 3.Z == 1);
1752 
1753     assert( 28.Z  % -3 == -1); // negative divisor gives negative remainder according to https://en.wikipedia.org/wiki/Remainder
1754     assert(-28.Z  %  3 == 1);  // dividend sign doesn't affect remainder
1755 
1756     //
1757     assert( 28.Z  % -3.Z == 1);  // TODO should be -1
1758     assert(-28.Z  %  3.Z == -1);  // TODO should be 1
1759     assert( 28  % -3.Z == 1);      // TODO should be -1
1760     assert(-28  %  3.Z == -1);     // TODO should be 1
1761 
1762     // modulo/remainder
1763 
1764     immutable one = 1.Z;
1765     const two = 2.Z;
1766     immutable three = 3.Z;
1767     const four = 4.Z;
1768     immutable five = 5.Z;
1769     const six = 6.Z;
1770     assert(six % one == 0);
1771     assert(six % two == 0);
1772     assert(six % three == 0);
1773     assert(six % four == 2);
1774     assert(six % five == 1);
1775     assert(six % six == 0);
1776 
1777     // subtraction
1778 
1779     assert(six - one == 5);
1780     assert(six - 1UL == 5);
1781     assert(six - 1 == 5);
1782     assert(1 - six == -5);
1783     assert(1L - six == -5);
1784     assert(1UL - six == -5);
1785 
1786     // exponentiation
1787 
1788     assert(0.Z^^0 == 1);
1789     assert(3.Z^^3 == 27);
1790     assert(3.Z^^3L == 27);
1791     assert(2.Z^^8 == 256);
1792     assert(2.Z^^8L == 256);
1793     assert(2.Z^^8UL == 256);
1794 
1795     assert(Z.pow(2UL, 8UL) == 256);
1796     assert(Z.pow(2UL, 8) == 256);
1797     assert(Z.pow(2UL, 8) == 256);
1798     assert(Z.pow(2, 8) == 256);
1799     assert(Z.pow(-2, 8) == 256);
1800     assert(Z.pow(-2, 7) == -128);
1801 
1802     // disallow power exponent to be an `MpZ`
1803     assert(!__traits(compiles, 2^^8.Z == 256));
1804     assert(!__traits(compiles, 2L^^8.Z == 256));
1805     assert(!__traits(compiles, 2UL^^8.Z == 256));
1806 
1807     // exponentiation plus modulus
1808 
1809     assert(2.Z.powm(8.Z, 8.Z) == 0.Z);
1810     assert(2.Z.powm(3.Z, 16.Z) == 8.Z);
1811     assert(3.Z.powm(3.Z, 16.Z) == 11.Z);
1812 
1813     assert(2.Z.powm(8, 8.Z) == 0.Z);
1814     assert(2.Z.powm(3, 16.Z) == 8.Z);
1815     assert(3.Z.powm(3, 16.Z) == 11.Z);
1816 
1817     // bitwise and, or and xor
1818 
1819     {
1820         foreach (immutable i; 0 .. 10)
1821         {
1822             foreach (immutable j; 0 .. 10)
1823             {
1824                 assert((i.Z & j.Z) == (i & j));
1825                 assert((i.Z | j.Z) == (i | j));
1826                 assert((i.Z ^ j.Z) == (i ^ j));
1827 
1828                 Z x = null;
1829 
1830                 x = i.Z;
1831                 x &= j.Z;
1832                 assert(x == (i & j));
1833 
1834                 x = i.Z;
1835                 x |= j.Z;
1836                 assert(x == (i | j));
1837 
1838                 x = i.Z;
1839                 x ^= j.Z;
1840                 assert(x == (i ^ j));
1841             }
1842         }
1843     }
1844 
1845     // swap
1846 
1847     auto x = 42.Z;
1848     auto y = 43.Z;
1849 
1850     assert(x == 42);
1851     assert(y == 43);
1852 
1853     x.swap(y);
1854 
1855     assert(y == 42);
1856     assert(x == 43);
1857 
1858     swap(x, y);
1859 
1860     assert(x == 42);
1861     assert(y == 43);
1862 
1863     assert(null.Z.fromString("42") == 42.Z);
1864     assert(null.Z.fromString("42") < 43.Z);
1865     assert(null.Z.fromString("42") > 41.Z);
1866     assert(null.Z.fromString("42") == 42);
1867     assert(null.Z.fromString("11", 2) == 3);
1868     assert(null.Z.fromString("7", 8) == 7);
1869     assert(null.Z.fromString("7") == 7);
1870     assert(null.Z.fromString("e", 16) == 14);
1871     assert(null.Z.fromString("f", 16) == 15);
1872     assert(null.Z.fromString("0xe") == 14);
1873     assert(null.Z.fromString("0xf") == 15);
1874     assert(null.Z.fromString("10", 16) == 16);
1875     assert(null.Z.fromString("10", 32) == 32);
1876 
1877     // odd and even
1878 
1879     assert(0.Z.isEven);
1880     assert(1.Z.isOdd);
1881     assert(2.Z.isEven);
1882     assert(3.Z.isOdd);
1883 
1884     assert((-1).Z.isOdd);
1885     assert((-2).Z.isEven);
1886     assert((-3).Z.isOdd);
1887 
1888     assert("300000000000000000000000000000000000000".Z.isEven);
1889     assert("300000000000000000000000000000000000001".Z.isOdd);
1890     assert("300000000000000000000000000000000000002".Z.isEven);
1891     assert("300000000000000000000000000000000000003".Z.isOdd);
1892 
1893     // negative and positive
1894 
1895     assert(0.Z.isPositive);
1896     assert(1.Z.isPositive);
1897     assert(2.Z.isPositive);
1898     assert(3.Z.isPositive);
1899 
1900     assert((-1).Z.isNegative);
1901     assert((-2).Z.isNegative);
1902     assert((-3).Z.isNegative);
1903 
1904     // sign function (sgn)
1905 
1906     assert(long.min.Z.sgn == -1);
1907     assert(int.min.Z.sgn  == -1);
1908     assert(-2.Z.sgn == -1);
1909     assert(-1.Z.sgn == -1);
1910     assert( 0.Z.sgn ==  0);
1911     assert( 1.Z.sgn ==  1);
1912     assert( 2.Z.sgn ==  1);
1913     assert(int.max.Z.sgn  == 1);
1914     assert(long.max.Z.sgn == 1);
1915 
1916     assert(!long.min.Z.isZero);
1917     assert(!int.min.Z.isZero);
1918     assert(!(-2).Z.isZero);
1919     assert(!(-1).Z.isZero);
1920     assert( 0.Z.isZero);
1921     assert(! 1.Z.isZero);
1922     assert(! 2.Z.isZero);
1923     assert(!int.max.Z.isZero);
1924     assert(!long.max.Z.isZero);
1925 
1926     // fits in type
1927 
1928     foreach (Integral; AliasSeq!(short, int, long,
1929                                  ushort, uint, ulong))
1930     {
1931         assert(Integral.min.Z.fitsIn!Integral);
1932         assert(Integral.max.Z.fitsIn!Integral);
1933     }
1934 
1935     // TODO
1936     // assert(short.min.Z.fitsIn!short);
1937     // assert(short.max.Z.fitsIn!short);
1938     // assert(ushort.min.Z.fitsIn!ushort);
1939     // assert(ushort.max.Z.fitsIn!ushort);
1940 
1941     // internal limb count
1942 
1943     assert(0.Z._limbCount == 0);
1944     assert(1.Z._limbCount == 1);
1945     assert(2.Z._limbCount == 1);
1946 
1947     assert(Z.pow(2UL, 32UL)._limbCount == 1);
1948 
1949     assert(Z.pow(2UL, 63UL)._limbCount == 1);
1950     assert(Z.pow(2UL, 63UL + 1)._limbCount == 2);
1951 
1952     assert(Z.pow(2UL, 127UL)._limbCount == 2);
1953     assert(Z.pow(2UL, 127UL + 1)._limbCount == 3);
1954 
1955     assert(Z.pow(2UL, 255UL)._limbCount == 4);
1956     assert(Z.pow(2UL, 255UL + 1)._limbCount == 5);
1957 }
1958 
1959 /// generators
1960 @safe @nogc unittest
1961 {
1962     assert(Z.mersennePrime(15) == 2^^15 - 1);
1963     assert(Z.mersennePrime(15UL) == 2^^15 - 1);
1964 }
1965 
1966 /// Phobos unittests
1967 @safe @nogc unittest
1968 {
1969     alias bigInt = mpz;
1970     alias BigInt = Z;     // Phobos naming convention
1971 
1972     {
1973         const BigInt a = "9588669891916142";
1974         const BigInt b = "7452469135154800";
1975         const c = a * b;
1976         assert(c == BigInt("71459266416693160362545788781600"));
1977         auto d = b * a;
1978         assert(d == BigInt("71459266416693160362545788781600"));
1979         assert(d == c);
1980         d = c * BigInt("794628672112");
1981         assert(d == BigInt("56783581982794522489042432639320434378739200"));
1982         auto e = c + d;
1983         assert(e == BigInt("56783581982865981755459125799682980167520800"));
1984         const f = d + c;
1985         assert(f == e);
1986         auto g = f - c;
1987         assert(g == d);
1988         g = f - d;
1989         assert(g == c);
1990         e = 12_345_678;
1991         g = c + e;
1992         immutable h = g / b;
1993         const i = g % b;
1994         assert(h == a);
1995         assert(i == e);
1996         BigInt j = "-0x9A56_57f4_7B83_AB78";
1997         j ^^= 11UL;
1998         j ^^= 2L;
1999         j ^^= 2;
2000     }
2001 
2002     {
2003         auto b = BigInt("1_000_000_000");
2004 
2005         b += 12_345;
2006         assert(b == 1_000_012_345);
2007         b += -12_345;
2008         assert(b == 1_000_000_000);
2009 
2010         b -= -12_345;
2011         assert(b == 1_000_012_345);
2012         b -= +12_345;
2013         assert(b == 1_000_000_000);
2014 
2015         b += 12_345;
2016         assert(b == 1_000_012_345);
2017 
2018         b /= 5UL;
2019         assert(b == 200_002_469);
2020     }
2021 
2022     {
2023         auto x = BigInt("123");
2024         const y = BigInt("321");
2025         x += y;
2026         assert(x == 444);
2027     }
2028 
2029     {
2030         const x = BigInt("123");
2031         const y = BigInt("456");
2032         const BigInt z = x * y;
2033         assert(z == 56_088);
2034     }
2035 
2036     {
2037         auto x = BigInt("123");
2038         x *= 300;
2039         assert(x == 36_900);
2040     }
2041 
2042     {
2043         auto x = BigInt("123");
2044         x *= -1;
2045         assert(x == -123);
2046     }
2047 
2048     {
2049         const x  = BigInt("1_000_000_500");
2050 
2051         immutable ulong ul  = 2_000_000UL;
2052         immutable uint ui   = 500_000;
2053         immutable ushort us = 30_000;
2054         immutable ubyte ub  = 50;
2055 
2056         immutable long  l = 1_000_000L;
2057         immutable int   i = 500_000;
2058         immutable short s = 30_000;
2059         immutable byte b  = 50;
2060 
2061         static assert(is(typeof(x % ul)  == ulong));
2062         static assert(is(typeof(x % ui)  == uint));
2063         static assert(is(typeof(x % us)  == ushort));
2064         static assert(is(typeof(x % ub)  == ubyte));
2065 
2066         static assert(is(typeof(x % l)  == long));
2067         static assert(is(typeof(x % i)  == int));
2068         static assert(is(typeof(x % s)  == short));
2069         static assert(is(typeof(x % b)  == byte));
2070 
2071         assert(x % ul == 500);
2072         assert(x % ui == 500);
2073         assert(x % us  == 10_500);
2074         assert(x % ub == 0);
2075 
2076         assert(x % l  == 500L);
2077         assert(x % i  == 500);
2078         assert(x % s  == 10_500);
2079         assert(x % b == 0);
2080     }
2081 
2082     {
2083         const x = BigInt("100");
2084         const BigInt y = 123 + x;
2085         assert(y == BigInt("223"));
2086 
2087         const BigInt z = 123 - x;
2088         assert(z == BigInt("23"));
2089 
2090         // Dividing a built-in integer type by BigInt always results in
2091         // something that fits in a built-in type, so the built-in type is
2092         // returned, not BigInt.
2093         static assert(is(typeof(1000 / x) == int));
2094         assert(1000 / x == 10);
2095     }
2096 
2097     {
2098         auto x = BigInt("1234");
2099         assert(+x  == BigInt(" 1234"));
2100         assert(-x  == BigInt("-1234"));
2101         ++x;
2102         assert(x == BigInt("1235"));
2103         --x;
2104         assert(x == BigInt("1234"));
2105     }
2106 
2107     {
2108         const x = BigInt("12345");
2109         const y = BigInt("12340");
2110         immutable int z = 12_345;
2111         immutable int w = 54_321;
2112         assert(x == x);
2113         assert(x != y);
2114         assert(x == y + 5);
2115         assert(x == z);
2116         assert(x != w);
2117     }
2118 
2119     {
2120         // non-zero values are regarded as `true`
2121         const x = BigInt("1");
2122         const y = BigInt("10");
2123         assert(x);
2124         assert(y);
2125 
2126         // zero value is regarded as `false`
2127         const z = BigInt("0");
2128         assert(!z);
2129     }
2130 
2131     {
2132         assert(cast(int)BigInt("0") == 0);
2133         assert(cast(ubyte)BigInt("0") == 0);
2134 
2135         assert(cast(ubyte)BigInt(255) == 255);
2136         assert(cast(ushort)BigInt(65_535) == 65_535);
2137         assert(cast(uint)BigInt(uint.max) == uint.max);
2138         assert(cast(ulong)BigInt(ulong.max) == ulong.max);
2139 
2140         assert(cast(byte)BigInt(-128) == -128);
2141         assert(cast(short)BigInt(-32_768) == -32_768);
2142         assert(cast(int)BigInt(int.min) == int.min);
2143         assert(cast(long)BigInt(long.min) == long.min);
2144 
2145         assert(cast(byte)BigInt(127) == 127);
2146         assert(cast(short)BigInt(32_767) == 32_767);
2147         assert(cast(int)BigInt(int.max) == int.max);
2148         assert(cast(long)BigInt(long.max) == long.max);
2149 
2150         // TODO:
2151         // import std.conv : to, ConvOverflowException;
2152         // import std.exception : assertThrown;
2153         // assertThrown!ConvOverflowException(BigInt("256").to!ubyte);
2154         // assertThrown!ConvOverflowException(BigInt("-1").to!ubyte);
2155     }
2156 
2157     {
2158         // TODO
2159         // segfaults because with double free
2160         // const(BigInt) x = BigInt("123");
2161         // BigInt y = cast()x;    // cast away const
2162         // assert(y == x);
2163     }
2164 
2165     {
2166         const x = BigInt("100");
2167         const y = BigInt("10");
2168         const int z = 50;
2169         const int w = 200;
2170         assert(y < x);
2171         assert(x > z);
2172         assert(z > y);
2173         assert(x < w);
2174     }
2175 
2176     {
2177         assert(BigInt("12345").toLong() == 12_345);
2178         assert(BigInt("-123450000000000000000000000000").toLong() == long.min);
2179         assert(BigInt("12345000000000000000000000000000").toLong() == long.max);
2180     }
2181 
2182     {
2183         assert(BigInt("12345").toInt() == 12_345);
2184         assert(BigInt("-123450000000000000000000000000").toInt() == int.min);
2185         assert(BigInt("12345000000000000000000000000000").toInt() == int.max);
2186     }
2187 }
2188 
2189 // Fermats Little Theorem
2190 @safe @nogc unittest
2191 {
2192     if (unittestLong) // compile but not run unless flagged for because running is slow
2193     {
2194         /*
2195           Fermats little theorem: a ^ p ≡ a (mod p) ∀ prime p check Fermats
2196           little theorem for a ≤ 100000 and all mersene primes M(p) : p ≤ 127
2197         */
2198         foreach (immutable ulong i; [2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127])
2199         {
2200             foreach (immutable ulong j; 2 .. 100_000)
2201             {
2202                 const p = Z.mersennePrime(i); // power
2203                 const a = Z(j); // base
2204                 const amp = a % p;
2205                 const b = a.powm(p, p); // result
2206                 assert(b == amp);
2207             }
2208         }
2209     }
2210 }
2211 
2212 // Euler's Sum of Powers Conjecture counter example
2213 pure @nogc unittest
2214 {
2215     /*
2216       a^5 + b^5 + c^5 + d^5 = e^5 Lander & Parkin, 1966 found the first counter
2217       example: 27^5 + 84^5 + 110^5 + 133^5 = 144^5. This test searches for this
2218       counter example by brute force for all positive a, b, c, d ≤ 144
2219     */
2220 
2221     enum LIMIT = 144 + 1;
2222     enum POWER = 5;
2223 
2224     if (unittestLong) // compile but not run unless flagged for because running is slow
2225     {
2226         bool found = false;
2227         Z r1 = 0;
2228     outermost:
2229         foreach (immutable ulong a; 1 .. LIMIT)
2230         {
2231             foreach (immutable ulong b; a .. LIMIT)
2232             {
2233                 foreach (immutable ulong c; b .. LIMIT)
2234                 {
2235                     foreach (immutable ulong d; c .. LIMIT)
2236                     {
2237                         r1 = ((Z(a) ^^ POWER) +
2238                               (Z(b) ^^ POWER) +
2239                               (Z(c) ^^ POWER) +
2240                               (Z(d) ^^ POWER));
2241                         Z rem = 0;
2242                         __gmpz_rootrem(r1._ptr, rem._ptr, r1._ptr, POWER); ++r1._ccc;
2243                         if (rem == 0UL)
2244                         {
2245                             debug printf("Counter Example Found: %lu^5 + %lu^5 + %lu^5 + %lu^5 = %lu^5\n",
2246                                          a, b, c, d, cast(ulong)r1);
2247                             found = true;
2248                             break outermost;
2249                         }
2250                     }
2251                 }
2252             }
2253         }
2254         assert(found);
2255     }
2256 }
2257 
2258 /// expression template types
2259 
2260 /// `MpZ`-`MpZ` adding expression.
2261 struct MpzAddExpr(T1, T2)
2262     if (isMpZExpr!T1 &&
2263         isMpZExpr!T2)
2264 {
2265     T1 e1;                      // first term
2266     T2 e2;                      // second term
2267     pragma(inline, true)
2268     MpZ eval() const @trusted
2269     {
2270         typeof(return) y = null;
2271         __gmpz_add(y._ptr, e1.eval()._ptr, e2.eval()._ptr); ++y._ccc;
2272         return y;
2273     }
2274 }
2275 version(unittest) static assert(isMpZExpr!(MpzAddExpr!(MpZ, MpZ)));
2276 
2277 /// Instantiator for `MpzAddExpr`.
2278 MpzAddExpr!(T1, T2) mpzAddExpr(T1, T2)(T1 t1, T2 t2)
2279     if (isMpZExpr!T1 &&
2280         isMpZExpr!T2)
2281 {
2282     // TODO don't eval certain type combinations of t1 and t2
2283     return MpzAddExpr!(T1, T2)(t1.eval(), t2.eval());
2284 }
2285 
2286 @safe @nogc unittest
2287 {
2288     assert(MpzAddExpr!(Z, Z)(3.Z, 4.Z).eval() == 3 + 4);
2289 
2290     const Z x = MpzAddExpr!(Z, Z)(3.Z, 4.Z);
2291     assert(x.mutatingCallCount == 1); // lower to `mpz_add`
2292     assert(x == 7);
2293 
2294     Z y = null;
2295     y = MpzAddExpr!(Z, Z)(3.Z, 4.Z);
2296     assert(y.mutatingCallCount == 2); // lowers to `mpz_init`, `mpz_add`
2297     assert(y == 7);
2298 
2299 }
2300 
2301 /// `MpZ`-`MpZ` subtraction expression.
2302 struct MpzSubExpr(T1, T2)
2303     if (isMpZExpr!T1 &&
2304         isMpZExpr!T2)
2305 {
2306     T1 e1;                      // first term
2307     T2 e2;                      // second term
2308     pragma(inline, true)
2309     MpZ eval() const @trusted
2310     {
2311         typeof(return) y = null;
2312         __gmpz_sub(y._ptr, e1.eval()._ptr, e2.eval()._ptr); ++y._ccc;
2313         return y;
2314     }
2315 }
2316 version(unittest) static assert(isMpZExpr!(MpzSubExpr!(MpZ, MpZ)));
2317 
2318 @safe @nogc unittest
2319 {
2320     assert(MpzSubExpr!(Z, Z)(3.Z, 4.Z).eval() == 3 - 4);
2321     const Z x = MpzSubExpr!(Z, Z)(3.Z, 4.Z);
2322     assert(x.mutatingCallCount == 1); // lowers to `mpz_sub`
2323     assert(x == -1);
2324 }
2325 
2326 /// `MpZ`-`MpZ` multiplication expression.
2327 struct MpzMulExpr(F1, F2)
2328     if (isMpZExpr!F1 &&
2329         isMpZExpr!F2)
2330 {
2331     F1 e1;                      // first factor
2332     F2 e2;                      // second factor
2333     pragma(inline, true)
2334     MpZ eval() const @trusted
2335     {
2336         typeof(return) y = null;
2337         __gmpz_mul(y._ptr, e1.eval()._ptr, e2.eval()._ptr); ++y._ccc;
2338         return y;
2339     }
2340 }
2341 version(unittest) static assert(isMpZExpr!(MpzMulExpr!(MpZ, MpZ)));
2342 
2343 @safe @nogc unittest
2344 {
2345     assert(MpzMulExpr!(Z, Z)(3.Z, 4.Z).eval() == 3 * 4);
2346 
2347     const Z x = MpzMulExpr!(Z, Z)(3.Z, 4.Z);
2348     assert(x == 12);
2349     assert(x.mutatingCallCount == 1); // lowers to `mpz_mul`
2350 }
2351 
2352 /// `MpZ`-`MpZ` division expression.
2353 struct MpzDivExpr(P, Q)
2354     if (isMpZExpr!P &&
2355         isMpZExpr!Q)
2356 {
2357     P e1;                       // divisor
2358     Q e2;                       // dividend
2359     pragma(inline, true)
2360     MpZ eval() const @trusted
2361     {
2362         typeof(return) y = null;
2363         __gmpz_tdiv_q(y._ptr, e1.eval()._ptr, e2.eval()._ptr); ++y._ccc;
2364         return y;
2365     }
2366 }
2367 version(unittest) static assert(isMpZExpr!(MpzDivExpr!(MpZ, MpZ)));
2368 
2369 @safe @nogc unittest
2370 {
2371     assert(MpzDivExpr!(Z, Z)(27.Z, 3.Z).eval() == 27 / 3);
2372     assert(MpzDivExpr!(Z, Z)(28.Z, 3.Z).eval() == 28 / 3);
2373     assert(MpzDivExpr!(Z, Z)(29.Z, 3.Z).eval() == 29 / 3);
2374     assert(MpzDivExpr!(Z, Z)(30.Z, 3.Z).eval() == 30 / 3);
2375 
2376     const Z x = MpzDivExpr!(Z, Z)(28.Z, 3.Z);
2377     assert(x == 9);
2378     assert(x.mutatingCallCount == 1); // lowers to `mpz_tdiv_q`
2379 }
2380 
2381 /// `MpZ`-`MpZ` modulus expression.
2382 struct MpzModExpr(P, Q)
2383     if (isMpZExpr!P &&
2384         isMpZExpr!Q)
2385 {
2386     P e1;                       // divisor
2387     Q e2;                       // dividend
2388     pragma(inline, true)
2389     MpZ eval() const @trusted
2390     {
2391         typeof(return) y = null;
2392         __gmpz_tdiv_r(y._ptr, e1.eval()._ptr, e2.eval()._ptr); ++y._ccc;
2393         return y;
2394     }
2395 }
2396 version(unittest) static assert(isMpZExpr!(MpzModExpr!(MpZ, MpZ)));
2397 
2398 @safe @nogc unittest
2399 {
2400     assert(MpzModExpr!(Z, Z)(27.Z, 3.Z).eval() == 27 % 3);
2401     assert(MpzModExpr!(Z, Z)(28.Z, 3.Z).eval() == 28 % 3);
2402     assert(MpzModExpr!(Z, Z)(29.Z, 3.Z).eval() == 29 % 3);
2403     assert(MpzModExpr!(Z, Z)(30.Z, 3.Z).eval() == 30 % 3);
2404 
2405     const Z x = MpzModExpr!(Z, Z)(29.Z, 3.Z);
2406     assert(x == 2);
2407     assert(x.mutatingCallCount == 1); // lowers to `mpz_tdiv_r`
2408 }
2409 
2410 /// `MpZ`-`ulong` power expression.
2411 struct MpzPowUExpr(P, Q)
2412     if (isMpZExpr!P &&
2413         isUnsigned!Q)
2414 {
2415     P e1;                       // base
2416     Q e2;                       // exponent
2417     pragma(inline, true)
2418     MpZ eval() const @trusted
2419     {
2420         typeof(return) y = null;
2421         __gmpz_pow_ui(y._ptr, e1.eval()._ptr, e2); ++y._ccc;
2422         return y;
2423     }
2424 }
2425 version(unittest) static assert(isMpZExpr!(MpzPowUExpr!(MpZ, ulong)));
2426 
2427 @safe @nogc unittest
2428 {
2429     assert(MpzPowUExpr!(Z, ulong)(3.Z, 3).eval() == 3^^3);
2430 }
2431 
2432 /// `MpZ`-`ulong`-`MpZ` power-modulo expression.
2433 struct MpzPowMUExpr(P, Q, M)
2434     if (isMpZExpr!P &&
2435         isUnsigned!Q &&
2436         isMpZExpr!M)
2437 {
2438     P base;                     // base
2439     Q exp;                      // exponent
2440     M mod;                      // modulo
2441     pragma(inline, true)
2442     MpZ eval() const @trusted
2443     {
2444         typeof(return) y = null;
2445         __gmpz_powm_ui(y._ptr, base.eval()._ptr, exp, mod._ptr); ++y._ccc;
2446         return y;
2447     }
2448 }
2449 version(unittest) static assert(isMpZExpr!(MpzPowMUExpr!(MpZ, ulong, MpZ)));
2450 
2451 @safe @nogc unittest
2452 {
2453     assert(MpzPowMUExpr!(Z, ulong, Z)(3.Z, 3, 20.Z).eval() == 3^^3 % 20);
2454 }
2455 
2456 /// `MpZ` negation expression.
2457 struct MpzNegExpr(A)
2458     if (isMpZExpr!A)
2459 {
2460     A e1;
2461     pragma(inline, true)
2462     MpZ eval() const @trusted
2463     {
2464         typeof(return) y = null;
2465         __gmpz_neg(y._ptr, e1.eval()._ptr); ++y._ccc;
2466         return y;
2467     }
2468 }
2469 version(unittest) static assert(isMpZExpr!(MpzNegExpr!(MpZ)));
2470 
2471 @safe @nogc unittest
2472 {
2473     assert(MpzNegExpr!(Z)(27.Z).eval() == -27);
2474 
2475     const Z x = MpzNegExpr!(Z)(27.Z);
2476     assert(x.mutatingCallCount == 1); // lowers to `mpz_neg`
2477     assert(x == -27);
2478 }
2479 
2480 /// Copied from `std.numeric` to prevent unnecessary Phobos deps.
2481 T _integralAbs(T)(T x)
2482     if (isIntegral!T)
2483 {
2484     return x >= 0 ? x : -x;
2485 }
2486 
2487 /** Faster than `std.traits`.
2488     See https://github.com/dlang/phobos/pull/5038
2489 */
2490 enum isArithmetic(T) = __traits(isArithmetic, T);
2491 enum isFloating(T) = __traits(isFloating, T);
2492 enum isFloatingPoint(T) = __traits(isFloating, T);
2493 enum isScalar(T) = __traits(isScalar, T);
2494 enum isScalarType(T) = __traits(isScalar, T);
2495 enum isIntegral(T) = __traits(isIntegral, T);
2496 enum isUnsigned(T) = __traits(isUnsigned, T);
2497 enum isSigned(T) = __traits(isArithmetic, T) && !__traits(isUnsigned, T);
2498 enum isStaticArray(T) = __traits(isStaticArray, T);
2499 enum isAssociativeArray(T) = __traits(isAssociativeArray, T);
2500 enum isAbstractClass(T) = __traits(isAbstractClass, T);
2501 enum isFinalClass(T) = __traits(isFinalClass, T);
2502 enum isPOD(T) = __traits(isPOD, T);
2503 enum isNested(T) = __traits(isNested, T);
2504 enum isVirtualFunction(alias fn) = __traits(isVirtualFunction, fn);
2505 enum isVirtualMethod(alias m) = __traits(isVirtualMethod, m);
2506 enum isAbstractFunction(alias fn) = __traits(isAbstractFunction, fn);
2507 enum isFinalFunction(alias fn) = __traits(isFinalFunction, fn);
2508 enum isOverrideFunction(alias fn) = __traits(isOverrideFunction, fn);
2509 enum isStaticFunction(alias fn) = __traits(isStaticFunction, fn);
2510 enum isOut(alias fn) = __traits(isOut, fn);
2511 enum isLazy(alias fn) = __traits(isLazy, fn);
2512 enum isTemplate(alias sym) = __traits(isTemplate, sym);
2513 enum hasMember(T, string member) = __traits(hasMember, T, member);
2514 enum IdentifierStringOfSymbol(alias sym) = __traits(identifier, sym);
2515 
2516 // TODO WARNING disabled because this cannot be wrapped in a template
2517 enum isRef(alias fn) = __traits(isRef, fn);
2518 
2519 /// http://forum.dlang.org/post/llwrbirvlqxawifyytqq@forum.dlang.org
2520 @safe pure nothrow @nogc unittest
2521 {
2522     struct S { int x, y; }
2523 
2524     static void f()(auto ref const S s)
2525     {
2526         static assert(__traits(isRef, s));
2527     }
2528 
2529     static void g()(auto ref const S s)
2530     {
2531         static assert(!__traits(isRef, s));
2532     }
2533 
2534     S s;
2535     static assert(!__traits(isRef, s));
2536 
2537     f(s);
2538 
2539     g(S.init);
2540 }
2541 
2542 /// as hash table key
2543 @safe unittest
2544 {
2545     // TODO disabled until non-copyable types work in AA's
2546     // string[Z] aa;
2547     // aa[123.Z] = "abc";
2548     // aa[456.Z] = "def";
2549     // assert(aa[123.Z] == "abc");
2550     // assert(aa[456.Z] == "def");
2551 }
2552 
2553 version(unittest)
2554 {
2555     import dbgio;
2556     alias Z = MpZ;
2557 }
2558 
2559 // C API
2560 extern(C)
2561 {
2562     alias __mp_limb_t = ulong;    // see `mp_limb_t` gmp.h. TODO detect when it is `uint` instead
2563     struct __mpz_struct
2564     {
2565         int _mp_alloc;		/* Number of *limbs* allocated and pointed to by
2566                                    the _mp_d field.  */
2567         int _mp_size;           /* abs(_mp_size) is the number of limbs the last
2568                                    field points to.  If _mp_size is negative
2569                                    this is a negative number.  */
2570         __mp_limb_t* _mp_d;       /* Pointer to the limbs. */
2571     }
2572     static assert(__mpz_struct.sizeof == 16); // fits in two 64-bit words
2573 
2574     alias mpz_srcptr = const(__mpz_struct)*;
2575     alias mpz_ptr = __mpz_struct*;
2576     alias mp_bitcnt_t = ulong;
2577 
2578     pure nothrow @nogc:
2579 
2580     void __gmpz_init (mpz_ptr);
2581     void __gmpz_init_set (mpz_ptr, mpz_srcptr);
2582 
2583     void __gmpz_init_set_d (mpz_ptr, double);
2584     void __gmpz_init_set_si (mpz_ptr, long);
2585     void __gmpz_init_set_ui (mpz_ptr, ulong);
2586     int __gmpz_init_set_str (mpz_ptr, const char*, int);
2587 
2588     void __gmpz_set (mpz_ptr rop, mpz_srcptr op);
2589     void __gmpz_set_ui (mpz_ptr rop, ulong op);
2590     void __gmpz_set_si (mpz_ptr rop, long op);
2591     void __gmpz_set_d (mpz_ptr rop, double op);
2592     int __gmpz_set_str (mpz_ptr, const char*, int);
2593 
2594     void __gmpz_clear (mpz_ptr);
2595 
2596     void __gmpz_abs (mpz_ptr, mpz_srcptr);
2597     void __gmpz_neg (mpz_ptr, mpz_srcptr);
2598 
2599     void __gmpz_add (mpz_ptr, mpz_srcptr, mpz_srcptr);
2600     void __gmpz_add_ui (mpz_ptr, mpz_srcptr, ulong);
2601     void __gmpz_addmul (mpz_ptr, mpz_srcptr, mpz_srcptr);
2602     void __gmpz_addmul_ui (mpz_ptr, mpz_srcptr, ulong);
2603 
2604     void __gmpz_sub (mpz_ptr, mpz_srcptr, mpz_srcptr);
2605     void __gmpz_sub_ui (mpz_ptr, mpz_srcptr, ulong);
2606     void __gmpz_ui_sub (mpz_ptr, ulong, mpz_srcptr);
2607 
2608     void __gmpz_mul (mpz_ptr, mpz_srcptr, mpz_srcptr);
2609     void __gmpz_mul_2exp (mpz_ptr, mpz_srcptr, mp_bitcnt_t);
2610     void __gmpz_mul_si (mpz_ptr, mpz_srcptr, long);
2611     void __gmpz_mul_ui (mpz_ptr, mpz_srcptr, ulong);
2612 
2613     void __gmpz_cdiv_q_ui (mpz_ptr, mpz_srcptr, ulong);
2614     ulong __gmpz_cdiv_r_ui (mpz_ptr, mpz_srcptr, ulong);
2615     void __gmpz_cdiv_qr_ui (mpz_ptr, mpz_ptr, mpz_srcptr, ulong);
2616     void __gmpz_cdiv_ui (mpz_srcptr, ulong);
2617 
2618     void __gmpz_fdiv_q_ui (mpz_ptr, mpz_srcptr, ulong);
2619     ulong __gmpz_fdiv_r_ui (mpz_ptr, mpz_srcptr, ulong);
2620     void __gmpz_fdiv_qr_ui (mpz_ptr, mpz_ptr, mpz_srcptr, ulong);
2621     void __gmpz_fdiv_ui (mpz_srcptr, ulong);
2622 
2623     void __gmpz_tdiv_q (mpz_ptr, mpz_srcptr, mpz_srcptr);
2624     void __gmpz_tdiv_r (mpz_ptr, mpz_srcptr, mpz_srcptr);
2625     void __gmpz_tdiv_q_ui (mpz_ptr, mpz_srcptr, ulong);
2626     ulong __gmpz_tdiv_r_ui (mpz_ptr, mpz_srcptr, ulong);
2627     void __gmpz_tdiv_qr_ui (mpz_ptr, mpz_ptr, mpz_srcptr, ulong);
2628     void __gmpz_tdiv_ui (mpz_srcptr, ulong);
2629 
2630     void __gmpz_mod (mpz_ptr, mpz_srcptr, mpz_srcptr);
2631 
2632     void __gmpz_pow_ui (mpz_ptr, mpz_srcptr, ulong);
2633     void __gmpz_ui_pow_ui (mpz_ptr, ulong, ulong);
2634 
2635     void __gmpz_swap (mpz_ptr, mpz_ptr);
2636 
2637     int __gmpz_cmp (mpz_srcptr, mpz_srcptr);
2638     int __gmpz_cmp_d (mpz_srcptr, double);
2639     int __gmpz_cmp_si (mpz_srcptr, long);
2640     int __gmpz_cmp_ui (mpz_srcptr, ulong);
2641 
2642     int __gmpz_cmpabs (mpz_srcptr, mpz_srcptr);
2643     int __gmpz_cmpabs_d (mpz_srcptr, double);
2644     int __gmpz_cmpabs_ui (mpz_srcptr, ulong);
2645 
2646     void __gmpz_nextprime (mpz_ptr, mpz_srcptr);
2647 
2648     void __gmpz_gcd (mpz_ptr, mpz_srcptr, mpz_srcptr);
2649     ulong __gmpz_gcd_ui (mpz_ptr, mpz_srcptr, ulong);
2650 
2651     void __gmpz_lcm (mpz_ptr, mpz_srcptr, mpz_srcptr);
2652     void __gmpz_lcm_ui (mpz_ptr, mpz_srcptr, ulong);
2653 
2654     char *__gmpz_get_str (char*, int, mpz_srcptr);
2655     size_t __gmpz_sizeinbase (mpz_srcptr, int);
2656 
2657     void __gmpz_powm (mpz_ptr, mpz_srcptr, mpz_srcptr, mpz_srcptr);
2658     void __gmpz_powm_ui (mpz_ptr, mpz_srcptr, ulong, mpz_srcptr);
2659 
2660     ulong __gmpz_get_ui (mpz_srcptr);
2661     long __gmpz_get_si (mpz_srcptr);
2662     double __gmpz_get_d (mpz_srcptr);
2663 
2664     int __gmpz_fits_ulong_p(mpz_srcptr);
2665     int __gmpz_fits_slong_p(mpz_srcptr);
2666     int __gmpz_fits_uint_p(mpz_srcptr);
2667     int __gmpz_fits_sint_p(mpz_srcptr);
2668     int __gmpz_fits_ushort_p(mpz_srcptr);
2669     int __gmpz_fits_sshort_p(mpz_srcptr);
2670 
2671     void __gmpz_and (mpz_ptr, mpz_srcptr, mpz_srcptr);
2672     void __gmpz_ior (mpz_ptr, mpz_srcptr, mpz_srcptr);
2673     void __gmpz_xor (mpz_ptr, mpz_srcptr, mpz_srcptr);
2674     void __gmpz_com (mpz_ptr, mpz_srcptr);
2675 
2676     // TODO wrap:
2677     void __gmpz_root (mpz_ptr, mpz_srcptr, ulong);
2678     void __gmpz_rootrem (mpz_ptr, mpz_ptr, mpz_srcptr, ulong);
2679 
2680     void __gmpz_sqrt (mpz_ptr, mpz_srcptr);
2681     void __gmpz_sqrtrem (mpz_ptr, mpz_ptr, mpz_srcptr);
2682 
2683     void __gmpz_perfect_power_p (mpz_ptr, mpz_srcptr);
2684 }
2685 
2686 // link with C library GNU MP
2687 pragma(lib, "gmp");

Page generated by adrdox