c - 如何使用 Arb 库获得更高精度的正弦波？

Question

我找到了 Arb library ，这should如果有足够的时间，能够计算非常高精度的正弦值。但是，我做不到。

尝试 sine example , 我可以获得预测的输出。

但是，当我试图通过将位数从 4096 增加到 32768 来提高精度时，我无法:

Using    64 bits, sin(x) = [+/- 2.67e+859]
Using   128 bits, sin(x) = [+/- 1.30e+840]
Using   256 bits, sin(x) = [+/- 3.60e+801]
Using   512 bits, sin(x) = [+/- 3.01e+724]
Using  1024 bits, sin(x) = [+/- 2.18e+570]
Using  2048 bits, sin(x) = [+/- 1.22e+262]
Using  4096 bits, sin(x) = [-0.7190842207 +/- 1.20e-11]
Using  8192 bits, sin(x) = [-0.7190842207 +/- 1.20e-11]
Using 16384 bits, sin(x) = [-0.7190842207 +/- 1.20e-11]
Using 32768 bits, sin(x) = [-0.7190842207 +/- 1.20e-11]

给定的示例有 x = 2016.1。

使用x = 0.1，我们得到以下输出:

Using    64 bits, sin(x) = [0.09983341665 +/- 3.18e-12]
Using   128 bits, sin(x) = [0.09983341665 +/- 3.18e-12]
Using   256 bits, sin(x) = [0.09983341665 +/- 3.18e-12]
Using   512 bits, sin(x) = [0.09983341665 +/- 3.18e-12]
Using  1024 bits, sin(x) = [0.09983341665 +/- 3.18e-12]
Using  2048 bits, sin(x) = [0.09983341665 +/- 3.18e-12]
Using  4096 bits, sin(x) = [0.09983341665 +/- 3.18e-12]
Using  8192 bits, sin(x) = [0.09983341665 +/- 3.18e-12]
Using 16384 bits, sin(x) = [0.09983341665 +/- 3.18e-12]
Using 32768 bits, sin(x) = [0.09983341665 +/- 3.18e-12]

这个精度似乎甚至低于 math.h 的 sin 函数。

我想提高说 e-40 的精度(或任何其他精度)。如果有人能指导我，我将不胜感激🙏🏼🙏🏼。

代码:

#include "arb.h"

void arb_sin_naive(arb_t res, const arb_t x, slong prec)
{
    arb_t s, t, u, tol;
    slong k;
    arb_init(s); arb_init(t); arb_init(u); arb_init(tol);

    arb_one(tol);
    arb_mul_2exp_si(tol, tol, -prec);  /* tol = 2^-prec */

    for (k = 0; ; k++)
    {
        arb_pow_ui(t, x, 2 * k + 1, prec);
        arb_fac_ui(u, 2 * k + 1, prec);
        arb_div(t, t, u, prec);  /* t = x^(2k+1) / (2k+1)! */

        arb_abs(u, t);
        if (arb_le(u, tol))   /* if |t| <= 2^-prec */
        {
            arb_add_error(s, u);    /* add |t| to the radius and stop */
            break;
        }

        if (k % 2 == 0)
            arb_add(s, s, t, prec);
        else
            arb_sub(s, s, t, prec);

    }

    arb_set(res, s);
    arb_clear(s); arb_clear(t); arb_clear(u); arb_clear(tol);
}

void main()
{
    arb_t x, y;
    slong prec;
    arb_init(x); arb_init(y);

    for (prec = 64; prec <= 32768 ; prec *= 2)
    {
        arb_set_str(x, "0.1", prec);
        arb_sin_naive(y, x, prec);
        printf("Using %5ld bits, sin(x) = ", prec);
        arb_printn(y, 10, 0); printf("\n");
    }

    arb_clear(x); arb_clear(y);
}

Answer 1

使用 x*(1 - x^2/3! + x^4/5! - x^6/7! ...) 实现更好的初始加法和更清晰的循环终止条件。

---

通常是正弦泰勒级数:sine(x) 是 x - x^3/3! + x^5/5! - x^6/7! ... 并且是 OP 使用的形式。

预计正弦泰勒级数的每一项都具有大约相同的相对进动，只是在后面的项中失去了一点精度。

然而，通过添加项(和跟踪容差)，总和并不比最大的 2 个项的绝对精度更精确。

通过将 sine(x) 形成为 x*(1 - x^2/3! + x^4/5! - x^6/7! ...)，我们在第一项 1.0 中具有无限精度，因此对于小的 x，精度受第二项限制，并且在向 1.0 添加一项时循环可以停止没有区别。

---

这并不能很好地解释为什么 OP 的结果停留在 0.09983341665 +/- 3.18e-12。

然而，注意对最大项求和(通过使其中一项具有无限精度的 1.0)会有所帮助。