1 /// Multiple precision rational numbers (Q).
2 module gmp.q;
3 
4 import gmp.traits;
5 import gmp.z;
6 
7 import std.algorithm.mutation : move, moveEmplace;
8 
9 /** Arbitrary (multi) precision rational number (Q).
10     Wrapper for GNU MP (GMP)'s type `mpq_t` and functions `__gmpq_.*`.
11  */
12 struct MpQ
13 {
14     pure nothrow:
15 
16     /// Convert to `string` in base `base`.
17     string toString(uint base = defaultBase,
18                     bool upperCaseDigits = false) const @trusted
19     {
20         assert((base >= -2 && base <= -36) ||
21                (base >= 2 && base <= 62));
22 
23         assert(false, "TODO");
24     }
25 
26 pragma(inline, true):
27 
28     // TODO toRCString wrapped in UniqueRange
29 
30     /// Returns: A unique hash of the `MpQ` value suitable for use in a hash table.
31     size_t toHash() const
32     {
33         assert(false, "TODO");
34     }
35 
36     @nogc:
37 
38     /** No default construction for now, because `mpq_init` initialize
39         `__mpq_struct`-fields to non-zero values.
40 
41         TODO Allow default construction by delaying call to initialize().
42     */
43     @disable this();
44 
45     /// No copy construction.
46     @disable this(this);
47 
48     /// Construct empty (undefined) from explicit `null`.
49     this(typeof(null)) @safe
50     {
51         initialize();
52     }
53 
54     /** Construct from floating-point `value`.
55      */
56     this(P)(P value) @safe
57         if (isFloating!P)
58     {
59         initialize();
60         this = value;           // reuse opAssign
61     }
62 
63     /** Construct from `pValue` / `qValue`.
64 
65         Note that `qValue` must be explicitly given, to prevent accidental
66         storage of integers as rations with denominator being 1.
67      */
68     this(P, Q)(P pValue, Q qValue,
69                bool canonicalizeFlag = false) @trusted
70         if (isIntegral!P &&
71             isIntegral!Q)
72     {
73         initialize();
74 
75         version(ccc) ++_ccc;
76 
77         static if (isSigned!Q)
78         {
79             assert(qValue >= 1, "Negative denominator");
80         }
81 
82         // dln("qValue:", qValue);
83 
84         static      if (isUnsigned!P)
85         {
86             // dln("unsigned pValue:", pValue);
87             __gmpq_set_ui(_ptr, pValue, qValue);
88         }
89         else                    // signed integral
90         {
91             // dln("signed pValue:", pValue);
92             __gmpq_set_si(_ptr, pValue, qValue);
93         }
94 
95         if (canonicalizeFlag) { canonicalize(); }
96     }
97 
98     /** Construct from floating-point `value`.
99      */
100     ref MpQ opAssign(P)(P value) @trusted return scope
101         if (isFloating!P)
102     {
103         version(ccc) ++_ccc;
104         __gmpq_set_d(_ptr, value);
105         return this;
106     }
107 
108     /** Assign from integer `value`. */
109     ref MpQ opAssign(P)(P value) @trusted return scope
110         if (isIntegral!P)
111     {
112         version(ccc) ++_ccc;
113 
114         static      if (isUnsigned!P)
115             __gmpq_set_ui(_ptr, value, 1);
116         else                    // signed integral
117             __gmpq_set_si(_ptr, value, 1);
118 
119         return this;
120     }
121 
122     /** Canonicalize `this`. */
123     void canonicalize() @trusted
124     {
125         __gmpq_canonicalize(_ptr); version(ccc) ++_ccc;
126     }
127 
128     /// Destruct `this`.
129     ~this() @trusted
130     {
131         assert(_ptr, "Pointer is null");
132         __gmpq_clear(_ptr); version(ccc) ++_ccc;
133     }
134 
135     /// Returns: `true` iff `this` equals `rhs`.
136     bool opEquals()(auto ref const MpQ rhs) const @trusted
137     {
138         if (_ptr == rhs._ptr)   // fast equality
139         {
140             return true;        // fast bailout
141         }
142         return __gmpq_equal(_ptr, rhs._ptr) != 0;
143     }
144     /// ditto
145     int opEquals(T)(T rhs) const @safe
146         if (isIntegral!T)
147     {
148         if (rhs == 0)
149         {
150             return numerator.isZero; // optimization
151         }
152         return numerator == rhs && denominator == 1;
153     }
154 
155     /// Compare `this` to `rhs`.
156     int opCmp()(auto ref const MpQ rhs) const @trusted
157     {
158         if (rhs.numerator == 0)
159         {
160             return sgn;         // optimization
161         }
162         return __gmpq_cmp(_ptr, rhs._ptr);
163     }
164     /// Compare `this` to `rhs`.
165     int opCmp()(auto ref const MpZ rhs) const @trusted
166     {
167         if (rhs == 0)
168         {
169             return sgn;         // optimization
170         }
171         return __gmpq_cmp_z(_ptr,
172                             cast(const(__mpz_struct)*)&rhs); // TODO wrap cast?
173     }
174     /// ditto
175     int opCmp(T)(T rhs) const @trusted
176         if (isIntegral!T)
177     {
178         if (rhs == 0)
179         {
180             return sgn;         // optimization
181         }
182         static if (isUnsigned!T)
183         {
184             return __gmpq_cmp_ui(_ptr, rhs, 1UL);
185         }
186         else                    // isSigned integral
187         {
188             return __gmpq_cmp_si(_ptr, rhs, 1UL);
189         }
190     }
191 
192     /// Returns: numerator reference of `this`.
193     @property ref inout(MpZ) numerator() @trusted inout return scope
194     {
195         return *(cast(inout(MpZ)*)_num_ptr);
196     }
197 
198     /// Returns: denominator reference of `this`.
199     @property ref inout(MpZ) denominator() @trusted inout return scope
200     {
201         return *(cast(inout(MpZ)*)_den_ptr);
202     }
203 
204     /// Returns: the integer part of `this`, with any remainder truncated.
205     @property MpZ integerPart() @safe
206     {
207         return numerator / denominator;
208     }
209 
210     /// Returns: the fractional part of `this`.
211     // TODO activate when sub(MpQ, MpZ) has been added
212     // @property MpQ fractionPart()
213     // {
214     //     return this - integerPart;
215     // }
216 
217     /// Cast to arithmetic type `T`.
218     T opCast(T)() const @trusted /*TODO scope*/
219         if (isFloating!T)
220     {
221         return cast(T)__gmpq_get_d(_ptr);
222     }
223 
224     /** Invert `this` in-place.
225         Returns: `void` to make it obvious that `this` is mutated.
226     */
227     void invert() @trusted
228     {
229         import std.algorithm.mutation : swap;
230         const bool negative = numerator < 0;
231         if (negative)
232         {
233             numerator.absolute();         // fast inline
234             swap(numerator, denominator); // fast inline
235             numerator.negate();           // fast inline
236         }
237         else
238         {
239             swap(numerator, denominator); // fast inline
240         }
241     }
242 
243     /** Returns: sign as either
244         - -1 (`this` < 0),
245         -  0 (`this` == 0), or
246         - +1 (`this` > 0).
247     */
248     @property int sgn() const @safe
249     {
250         assert(denominator >= 1);
251         return numerator.sgn;   // sign always stored in numerator so reuse fast
252     }
253 
254     /** Make `this` the absolute value of itself in-place.
255         Returns: `void` to make it obvious that `this` is mutated.
256     */
257     void absolute() @safe
258     {
259         numerator.absolute();
260     }
261 
262     MpQ opBinary(string s)(auto ref const MpQ rhs) const @trusted // direct value
263         if ((s == "+" || s == "-" ||
264              s == "*" || s == "/"))
265     {
266         static if (!__traits(isRef, rhs)) // r-value `rhs`
267         {
268             MpQ* mut_rhs = (cast(MpQ*)(&rhs)); // @trusted because `MpQ` has no aliased indirections
269             static      if (s == "+")
270             {
271                 __gmpq_add(mut_rhs._ptr, _ptr, rhs._ptr); version(ccc) ++mut_rhs._ccc;
272             }
273             else static if (s == "-")
274             {
275                 __gmpq_sub(mut_rhs._ptr, _ptr, rhs._ptr); version(ccc) ++mut_rhs._ccc;
276             }
277             else static if (s == "*")
278             {
279                 __gmpq_mul(mut_rhs._ptr, _ptr, rhs._ptr); version(ccc) ++mut_rhs._ccc;
280             }
281             else static if (s == "/")
282             {
283                 assert(rhs != 0, "Divison by zero");
284                 __gmpq_div(mut_rhs._ptr, _ptr, rhs._ptr); version(ccc) ++mut_rhs._ccc;
285             }
286             else
287             {
288                 static assert(false);
289             }
290             return move(*mut_rhs); // TODO shouldn't have to call `move` here
291         }
292         else
293         {
294             typeof(return) y = null;
295             static      if (s == "+")
296             {
297                 __gmpq_add(y._ptr, _ptr, rhs._ptr); version(ccc) ++y._ccc;
298             }
299             else static if (s == "-")
300             {
301                 __gmpq_sub(y._ptr, _ptr, rhs._ptr); version(ccc) ++y._ccc;
302             }
303             else static if (s == "*")
304             {
305                 __gmpq_mul(y._ptr, _ptr, rhs._ptr); version(ccc) ++y._ccc;
306             }
307             else static if (s == "/")
308             {
309                 assert(rhs != 0, "Divison by zero");
310                 __gmpq_div(y._ptr, _ptr, rhs._ptr); version(ccc) ++y._ccc;
311             }
312             else
313             {
314                 static assert(false);
315             }
316             return y;
317         }
318     }
319 
320     // /// Divided integral with `this`.
321     // Unqual!Lhs opBinaryRight(string s, Lhs)(Lhs lhs) const @trusted
322     //     if ((s == "/") &&
323     //         isIntegral!Lhs)
324     // {
325     //     if (lhs == 1)
326     //     {
327     //         return inv(this);
328     //     }
329     //     else
330     //     {
331     //         MpQ y = null; // TODO avoid if !__traits(isRef, this)
332     //         version(ccc) ++y._ccc;
333     //         assert(this != 0, "Divison by zero");
334     //         denominator *= lhs;
335     //         __gmpq_div(y._ptr, MpZ(lhs)._ptr, _ptr);
336     //         return y;
337     //     }
338     // }
339 
340 private:
341 
342     /// Default conversion base.
343     enum defaultBase = 10;
344 
345     /** Initialize internal struct. */
346     private void initialize() @trusted // cannot be called `init` as that will override builtin type property
347     {
348         __gmpq_init(_ptr); version(ccc) ++_ccc;
349     }
350 
351     /// Returns: pointer to internal rational C struct.
352     inout(__mpq_struct)* _ptr() inout return @system
353     {
354         return &_q;
355     }
356 
357     /// Returns: pointer to internal numerator C struct.
358     inout(__mpz_struct)* _num_ptr() inout return @system
359     {
360         return cast(typeof(return))&_q._mp_num;
361     }
362 
363     /// Returns: pointer to internal denominator C struct.
364     inout(__mpz_struct)* _den_ptr() inout return @system
365     {
366         return cast(typeof(return))&_q._mp_den;
367     }
368 
369     __mpq_struct _q;            // internal libgmp C struct
370 
371     version(ccc)
372     {
373 
374         /** Number of calls made to `__gmpq`--functions that construct or changes
375             this value. Used to verify correct lowering and evaluation of template
376             expressions.
377 
378             For instance the `x` in `x = y + z` should be assigned only once inside
379             a call to `mpq_add`.
380         */
381         @property size_t mutatingCallCount() const @safe { return _ccc; }
382         size_t _ccc;  // C mutation call count. number of calls to C GMP function calls that mutate this object
383     }
384 }
385 
386 pure nothrow pragma(inline, true):
387 
388 /// Swap contents of `x` with contents of `y`.
389 void swap()(ref MpQ x,
390             ref MpQ y)
391 {
392     import std.algorithm.mutation : swap;
393     swap(x, y); // x.swap(y);
394 }
395 
396 /// Returns: absolute value of `x`.
397 MpQ abs()(auto ref const MpQ x) @trusted
398 {
399     static if (__traits(isRef, x)) // l-value `x`
400     {
401         MpQ y = null;
402         __gmpq_abs(y._ptr, x._ptr);
403         return y;
404     }
405     else                        // r-value `x`
406     {
407         MpQ* mut_x = (cast(MpQ*)(&x)); // @trusted because `MpQ` has no aliased indirections
408         mut_x.absolute();
409         return move(*mut_x);    // TODO shouldn't have to call `move` here
410     }
411 }
412 
413 /// Returns: inverse of `x`.
414 MpQ inverse()(auto ref const MpQ x) @trusted
415 {
416     static if (__traits(isRef, x)) // l-value `x`
417     {
418         MpQ y = null;
419         __gmpq_inv(y._ptr, x._ptr);
420         return y;
421     }
422     else                        // r-value `x`
423     {
424         MpQ* mut_x = (cast(MpQ*)(&x)); // @trusted because `MpQ` has no aliased indirections
425         mut_x.invert();
426         return move(*mut_x);    // TODO shouldn't have to call `move` here
427     }
428 }
429 alias inv = inverse;
430 
431 /// construction and assignment
432 @safe @nogc unittest
433 {
434     scope Q x = null;
435     assert(x.numerator == 0);
436     assert(x.denominator == 1);
437 
438     scope const Q y = Q(11, 13UL);
439     assert(y.numerator == 11);
440     assert(y.denominator == 13);
441 
442     scope immutable Q z = Q(7UL, 13UL);
443     assert(z.numerator == 7);
444     assert(z.denominator == 13);
445 
446     scope Q w = 0.25;           // construct from `double`
447     assert(w.numerator == 1);
448     assert(w.denominator == 4);
449 
450     w = 0.125f;                 // assign from `float`
451     assert(w.numerator == 1);
452     assert(w.denominator == 8);
453 
454     w = 2;                      // assign from `int`
455     assert(w.numerator == 2);
456     assert(w.denominator == 1);
457 
458     w = 3;                      // assign from `int`
459     assert(w.numerator == 3);
460     assert(w.denominator == 1);
461 }
462 
463 /// canonicalization
464 @safe @nogc unittest
465 {
466     Q x = Q(2, 4);
467     assert(x.numerator == 2);
468     assert(x.denominator == 4);
469     x.canonicalize();
470     assert(x.numerator == 1);
471     assert(x.denominator == 2);
472 }
473 
474 /// negative numerator canonicalization
475 @safe @nogc unittest
476 {
477     Q x = Q(-2, 4);
478     assert(x.numerator == -2);
479     assert(x.denominator == 4);
480     x.canonicalize();
481     assert(x.numerator == -1);
482     assert(x.denominator == 2);
483 }
484 
485 /// swap
486 @safe @nogc unittest
487 {
488     Q x = Q(1, 2);
489     Q y = Q(1, 3);
490     swap(x, y);
491     assert(x == Q(1, 3));
492     assert(y == Q(1, 2));
493 }
494 
495 /// invert
496 @safe unittest
497 {
498     Q x = Q(1, 2);
499     assert(x.numerator == 1);
500     assert(x.denominator == 2);
501 
502     x.invert();
503     assert(x.numerator == 2);
504     assert(x.denominator == 1);
505 
506     Q y = Q(-1, 2);
507 
508     // TODO:
509     // assert(x.numerator == -1);
510     // assert(x.denominator == 2);
511 
512     // x.invert();
513     // assert(x.numerator == -2);
514     // assert(x.denominator == 1);
515 }
516 
517 /// inversion
518 @safe unittest
519 {
520     const Q q = Q(-2, 3);
521     assert(inverse(q) == Q(-3, 2));
522 
523     assert(inverse(Q(2, 3)) == Q(3, 2));
524     assert(inverse(Q(1, 10)) == 10);
525     assert(inverse(Q(10, 1)) == Q(1, 10));
526 }
527 
528 /// absolute value
529 @safe unittest
530 {
531     const Q q = Q(-2, 3);
532     assert(abs(q) == Q(2, 3));
533     assert(abs(Q(-2, 3)) == Q(2, 3));
534 }
535 
536 /// integer and fractional part
537 @safe unittest
538 {
539     Q x = Q(5, 2);
540 
541     assert(x.integerPart == 2);
542     // TODO assert(x.fractionalPart == Q(1, 2));
543 
544     x = Q(7, 2);
545     assert(x.integerPart == 3);
546     // TODO assert(x.fractionalPart == Q(1, 2));
547 
548     x = Q(10, 2);
549     assert(x.integerPart == 5);
550     // TODO assert(x.fractionalPart == 0);
551 
552     x = Q(11, 3);
553     assert(x.integerPart == 3);
554     // TODO assert(x.fractionalPart == Q(2, 3));
555 
556     x = Q(12, 2);
557     assert(x.integerPart == 6);
558     // TODO assert(x.fractionalPart == 0);
559 }
560 
561 /// casting
562 @safe @nogc unittest
563 {
564     assert(cast(double)Q(1, 2) == 0.5f);
565     assert(cast(double)Q(1, 2) == 0.5);
566     assert(cast(double)Q(2, 4) == 0.5);
567     assert(cast(double)Q(1, 8) == 1.0/8);
568 }
569 
570 /// equality
571 @safe unittest
572 {
573     assert(Q(1, 1) == 1);
574     assert(Q(2, 1) == 2);
575     assert(Q(1, 1) == Q(1, 1));
576     assert(Q(1, 1) != Q(1, 2));
577     const x = Q(1, 3);
578     assert(x == x);             // same
579 }
580 
581 /// sign
582 @safe unittest
583 {
584     assert(Q(-1, 3).sgn == -1);
585     assert(Q( 0, 3).sgn ==  0);
586     assert(Q( 1, 3).sgn ==  1);
587 }
588 
589 /// comparison
590 @safe unittest
591 {
592     assert(Q( 1, 3) < Q(1, 2));
593     assert(Q( 1, 2) > Q(1, 3));
594     assert(Q( 1, 2) > Q(0, 1));
595     assert(Q( 0, 1) == Q(0, 1));
596     assert(Q( 0, 2) == Q(0, 1));
597     assert(Q(-1, 2) < Q(0, 1));
598 
599     assert(Q( 1, 3) < 1);
600     assert(Q( 1, 3) > 0);
601     assert(Q(-1, 3) < 0);
602 
603     assert(Q( 1, 3) < 1UL);
604     assert(Q( 1, 3) > 0UL);
605     assert(Q(-1, 3) < 0UL);
606 
607     assert(Q( 1, 3) < 1.Z);
608     assert(Q( 1, 3) > 0.Z);
609     assert(Q(-1, 3) < 0.Z);
610 }
611 
612 /// addition
613 @safe unittest
614 {
615     assert(Q(1, 2) + Q(1, 2) == Q(1, 1));
616     assert(Q(1, 3) + Q(1, 3) == Q(2, 3));
617     assert(Q(1, 2) + Q(1, 3) == Q(5, 6));
618 }
619 
620 /// subtraction
621 @safe unittest
622 {
623     assert(Q(1, 2) - Q(1, 2) == Q( 0, 1));
624     assert(Q(1, 2) - Q(1, 3) == Q (1, 6));
625     assert(Q(1, 3) - Q(1, 2) == Q(-1, 6));
626 }
627 
628 /// multiplication
629 @safe unittest
630 {
631     assert(Q(1, 2) * Q(1, 2) == Q(1, 4));
632     assert(Q(2, 3) * Q(2, 3) == Q(4, 9));
633     assert(Q(1, 2) * Q(1, 3) == Q(1, 6));
634 }
635 
636 /// division
637 @safe unittest
638 {
639     assert(Q(2, 3) / Q(2, 3) == Q(1, 1));
640     assert(Q(2, 3) / Q(2, 3) == 1);
641     assert(Q(2, 3) / Q(3, 2) == Q(4, 9));
642     assert(Q(3, 2) / Q(2, 3) == Q(9, 4));
643     // TODO assert(1 / Q(2, 3) == Q(3, 2));
644 }
645 
646 version(unittest)
647 {
648     // import dbgio : dln;
649     alias Z = MpZ;
650     alias Q = MpQ;
651     version = ccc;              // do C mutation call count
652     static assert(!isMpZExpr!int);
653     import std.meta : AliasSeq;
654 }
655 
656 // C API
657 extern(C) pragma(inline, false)
658 {
659     struct __mpq_struct
660     {
661         __mpz_struct _mp_num;
662         __mpz_struct _mp_den;
663     }
664     static assert(__mpq_struct.sizeof == 32); // fits in four 64-bit words
665 
666     alias mpq_srcptr = const(__mpq_struct)*;
667     alias mpq_ptr = __mpq_struct*;
668     alias mp_bitcnt_t = ulong;
669 
670     pure nothrow @nogc:
671 
672     void __gmpq_init (mpq_ptr);
673     void __gmpq_clear (mpq_ptr);
674 
675     void __gmpq_set (mpq_ptr, mpq_srcptr);
676     void __gmpq_set_z (mpq_ptr, mpz_srcptr);
677 
678     void __gmpq_set_ui (mpq_ptr, ulong, ulong);
679     void __gmpq_set_si (mpq_ptr, long, ulong);
680     void __gmpq_set_d (mpq_ptr, double);
681 
682     double __gmpq_get_d (mpq_srcptr);
683 
684     void __gmpq_canonicalize (mpq_ptr);
685 
686     int __gmpq_equal (mpq_srcptr, mpq_srcptr);
687     int __gmpq_cmp (mpq_srcptr, mpq_srcptr);
688     int __gmpq_cmp_z (mpq_srcptr, mpz_srcptr);
689 
690     int __gmpq_cmp_ui (mpq_srcptr, ulong, ulong);
691     int __gmpq_cmp_si (mpq_srcptr, long, ulong);
692 
693     void __gmpq_add (mpq_ptr, mpq_srcptr, mpq_srcptr);
694     void __gmpq_sub (mpq_ptr, mpq_srcptr, mpq_srcptr);
695     void __gmpq_mul (mpq_ptr, mpq_srcptr, mpq_srcptr);
696     void __gmpq_div (mpq_ptr, mpq_srcptr, mpq_srcptr);
697 
698     void __gmpq_abs (mpq_ptr, mpq_srcptr);
699     void __gmpq_inv (mpq_ptr, mpq_srcptr);
700     }
701 
702 pragma(lib, "gmp");