illinois.c

#include "acb_dirichlet.h"

/* todo: separate prec, eval_prec... */
void
acb_dirichlet_zeta_zero_refine_illinois(arb_t res, const arf_t ra, const arf_t rb, slong prec)
{
    arf_t a, b, fa, fb, c, fc, t;
    acb_t z;
    slong k;
    int asign, bsign, csign;

    arf_init(a);
    arf_init(b);
    arf_init(c);
    arf_init(fa);
    arf_init(fb);
    arf_init(fc);
    arf_init(t);
    acb_init(z);

    arf_set(a, ra);
    arf_set(b, rb);

    acb_zero(z);
    arf_set(arb_midref(acb_realref(z)), a);
    acb_dirichlet_hardy_z(z, z, NULL, NULL, 1, prec);
    arf_set(fa, arb_midref(acb_realref(z)));
    asign = arb_sgn_nonzero(acb_realref(z));
    if (asign == 0)
        goto cleanup;

    acb_zero(z);
    arf_set(arb_midref(acb_realref(z)), b);
    acb_dirichlet_hardy_z(z, z, NULL, NULL, 1, prec);
    arf_set(fb, arb_midref(acb_realref(z)));
    bsign = arb_sgn_nonzero(acb_realref(z));
    if (bsign == 0)
        goto cleanup;

    /* todo: verify sign change here? */

    for (k = 0; k < 30; k++)
    {
        /* c = a - fa * (b - a) / (fb - fa) */
        arf_sub(c, b, a, prec, ARF_RND_DOWN);
        arf_sub(t, fb, fa, prec, ARF_RND_DOWN);
        arf_div(c, c, t, prec, ARF_RND_DOWN);
        arf_mul(c, c, fa, prec, ARF_RND_DOWN);
        arf_sub(c, a, c, prec, ARF_RND_DOWN);

        /* verify that c is sandwiched between ra and rb */
        if (!arf_is_finite(c) || 
            !(arf_cmp(ra, c) < 0 && arf_cmp(c, rb) < 0))
        {
            printf("no sandwich\n");
            arf_printd(a, 30); printf("\n");
            arf_printd(b, 30); printf("\n");
            arf_printd(c, 30); printf("\n");

            goto cleanup;
        }

        acb_zero(z);
        arf_set(arb_midref(acb_realref(z)), c);
        acb_dirichlet_hardy_z(z, z, NULL, NULL, 1, prec);
        arf_set(fc, arb_midref(acb_realref(z)));
        csign = arb_sgn_nonzero(acb_realref(z));
        if (csign == 0)
            goto cleanup;

        if (csign != bsign)
        {
            arf_set(a, b);
            arf_set(fa, fb);
            asign = bsign;

            arf_set(b, c);
            arf_set(fb, fc);
            bsign = csign;
        }
        else
        {
            arf_set(b, c);
            arf_set(fb, fc);
            bsign = csign;

            arf_mul_2exp_si(fa, fa, -1);
        }

        if (asign == bsign)
            flint_abort();

        arf_sub(t, a, b, prec, ARF_RND_DOWN);
        arf_abs(t, t);

        arf_printd(a, 30); printf("    ");
        arf_printd(b, 30); printf("    ");
        arf_printd(t, 6); printf("\n");

        if (arf_cmpabs_2exp_si(t, -prec + arf_abs_bound_lt_2exp_si(rb) + 8) < 0)
            break;
    }

cleanup:

    if (arf_cmp(a, b) > 0)
        arf_swap(a, b);

    arb_set_interval_arf(res, a, b, prec);

    arf_clear(a);
    arf_clear(b);
    arf_clear(c);
    arf_clear(fa);
    arf_clear(fb);
    arf_clear(fc);
    arf_clear(t);
    acb_clear(z);
}


int main()
{
    arf_t a, b;
    arb_t r;

    arf_init(a);
    arf_init(b);
    arb_init(r);

    arf_set_d(a, 14.0);
    arf_set_d(b, 15.0);
    acb_dirichlet_zeta_zero_refine_illinois(r, a, b, 53);
    arb_printn(r, 30, 0); printf("\n\n");

    arf_set_d(a, 600269.01);
    arf_set_d(b, 600270.30);
    acb_dirichlet_zeta_zero_refine_illinois(r, a, b, 53);
    arb_printn(r, 30, 0); printf("\n\n");

    arf_clear(a);
    arf_clear(b);
    arb_clear(r);
}