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