在 C++ 中使用 int 分子和分母对 double 进行有理逼近答案

【问题标题】：Rational approximation of double using int numerator and denominator in C++在 C++ 中使用 int 分子和分母对 double 进行有理逼近
【发布时间】：2022-01-05 14:16:03
【问题描述】：

真实世界的第三方 API 采用 fraction 类型的参数，它是 struct 的 int 分子和分母。我需要传递的值对我来说是一个转换为双精度的十进制字符串。

可能值的范围是，比如 10K 到 300M，但如果小数点后有小数部分，则很重要。

我这里有两种近似方法的代码，一种使用扩展欧几里得算法，另一种是蛮力。对于给定的double，这两种方法都使用int 类型找到有理逼近。

蛮力当然是两者中更准确的，并且当转换的数字很大时实际上更快。我的问题是，我能对使用欧几里得算法的近似质量说些什么吗？更正式地说，我可以对使用欧几里得算法的近似值与蛮力算法的近似值（我认为是最优的）进行限制吗？

绑定示例：
如果最优逼近的误差是r，那么欧几里得算法逼近会产生小于2*r的误差。
（我并不是说这是界限，我当然也无法证明这一点，这只是一个良好界限的示例）。

这是一个测试程序的代码：

#include <iostream>
#include <iomanip>
#include <cmath>
#include <limits>
#include <chrono>
#include <random>

// extended euclidian algorithm
// finds the coefficients that produce the gcd
// in u, we store m,n the coefficients that produce m*a - n*b == gcd.
// in v, we store m,n the coefficients that produce m*a - n*b == 0.
// breaks early if the coefficients become larger than INT_MAX
int gcd_e(uint64_t a, int b, int u[2], int v[2])
{
    auto w = lldiv(a, b);

    // u[0] * a' - u[1] * b' == a
    // v[0] * a' - v[1] * b' == b

    // a - w.quot * b == w.rem
    // (u[0] * a' - u[1] * b') - w.quot * (v[0] * a' - v[1] * b') == w.rem
    // (u[0] - w.quot * v[0]) * a' - u[1] * b' + w.quot * v[1] * b' == w.rem
    // (u[0] - w.quot * v[0]) * a' + (w.quot * v[1] - u[1]) * b' == w.rem
    // (u[0] - w.quot * v[0]) * a' - (u[1] - w.quot * v[1]) * b' == w.rem


    auto m = u[0] - w.quot * v[0];
    auto n = u[1] - w.quot * v[1];
    u[0] = v[0];
    u[1] = v[1];

    constexpr auto L = std::numeric_limits<int>::max();
    if (m > L || n > L)
        throw 0;  // break early
    if (m < -L || n < -L)
        throw 0;  // break early

    v[0] = int(m);
    v[1] = int(n);

    if (w.rem == 0)
        return b;

    return gcd_e(b, int(w.rem), u, v);
}


inline double helper_pre(double d, bool* negative, bool* inverse)
{
    bool v = (d < 0);
    *negative = v;
    if (v)
        d = -d;

    v = (d < 1);
    *inverse = v;
    if (v)
        d = 1 / d;

    return d;
}

inline void helper_post(int* m, int* n, bool negative, bool inverse)
{
    if (inverse)
        std::swap(*n, *m);

    if (negative)
        *n = -(*n);
}

// gets a rational approximation for double d
// numerator is stored in n
// denominator is stored in m
void approx(double d, int* n, int *m)
{
    int u[] = { 1, 0 };  // 1*a - 0*b == a
    int v[] = { 0, -1 }; // 0*a - (-1)*b == b

    bool negative, inverse;
    d = helper_pre(d, &negative, &inverse);

    constexpr int q = 1 << 30;
    auto round_d = std::round(d);
    if (d == round_d)
    {
        // nothing to do, it's an integer.
        v[1] = int(d);
        v[0] = 1;
    }
    else try
    {
        uint64_t k = uint64_t(std::round(d*q));
        gcd_e(k, q, u, v);
    }
    catch (...)
    {
        // OK if we got here.
        // int limits
    }

    // get the approximate numerator and denominator
    auto nn = v[1];
    auto mm = v[0];

    // make them positive
    if (mm < 0)
    {
        mm = -mm;
        nn = -nn;
    }

    helper_post(&mm, &nn, negative, inverse);

    *m = mm;
    *n = nn;
}


// helper to test a denominator
// returns the magnitude of the error
double helper_rattest(double x, int tryDenom, int* numerator)
{
    double r = x * tryDenom;
    double rr = std::round(r);
    auto num = int(rr);
    auto err = std::abs(r - rr) / tryDenom;
    *numerator = num;
    return err;
}

// helper to reduce the rational number
int gcd(int a, int b)
{
    auto c = a % b;
    if (c == 0)
        return b;
    return gcd(b, int(c));
}

// gets a rational approximation for double d
// numerator is stored in n
// denominator is stored in m
// uses brute force by scanning denominator range
void approx_brute(double d, int* n, int* m)
{
    bool negative, inverse;
    d = helper_pre(d, &negative, &inverse);

    int upto = int(std::numeric_limits<int>::max() / d);
    int bestNumerator;
    int bestDenominator = 1;
    auto bestErr = helper_rattest(d, 1, &bestNumerator);
    for (int kk = 2; kk < upto; ++kk)
    {
        int n;
        auto e = helper_rattest(d, kk, &n);
        if (e < bestErr)
        {
            bestErr = e;
            bestNumerator = n;
            bestDenominator = kk;
        }
        if (bestErr == 0)
            break;
    }

    // reduce, just in case
    auto g = gcd(bestNumerator, bestDenominator);
    bestNumerator /= g;
    bestDenominator /= g;

    helper_post(&bestDenominator, &bestNumerator, negative, inverse);

    *n = bestNumerator;
    *m = bestDenominator;
}

int main()
{
    int n, m;

    auto re = std::default_random_engine();
    std::random_device rd;
    re.seed(rd());

    for (auto& u : {
        std::uniform_real_distribution<double>(10000,    15000),
        std::uniform_real_distribution<double>(100000,   150000),
        std::uniform_real_distribution<double>(200000,   250000),
        std::uniform_real_distribution<double>(400000,   450000),
        std::uniform_real_distribution<double>(800000,   850000),
        std::uniform_real_distribution<double>(1000000,  1500000),
        std::uniform_real_distribution<double>(2000000,  2500000),
        std::uniform_real_distribution<double>(4000000,  4500000),
        std::uniform_real_distribution<double>(8000000,  8500000),
        std::uniform_real_distribution<double>(10000000, 15000000)
        })
    {
        auto dd = u(re);
        std::cout << "approx: " << std::setprecision(14) << dd << std::endl;

        auto before = std::chrono::steady_clock::now();
        approx_brute(dd, &n, &m);
        auto after = std::chrono::steady_clock::now();
        std::cout << n << " / " << m << "  dur: " << (after - before).count() << std::endl;
        before = std::chrono::steady_clock::now();
        approx(dd, &n, &m);
        after = std::chrono::steady_clock::now();
        std::cout << n << " / " << m << "  dur: " << (after - before).count()
            << std::endl
            << std::endl;
    }
}

这是一些示例输出：

approx: 13581.807792679
374722077 / 27590  dur: 3131300
374722077 / 27590  dur: 15000

approx: 103190.31976517
263651267 / 2555  dur: 418700
263651267 / 2555  dur: 6300

approx: 223753.78683426
1726707973 / 7717  dur: 190100
1726707973 / 7717  dur: 5800

approx: 416934.79214075
1941665327 / 4657  dur: 102100
403175944 / 967  dur: 5700

approx: 824300.61241502
1088901109 / 1321  dur: 51900
1088901109 / 1321  dur: 5900

approx: 1077460.29557
1483662827 / 1377  dur: 39600
1483662827 / 1377  dur: 5600

approx: 2414781.364653
1079407270 / 447  dur: 17900
1079407270 / 447  dur: 7300

approx: 4189869.294816
1776504581 / 424  dur: 10600
1051657193 / 251  dur: 9900

approx: 8330270.2432111
308219999 / 37  dur: 5400
308219999 / 37  dur: 10300

approx: 11809264.006453
1830435921 / 155  dur: 4000
1830435921 / 155  dur: 10500

【问题讨论】：

因为计算机只能容纳有限数量的数字，所以任何双精度数在技术上都已经是 x / (2^y) 形式的有理数。
谢谢@unddoch，这有助于找到两个适合 int 变量的数字并且其比率接近双精度数中的数字吗？
只是总结一下管道：你从十进制字符串开始，然后用二进制双精度近似它，然后用有理数近似？
是的，因为 x 和 2^y 是两个整数，其比率是双精度数中的数字。（除非 y 是负数，但双精度数是整数）。您仍然会遇到范围问题，但无论如何您都不能用整数近似 1e30 或 NaN。
@unddoch，非暴力算法几乎完全建立在此之上，并解决了范围问题。然而，它显然是次优的，问题主要是一个幻想“它有多糟糕？”

标签： c++ algorithm floating-point approximation rational-number

【解决方案1】：

感谢所有评论并提请我注意连分数概念的人。根据this paper by (William F. Hammond) 欧几里得算法和连分数法是等价的。

次优结果是由于分子和分母都受到约束，因此如果非暴力算法仅产生“收敛”，则意味着它忽略了第一个收敛之间的分母范围违反约束和它之前的约束。

返回的收敛后的分母和后面的分母可能近似于后一个收敛，后续收敛之间的差异可以表示为：

所以我想这将是蛮力和欧几里得算法之间差异的界限。它们之间的误差比率几乎可以是任何值。
（可以轻松找到错误率超过 100 的示例）

我希望我能正确阅读所有内容。我没有这方面的权威。

【讨论】：

一系列分母的最佳有理近似可以被证明是 c-frac 展开中（包括）收敛之间的semiconvergent 项。
另外 - 使用无符号输入的扩展 gcd 的一个非常简洁的实现，具有边界证明等，是 here。
如果有任何帮助，您可以在此处查看 Python 的标准库代码以找到最佳有理逼近：github.com/python/cpython/blob/…。正如 Brett Hale 所说，它使用连分数，但确保检查半收敛和收敛。