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");