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

Question

I found the Arb library, which should be able to compute very high precision values of sine, given enough time. However, I am unable to.

Trying the sine example, I could get the predicted output.

However, when I tried to increase the precision by increasing the number of bits from 4096 to 32768, I was unable to:

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]

The given example has x = 2016.1.

With x = 0.1, we get the following output:

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]

This precision seems to be less than even the sin function of math.h.

I would like to increase the precision to say e-40 (or any other precision). I would be most grateful if someone could guide me 🙏🏼🙏🏼.

The code:

#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

Use x*(1 - x^2/3! + x^4/5! - x^6/7! ...) to effect a better initial addition and clearer loop terminating condition.

---

The usually sine Taylor's series: sine(x) is x - x^3/3! + x^5/5! - x^6/7! ... and is the form used by OP.

Each term of the sine Taylor's series is expected to have about the same relative precession, only losing a little precision in later terms.

Yet by adding terms (and tracking tolerance), the sum is no more precise than the absolute precision of the largest 2 terms.

By forming sine(x) as x*(1 - x^2/3! + x^4/5! - x^6/7! ...), we have unlimited precision in the first term 1.0, thus the precision, for small x, is limited by the 2nd term and the loop can stop when adding a term to 1.0 makes no difference.

---

This does not well explain why OP's result was stuck at 0.09983341665 +/- 3.18e-12.

Yet taking care of summing the largest terms (by having one of them 1.0 with unlimited precision) helps.