1 /** This module contains an implementation of rational numbers that is templated
2  * on the underlying integer type.  It can be used with either builtin fixed
3  * width integers or arbitrary precision integers.  All relevant operators are
4  * overloaded for both rational-rational and rational-integer operations.
5  *
6  * Synopsis:
7  * ---
8  * // Compute pi using the generalized continued fraction approximation.
9  * import std.bigint, std.rational, std.stdio;
10  *
11  * enum maxTerm = 30;
12  *
13  * Rational!(BigInt) getTerm(int termNumber)
14  * {
15  *	 auto addFactor = 2 * termNumber - 1;
16  *
17  *	 if (termNumber == maxTerm)
18  *	 {
19  *		 return rational(BigInt(addFactor));
20  *	 }
21  *
22  *	 auto termNumberSquared = BigInt(termNumber * termNumber);
23  *	 auto continued = termNumberSquared / getTerm(termNumber + 1);
24  *
25  *	 continued += addFactor;
26  *	 return continued;
27  * }
28  *
29  * void main()
30  * {
31  *	 auto pi = rational(BigInt(4)) / getTerm(1);
32  *
33  *	 // Display the result in rational form.
34  *	 writeln(pi);
35  *
36  *	 // Display the decimal equivalent, which is accurate to 18 decimal places.
37  *	 writefln("%.18f", cast(real) pi);
38  * }
39  * ---
40  *
41  *
42  * Author:  David Simcha
43  * Copyright:  Copyright (c) 2009-2011, David Simcha.
44  * License:	$(LINK2 http://www.boost.org/LICENSE_1_0.txt, Boost License 1.0)
45  */
46 module nxt.rational;
47 
48 import std.conv : to;
49 import std.math : abs;
50 
51 /** Checks whether $(D T) is structurally an integer, i.e. whether it supports
52  * all of the operations an integer type should support.  Does not check the
53  * nominal type of $(D T).  In particular, for a mutable type $(D T) the
54  * following must compile:
55  *
56  * ---
57  * T n;
58  * n = 2;
59  * n <<= 1;
60  * n >>= 1;
61  * n += n;
62  * n += 2;
63  * n *= n;
64  * n *= 2;
65  * n /= n;
66  * n /= 2;
67  * n -= n;
68  * n -= 2;
69  * n %= 2;
70  * n %= n;
71  * bool foo = n < 2;
72  * bool bar = n == 2;
73  * bool goo = n < n + 1;
74  * bool tar = n == n;
75  * ---
76  *
77  * while for a non-mutable type, the above must compile for its unqualified,
78  * mutable variant.
79  *
80  * All built-in D integers and character types and $(XREF bigint, BigInt) are
81  * integer-like by this definition.
82  */
83 template isIntegerLike(T)
84 {
85 	import std.traits : isMutable;
86 	static if (isMutable!T)
87 	{
88 		import std.traits : isIntegral, isSomeChar, isArray, isFloatingPoint;
89 		static if (isIntegral!T ||
90 				   isSomeChar!T)
91 		{
92 			enum isIntegerLike = true;
93 		}
94 		else static if (isFloatingPoint!T ||
95 						is(T == bool) ||
96 						isArray!T)
97 		{
98 			enum isIntegerLike = false;
99 		}
100 		else
101 		{
102 			enum isIntegerLike = is(typeof(
103 			{
104 				T n;
105 				n = 2;
106 				n = n;
107 				n <<= 1;
108 				n >>= 1;
109 				n += n;
110 				n += 2;
111 				n *= n;
112 				n *= 2;
113 				n /= n;
114 				n /= 2;
115 				n -= n;
116 				n -= 2;
117 				n %= 2;
118 				n %= n;
119 				/+ TODO: what about ^^= ? +/
120 				bool lt = n < 2; // less than
121 				bool eq = n == 2; // equal to literal
122 				bool ltg = n < n + 1;
123 				bool seq = n == n; // reflexive equal
124 				return n;
125 			}));
126 		}
127 	}
128 	else
129 	{
130 		import core.internal.traits : Unqual;
131 		alias isIntegerLike = isIntegerLike!(Unqual!T);
132 	}
133 }
134 
135 pure nothrow @safe @nogc unittest {
136 	import std.bigint : BigInt;
137 	import std.meta : AliasSeq;
138 	foreach (T; AliasSeq!(BigInt,
139 						  long, ulong, int, uint,
140 						  short, ushort, byte, ubyte,
141 						  char, wchar, dchar))
142 	{
143 		static assert(isIntegerLike!T);
144 		static assert(isIntegerLike!(const(T)));
145 		static assert(isIntegerLike!(immutable(T)));
146 	}
147 
148 	foreach (T; AliasSeq!(real, double, float, bool))
149 	{
150 		static assert(!isIntegerLike!T);
151 		static assert(!isIntegerLike!(const(T)));
152 		static assert(!isIntegerLike!(immutable(T)));
153 	}
154 }
155 
156 /** Checks if $(D T) has the basic properties of a rational type, i.e.  it has a
157  * numerator and a denominator.
158  */
159 enum isRational(T) = (is(typeof(T.init.numerator)) && /+ TODO: faster to use hasMember? TODO: check that member `isIntegerLike`? +/
160 					  is(typeof(T.init.denominator))); /+ TODO: faster to use hasMember? TODO: check that member `isIntegerLike`? +/
161 
162 /** Returns a Common Integral Type between $(D I1) and $(D I2).  This is defined
163  * as the type returned by I1.init * I2.init.
164  */
165 private template CommonInteger(I1, I2)
166 if (isIntegerLike!I1 &&
167 	isIntegerLike!I2)
168 {
169 	import core.internal.traits : Unqual;
170 	alias CommonInteger = typeof(Unqual!(I1).init *
171 								 Unqual!(I2).init);
172 }
173 
174 pure nothrow @safe @nogc unittest {
175 	import std.bigint : BigInt;
176 	static assert(is(CommonInteger!(BigInt, int) == BigInt));
177 	static assert(is(CommonInteger!(byte, int) == int));
178 }
179 
180 /** Returns a Common Rational Type between $(D R1) and $(D R2), which will be a
181  * Rational based on the CommonInteger of their underlying integer types (or
182  * just on the CommonInteger of ($D R1) and $(D R2), if they themselves are
183  * integers).
184  */
185 private template CommonRational(R1, R2) /+ TODO: avoid recursions below +/
186 {
187 	static if (isRational!R1)
188 	{
189 		alias CommonRational = CommonRational!(typeof(R1.numerator), R2); // recurse. TODO: avoid
190 	}
191 	else static if (isRational!R2)
192 	{
193 		alias CommonRational = CommonRational!(R1, typeof(R2.numerator)); // recurse. TODO: avoid
194 	}
195 	else static if (is(CommonInteger!(R1, R2)))
196 	{
197 		alias CommonRational = Rational!(CommonInteger!(R1, R2));
198 	}
199 }
200 
201 /** Implements rational numbers on top of whatever integer type is specified by
202  * the user.  The integer type used may be any type that behaves as an integer.
203  * Specifically, $(D isIntegerLike) must return true, the integer type must have
204  * value semantics, and the semantics of all integer operations must follow the
205  * normal rules of integer arithmetic.
206  *
207  * A regular integer can be converted to rational type simply by passing it as
208  * a single argument.  In this case the denominator will simply be set to 1.
209  *
210  * Examples:
211  * ---
212  * auto r1 = rational(BigInt("314159265"), BigInt("27182818"));
213  * auto r2 = rational(BigInt("8675309"), BigInt("362436"));
214  * r1 += r2;
215  * assert(r1 == rational(BigInt("174840986505151"),
216  *					   BigInt("4926015912324")));
217  *
218  * // Print result.  Prints:
219  * // "174840986505151/4926015912324"
220  * writeln(r1);
221  *
222  * // Print result in decimal form.  Prints:
223  * // "35.4934"
224  * writeln(cast(real) r1);
225  *
226  * auto r3 = rational(10);
227  * assert(r3.numerator == 10);
228  * assert(r3.denominator == 1);
229  * assert(r3 == 10);
230  * ---
231  */
232 Rational!(CommonInteger!(I1, I2)) rational(I1, I2)(I1 i1, I2 i2)
233 if (isIntegerLike!I1 &&
234 	isIntegerLike!I2)
235 {
236 	static if (is(typeof(typeof(return)(i1, i2))))
237 	{
238 		// Avoid initializing and then reassigning.
239 		auto ret = typeof(return)(i1, i2);
240 	}
241 	else
242 	{
243 		/* Don't want to use void initialization b/c BigInts probably use
244 		 * assignment operator, copy c'tor, etc.
245 		 */
246 		typeof(return) ret;
247 		ret._num = i1;
248 		ret._den = i2;
249 	}
250 	ret.simplify();
251 	return ret;
252 }
253 /// ditto
254 Rational!(I) rational(I)(I val)
255 if (isIntegerLike!I)
256 {
257 	return rational(val, 1);
258 }
259 
260 /** The struct that implements rational numbers.  All relevant operators
261  * (addition, subtraction, multiplication, division, exponentiation by a
262  * non-negative integer, equality and comparison) are overloaded.  The second
263  * operand for all binary operators except exponentiation may be either another
264  * $(D Rational) or another integer type.
265  */
266 struct Rational(SomeIntegral)
267 if (isIntegerLike!SomeIntegral)
268 {
269 public:
270 
271 	// ----------------Multiplication operators----------------------------------
272 	Rational!(SomeIntegral) opBinary(string op, Rhs)(const scope Rhs that) const
273 	if (op == "*" && is(CommonRational!(SomeIntegral, Rhs)) && isRational!Rhs)
274 	{
275 		auto ret = CommonRational!(SomeIntegral, Rhs)(this.numerator, this.denominator);
276 		return ret *= that;
277 	}
278 	/// ditto
279 	Rational!(SomeIntegral) opBinary(string op, Rhs)(const scope Rhs that) const
280 	if (op == "*" && is(CommonRational!(SomeIntegral, Rhs)) && isIntegerLike!Rhs)
281 	{
282 		const factor = gcf(this._den, that);
283 		const adjusted_den = cast(SomeIntegral)(this._den / factor);
284 		const adjusted_rhs = that / factor;
285 		const long new_den = this._num * adjusted_rhs;
286 		if (new_den > SomeIntegral.max)
287 			throw new Exception("");
288 		else
289 			return typeof(this)(cast(SomeIntegral)new_den,
290 								adjusted_den);
291 	}
292 
293 	auto opBinaryRight(string op, Rhs)(const scope Rhs that) const
294 	if (op == "*" && is(CommonRational!(SomeIntegral, Rhs)) && isIntegerLike!Rhs)
295 	{
296 		return opBinary!(op, Rhs)(that); // commutative
297 	}
298 
299 	typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
300 	if (op == "*" && isRational!Rhs)
301 	{
302 		/* Cancel common factors first, then multiply.  This prevents
303 		 * overflows and is much more efficient when using BigInts. */
304 		auto divisor = gcf(this._num, that._den);
305 		this._num /= divisor;
306 		const rhs_den = that._den / divisor;
307 
308 		divisor = gcf(this._den, that._num); // reuse divisor
309 		this._den /= divisor;
310 		const rhs_num = that._num / divisor;
311 
312 		this._num *= rhs_num;
313 		this._den *= rhs_den;
314 
315 		/* no simplify. already cancelled common factors before multiplying. */
316 		fixSigns();
317 
318 		return this;
319 	}
320 	/// ditto
321 	typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
322 	if (op == "*" && isIntegerLike!Rhs)
323 	{
324 		const divisor = gcf(this._den, that);
325 		this._den /= divisor;
326 		this._num *= that / divisor;
327 		/* no to simplify. already cancelled common factors before multiplying. */
328 		fixSigns();
329 		return this;
330 	}
331 
332 	typeof(this) opBinary(string op, Rhs)(const scope Rhs that) const
333 	if (op == "/" &&
334 		is(CommonRational!(int, Rhs)) &&
335 		isRational!Rhs)
336 	{
337 		return typeof(return)(_num * that._den,
338 							  _den * that._num);
339 	}
340 	// ditto
341 	typeof(this) opBinary(string op, Rhs)(const scope Rhs that) const
342 	if (op == "/" &&
343 		is(CommonRational!(int, Rhs)) &&
344 		isIntegerLike!(Rhs))
345 	{
346 		auto ret = CommonRational!(int, Rhs)(this.numerator, this.denominator);
347 		return ret /= that;
348 	}
349 
350 	typeof(this) opBinaryRight(string op, Rhs)(const scope Rhs that) const
351 	if (op == "/" &&
352 		is(CommonRational!(int, Rhs)) &&
353 		isIntegerLike!Rhs)
354 	{
355 		auto ret = CommonRational!(int, Rhs)(this.denominator, this.numerator);
356 		return ret *= that;
357 	}
358 
359 	typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
360 	if (op == "/" &&
361 		isIntegerLike!Rhs)
362 	{
363 		const divisor = gcf(this._num, that);
364 		this._num /= divisor;
365 		this._den *= that / divisor;
366 
367 		/* no to simplify. already cancelled common factors before multiplying. */
368 		fixSigns();
369 		return this;
370 	}
371 
372 	typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
373 	if (op == "/" &&
374 		isRational!Rhs)
375 	{
376 		return this *= that.inverse();
377 	}
378 
379 	// ---------------------addition operators-------------------------------------
380 
381 	auto opBinary(string op, Rhs)(const scope Rhs that) const
382 	if (op == "+" &&
383 		(isRational!Rhs ||
384 		 isIntegerLike!Rhs))
385 	{
386 		auto ret = CommonRational!(typeof(this), Rhs)(this.numerator, this.denominator);
387 		return ret += that;
388 	}
389 
390 	auto opBinaryRight(string op, Rhs)(const scope Rhs that) const
391 	if (op == "+" &&
392 		is(CommonRational!(SomeIntegral, Rhs)) &&
393 		isIntegerLike!Rhs)
394 	{
395 		return opBinary!(op, Rhs)(that); // commutative
396 	}
397 
398 	typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
399 	if (op == "+" &&
400 		isRational!Rhs)
401 	{
402 		if (this._den == that._den)
403 		{
404 			this._num += that._num;
405 			simplify();
406 			return this;
407 		}
408 
409 		SomeIntegral commonDenom = lcm(this._den, that._den);
410 		this._num *= commonDenom / this._den;
411 		this._num += (commonDenom / that._den) * that._num;
412 		this._den = commonDenom;
413 
414 		simplify();
415 		return this;
416 	}
417 
418 	typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
419 	if (op == "+" &&
420 		isIntegerLike!Rhs)
421 	{
422 		this._num += that * this._den;
423 
424 		simplify();
425 		return this;
426 	}
427 
428 	// -----------------------Subtraction operators-------------------------------
429 	auto opBinary(string op, Rhs)(const scope Rhs that) const
430 	if (op == "-" &&
431 		is(CommonRational!(SomeIntegral, Rhs)))
432 	{
433 		auto ret = CommonRational!(typeof(this), Rhs)(this.numerator,
434 													  this.denominator);
435 		return ret -= that;
436 	}
437 
438 	typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
439 	if (op == "-" &&
440 		isRational!Rhs)
441 	{
442 		if (this._den == that._den)
443 		{
444 			this._num -= that._num;
445 			simplify();
446 			return this;
447 		}
448 
449 		auto commonDenom = lcm(this._den, that._den);
450 		this._num *= commonDenom / this._den;
451 		this._num -= (commonDenom / that._den) * that._num;
452 		this._den = commonDenom;
453 
454 		simplify();
455 		return this;
456 	}
457 
458 	typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
459 	if (op == "-" &&
460 		isIntegerLike!Rhs)
461 	{
462 		this._num -= that * this._den;
463 
464 		simplify();
465 		return this;
466 	}
467 
468 	typeof(this) opBinaryRight(string op, Rhs)(const scope Rhs that) const
469 	if (op == "-" &&
470 		is(CommonInteger!(SomeIntegral, Rhs)) &&
471 		isIntegerLike!Rhs)
472 	{
473 		Rational!(SomeIntegral) ret;
474 		ret._den = this._den;
475 		ret._num = (that * this._den) - this._num;
476 
477 		ret.simplify();
478 		return ret;
479 	}
480 
481 	// ----------------------Unary operators---------------------------------------
482 	typeof(this) opUnary(string op)() const
483 	if (op == "-" || op == "+")
484 	{
485 		mixin("return typeof(this)(" ~ op ~ "_num, _den);");
486 	}
487 
488 	// Can only handle integer powers if the result has to also be rational.
489 	typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
490 	if (op == "^^" &&
491 		isIntegerLike!Rhs)
492 	{
493 		if (that < 0)
494 		{
495 			this.invert();
496 			_num ^^= that * -1;
497 			_den ^^= that * -1;
498 		}
499 		else
500 		{
501 			_num ^^= that;
502 			_den ^^= that;
503 		}
504 
505 		/* don't need to simplify here.  this is already simplified, meaning
506 		 * the numerator and denominator don't have any common factors.  raising
507 		 * both to a positive integer power won't create any.
508 		 */
509 		return this;
510 	}
511 	/// ditto
512 	auto opBinary(string op, Rhs)(const scope Rhs that) const
513 	if (op == "^^" &&
514 		isIntegerLike!Rhs &&
515 		is(CommonRational!(SomeIntegral, Rhs)))
516 	{
517 		auto ret = CommonRational!(SomeIntegral, Rhs)(this.numerator, this.denominator);
518 		ret ^^= that;
519 		return ret;
520 	}
521 
522 	import std.traits : isAssignable;
523 
524 	typeof(this) opAssign(Rhs)(const scope Rhs that)
525 	if (isIntegerLike!Rhs &&
526 		isAssignable!(SomeIntegral, Rhs))
527 	{
528 		this._num = that;
529 		this._den = 1;
530 		return this;
531 	}
532 
533 	typeof(this) opAssign(Rhs)(const scope Rhs that)
534 	if (isRational!Rhs &&
535 		isAssignable!(SomeIntegral, typeof(Rhs.numerator)))
536 	{
537 		this._num = that.numerator;
538 		this._den = that.denominator;
539 		return this;
540 	}
541 
542 	bool opEquals(Rhs)(const scope Rhs that) const
543 	if (isRational!Rhs ||
544 		isIntegerLike!Rhs)
545 	{
546 		static if (isRational!Rhs)
547 		{
548 			return (that._num == this._num &&
549 					that._den == this._den);
550 		}
551 		else
552 		{
553 			static assert(isIntegerLike!Rhs);
554 			return (that == this._num &&
555 					this._den == 1);
556 		}
557 	}
558 
559 	int opCmp(Rhs)(const scope Rhs that) const
560 	if (isRational!Rhs)
561 	{
562 		if (opEquals(that))
563 			return 0;
564 
565 		/* Check a few obvious cases first, see if we can avoid having to use a
566 		 * common denominator.  These are basically speed hacks.
567 		 *
568 		 * Assumption:  When simplify() is called, rational will be written in
569 		 * canonical form, with any negative signs being only in the numerator.
570 		 */
571 		if (this._num < 0 &&
572 			that._num > 0)
573 		{
574 			return -1;
575 		}
576 		else if (this._num > 0 &&
577 				 that._num < 0)
578 		{
579 			return 1;
580 		}
581 		else if (this._num >= that._num &&
582 				 this._den <= that._den)
583 		{
584 			// We've already ruled out equality, so this must be > that.
585 			return 1;
586 		}
587 		else if (that._num >= this._num &&
588 				 that._den <= this._den)
589 		{
590 			return -1;
591 		}
592 
593 		// Can't do it without common denominator.  Argh.
594 		auto commonDenom = lcm(this._den, that._den);
595 		auto lhsNum = this._num * (commonDenom / this._den);
596 		auto rhsNum = that._num * (commonDenom / that._den);
597 
598 		if (lhsNum > rhsNum)
599 			return 1;
600 		else if (lhsNum < rhsNum)
601 			return -1;
602 
603 		/* We've checked for equality already.  If we get to this point,
604 		 * there's clearly something wrong.
605 		 */
606 		assert(0);
607 	}
608 
609 	int opCmp(Rhs)(const scope Rhs that) const
610 	if (isIntegerLike!Rhs)
611 	{
612 		if (opEquals(that))
613 			return 0;
614 
615 		// Again, check the obvious cases first.
616 		if (that >= this._num)
617 			return -1;
618 
619 		const that_ = that * this._den;
620 		if (that_ > this._num)
621 			return -1;
622 		else if (that_ < this._num)
623 			return 1;
624 
625 		// Already checked for equality.  If we get here, something's wrong.
626 		assert(0);
627 	}
628 
629 	///////////////////////////////////////////////////////////////////////////////
630 
631 	/// Get inverse.
632 	typeof(this) inverse() const
633 	{
634 		return typeof(return)(_den, _num);
635 	}
636 
637 	/// Invert `this`.
638 	void invert()
639 	{
640 		import std.algorithm.mutation : swap;
641 		swap(_den, _num);
642 	}
643 
644 	import std.traits : isFloatingPoint;
645 
646 	///Convert to floating point representation.
647 	F opCast(F)() const
648 	if (isFloatingPoint!F)
649 	{
650 		import std.traits : isIntegral;
651 		// Do everything in real precision, then convert to F at the end.
652 		static if (isIntegral!(SomeIntegral))
653 		{
654 			return cast(real) _num / _den;
655 		}
656 		else
657 		{
658 			Rational!SomeIntegral temp = this;
659 			real expon = 1.0;
660 			real ans = 0;
661 			byte sign = 1;
662 			if (temp._num < 0)
663 			{
664 				temp._num *= -1;
665 				sign = -1;
666 			}
667 
668 			while (temp._num > 0)
669 			{
670 				while (temp._num < temp._den)
671 				{
672 
673 					assert(temp._den > 0);
674 
675 					static if (is(typeof(temp._den & 1)))
676 					{
677 						// Try to make numbers smaller instead of bigger.
678 						if ((temp._den & 1) == 0)
679 							temp._den >>= 1;
680 						else
681 							temp._num <<= 1;
682 					}
683 					else
684 					{
685 						temp._num <<= 1;
686 					}
687 
688 					expon *= 0.5;
689 
690 					/* This checks for overflow in case we're working with a
691 					 * user-defined fixed-precision integer.
692 					 */
693 					import std.exception : enforce;
694 					import std.conv : text;
695 					enforce(temp._num > 0, text(
696 						"Overflow while converting ", typeof(this).stringof,
697 						" to ", F.stringof, "."));
698 
699 				}
700 
701 				auto intPart = temp._num / temp._den;
702 
703 				static if (is(SomeIntegral == struct) ||
704 						   is(SomeIntegral == class))
705 				{
706 					import std.bigint : BigInt;
707 				}
708 				static if (is(SomeIntegral == BigInt))
709 				{
710 					/* This should really be a cast, but BigInt still has a few
711 					 * issues.
712 					 */
713 					long lIntPart = intPart.toLong();
714 				}
715 				else
716 				{
717 					long lIntPart = cast(long)intPart;
718 				}
719 
720 				// Test for changes.
721 				real oldAns = ans;
722 				ans += lIntPart * expon;
723 				if (ans == oldAns)  // Smaller than epsilon.
724 					return ans * sign;
725 
726 				// Subtract out int part.
727 				temp._num -= intPart * temp._den;
728 			}
729 
730 			return ans * sign;
731 		}
732 	}
733 
734 	/** Casts $(D this) to an integer by truncating the fractional part.
735 	 * Equivalent to $(D integerPart), and then casting it to type $(D I).
736 	 */
737 	I opCast(I)() const
738 	if (isIntegerLike!I &&
739 		is(typeof(cast(I) SomeIntegral.init)))
740 	{
741 		return cast(I) integerPart;
742 	}
743 
744 	///Returns the numerator.
745 	@property inout(SomeIntegral) numerator() inout
746 	{
747 		return _num;
748 	}
749 
750 	///Returns the denominator.
751 	@property inout(SomeIntegral) denominator() inout
752 	{
753 		return _den;
754 	}
755 
756 	/// Returns the integer part of this rational, with any remainder truncated.
757 	@property SomeIntegral integerPart() const
758 	{
759 		return this.numerator / this.denominator;
760 	}
761 
762 	/// Returns the fractional part of this rational.
763 	@property typeof(this) fractionPart() const
764 	{
765 		return this - integerPart;
766 	}
767 
768 	/// Returns a string representation of $(D this) in the form a/b.
769 	string toString() const
770 	{
771 		import std.bigint : BigInt, toDecimalString;
772 		static if (is(SomeIntegral == BigInt))
773 		{
774 			// Special case it for now.  This should be fixed later.
775 			return toDecimalString(_num) ~ "/" ~ toDecimalString(_den);
776 		}
777 		else
778 		{
779 			return to!string(_num) ~ "/" ~ to!string(_den);
780 		}
781 	}
782 
783 private:
784 	SomeIntegral _num;					///< Numerator.
785 	SomeIntegral _den;					///< Denominator.
786 
787 	void simplify()
788 	{
789 		if (_num == 0)
790 		{
791 			_den = 1;
792 			return;
793 		}
794 
795 		const divisor = gcf(_num, _den);
796 		_num /= divisor;
797 		_den /= divisor;
798 
799 		fixSigns();
800 	}
801 
802 	void fixSigns() scope
803 	{
804 		static if (!is(SomeIntegral == ulong) &&
805 				   !is(SomeIntegral == uint) &&
806 				   !is(SomeIntegral == ushort) &&
807 				   !is(SomeIntegral == ubyte))
808 		{
809 			// Write in canonical form w.r.t. signs.
810 			if (_den < 0)
811 			{
812 				_den *= -1;
813 				_num *= -1;
814 			}
815 		}
816 	}
817 }
818 
819 class OverflowException : Exception
820 {
821 pure nothrow @safe @nogc:
822 	this(string msg) { super(msg); }
823 }
824 
825 pure unittest {
826 	import std.bigint : BigInt;
827 	import std.math : isClose;
828 
829 	// All reference values from the Maxima computer algebra system.
830 	// Test c'tor and simplification first.
831 	auto _num = BigInt("295147905179352825852");
832 	auto _den = BigInt("147573952589676412920");
833 	auto simpNum = BigInt("24595658764946068821");
834 	auto simpDen = BigInt("12297829382473034410");
835 	auto f1 = rational(_num, _den);
836 	auto f2 = rational(simpNum, simpDen);
837 	assert(f1 == f2);
838 	// Check that signs of numerator/denominator are corrected
839 	assert(rational(10, -3).numerator == -10);
840 	assert(rational(7, -5).denominator == 5);
841 
842 	// Test multiplication.
843 	assert((rational(0, 1) * rational(1, 1)) == 0);
844 	assert(rational(8, 42) * rational(cast(byte) 7, cast(byte) 68)
845 		   == rational(1, 51));
846 	assert(rational(20_000L, 3_486_784_401U) * rational(3_486_784_401U, 1_000U)
847 		   == rational(20, 1));
848 	auto f3 = rational(7, 57);
849 	f3 *= rational(2, 78);
850 	assert(f3 == rational(7, 2223));
851 	f3 = 5 * f3;
852 	assert(f3 == rational(35, 2223));
853 	assert(f3 * 5UL == 5 * f3);
854 
855 	/* Test division.  Since it's implemented in terms of multiplication,
856 	 * quick and dirty tests should be good enough.
857 	 */
858 	assert(rational(7, 38) / rational(8, 79) == rational(553, 304));
859 	assert(rational(7, 38) / rational(8, 79) == rational(553, 304));
860 	auto f4 = rational(7, 38);
861 	f4 /= rational(8UL, 79);
862 	assert(f4 == rational(553, 304));
863 	f4 = f4 / 2;
864 	assert(f4 == rational(553, 608));
865 	f4 = 2 / f4;
866 	assert(f4 == rational(1216, 553));
867 	assert(f4 * 2 == f4 * rational(2));
868 	f4 = 2;
869 	assert(f4 == 2);
870 
871 	// Test addition.
872 	assert(rational(1, 3) + rational(cast(byte) 2, cast(byte) 3) == rational(1, 1));
873 	assert(rational(1, 3) + rational(1, 2L) == rational(5, 6));
874 	auto f5 = rational( BigInt("314159265"), BigInt("27182818"));
875 	auto f6 = rational( BigInt("8675309"), BigInt("362436"));
876 	f5 += f6;
877 	assert(f5 == rational( BigInt("174840986505151"), BigInt("4926015912324")));
878 	assert(rational(1, 3) + 2UL == rational(7, 3));
879 	assert(5UL + rational(1, 5) == rational(26, 5));
880 
881 	// Test subtraction.
882 	assert(rational(2, 3) - rational(1, 3) == rational(1, 3UL));
883 	assert(rational(1UL, 2) - rational(1, 3) == rational(1, 6));
884 	f5 = rational( BigInt("314159265"), BigInt("27182818"));
885 	f5 -= f6;
886 	assert(f5 == rational( BigInt("-60978359135611"), BigInt("4926015912324")));
887 	assert(rational(4, 3) - 1 == rational(1, 3));
888 	assert(1 - rational(1, 4) == rational(3, 4));
889 
890 	// Test unary operators.
891 	auto fExp = rational(2, 5);
892 	assert(-fExp == rational(-2, 5));
893 	assert(+fExp == rational(2, 5));
894 
895 	// Test exponentiation.
896 	fExp ^^= 3;
897 	assert(fExp == rational(8, 125));
898 	fExp = fExp ^^ 2;
899 	assert(fExp == rational(64, 125 * 125));
900 	assert(rational(2, 5) ^^ -2 == rational(25, 4));
901 
902 	// Test decimal conversion.
903 	assert(isClose(cast(real) f5, -12.37883925284411L));
904 
905 	// Test comparison.
906 	assert(rational(1UL, 6) < rational(1, 2));
907 	assert(rational(cast(byte) 1, cast(byte) 2) > rational(1, 6));
908 	assert(rational(-1, 7) < rational(7, 2));
909 	assert(rational(7, 2) > rational(-1, 7));
910 	assert(rational(7, 9) > rational(8, 11));
911 	assert(rational(8, 11) < rational(7, 9));
912 
913 	assert(rational(9, 10) < 1UL);
914 	assert(1UL > rational(9, 10));
915 	assert(10 > rational(9L, 10));
916 	assert(2 > rational(5, 4));
917 	assert(1 < rational(5U, 4));
918 
919 	// Test creating rationals of value zero.
920 	auto zero = rational(0, 8);
921 	assert(zero == 0);
922 	assert(zero == rational(0, 16));
923 	assert(zero.numerator == 0);
924 	assert(zero.denominator == 1);
925 	auto one = zero + 1;
926 	one -= one;
927 	assert(one == zero);
928 
929 	// Test integerPart, fraction part.
930 	auto intFract = rational(5, 4);
931 	assert(intFract.integerPart == 1);
932 	assert(intFract.fractionPart == rational(1, 4));
933 	assert(cast(long) intFract == 1);
934 
935 	// Test whether CTFE works for primitive types.  Doesn't work yet.
936 	version (none)
937 	{
938 		enum myRational = (((rational(1, 2) + rational(1, 4)) * 2 - rational(1, 4))
939 						   / 2 + 1 * rational(1, 2) - 1) / rational(2, 5);
940 		import std.stdio;
941 		writeln(myRational);
942 		static assert(myRational == rational(-15, 32));
943 	}
944 }
945 
946 /** Convert a floating point number to a $(D Rational) based on integer type $(D
947  * SomeIntegral).  Allows an error tolerance of $(D epsilon).  (Default $(D epsilon) =
948  * 1e-8.)
949  *
950  * $(D epsilon) must be greater than 1.0L / long.max.
951  *
952  * Throws:  Exception on infinities, NaNs, numbers with absolute value
953  * larger than long.max and epsilons smaller than 1.0L / long.max.
954  *
955  * Examples:
956  * ---
957  * // Prints "22/7".
958  * writeln(toRational!int(PI, 1e-1));
959  * ---
960  */
961 Rational!(SomeIntegral) toRational(SomeIntegral)(real floatNum, real epsilon = 1e-8)
962 {
963 	import std.math: isNaN;
964 	import std.exception : enforce;
965 	enforce(floatNum != real.infinity &&
966 			floatNum != -real.infinity &&
967 			!isNaN(floatNum),
968 			"Can't convert NaNs and infinities to rational.");
969 	enforce(floatNum < long.max &&
970 			floatNum > -long.max,
971 			"Rational conversions of very large numbers not yet implemented.");
972 	enforce(1.0L / epsilon < long.max,
973 			"Can't handle very small epsilons < long.max in toRational.");
974 
975 	/* Handle this as a special case to make the rest of the code less
976 	 * complicated:
977 	 */
978 	if (abs(floatNum) < epsilon)
979 	{
980 		Rational!SomeIntegral ret;
981 		ret._num = 0;
982 		ret._den = 1;
983 		return ret;
984 	}
985 
986 	return toRationalImpl!(SomeIntegral)(floatNum, epsilon);
987 }
988 
989 private Rational!SomeIntegral toRationalImpl(SomeIntegral)(real floatNum, real epsilon)
990 {
991 	import std.conv : roundTo;
992 	import std.traits : isIntegral;
993 
994 	real actualEpsilon;
995 	Rational!SomeIntegral ret;
996 
997 	if (abs(floatNum) < 1)
998 	{
999 		real invFloatNum = 1.0L / floatNum;
1000 		long intPart = roundTo!long(invFloatNum);
1001 		actualEpsilon = floatNum - 1.0L / intPart;
1002 
1003 		static if (isIntegral!(SomeIntegral))
1004 		{
1005 			ret._den = cast(SomeIntegral) intPart;
1006 			ret._num = cast(SomeIntegral) 1;
1007 		}
1008 		else
1009 		{
1010 			ret._den = intPart;
1011 			ret._num = 1;
1012 		}
1013 	}
1014 	else
1015 	{
1016 		long intPart = roundTo!long(floatNum);
1017 		actualEpsilon = floatNum - intPart;
1018 
1019 		static if (isIntegral!(SomeIntegral))
1020 		{
1021 			ret._den = cast(SomeIntegral) 1;
1022 			ret._num = cast(SomeIntegral) intPart;
1023 		}
1024 		else
1025 		{
1026 			ret._den = 1;
1027 			ret._num = intPart;
1028 		}
1029 	}
1030 
1031 	if (abs(actualEpsilon) <= epsilon)
1032 		return ret;
1033 
1034 	// Else get results from downstream recursions, add them to this result.
1035 	return ret + toRationalImpl!(SomeIntegral)(actualEpsilon, epsilon);
1036 }
1037 
1038 unittest {
1039 	import std.bigint : BigInt;
1040 	import std.math : PI, E;
1041 	// Start with simple cases.
1042 	assert(toRational!int(0.5) == rational(1, 2));
1043 	assert(toRational!BigInt(0.333333333333333L) == rational(BigInt(1), BigInt(3)));
1044 	assert(toRational!int(2.470588235294118) == rational(cast(int) 42, cast(int) 17));
1045 	assert(toRational!long(2.007874015748032) == rational(255L, 127L));
1046 	assert(toRational!int( 3.0L / 7.0L) == rational(3, 7));
1047 	assert(toRational!int( 7.0L / 3.0L) == rational(7, 3));
1048 
1049 	// Now for some fun.
1050 	real myEpsilon = 1e-8;
1051 	auto piRational = toRational!long(PI, myEpsilon);
1052 	assert(abs(cast(real) piRational - PI) < myEpsilon);
1053 
1054 	auto eRational = toRational!long(E, myEpsilon);
1055 	assert(abs(cast(real) eRational - E) < myEpsilon);
1056 }
1057 
1058 /** Find the Greatest Common Factor (GCF), aka Greatest Common Divisor (GCD), of
1059  * $(D m) and $(D n).
1060  */
1061 CommonInteger!(I1, I2) gcf(I1, I2)(I1 m, I2 n)
1062 if (isIntegerLike!I1 &&
1063 	isIntegerLike!I2)
1064 {
1065 	static if (is(I1 == const) || is(I1 == immutable) ||
1066 			   is(I2 == const) || is(I2 == immutable))
1067 	{
1068 		// Doesn't work with immutable(BigInt).
1069 		import core.internal.traits : Unqual;
1070 		return gcf!(Unqual!I1,
1071 					Unqual!I2)(m, n);
1072 	}
1073 	else
1074 	{
1075 		typeof(return) a = abs(m);
1076 		typeof(return) b = abs(n);
1077 
1078 		while (b)
1079 		{
1080 			auto t = b;
1081 			b = a % b;
1082 			a = t;
1083 		}
1084 
1085 		return a;
1086 	}
1087 }
1088 
1089 pure unittest {
1090 	import std.bigint : BigInt;
1091 	assert(gcf(0, 0) == 0);
1092 	assert(gcf(0, 1) == 1);
1093 	assert(gcf(999, 0) == 999);
1094 	assert(gcf(to!(immutable(int))(8), to!(const(int))(12)) == 4);
1095 
1096 	// Values from the Maxima computer algebra system.
1097 	assert(gcf(BigInt(314_156_535UL), BigInt(27_182_818_284UL)) == BigInt(3));
1098 	assert(gcf(8675309, 362436) == 1);
1099 	assert(gcf(BigInt("8589934596"), BigInt("295147905179352825852")) == 12);
1100 }
1101 
1102 /// Find the Least Common Multiple (LCM) of $(D a) and $(D b).
1103 CommonInteger!(I1, I2) lcm(I1, I2)(const scope I1 a, const scope I2 b)
1104 if (isIntegerLike!I1 &&
1105 	isIntegerLike!I2)
1106 {
1107 	const n1 = abs(a);
1108 	const n2 = abs(b);
1109 	if (n1 == n2)
1110 		return n1;
1111 	return (n1 / gcf(n1, n2)) * n2;
1112 }
1113 
1114 /// Returns the largest integer less than or equal to $(D r).
1115 SomeIntegral floor(SomeIntegral)(const scope Rational!SomeIntegral r)
1116 {
1117 	SomeIntegral intPart = r.integerPart;
1118 	if (r > 0 || intPart == r)
1119 		return intPart;
1120 	else
1121 	{
1122 		intPart -= 1;
1123 		return intPart;
1124 	}
1125 }
1126 
1127 pure nothrow @safe @nogc unittest {
1128 	assert(floor(rational(1, 2)) == 0);
1129 	assert(floor(rational(-1, 2)) == -1);
1130 	assert(floor(rational(2)) == 2);
1131 	assert(floor(rational(-2)) == -2);
1132 	assert(floor(rational(-1, 2)) == -1);
1133 }
1134 
1135 /// Returns the smallest integer greater than or equal to $(D r).
1136 SomeIntegral ceil(SomeIntegral)(const scope Rational!SomeIntegral r)
1137 {
1138 	SomeIntegral intPart = r.integerPart;
1139 	if (intPart == r || r < 0)
1140 		return intPart;
1141 	else
1142 	{
1143 		intPart += 1;
1144 		return intPart;
1145 	}
1146 }
1147 
1148 pure nothrow @safe @nogc unittest {
1149 	assert(ceil(rational(1, 2)) == 1);
1150 	assert(ceil(rational(0)) == 0);
1151 	assert(ceil(rational(-1, 2)) == 0);
1152 	assert(ceil(rational(1)) == 1);
1153 	assert(ceil(rational(-2)) == -2);
1154 }
1155 
1156 /** Round $(D r) to the nearest integer.  If the fractional part is exactly 1/2,
1157  * $(D r) will be rounded such that the absolute value is increased by rounding.
1158  */
1159 SomeIntegral round(SomeIntegral)(const scope Rational!SomeIntegral r)
1160 {
1161 	auto intPart = r.integerPart;
1162 	const fractPart = r.fractionPart;
1163 
1164 	bool added;
1165 	if (fractPart >= rational(1, 2))
1166 	{
1167 		added = true;
1168 		intPart += 1;
1169 	}
1170 
1171 	import std.traits : isUnsigned;
1172 
1173 	static if (!isUnsigned!SomeIntegral)
1174 	{
1175 		if (!added && fractPart <= rational(-1, 2))
1176 			intPart -= 1;
1177 	}
1178 
1179 	return intPart;
1180 }
1181 
1182 pure nothrow @safe @nogc unittest {
1183 	assert(round(rational(1, 3)) == 0);
1184 	assert(round(rational(7, 2)) == 4);
1185 	assert(round(rational(-3, 4)) == -1);
1186 	assert(round(rational(8U, 15U)) == 1);
1187 }