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

Question

我发现 Arb library, which 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）会有所帮助。

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

How to get higher precision for sine using the Arb library?

c

precision

trigonometry

arbitrary-precision