view release on metacpan or  search on metacpan

quickjs/libbf.c  view on Meta::CPAN

#define CHUD_C 640320
#define CHUD_BITS_PER_TERM 47

static void chud_bs(bf_t *P, bf_t *Q, bf_t *G, int64_t a, int64_t b, int need_g,
                    limb_t prec)
{
    bf_context_t *s = P->ctx;
    int64_t c;

    if (a == (b - 1)) {
        bf_t T0, T1;
        
        bf_init(s, &T0);
        bf_init(s, &T1);
        bf_set_ui(G, 2 * b - 1);
        bf_mul_ui(G, G, 6 * b - 1, prec, BF_RNDN);
        bf_mul_ui(G, G, 6 * b - 5, prec, BF_RNDN);
        bf_set_ui(&T0, CHUD_B);
        bf_mul_ui(&T0, &T0, b, prec, BF_RNDN);
        bf_set_ui(&T1, CHUD_A);
        bf_add(&T0, &T0, &T1, prec, BF_RNDN);
        bf_mul(P, G, &T0, prec, BF_RNDN);
        P->sign = b & 1;

        bf_set_ui(Q, b);
        bf_mul_ui(Q, Q, b, prec, BF_RNDN);
        bf_mul_ui(Q, Q, b, prec, BF_RNDN);
        bf_mul_ui(Q, Q, (uint64_t)CHUD_C * CHUD_C * CHUD_C / 24, prec, BF_RNDN);
        bf_delete(&T0);
        bf_delete(&T1);
    } else {
        bf_t P2, Q2, G2;
        
        bf_init(s, &P2);
        bf_init(s, &Q2);
        bf_init(s, &G2);

        c = (a + b) / 2;
        chud_bs(P, Q, G, a, c, 1, prec);
        chud_bs(&P2, &Q2, &G2, c, b, need_g, prec);
        
        /* Q = Q1 * Q2 */
        /* G = G1 * G2 */
        /* P = P1 * Q2 + P2 * G1 */
        bf_mul(&P2, &P2, G, prec, BF_RNDN);
        if (!need_g)
            bf_set_ui(G, 0);
        bf_mul(P, P, &Q2, prec, BF_RNDN);
        bf_add(P, P, &P2, prec, BF_RNDN);
        bf_delete(&P2);

        bf_mul(Q, Q, &Q2, prec, BF_RNDN);
        bf_delete(&Q2);
        if (need_g)
            bf_mul(G, G, &G2, prec, BF_RNDN);
        bf_delete(&G2);
    }
}

/* compute Pi with faithful rounding at precision 'prec' using the
   Chudnovsky formula */
static void bf_const_pi_internal(bf_t *Q, limb_t prec)
{
    bf_context_t *s = Q->ctx;
    int64_t n, prec1;
    bf_t P, G;

    /* number of serie terms */
    n = prec / CHUD_BITS_PER_TERM + 1;
    /* XXX: precision analysis */
    prec1 = prec + 32;

    bf_init(s, &P);
    bf_init(s, &G);

    chud_bs(&P, Q, &G, 0, n, 0, BF_PREC_INF);
    
    bf_mul_ui(&G, Q, CHUD_A, prec1, BF_RNDN);
    bf_add(&P, &G, &P, prec1, BF_RNDN);
    bf_div(Q, Q, &P, prec1, BF_RNDF);
 
    bf_set_ui(&P, CHUD_C);
    bf_sqrt(&G, &P, prec1, BF_RNDF);
    bf_mul_ui(&G, &G, (uint64_t)CHUD_C / 12, prec1, BF_RNDF);
    bf_mul(Q, Q, &G, prec, BF_RNDN);
    bf_delete(&P);
    bf_delete(&G);
}

static int bf_const_get(bf_t *T, limb_t prec, bf_flags_t flags,
                        BFConstCache *c,
                        void (*func)(bf_t *res, limb_t prec), int sign)
{
    limb_t ziv_extra_bits, prec1;

    ziv_extra_bits = 32;
    for(;;) {
        prec1 = prec + ziv_extra_bits;
        if (c->prec < prec1) {
            if (c->val.len == 0)
                bf_init(T->ctx, &c->val);
            func(&c->val, prec1);
            c->prec = prec1;
        } else {
            prec1 = c->prec;
        }
        bf_set(T, &c->val);
        T->sign = sign;
        if (!bf_can_round(T, prec, flags & BF_RND_MASK, prec1)) {
            /* and more precision and retry */
            ziv_extra_bits = ziv_extra_bits  + (ziv_extra_bits / 2);
        } else {
            break;
        }
    }
    return bf_round(T, prec, flags);
}

static void bf_const_free(BFConstCache *c)
{
    bf_delete(&c->val);