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

Posted 2023-02-22

【中文标题】计算双精度分子和分母的最准确方法

【英文标题】：The most accurate way to calculate numerator and denominator of a double

【发布时间】：2014-01-18 23:38:16

【问题描述】：

我已经实现了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_”开头的变量是成员变量。

谢谢

【问题讨论】：

将浮点数分解为符号、尾数和指数，然后你就有足够的信息来准确地表示它了。

我不会冒着在不同编译器（或 CPU 架构）上使用不同 FP 表示的风险吗？

@barakmanos 您是为 IBM 大型机还是 Cray 超级计算机编程？如果您不是，那么请忘记非 IEEE 754 表示。

@Pascal Cuoq：谢谢。如您所见，我不依赖任何 FP 表示。

m_numerator 和 m_denominator 有哪些类型？

【参考方案1】：

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

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

【讨论】：

谢谢；问与上述 cmets 之一相同的问题：我不会冒着在不同编译器（或 CPU 架构）上出现不同 FP 表示的风险吗？或者frexp在任何情况下都能保证相应地实现？

有效数字（不是尾数；有效数字是线性的，但尾数是对数的）有 53 个有效位，而不是 52。52 只是显式编码的位数。

@EricPostpischil：谢谢，已更新（出于某种原因，我认为指数是 12，而不是 11 位）。

@EricPostpischil：我已经更改了术语，但请注意，wikipedia 将“尾数”作为有效数字的同义词。

@JanHudec：***页面上说“尾数”是 IEEE 浮点委员会和其他人不鼓励的。 IEEE 754 标准根本不使用“尾数”。我们中的一些人记得何时必须在对数表中查找尾数。当您在纸上手工进行浮点数学运算时，尾数和有效数字之间的区别将非常明显。

【参考方案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 功能可能不起作用。

【讨论】：

所以？现在我被另一个double 值困住了……还是这个值总是0.5？

@barakmanos：有效数字以double 的形式返回，但很容易转换为其他形式，因为它的比例是已知的。如果您有 64 位整数可用，只需乘以 2**53 并转换为整数。

事实证明这比我想象的要复杂得多，我开始实现自己的解决方案，但它很难看，我不喜欢它。如果我不得不这样做，我会使用 Boost 有理数库：boost.org/doc/libs/1_55_0/libs/rational/rational.html

@vmrob,@Eric Postpischil：请记住，我的主要目标是计算 double 值的分子和分母。我仍然看不到上面的解决方案如何帮助我实现这一目标（更不用说，以比我的问题中实现的更准确的方式）。

@barakmanos 您的方法和此答案中描述的方法都适用于 IEEE 754 双精度（您的方法仅在它完全有效时：它在 (unsigned int)value 溢出时调用未定义的行为）。在源代码中写0.1得到的双精度数是3602879701896397乘以2的幂。不涉及“证据”。