2

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

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

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

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

一个界限的例子：
如果最优逼近的误差是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
4

1 回答 1

0

感谢所有评论并提请我注意连分数概念的人。根据(William F. Hammond) 的这篇论文， 欧几里得算法和连分数法之间存在等价性。 定理 5。

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

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

收敛改进

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

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

于 2021-11-29T16:44:08.187 回答