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

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

我对实际计算的数学没有百分之百的信心，只是因为我还没有真正检查它，但我认为下面的方法消除了使用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函数可能无法工作。