/** This module contains an implementation of rational numbers that is templated
 * on the underlying integer type.  It can be used with either builtin fixed
 * width integers or arbitrary precision integers.  All relevant operators are
 * overloaded for both rational-rational and rational-integer operations.
 *
 * Synopsis:
 * ---
 * // Compute pi using the generalized continued fraction approximation.
 * import std.bigint, std.rational, std.stdio;
 *
 * enum maxTerm = 30;
 *
 * Rational!(BigInt) getTerm(int termNumber)
 * {
 *     auto addFactor = 2 * termNumber - 1;
 *
 *     if (termNumber == maxTerm)
 *     {
 *         return rational(BigInt(addFactor));
 *     }
 *
 *     auto termNumberSquared = BigInt(termNumber * termNumber);
 *     auto continued = termNumberSquared / getTerm(termNumber + 1);
 *
 *     continued += addFactor;
 *     return continued;
 * }
 *
 * void main()
 * {
 *     auto pi = rational(BigInt(4)) / getTerm(1);
 *
 *     // Display the result in rational form.
 *     writeln(pi);
 *
 *     // Display the decimal equivalent, which is accurate to 18 decimal places.
 *     writefln("%.18f", cast(real) pi);
 * }
 * ---
 *
 *
 * Author:  David Simcha
 * Copyright:  Copyright (c) 2009-2011, David Simcha.
 * License:    $(LINK2 http://www.boost.org/LICENSE_1_0.txt, Boost License 1.0)
 */
module nxt.rational;

import std.conv : to;
import std.math : abs;

/** Checks whether $(D T) is structurally an integer, i.e. whether it supports
 * all of the operations an integer type should support.  Does not check the
 * nominal type of $(D T).  In particular, for a mutable type $(D T) the
 * following must compile:
 *
 * ---
 * T n;
 * n = 2;
 * n <<= 1;
 * n >>= 1;
 * n += n;
 * n += 2;
 * n *= n;
 * n *= 2;
 * n /= n;
 * n /= 2;
 * n -= n;
 * n -= 2;
 * n %= 2;
 * n %= n;
 * bool foo = n < 2;
 * bool bar = n == 2;
 * bool goo = n < n + 1;
 * bool tar = n == n;
 * ---
 *
 * while for a non-mutable type, the above must compile for its unqualified,
 * mutable variant.
 *
 * All built-in D integers and character types and $(XREF bigint, BigInt) are
 * integer-like by this definition.
 */
template isIntegerLike(T)
{
    import std.traits : isMutable;
    static if (isMutable!T)
    {
        import std.traits : isIntegral, isSomeChar, isArray, isFloatingPoint;
        static if (isIntegral!T ||
                   isSomeChar!T)
        {
            enum isIntegerLike = true;
        }
        else static if (isFloatingPoint!T ||
                        is(T == bool) ||
                        isArray!T)
        {
            enum isIntegerLike = false;
        }
        else
        {
            enum isIntegerLike = is(typeof(
            {
                T n;
                n = 2;
                n = n;
                n <<= 1;
                n >>= 1;
                n += n;
                n += 2;
                n *= n;
                n *= 2;
                n /= n;
                n /= 2;
                n -= n;
                n -= 2;
                n %= 2;
                n %= n;
                // TODO: what about ^^= ?
                bool lt = n < 2; // less than
                bool eq = n == 2; // equal to literal
                bool ltg = n < n + 1;
                bool seq = n == n; // reflexive equal
                return n;
            }));
        }
    }
    else
    {
        import core.internal.traits : Unqual;
        alias isIntegerLike = isIntegerLike!(Unqual!T);
    }
}

@safe pure nothrow @nogc unittest
{
    import std.bigint : BigInt;
    import std.meta : AliasSeq;
    foreach (T; AliasSeq!(BigInt,
                          long, ulong, int, uint,
                          short, ushort, byte, ubyte,
                          char, wchar, dchar))
    {
        static assert(isIntegerLike!T);
        static assert(isIntegerLike!(const(T)));
        static assert(isIntegerLike!(immutable(T)));
    }

    foreach (T; AliasSeq!(real, double, float, bool))
    {
        static assert(!isIntegerLike!T);
        static assert(!isIntegerLike!(const(T)));
        static assert(!isIntegerLike!(immutable(T)));
    }
}

/** Checks if $(D T) has the basic properties of a rational type, i.e.  it has a
 * numerator and a denominator.
 */
enum isRational(T) = (is(typeof(T.init.numerator)) && // TODO: faster to use hasMember? TODO: check that member `isIntegerLike`?
                      is(typeof(T.init.denominator))); // TODO: faster to use hasMember? TODO: check that member `isIntegerLike`?

/** Returns a Common Integral Type between $(D I1) and $(D I2).  This is defined
 * as the type returned by I1.init * I2.init.
 */
private template CommonInteger(I1, I2)
if (isIntegerLike!I1 &&
    isIntegerLike!I2)
{
    import core.internal.traits : Unqual;
    alias CommonInteger = typeof(Unqual!(I1).init *
                                 Unqual!(I2).init);
}

@safe pure nothrow @nogc unittest
{
    import std.bigint : BigInt;
    static assert(is(CommonInteger!(BigInt, int) == BigInt));
    static assert(is(CommonInteger!(byte, int) == int));
}

/** Returns a Common Rational Type between $(D R1) and $(D R2), which will be a
 * Rational based on the CommonInteger of their underlying integer types (or
 * just on the CommonInteger of ($D R1) and $(D R2), if they themselves are
 * integers).
 */
private template CommonRational(R1, R2) // TODO: avoid recursions below
{
    static if (isRational!R1)
    {
        alias CommonRational = CommonRational!(typeof(R1.numerator), R2); // recurse. TODO: avoid
    }
    else static if (isRational!R2)
    {
        alias CommonRational = CommonRational!(R1, typeof(R2.numerator)); // recurse. TODO: avoid
    }
    else static if (is(CommonInteger!(R1, R2)))
    {
        alias CommonRational = Rational!(CommonInteger!(R1, R2));
    }
}

/** Implements rational numbers on top of whatever integer type is specified by
 * the user.  The integer type used may be any type that behaves as an integer.
 * Specifically, $(D isIntegerLike) must return true, the integer type must have
 * value semantics, and the semantics of all integer operations must follow the
 * normal rules of integer arithmetic.
 *
 * A regular integer can be converted to rational type simply by passing it as
 * a single argument.  In this case the denominator will simply be set to 1.
 *
 * Examples:
 * ---
 * auto r1 = rational(BigInt("314159265"), BigInt("27182818"));
 * auto r2 = rational(BigInt("8675309"), BigInt("362436"));
 * r1 += r2;
 * assert(r1 == rational(BigInt("174840986505151"),
 *                       BigInt("4926015912324")));
 *
 * // Print result.  Prints:
 * // "174840986505151/4926015912324"
 * writeln(r1);
 *
 * // Print result in decimal form.  Prints:
 * // "35.4934"
 * writeln(cast(real) r1);
 *
 * auto r3 = rational(10);
 * assert(r3.numerator == 10);
 * assert(r3.denominator == 1);
 * assert(r3 == 10);
 * ---
 */
Rational!(CommonInteger!(I1, I2)) rational(I1, I2)(I1 i1, I2 i2)
if (isIntegerLike!I1 &&
    isIntegerLike!I2)
{
    static if (is(typeof(typeof(return)(i1, i2))))
    {
        // Avoid initializing and then reassigning.
        auto ret = typeof(return)(i1, i2);
    }
    else
    {
        /* Don't want to use void initialization b/c BigInts probably use
         * assignment operator, copy c'tor, etc.
         */
        typeof(return) ret;
        ret._num = i1;
        ret._den = i2;
    }
    ret.simplify();
    return ret;
}
/// ditto
Rational!(I) rational(I)(I val)
if (isIntegerLike!I)
{
    return rational(val, 1);
}

/** The struct that implements rational numbers.  All relevant operators
 * (addition, subtraction, multiplication, division, exponentiation by a
 * non-negative integer, equality and comparison) are overloaded.  The second
 * operand for all binary operators except exponentiation may be either another
 * $(D Rational) or another integer type.
 */
struct Rational(SomeIntegral)
if (isIntegerLike!SomeIntegral)
{
public:

    // ----------------Multiplication operators----------------------------------
    Rational!(SomeIntegral) opBinary(string op, Rhs)(const scope Rhs that) const
    if (op == "*" && is(CommonRational!(SomeIntegral, Rhs)) && isRational!Rhs)
    {
        auto ret = CommonRational!(SomeIntegral, Rhs)(this.numerator, this.denominator);
        return ret *= that;
    }
    /// ditto
    Rational!(SomeIntegral) opBinary(string op, Rhs)(const scope Rhs that) const
    if (op == "*" && is(CommonRational!(SomeIntegral, Rhs)) && isIntegerLike!Rhs)
    {
        const factor = gcf(this._den, that);
        const adjusted_den = cast(SomeIntegral)(this._den / factor);
        const adjusted_rhs = that / factor;
        const long new_den = this._num * adjusted_rhs;
        if (new_den > SomeIntegral.max)
            throw new Exception("");
        else
            return typeof(this)(cast(SomeIntegral)new_den,
                                adjusted_den);
    }

    auto opBinaryRight(string op, Rhs)(const scope Rhs that) const
    if (op == "*" && is(CommonRational!(SomeIntegral, Rhs)) && isIntegerLike!Rhs)
    {
        return opBinary!(op, Rhs)(that); // commutative
    }

    typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
    if (op == "*" && isRational!Rhs)
    {
        /* Cancel common factors first, then multiply.  This prevents
         * overflows and is much more efficient when using BigInts. */
        auto divisor = gcf(this._num, that._den);
        this._num /= divisor;
        const rhs_den = that._den / divisor;

        divisor = gcf(this._den, that._num); // reuse divisor
        this._den /= divisor;
        const rhs_num = that._num / divisor;

        this._num *= rhs_num;
        this._den *= rhs_den;

        /* no simplify. already cancelled common factors before multiplying. */
        fixSigns();

        return this;
    }
    /// ditto
    typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
    if (op == "*" && isIntegerLike!Rhs)
    {
        const divisor = gcf(this._den, that);
        this._den /= divisor;
        this._num *= that / divisor;
        /* no to simplify. already cancelled common factors before multiplying. */
        fixSigns();
        return this;
    }

    typeof(this) opBinary(string op, Rhs)(const scope Rhs that) const
    if (op == "/" &&
        is(CommonRational!(int, Rhs)) &&
        isRational!Rhs)
    {
        return typeof(return)(_num * that._den,
                              _den * that._num);
    }
    // ditto
    typeof(this) opBinary(string op, Rhs)(const scope Rhs that) const
    if (op == "/" &&
        is(CommonRational!(int, Rhs)) &&
        isIntegerLike!(Rhs))
    {
        auto ret = CommonRational!(int, Rhs)(this.numerator, this.denominator);
        return ret /= that;
    }

    typeof(this) opBinaryRight(string op, Rhs)(const scope Rhs that) const
    if (op == "/" &&
        is(CommonRational!(int, Rhs)) &&
        isIntegerLike!Rhs)
    {
        auto ret = CommonRational!(int, Rhs)(this.denominator, this.numerator);
        return ret *= that;
    }

    typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
    if (op == "/" &&
        isIntegerLike!Rhs)
    {
        const divisor = gcf(this._num, that);
        this._num /= divisor;
        this._den *= that / divisor;

        /* no to simplify. already cancelled common factors before multiplying. */
        fixSigns();
        return this;
    }

    typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
    if (op == "/" &&
        isRational!Rhs)
    {
        return this *= that.inverse();
    }

    // ---------------------addition operators-------------------------------------

    auto opBinary(string op, Rhs)(const scope Rhs that) const
    if (op == "+" &&
        (isRational!Rhs ||
         isIntegerLike!Rhs))
    {
        auto ret = CommonRational!(typeof(this), Rhs)(this.numerator, this.denominator);
        return ret += that;
    }

    auto opBinaryRight(string op, Rhs)(const scope Rhs that) const
    if (op == "+" &&
        is(CommonRational!(SomeIntegral, Rhs)) &&
        isIntegerLike!Rhs)
    {
        return opBinary!(op, Rhs)(that); // commutative
    }

    typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
    if (op == "+" &&
        isRational!Rhs)
    {
        if (this._den == that._den)
        {
            this._num += that._num;
            simplify();
            return this;
        }

        SomeIntegral commonDenom = lcm(this._den, that._den);
        this._num *= commonDenom / this._den;
        this._num += (commonDenom / that._den) * that._num;
        this._den = commonDenom;

        simplify();
        return this;
    }

    typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
    if (op == "+" &&
        isIntegerLike!Rhs)
    {
        this._num += that * this._den;

        simplify();
        return this;
    }

    // -----------------------Subtraction operators-------------------------------
    auto opBinary(string op, Rhs)(const scope Rhs that) const
    if (op == "-" &&
        is(CommonRational!(SomeIntegral, Rhs)))
    {
        auto ret = CommonRational!(typeof(this), Rhs)(this.numerator,
                                                      this.denominator);
        return ret -= that;
    }

    typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
    if (op == "-" &&
        isRational!Rhs)
    {
        if (this._den == that._den)
        {
            this._num -= that._num;
            simplify();
            return this;
        }

        auto commonDenom = lcm(this._den, that._den);
        this._num *= commonDenom / this._den;
        this._num -= (commonDenom / that._den) * that._num;
        this._den = commonDenom;

        simplify();
        return this;
    }

    typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
    if (op == "-" &&
        isIntegerLike!Rhs)
    {
        this._num -= that * this._den;

        simplify();
        return this;
    }

    typeof(this) opBinaryRight(string op, Rhs)(const scope Rhs that) const
    if (op == "-" &&
        is(CommonInteger!(SomeIntegral, Rhs)) &&
        isIntegerLike!Rhs)
    {
        Rational!(SomeIntegral) ret;
        ret._den = this._den;
        ret._num = (that * this._den) - this._num;

        ret.simplify();
        return ret;
    }

    // ----------------------Unary operators---------------------------------------
    typeof(this) opUnary(string op)() const
    if (op == "-" || op == "+")
    {
        mixin("return typeof(this)(" ~ op ~ "_num, _den);");
    }

    // Can only handle integer powers if the result has to also be rational.
    typeof(this) opOpAssign(string op, Rhs)(const scope Rhs that)
    if (op == "^^" &&
        isIntegerLike!Rhs)
    {
        if (that < 0)
        {
            this.invert();
            _num ^^= that * -1;
            _den ^^= that * -1;
        }
        else
        {
            _num ^^= that;
            _den ^^= that;
        }

        /* don't need to simplify here.  this is already simplified, meaning
         * the numerator and denominator don't have any common factors.  raising
         * both to a positive integer power won't create any.
         */
        return this;
    }
    /// ditto
    auto opBinary(string op, Rhs)(const scope Rhs that) const
    if (op == "^^" &&
        isIntegerLike!Rhs &&
        is(CommonRational!(SomeIntegral, Rhs)))
    {
        auto ret = CommonRational!(SomeIntegral, Rhs)(this.numerator, this.denominator);
        ret ^^= that;
        return ret;
    }

    import std.traits : isAssignable;

    typeof(this) opAssign(Rhs)(const scope Rhs that)
    if (isIntegerLike!Rhs &&
        isAssignable!(SomeIntegral, Rhs))
    {
        this._num = that;
        this._den = 1;
        return this;
    }

    typeof(this) opAssign(Rhs)(const scope Rhs that)
    if (isRational!Rhs &&
        isAssignable!(SomeIntegral, typeof(Rhs.numerator)))
    {
        this._num = that.numerator;
        this._den = that.denominator;
        return this;
    }

    bool opEquals(Rhs)(const scope Rhs that) const
    if (isRational!Rhs ||
        isIntegerLike!Rhs)
    {
        static if (isRational!Rhs)
        {
            return (that._num == this._num &&
                    that._den == this._den);
        }
        else
        {
            static assert(isIntegerLike!Rhs);
            return (that == this._num &&
                    this._den == 1);
        }
    }

    int opCmp(Rhs)(const scope Rhs that) const
    if (isRational!Rhs)
    {
        if (opEquals(that))
            return 0;

        /* Check a few obvious cases first, see if we can avoid having to use a
         * common denominator.  These are basically speed hacks.
         *
         * Assumption:  When simplify() is called, rational will be written in
         * canonical form, with any negative signs being only in the numerator.
         */
        if (this._num < 0 &&
            that._num > 0)
        {
            return -1;
        }
        else if (this._num > 0 &&
                 that._num < 0)
        {
            return 1;
        }
        else if (this._num >= that._num &&
                 this._den <= that._den)
        {
            // We've already ruled out equality, so this must be > that.
            return 1;
        }
        else if (that._num >= this._num &&
                 that._den <= this._den)
        {
            return -1;
        }

        // Can't do it without common denominator.  Argh.
        auto commonDenom = lcm(this._den, that._den);
        auto lhsNum = this._num * (commonDenom / this._den);
        auto rhsNum = that._num * (commonDenom / that._den);

        if (lhsNum > rhsNum)
            return 1;
        else if (lhsNum < rhsNum)
            return -1;

        /* We've checked for equality already.  If we get to this point,
         * there's clearly something wrong.
         */
        assert(0);
    }

    int opCmp(Rhs)(const scope Rhs that) const
    if (isIntegerLike!Rhs)
    {
        if (opEquals(that))
            return 0;

        // Again, check the obvious cases first.
        if (that >= this._num)
            return -1;

        const that_ = that * this._den;
        if (that_ > this._num)
            return -1;
        else if (that_ < this._num)
            return 1;

        // Already checked for equality.  If we get here, something's wrong.
        assert(0);
    }

    ///////////////////////////////////////////////////////////////////////////////

    /// Get inverse.
    typeof(this) inverse() const
    {
        return typeof(return)(_den, _num);
    }

    /// Invert `this`.
    void invert()
    {
        import std.algorithm.mutation : swap;
        swap(_den, _num);
    }

    import std.traits : isFloatingPoint;

    ///Convert to floating point representation.
    F opCast(F)() const
    if (isFloatingPoint!F)
    {
        import std.traits : isIntegral;
        // Do everything in real precision, then convert to F at the end.
        static if (isIntegral!(SomeIntegral))
        {
            return cast(real) _num / _den;
        }
        else
        {
            Rational!SomeIntegral temp = this;
            real expon = 1.0;
            real ans = 0;
            byte sign = 1;
            if (temp._num < 0)
            {
                temp._num *= -1;
                sign = -1;
            }

            while (temp._num > 0)
            {
                while (temp._num < temp._den)
                {

                    assert(temp._den > 0);

                    static if (is(typeof(temp._den & 1)))
                    {
                        // Try to make numbers smaller instead of bigger.
                        if ((temp._den & 1) == 0)
                            temp._den >>= 1;
                        else
                            temp._num <<= 1;
                    }
                    else
                    {
                        temp._num <<= 1;
                    }

                    expon *= 0.5;

                    /* This checks for overflow in case we're working with a
                     * user-defined fixed-precision integer.
                     */
                    import std.exception : enforce;
                    import std.conv : text;
                    enforce(temp._num > 0, text(
                        "Overflow while converting ", typeof(this).stringof,
                        " to ", F.stringof, "."));

                }

                auto intPart = temp._num / temp._den;

                static if (is(SomeIntegral == struct) ||
                           is(SomeIntegral == class))
                {
                    import std.bigint : BigInt;
                }
                static if (is(SomeIntegral == BigInt))
                {
                    /* This should really be a cast, but BigInt still has a few
                     * issues.
                     */
                    long lIntPart = intPart.toLong();
                }
                else
                {
                    long lIntPart = cast(long)intPart;
                }

                // Test for changes.
                real oldAns = ans;
                ans += lIntPart * expon;
                if (ans == oldAns)  // Smaller than epsilon.
                    return ans * sign;

                // Subtract out int part.
                temp._num -= intPart * temp._den;
            }

            return ans * sign;
        }
    }

    /** Casts $(D this) to an integer by truncating the fractional part.
     * Equivalent to $(D integerPart), and then casting it to type $(D I).
     */
    I opCast(I)() const
    if (isIntegerLike!I &&
        is(typeof(cast(I) SomeIntegral.init)))
    {
        return cast(I) integerPart;
    }

    ///Returns the numerator.
    @property inout(SomeIntegral) numerator() inout
    {
        return _num;
    }

    ///Returns the denominator.
    @property inout(SomeIntegral) denominator() inout
    {
        return _den;
    }

    /// Returns the integer part of this rational, with any remainder truncated.
    @property SomeIntegral integerPart() const
    {
        return this.numerator / this.denominator;
    }

    /// Returns the fractional part of this rational.
    @property typeof(this) fractionPart() const
    {
        return this - integerPart;
    }

    /// Returns a string representation of $(D this) in the form a/b.
    string toString() const
    {
        import std.bigint : BigInt, toDecimalString;
        static if (is(SomeIntegral == BigInt))
        {
            // Special case it for now.  This should be fixed later.
            return toDecimalString(_num) ~ "/" ~ toDecimalString(_den);
        }
        else
        {
            return to!string(_num) ~ "/" ~ to!string(_den);
        }
    }

private:
    SomeIntegral _num;                    ///< Numerator.
    SomeIntegral _den;                    ///< Denominator.

    void simplify()
    {
        if (_num == 0)
        {
            _den = 1;
            return;
        }

        const divisor = gcf(_num, _den);
        _num /= divisor;
        _den /= divisor;

        fixSigns();
    }

    void fixSigns() scope
    {
        static if (!is(SomeIntegral == ulong) &&
                   !is(SomeIntegral == uint) &&
                   !is(SomeIntegral == ushort) &&
                   !is(SomeIntegral == ubyte))
        {
            // Write in canonical form w.r.t. signs.
            if (_den < 0)
            {
                _den *= -1;
                _num *= -1;
            }
        }
    }
}

class OverflowException : Exception
{
@safe pure nothrow @nogc:
    this(string msg) { super(msg); }
}

pure unittest
{
    import std.bigint : BigInt;
    import std.math : isClose;

    // All reference values from the Maxima computer algebra system.
    // Test c'tor and simplification first.
    auto _num = BigInt("295147905179352825852");
    auto _den = BigInt("147573952589676412920");
    auto simpNum = BigInt("24595658764946068821");
    auto simpDen = BigInt("12297829382473034410");
    auto f1 = rational(_num, _den);
    auto f2 = rational(simpNum, simpDen);
    assert(f1 == f2);
    // Check that signs of numerator/denominator are corrected
    assert(rational(10, -3).numerator == -10);
    assert(rational(7, -5).denominator == 5);

    // Test multiplication.
    assert((rational(0, 1) * rational(1, 1)) == 0);
    assert(rational(8, 42) * rational(cast(byte) 7, cast(byte) 68)
           == rational(1, 51));
    assert(rational(20_000L, 3_486_784_401U) * rational(3_486_784_401U, 1_000U)
           == rational(20, 1));
    auto f3 = rational(7, 57);
    f3 *= rational(2, 78);
    assert(f3 == rational(7, 2223));
    f3 = 5 * f3;
    assert(f3 == rational(35, 2223));
    assert(f3 * 5UL == 5 * f3);

    /* Test division.  Since it's implemented in terms of multiplication,
     * quick and dirty tests should be good enough.
     */
    assert(rational(7, 38) / rational(8, 79) == rational(553, 304));
    assert(rational(7, 38) / rational(8, 79) == rational(553, 304));
    auto f4 = rational(7, 38);
    f4 /= rational(8UL, 79);
    assert(f4 == rational(553, 304));
    f4 = f4 / 2;
    assert(f4 == rational(553, 608));
    f4 = 2 / f4;
    assert(f4 == rational(1216, 553));
    assert(f4 * 2 == f4 * rational(2));
    f4 = 2;
    assert(f4 == 2);

    // Test addition.
    assert(rational(1, 3) + rational(cast(byte) 2, cast(byte) 3) == rational(1, 1));
    assert(rational(1, 3) + rational(1, 2L) == rational(5, 6));
    auto f5 = rational( BigInt("314159265"), BigInt("27182818"));
    auto f6 = rational( BigInt("8675309"), BigInt("362436"));
    f5 += f6;
    assert(f5 == rational( BigInt("174840986505151"), BigInt("4926015912324")));
    assert(rational(1, 3) + 2UL == rational(7, 3));
    assert(5UL + rational(1, 5) == rational(26, 5));

    // Test subtraction.
    assert(rational(2, 3) - rational(1, 3) == rational(1, 3UL));
    assert(rational(1UL, 2) - rational(1, 3) == rational(1, 6));
    f5 = rational( BigInt("314159265"), BigInt("27182818"));
    f5 -= f6;
    assert(f5 == rational( BigInt("-60978359135611"), BigInt("4926015912324")));
    assert(rational(4, 3) - 1 == rational(1, 3));
    assert(1 - rational(1, 4) == rational(3, 4));

    // Test unary operators.
    auto fExp = rational(2, 5);
    assert(-fExp == rational(-2, 5));
    assert(+fExp == rational(2, 5));

    // Test exponentiation.
    fExp ^^= 3;
    assert(fExp == rational(8, 125));
    fExp = fExp ^^ 2;
    assert(fExp == rational(64, 125 * 125));
    assert(rational(2, 5) ^^ -2 == rational(25, 4));

    // Test decimal conversion.
    assert(isClose(cast(real) f5, -12.37883925284411L));

    // Test comparison.
    assert(rational(1UL, 6) < rational(1, 2));
    assert(rational(cast(byte) 1, cast(byte) 2) > rational(1, 6));
    assert(rational(-1, 7) < rational(7, 2));
    assert(rational(7, 2) > rational(-1, 7));
    assert(rational(7, 9) > rational(8, 11));
    assert(rational(8, 11) < rational(7, 9));

    assert(rational(9, 10) < 1UL);
    assert(1UL > rational(9, 10));
    assert(10 > rational(9L, 10));
    assert(2 > rational(5, 4));
    assert(1 < rational(5U, 4));

    // Test creating rationals of value zero.
    auto zero = rational(0, 8);
    assert(zero == 0);
    assert(zero == rational(0, 16));
    assert(zero.numerator == 0);
    assert(zero.denominator == 1);
    auto one = zero + 1;
    one -= one;
    assert(one == zero);

    // Test integerPart, fraction part.
    auto intFract = rational(5, 4);
    assert(intFract.integerPart == 1);
    assert(intFract.fractionPart == rational(1, 4));
    assert(cast(long) intFract == 1);

    // Test whether CTFE works for primitive types.  Doesn't work yet.
    version(none)
    {
        enum myRational = (((rational(1, 2) + rational(1, 4)) * 2 - rational(1, 4))
                           / 2 + 1 * rational(1, 2) - 1) / rational(2, 5);
        import std.stdio;
        writeln(myRational);
        static assert(myRational == rational(-15, 32));
    }
}

/** Convert a floating point number to a $(D Rational) based on integer type $(D
 * SomeIntegral).  Allows an error tolerance of $(D epsilon).  (Default $(D epsilon) =
 * 1e-8.)
 *
 * $(D epsilon) must be greater than 1.0L / long.max.
 *
 * Throws:  Exception on infinities, NaNs, numbers with absolute value
 * larger than long.max and epsilons smaller than 1.0L / long.max.
 *
 * Examples:
 * ---
 * // Prints "22/7".
 * writeln(toRational!int(PI, 1e-1));
 * ---
 */
Rational!(SomeIntegral) toRational(SomeIntegral)(real floatNum, real epsilon = 1e-8)
{
    import std.math: isNaN;
    import std.exception : enforce;
    enforce(floatNum != real.infinity &&
            floatNum != -real.infinity &&
            !isNaN(floatNum),
            "Can't convert NaNs and infinities to rational.");
    enforce(floatNum < long.max &&
            floatNum > -long.max,
            "Rational conversions of very large numbers not yet implemented.");
    enforce(1.0L / epsilon < long.max,
            "Can't handle very small epsilons < long.max in toRational.");

    /* Handle this as a special case to make the rest of the code less
     * complicated:
     */
    if (abs(floatNum) < epsilon)
    {
        Rational!SomeIntegral ret;
        ret._num = 0;
        ret._den = 1;
        return ret;
    }

    return toRationalImpl!(SomeIntegral)(floatNum, epsilon);
}

private Rational!SomeIntegral toRationalImpl(SomeIntegral)(real floatNum, real epsilon)
{
    import std.conv : roundTo;
    import std.traits : isIntegral;

    real actualEpsilon;
    Rational!SomeIntegral ret;

    if (abs(floatNum) < 1)
    {
        real invFloatNum = 1.0L / floatNum;
        long intPart = roundTo!long(invFloatNum);
        actualEpsilon = floatNum - 1.0L / intPart;

        static if (isIntegral!(SomeIntegral))
        {
            ret._den = cast(SomeIntegral) intPart;
            ret._num = cast(SomeIntegral) 1;
        }
        else
        {
            ret._den = intPart;
            ret._num = 1;
        }
    }
    else
    {
        long intPart = roundTo!long(floatNum);
        actualEpsilon = floatNum - intPart;

        static if (isIntegral!(SomeIntegral))
        {
            ret._den = cast(SomeIntegral) 1;
            ret._num = cast(SomeIntegral) intPart;
        }
        else
        {
            ret._den = 1;
            ret._num = intPart;
        }
    }

    if (abs(actualEpsilon) <= epsilon)
        return ret;

    // Else get results from downstream recursions, add them to this result.
    return ret + toRationalImpl!(SomeIntegral)(actualEpsilon, epsilon);
}

unittest
{
    import std.bigint : BigInt;
    import std.math : PI, E;
    // Start with simple cases.
    assert(toRational!int(0.5) == rational(1, 2));
    assert(toRational!BigInt(0.333333333333333L) == rational(BigInt(1), BigInt(3)));
    assert(toRational!int(2.470588235294118) == rational(cast(int) 42, cast(int) 17));
    assert(toRational!long(2.007874015748032) == rational(255L, 127L));
    assert(toRational!int( 3.0L / 7.0L) == rational(3, 7));
    assert(toRational!int( 7.0L / 3.0L) == rational(7, 3));

    // Now for some fun.
    real myEpsilon = 1e-8;
    auto piRational = toRational!long(PI, myEpsilon);
    assert(abs(cast(real) piRational - PI) < myEpsilon);

    auto eRational = toRational!long(E, myEpsilon);
    assert(abs(cast(real) eRational - E) < myEpsilon);
}

/** Find the Greatest Common Factor (GCF), aka Greatest Common Divisor (GCD), of
 * $(D m) and $(D n).
 */
CommonInteger!(I1, I2) gcf(I1, I2)(I1 m, I2 n)
if (isIntegerLike!I1 &&
    isIntegerLike!I2)
{
    static if (is(I1 == const) || is(I1 == immutable) ||
               is(I2 == const) || is(I2 == immutable))
    {
        // Doesn't work with immutable(BigInt).
        import core.internal.traits : Unqual;
        return gcf!(Unqual!I1,
                    Unqual!I2)(m, n);
    }
    else
    {
        typeof(return) a = abs(m);
        typeof(return) b = abs(n);

        while (b)
        {
            auto t = b;
            b = a % b;
            a = t;
        }

        return a;
    }
}

pure unittest
{
    import std.bigint : BigInt;
    assert(gcf(0, 0) == 0);
    assert(gcf(0, 1) == 1);
    assert(gcf(999, 0) == 999);
    assert(gcf(to!(immutable(int))(8), to!(const(int))(12)) == 4);

    // Values from the Maxima computer algebra system.
    assert(gcf(BigInt(314_156_535UL), BigInt(27_182_818_284UL)) == BigInt(3));
    assert(gcf(8675309, 362436) == 1);
    assert(gcf(BigInt("8589934596"), BigInt("295147905179352825852")) == 12);
}

/// Find the Least Common Multiple (LCM) of $(D a) and $(D b).
CommonInteger!(I1, I2) lcm(I1, I2)(const scope I1 a, const scope I2 b)
if (isIntegerLike!I1 &&
    isIntegerLike!I2)
{
    const n1 = abs(a);
    const n2 = abs(b);
    if (n1 == n2)
        return n1;
    return (n1 / gcf(n1, n2)) * n2;
}

/// Returns the largest integer less than or equal to $(D r).
SomeIntegral floor(SomeIntegral)(const scope Rational!SomeIntegral r)
{
    SomeIntegral intPart = r.integerPart;
    if (r > 0 || intPart == r)
        return intPart;
    else
    {
        intPart -= 1;
        return intPart;
    }
}

@safe pure nothrow @nogc unittest
{
    assert(floor(rational(1, 2)) == 0);
    assert(floor(rational(-1, 2)) == -1);
    assert(floor(rational(2)) == 2);
    assert(floor(rational(-2)) == -2);
    assert(floor(rational(-1, 2)) == -1);
}

/// Returns the smallest integer greater than or equal to $(D r).
SomeIntegral ceil(SomeIntegral)(const scope Rational!SomeIntegral r)
{
    SomeIntegral intPart = r.integerPart;
    if (intPart == r || r < 0)
        return intPart;
    else
    {
        intPart += 1;
        return intPart;
    }
}

@safe pure nothrow @nogc unittest
{
    assert(ceil(rational(1, 2)) == 1);
    assert(ceil(rational(0)) == 0);
    assert(ceil(rational(-1, 2)) == 0);
    assert(ceil(rational(1)) == 1);
    assert(ceil(rational(-2)) == -2);
}

/** Round $(D r) to the nearest integer.  If the fractional part is exactly 1/2,
 * $(D r) will be rounded such that the absolute value is increased by rounding.
 */
SomeIntegral round(SomeIntegral)(const scope Rational!SomeIntegral r)
{
    auto intPart = r.integerPart;
    const fractPart = r.fractionPart;

    bool added;
    if (fractPart >= rational(1, 2))
    {
        added = true;
        intPart += 1;
    }

    import std.traits : isUnsigned;

    static if (!isUnsigned!SomeIntegral)
    {
        if (!added && fractPart <= rational(-1, 2))
            intPart -= 1;
    }

    return intPart;
}

@safe pure nothrow @nogc unittest
{
    assert(round(rational(1, 3)) == 0);
    assert(round(rational(7, 2)) == 4);
    assert(round(rational(-3, 4)) == -1);
    assert(round(rational(8U, 15U)) == 1);
}