计算双精度分子和分母的最准确方法

Question

我已经实现了class NaturalNum 来表示“无限”大小（最大 4GB）的自然数。

我还实现了class RationalNum 以无限准确地表示有理数。它存储有理数的分子和分母，两者都是NaturalNum 实例，并在执行用户发出的任何算术运算时依赖它们。

精度“下降一定程度”的唯一地方是在打印时，因为小数点（或非小数点）后出现的位数有限制（由用户提供）。

我的问题涉及class RationalNum 的构造函数之一。即，构造函数接受double 值，并计算相应的分子和分母。

我的代码如下，我想知道是否有人看到了更准确的计算方法：

RationalNum::RationalNum(double value)
{
    if (value == value+1)
        throw "Infinite Value";
    if (value != value)
        throw "Undefined Value";

    m_sign        = false;
    m_numerator   = 0;
    m_denominator = 1;

    if (value < 0)
    {
        m_sign = true;
        value  = -value;
    }

    // Here is the actual computation
    while (value > 0)
    {
        unsigned int floor = (unsigned int)value;
        value         -= floor;
        m_numerator   += floor;
        value         *= 2;
        m_numerator   *= 2;
        m_denominator *= 2;
    }

    NaturalNum gcd = GCD(m_numerator,m_denominator);
    m_numerator   /= gcd;
    m_denominator /= gcd;
}

注意：以“m_”开头的变量是成员变量。

谢谢

Answer 1

标准库包含一个获取有效数和指数的函数frexp。

只需乘以有效数字即可得到小数点前的所有位并设置适当的分母。只是不要忘记有效数字被归一化为 0.5 和 1 之间（我认为 1 和 2 之间更自然，但无论如何）并且它有 53 个 IEEE double 有效位（没有实际使用的平台会使用不同的浮点数点格式）。

Answer 2

我不是 100% 对您对实际计算的数学有信心，只是因为我没有真正检查过它，但我认为下面的方法不需要使用 GCD 函数，这可能会带来一些不必要的运行时间。

这是我想出的课程。我还没有完全测试过它，但是我产生了几十亿个随机双打并且断言从未被触发，所以我对它的可用性有相当的信心，但我仍然会测试 INT64_MAX 周围的边缘情况还有一点。

如果我没记错的话，这个算法的运行时间复杂度与输入的比特大小成线性关系。

#include <iostream>
#include <cmath>
#include <cassert>
#include <limits>

class Real;

namespace std {
    inline bool isnan(const Real& r);
    inline bool isinf(const Real& r);
}

class Real {
public:
    Real(double val)
        : _val(val)
    {
        if (std::isnan(val)) { return; }
        if (std::isinf(val)) { return; }

        double d;

        if (modf(val, &d) == 0) {
            // already a whole number
            _num = val;
            _den = 1.0;
            return;
        }

        int exponent;
        double significand = frexp(val, &exponent); // val = significand * 2^exponent
        double numerator = val;
        double denominator = 1;

        // 0.5 <= significand < 1.0
        // significand is a fraction, multiply it by two until it's a whole number
        // subtract exponent appropriately to maintain val = significand * 2^exponent
        do {
            significand *= 2;
            --exponent;
            assert(std::ldexp(significand, exponent) == val);
        } while (modf(significand, &d) != 0);

        assert(exponent <= 0);  

        // significand is now a whole number
        _num = significand;
        _den = 1.0 / std::ldexp(1.0, exponent);

        assert(_val == _num / _den);
    }

    friend std::ostream& operator<<(std::ostream &os, const Real& rhs);
    friend bool std::isnan(const Real& r);
    friend bool std::isinf(const Real& r);

private:
    double _val = 0;
    double _num = 0;
    double _den = 0;
};

std::ostream& operator<<(std::ostream &os, const Real& rhs) {
    if (std::isnan(rhs) || std::isinf(rhs)) {
        return os << rhs._val;
    }
    if (rhs._den == 1.0) {
        return os << rhs._num;
    }
    return os << rhs._num << " / " << rhs._den;
}

namespace std {
    inline bool isnan(const Real& r) { return std::isnan(r._val); }
    inline bool isinf(const Real& r) { return std::isinf(r._val); }
}

#include <iomanip>

int main () {

    #define PRINT_REAL(num) \
        std::cout << std::setprecision(100) << #num << " = " << num << " = " << Real(num) << std::endl

    PRINT_REAL(1.5);
    PRINT_REAL(123.875);
    PRINT_REAL(0.125);

    // double precision issues
    PRINT_REAL(-10000000000000023.219238745);
    PRINT_REAL(-100000000000000000000000000000000000000000.5);

    return 0;
}

再看你的代码，你的无限值测试至少有问题。请注意以下程序：

#include <numeric>
#include <cassert>
#include <cmath>

int main() {
    {
        double d = std::numeric_limits<double>::max(); // about 1.7976931348623e+308

        assert(!std::isnan(d));
        assert(!std::isinf(d));

        // assert(d != d + 1); // fires
    }

    {
        double d = std::ldexp(1.0, 500); // 2 ^ 700
        assert(!std::isnan(d));
        assert(!std::isinf(d));

        // assert(d != d + 1); // fires
    }
}

除此之外，如果您的 GCD 函数不支持双精度，那么您将限制自己可以导入为双精度的值。尝试任何数字 > INT64_MAX，GCD 功能可能不起作用。