计算二重分子和分母的最准确方法

The most accurate way to calculate numerator and denominator of a double

本文关键字：方法 二重 计算      更新时间：2023-10-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_"开头的变量是成员变量

感谢

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

只需将有效位相乘，即可获得小数点前的所有位，并设置适当的分母。别忘了，有效位被标准化为0.5到1之间(我认为1到2之间更自然，但无论如何)，并且它有53个IEEE双有效位(没有实际使用的平台会使用不同的浮点格式)。

我对实际计算的数学不是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函数不支持doubles，那么您将限制自己可以作为doubles导入的值。尝试任何数字>INT64_MAX，GCD功能可能不起作用。

相关文章：