root/trunk/rational/rational.d

/**A rational number implementation
 * By:  David Simcha
 *
 * License:
 * Boost Software License - Version 1.0 - August 17th, 2003
 *
 * Permission is hereby granted, free of charge, to any person or organization
 * obtaining a copy of the software and accompanying documentation covered by
 * this license (the "Software") to use, reproduce, display, distribute,
 * execute, and transmit the Software, and to prepare derivative works of the
 * Software, and to permit third-parties to whom the Software is furnished to
 * do so, all subject to the following:
 *
 * The copyright notices in the Software and this entire statement, including
 * the above license grant, this restriction and the following disclaimer,
 * must be included in all copies of the Software, in whole or in part, and
 * all derivative works of the Software, unless such copies or derivative
 * works are solely in the form of machine-executable object code generated by
 * a source language processor.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
 * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
 * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
 * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 */
module rational;


import std.algorithm, std.stdio, std.bigint, std.conv, std.math, std.exception,
       std.conv, std.traits;

alias std.math.abs abs;  // Allow cross-module overloading.

private template isIntegerLike(T) {
    enum bool isIntegerLike = is(typeof({
        T num;
        num = 2;
        num <<= 1;
        num >>= 1;
        num += num;
        num *= num;
        num /= num;
        num -= num;
        num %= 2;
        num %= num;
        bool foo = num < 2;
        bool bar = num == 2;

        return num;
    }));
}

unittest {
    static assert(isIntegerLike!BigInt);
    static assert(isIntegerLike!int);
    static assert(isIntegerLike!byte);
    static assert(!isIntegerLike!real);
}

private template isRational(T) {
    enum bool isRational =
        is(typeof(T.init.denom)) && is(typeof(T.init.num));
}

private template CommonRational(R1, R2) {
    static if(isRational!R1) {
        alias CommonRational!(typeof(R1.num), R2) CommonRational;
    } else static if(isRational!R2) {
        alias CommonRational!(R1, typeof(R2.num)) CommonRational;
    } else static if(is(CommonInteger!(R1, R2))) {
        alias Rational!(CommonInteger!(R1, R2)) CommonRational;
    }
}

/**
Returns $(D true) iff a value of type $(D U) can be assigned to a variable of
type $(D T).

Examples:
---
static assert(isAssignable!(long, int));
static assert(!isAssignable!(int, long));
static assert(isAssignable!(const(char)[], string));
static assert(!isAssignable!(string, char[]));
---
*/
template isAssignable(T, U) {
    enum bool isAssignable = is(typeof({
        T t;
        U u;
        t = u;
        return t;
    }));
}

unittest {
    static assert(isAssignable!(long, int));
    static assert(!isAssignable!(int, long));
    static assert(isAssignable!(const(char)[], string));
    static assert(!isAssignable!(string, char[]));
}

/**Returns a common integral type between I1, I2.  This is defined as the type
 * returned by I1.init * I2.init.
 */
template CommonInteger(I1, I2) if(isIntegerLike!I1 && isIntegerLike!I2) {
    alias typeof(I1.init * I2.init) CommonInteger;
}

unittest {
    static assert(is(CommonInteger!(BigInt, int) == BigInt));
    static assert(is(CommonInteger!(byte, int) == int));
}

/**Implements rational numbers on top of whatever integer type is specified
 * by the user.  If you use a fixed-size int, the onus is on you to make
 * sure nothing overflows.  Neither the denominator nor the numerator may
 * be bigger than the maximum value of the underlying integer.
 *
 * If you use an arbitrary precision type, it must have the relevant operators
 * overloaded and have value semantics.
 *
 * 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(f1);
 *
 * // Print result in decimal form.  Prints:
 * // "35.4934"
 * writeln(cast(real) result);
 * ---
 */
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.numerator = i1;
        ret.denominator = i2;
    }
    ret.simplify();
    return ret;
}

/**Overload for creating rational that initially has an integer value.*/
Rational!(I) rational(I)(I val)
if(isIntegerLike!I) {
     return rational(val, 1);
}

///
struct Rational(Int) if(isIntegerLike!Int) {
public:

// ----------------Multiplication operators----------------------------------
    auto opBinary(string op, Rhs)(Rhs rhs)
    if(op == "*" && is(CommonRational!(Int, Rhs)) && isRational!Rhs) {
        auto ret = CommonRational!(Int, Rhs)(this.num, this.denom);
        return ret *= rhs;
    }

    auto opBinary(string op, Rhs)(Rhs rhs)
    if(op == "*" && is(CommonRational!(Int, Rhs)) && isIntegerLike!Rhs) {
        auto ret = this;
        return ret *= rhs;
    }

    auto opBinaryRight(string op, Rhs)(Rhs rhs)
    if(op == "*" && is(CommonRational!(Int, Rhs)) && isIntegerLike!Rhs) {
        return opBinary!(op, Rhs)(rhs);
    }

    typeof(this) opOpAssign(string op, Rhs)(Rhs rhs)
    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.numerator, rhs.denominator);
        this.numerator /= divisor;
        rhs.denominator /= divisor;

        divisor = gcf(this.denominator, rhs.numerator);
        this.denominator /= divisor;
        rhs.numerator /= divisor;

        this.numerator *= rhs.numerator;
        this.denominator *= rhs.denominator;

        // Don't need to simplify.  Already cancelled common factors before
        // multiplying.
        fixSigns();
        return this;
    }

    typeof(this) opOpAssign(string op, Rhs)(Rhs rhs)
    if(op == "*" && isIntegerLike!Rhs) {
        auto divisor = gcf(this.denominator, rhs);
        this.denominator /= divisor;
        rhs /= divisor;
        this.numerator *= rhs;

        // Don't need to simplify.  Already cancelled common factors before
        // multiplying.
        fixSigns();
        return this;
    }

// --------------------Division operators--------------------------------------

    auto opBinary(string op, Rhs)(Rhs rhs)
    if(op == "/" && is(CommonRational!(Int, Rhs)) && isRational!Rhs) {
        // Division = multiply by inverse.
        swap(rhs.numerator, rhs.denominator);
        return this *= rhs;
    }

    typeof(this) opBinary(string op, Rhs)(Rhs rhs)
    if(op == "/" && is(CommonRational!(Int, Rhs)) && isIntegerLike!(Rhs)) {
        auto ret = CommonRational!(Int, Rhs)(this.num, this.denom);
        return ret /= rhs;
    }

    typeof(this) opBinaryRight(string op, Rhs)(Rhs rhs)
    if(op == "/" && is(CommonRational!(Int, Rhs)) && isIntegerLike!Rhs) {
        auto ret = CommonRational!(Int, Rhs)(this.denom, this.num);
        return ret *= rhs;
    }

    typeof(this) opOpAssign(string op, Rhs)(Rhs rhs)
    if(op == "/" && isIntegerLike!Rhs) {
        auto divisor = gcf(this.numerator, rhs);
        this.numerator /= divisor;
        rhs /= divisor;
        this.denominator *= rhs;

        // Don't need to simplify.  Already cancelled common factors before
        // multiplying.
        fixSigns();
        return this;
    }

    typeof(this) opOpAssign(string op, Rhs)(Rhs rhs)
    if(op == "/" && isRational!Rhs) {
        // Division = multiply by inverse.
        swap(rhs.numerator, rhs.denominator);
        return this *= rhs;
    }

// ---------------------Addition operators-------------------------------------
    auto opBinary(string op, Rhs)(Rhs rhs)
    if(op == "+" && (isRational!Rhs || isIntegerLike!Rhs)) {
        auto ret = CommonRational!(typeof(this), Rhs)(this.num, this.denom);
        return ret += rhs;
    }

    auto opBinaryRight(string op, Rhs)(Rhs rhs)
    if(op == "+" && is(CommonRational!(Int, Rhs)) && isIntegerLike!Rhs) {
        return opBinary!(op, Rhs)(rhs);
    }

    typeof(this) opOpAssign(string op, Rhs)(Rhs rhs)
    if(op == "+" && isRational!Rhs) {
        if(this.denominator == rhs.denominator) {
            this.numerator += rhs.numerator;
            simplify();
            return this;
        }

        Int commonDenom = lcm(this.denominator, rhs.denominator);
        this.numerator *= commonDenom / this.denominator;
        this.numerator += (commonDenom / rhs.denominator) * rhs.numerator;
        this.denominator = commonDenom;

        simplify();
        return this;
    }

    typeof(this) opOpAssign(string op, Rhs)(Rhs rhs)
    if(op == "+" && isIntegerLike!Rhs) {
        this.numerator += rhs * this.denominator;

        simplify();
        return this;
    }

// -----------------------Subtraction operators-------------------------------

    auto opBinary(string op, Rhs)(Rhs rhs)
    if(op == "-" && is(CommonRational!(Int, Rhs))) {
        auto ret = CommonRational!(typeof(this), Rhs)(this.num, this.denom);
        return ret -= rhs;
    }

    typeof(this) opOpAssign(string op, Rhs)(Rhs rhs)
    if(op == "-" && isRational!Rhs) {
        if(this.denominator == rhs.denominator) {
            this.numerator -= rhs.numerator;
            simplify();
            return this;
        }

        auto commonDenom = lcm(this.denominator, rhs.denominator);
        this.numerator *= commonDenom / this.denominator;
        this.numerator -= (commonDenom / rhs.denominator) * rhs.numerator;
        this.denominator = commonDenom;

        simplify();
        return this;
    }

    typeof(this) opOpAssign(string op, Rhs)(Rhs rhs)
    if(op == "-" && isIntegerLike!Rhs) {
        this.numerator -= rhs * this.denominator;

        simplify();
        return this;
    }

    typeof(this) opBinaryRight(string op, Rhs)(Rhs rhs)
    if(op == "-" && is(CommonInteger!(Int, Rhs)) && isIntegerLike!Rhs) {
        typeof(this) ret;
        ret.denominator = this.denominator;
        ret.numerator = (rhs * this.denominator) - this.numerator;

        simplify();
        return ret;
    }

// ----------------------Unary operators---------------------------------------

    typeof(this) opUnary(string op)() if(op == "-" || op == "+") {
        mixin("return typeof(this)(" ~ op ~ "numerator, denominator);");
    }

// ---------------------Exponentiation operator---------------------------------

    // Can only handle integer powers if the result has to also be
    // rational.

    typeof(this) opOpAssign(string op, Rhs)(Rhs rhs)
    if(op == "^^" && isIntegerLike!Rhs) {
        if(rhs < 0) {
            this.invert();
            rhs *= -1;
        }

        /* 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.
         */
         numerator ^^= rhs;
         denominator ^^= rhs;
         return this;
    }

    auto opBinary(string op, Rhs)(Rhs rhs)
    if(op == "^^" && isIntegerLike!Rhs && is(CommonRational!(Int, Rhs))) {
        auto ret = CommonRational!(Int, Rhs)(this.num, this.denom);
        ret ^^= rhs;
        return ret;
    }

// ---------------------Assignment operators------------------------------------
    typeof(this) opAssign(Rhs)(Rhs rhs)
    if(isIntegerLike!Rhs && isAssignable!(Int, Rhs)) {
        this.numerator = rhs;
        this.denominator = 1;
        return this;
    }

    typeof(this) opAssign(Rhs)(Rhs rhs)
    if(isRational!Rhs && isAssignable!(Int, typeof(Rhs.num))) {
        this.numerator = rhs.num;
        this.denominator = rhs.denom;
        return this;
    }

// --------------------Comparison/Equality Operators---------------------------

    bool opEquals(Rhs)(Rhs rhs) if(isRational!Rhs || isIntegerLike!Rhs) {
        static if(isRational!Rhs) {
            return rhs.numerator == this.numerator &&
                rhs.denominator == this.denominator;
        } else {
            static assert(isIntegerLike!Rhs);
            return rhs == this.numerator && this.denominator == 1;
        }
    }

    int opCmp(Rhs)(Rhs rhs) if(isRational!Rhs) {
        if( opEquals(rhs)) {
            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.numerator < 0 && rhs.numerator > 0) {
            return -1;
        } else if(this.numerator > 0 && rhs.numerator < 0) {
            return 1;
        } else if(this.numerator >= rhs.numerator &&
            this.denominator <= rhs.denominator) {
            // We've already ruled out equality, so this must be > rhs.
            return 1;
        } else if(rhs.numerator >= this.numerator &&
            rhs.denominator <= this.denominator) {
            return -1;
        }

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

        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)(Rhs rhs) if(isIntegerLike!Rhs) {
        if( opEquals(rhs)) {
            return 0;
        }

        // Again, check the obvious cases first.
        if(rhs >= this.numerator) {
            return -1;
        }

        rhs *= this.denominator;
        if(rhs > this.numerator) {
            return -1;
        } else if(rhs < this.numerator) {
            return 1;
        }

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

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

    /**Fast inversion, equivalent to 1 / rational.*/
    typeof(this) invert() {
        swap(numerator, denominator);
        return this;
    }

    /**Convert to floating point representation.*/
    F opCast(F)() if(isFloatingPoint!F) {
        // Do everything in real precision, then convert to F at the end.

        static if(isIntegral!(Int)) {
            return cast(real) numerator / denominator;
        } else {
            auto temp = this;
            real expon = 1.0;
            real ans = 0;
            byte sign = 1;
            if(temp.numerator < 0) {
                temp.numerator *= -1;
                sign = -1;
            }

            while(temp.numerator > 0) {
                while(temp.numerator < temp.denominator) {

                    assert(temp.denominator > 0);

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

                    expon *= 0.5;

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

                }

                auto intPart = temp.numerator / temp.denominator;

                static if(is(Int == std.bigint.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.numerator -= intPart * temp.denominator;
            }

            return ans * sign;
        }
    }

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

    /**Returns the numerator.*/
    @property Int num() {
        return numerator;
    }

    /**Returns the denominator.*/
    @property Int denom() {
        return denominator;
    }

    /**Returns the integer part of this rational, with any remainder truncated.
     */
    @property Int integerPart() {
        return num / denom;
    }

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

    string toString() {
        static if(is(Int == std.bigint.BigInt)) {
            // Special case it for now.  This should be fixed later.
            return toDecimalString(numerator) ~ " / " ~
                toDecimalString(denominator);
        } else {
            return to!string(numerator) ~ " / " ~ to!string(denominator);
        }
    }

private :
    Int numerator;
    Int denominator;

    void simplify() {
        if(numerator == 0) {
            denominator = 1;
            return;
        }

        auto divisor = gcf(numerator, denominator);
        numerator /= divisor;
        denominator /= divisor;

        fixSigns();
    }

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

unittest {
    // 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);

    // Test multiplication.
    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(approxEqual(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.num == 0);
    assert(zero.denom == 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.
//    enum myRational = (((rational(1, 2) + rational(1, 4)) * 2 - rational(1, 4))
//        / 2 + 1 * rational(1, 2) - 1) / rational(2, 5);
//    writeln(myRational);
//    static assert(myRational == rational(-15, 32));
}

/**Convert a floating point number to a Rational based on integer type Int.
 * Allows an error tolerance of epsilon.  (Default epsilon = 1e-8.)
 *
 * 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!(Int) toRational(Int)(real floatNum, real epsilon = 1e-8) {
    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!Int ret;
        ret.numerator = 0;
        ret.denominator = 1;
        return ret;
    }

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

private Rational!Int toRationalImpl(Int)(real floatNum, real epsilon) {
    real actualEpsilon;
    Rational!Int ret;

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

        static if(isIntegral!(Int)) {
            ret.denominator = cast(Int) intPart;
            ret.numerator = cast(Int) 1;
        } else {
            ret.denominator = intPart;
            ret.numerator = 1;
        }
    } else {
        long intPart = roundTo!long(floatNum);
        actualEpsilon = floatNum - intPart;

        static if(isIntegral!(Int)) {
            ret.denominator = cast(Int) 1;
            ret.numerator = cast(Int) intPart;
        } else {
            ret.denominator = 1;
            ret.numerator = intPart;
        }
    }

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

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

unittest {
    // 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 of num1 and num2 using Euclid's Algorithm.*/
CommonInteger!(I1, I2) gcf(I1, I2)(I1 num1, I2 num2)
if(isIntegerLike!I1 && isIntegerLike!I2) {
    num1 = iAbs(num1);
    num2 = iAbs(num2);
    if(num2 > num1) {
        return gcf(num2, num1);
    } else if(num2 == num1) {
        typeof(return) ret = num1;
        return ret;
    }

    // Work around Bug 4742.
    static if(is(I1 == I2)) {
        auto remainder = num1 % num2;
    } else {
        typeof(return) workaround1 = num1;
        typeof(return) workaround2 = num2;
        auto remainder = workaround1 % workaround2;
    }

    if(remainder == 0) {
        typeof(return) ret = num2;
        return ret;
    } else {
        return gcf(num2, remainder);
    }
    assert(0);
}

unittest {
    // 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 of num1, num2.*/
CommonInteger!(I1, I2) lcm(I1, I2)(I1 num1, I2 num2)
if(isIntegerLike!I1 && isIntegerLike!I2) {
    num1 = iAbs(num1);
    num2 = iAbs(num2);
    if(num1 == num2) {
        return num1;
    }
    return (num1 / gcf(num1, num2)) * num2;
}

/**Absolute value function that should gracefully handle any reasonable
 * BigInt implementation.*/
Int iAbs(Int)(Int num1)
if(isIntegerLike!Int) {
    static if(isUnsigned!Int) {
        return num1;
    } else {
        // For some reason DMD insists that a byte multipled by -1 is an int
        // not a byte.
        return cast(Int) ((num1 < 0) ? -1 * num1 : num1);
    }
}

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

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).*/
Int ceil(Int)(Rational!Int r) {
    auto intPart = r.integerPart;
    if(intPart == r || r < 0) {
        return intPart;
    } else {
        intPart += 1;
        return intPart;
    }
}

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.
 */
Int round(Int)(Rational!Int r) {
    auto intPart = r.integerPart;
    auto fractPart = r.fractionPart;

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

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

    return intPart;
}

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

version(unittest) {
    void main() {}
}