zeta zero script

gistfile1.txt

#include "fmprb_poly.h"
#include "fmpcb_poly.h"
#include "profiler.h"

/* theta(t) = im(loggamma((2it+1)/4)) - log(pi)/2 * t */

void
theta(fmprb_poly_t res, fmprb_poly_t t, long n, long prec)
{
    fmpcb_poly_t s;
    fmprb_t c, d;
    long i;

    fmpcb_poly_init(s);
    fmprb_init(c);
    fmprb_init(d);

    fmpcb_poly_fit_length(s, fmprb_poly_length(t));
    _fmpcb_poly_set_length(s, fmprb_poly_length(t));
    for (i = 0; i < fmprb_poly_length(t); i++)
    {
        fmprb_set_si(fmpcb_realref(fmpcb_poly_get_coeff_ptr(s, i)), i == 0);
        fmprb_mul_2exp_si(fmpcb_imagref(fmpcb_poly_get_coeff_ptr(s, i)),
            fmprb_poly_get_coeff_ptr(t, i), 1);
    }
    _fmpcb_poly_normalise(s);
    fmpcb_poly_scalar_mul_2exp_si(s, s, -2);

    fmpcb_poly_lgamma_series(s, s, n, prec);

    fmprb_const_pi(c, prec);
    fmprb_log(c, c, prec);
    fmprb_mul_2exp_si(c, c, -1);

    fmprb_poly_fit_length(res, fmpcb_poly_length(s));
    _fmprb_poly_set_length(res, fmpcb_poly_length(s));
    for (i = 0; i < fmpcb_poly_length(s); i++)
    {
        if (i < fmprb_poly_length(t))
            fmprb_mul(d, fmprb_poly_get_coeff_ptr(t, i), c, prec);
        else
            fmprb_zero(d);

        fmprb_sub(fmprb_poly_get_coeff_ptr(res, i),
            fmpcb_imagref(fmpcb_poly_get_coeff_ptr(s, i)), d, prec);
    }
    _fmprb_poly_normalise(res);

    fmpcb_poly_clear(s);
    fmprb_clear(c);
    fmprb_clear(d);
}

void
siegelz(fmprb_poly_t res, fmprb_poly_t t, long n, long prec)
{
    fmpcb_poly_t s, z;
    fmprb_poly_t u, re, im, zre, zim;
    fmpcb_t a;
    long i;

    fmpcb_poly_init(s);
    fmpcb_poly_init(z);
    fmprb_poly_init(u);
    fmprb_poly_init(re);
    fmprb_poly_init(im);
    fmprb_poly_init(zre);
    fmprb_poly_init(zim);
    fmpcb_init(a);

    /* build s = 1/2 + it */
    fmpcb_poly_fit_length(s, fmprb_poly_length(t));
    _fmpcb_poly_set_length(s, fmprb_poly_length(t));
    for (i = 0; i < fmprb_poly_length(t); i++)
    {
        fmprb_set_si(fmpcb_realref(fmpcb_poly_get_coeff_ptr(s, i)), i == 0);
        fmprb_mul_2exp_si(fmpcb_realref(fmpcb_poly_get_coeff_ptr(s, i)),
            fmpcb_realref(fmpcb_poly_get_coeff_ptr(s, i)), -1);

        fmprb_set(fmpcb_imagref(fmpcb_poly_get_coeff_ptr(s, i)),
            fmprb_poly_get_coeff_ptr(t, i));
    }
    _fmpcb_poly_normalise(s);

    /* evaluate zeta */
    fmpcb_one(a);
    fmpcb_poly_zeta_series(z, s, a, 0, n, prec);

    for (i = 0; i < fmpcb_poly_length(z); i++)
    {
        fmprb_poly_set_coeff_fmprb(zre, i, fmpcb_realref(fmpcb_poly_get_coeff_ptr(z, i)));
        fmprb_poly_set_coeff_fmprb(zim, i, fmpcb_imagref(fmpcb_poly_get_coeff_ptr(z, i)));
    }

    /* evaluate exp(i theta(t)) */
    theta(u, t, n, prec);

    fmprb_poly_sin_cos_series(im, re, u, n, prec);

    fmprb_poly_mullow(t, re, zre, n, prec);
    fmprb_poly_mullow(u, im, zim, n, prec);

    fmprb_poly_sub(res, t, u, prec);

    fmpcb_poly_clear(s);
    fmpcb_poly_clear(z);
    fmprb_poly_clear(u);
    fmprb_poly_clear(re);
    fmprb_poly_clear(im);
    fmprb_poly_clear(zre);
    fmprb_poly_clear(zim);
    fmpcb_clear(a);
}

int
newton_step(fmprb_t xnew, const fmprb_t x, const fmprb_t convergence_interval, const fmpr_t convergence_factor, long prec)
{
    fmpr_t err;
    fmprb_t t, u, v;
    fmprb_poly_t tmp;
    int result;

    fmpr_init(err);
    fmprb_init(t);
    fmprb_init(u);
    fmprb_init(v);
    fmprb_poly_init(tmp);

    fmpr_mul(err, fmprb_radref(x), fmprb_radref(x), FMPRB_RAD_PREC, FMPR_RND_UP);
    fmpr_mul(err, err, convergence_factor, FMPRB_RAD_PREC, FMPR_RND_UP);

    fmpr_set(fmprb_midref(t), fmprb_midref(x));
    fmpr_zero(fmprb_radref(t));

    fmprb_poly_set_coeff_fmprb(tmp, 0, t);
    fmprb_poly_set_coeff_si(tmp, 1, 1);

    siegelz(tmp, tmp, 2, prec);

    fmprb_poly_get_coeff_fmprb(u, tmp, 0);
    fmprb_poly_get_coeff_fmprb(v, tmp, 1);

    fmprb_div(u, u, v, prec);
    fmprb_sub(u, t, u, prec);

    fmprb_add_error_fmpr(u, err);

    if (fmprb_contains(convergence_interval, u) &&
        (fmpr_cmp(fmprb_radref(u), fmprb_radref(x)) < 0))
    {
        fmprb_swap(xnew, u);
        result = 1;
    }
    else
    {
        fmprb_set(xnew, x);
        result = 0;
    }

    fmprb_clear(t);
    fmprb_clear(u);
    fmprb_clear(v);
    fmpr_clear(err);
    fmprb_poly_clear(tmp);

    return result;
}

long _fmpr_mag(const fmpr_t c);

void
refine_root(fmprb_t r, const fmprb_t start, const fmprb_t convergence_interval, const fmpr_t convergence_factor, long eval_extra_prec, long prec)
{
    long precs[FLINT_BITS];
    long i, iters, wp, padding, start_prec;

    start_prec = fmprb_rel_accuracy_bits(start);

    printf("initial accuracy: %ld\n", start_prec);

    padding = 5 + _fmpr_mag(convergence_factor);
    precs[0] = prec + padding;
    iters = 1;
    while ((iters < FLINT_BITS) && (precs[iters-1] + padding > 2*start_prec))
    {
        precs[iters] = (precs[iters-1] / 2) + padding;
        iters++;

        if (iters == FLINT_BITS)
        {
            printf("newton_refine_root: initial value too imprecise\n");
            abort();
        }
    }

    fmprb_set(r, start);

    for (i = iters - 1; i >= 0; i--)
    {
        wp = precs[i] + eval_extra_prec;

        printf("newton step: wp = %ld + %ld = %ld\n", precs[i], eval_extra_prec, wp);

        if (!newton_step(r, r, convergence_interval, convergence_factor, wp))
        {
            printf("warning: newton_refine_root: improvement failed\n");
            break;
        }

    }
}

void
newton_convergence_factor(fmpr_t convergence_factor, const fmprb_t convergence_interval, long prec)
{
    fmprb_poly_t tmp;
    fmprb_t t, u;

    fmprb_init(t);
    fmprb_init(u);
    fmprb_poly_init(tmp);

    fmprb_poly_set_coeff_fmprb(tmp, 0, convergence_interval);
    fmprb_poly_set_coeff_si(tmp, 1, 1);

    siegelz(tmp, tmp, 3, prec);

    fmprb_poly_get_coeff_fmprb(t, tmp, 1);
    fmprb_poly_get_coeff_fmprb(u, tmp, 2);

    fmprb_div(t, u, t, prec);
    fmprb_mul_2exp_si(t, t, -1);

    fmprb_get_abs_ubound_fmpr(convergence_factor, t, prec);

    fmprb_poly_clear(tmp);
    fmprb_clear(t);
    fmprb_clear(u);
}


int main()
{
    fmprb_t interval;
    fmprb_t root, root2;
    fmprb_t initial;
    fmpr_t convergence_factor;
    timeit_t t0;

    fmprb_init(interval);
    fmprb_init(root);
    fmprb_init(root2);
    fmprb_init(initial);
    fmpr_init(convergence_factor);

    long digits = 1000;
    long prec = digits * 3.32192809488736 + 2;

    fmpr_set_d(fmprb_midref(initial), 14.1347251417347);
    fmpr_set_d(fmprb_radref(initial), 1e-12);

    fmpr_set_d(fmprb_midref(interval), 14.1347251417347);
    fmpr_set_d(fmprb_radref(interval), 1e-6);

    newton_convergence_factor(convergence_factor, interval, 53);

    printf("factor:\n");
    fmpr_printd(convergence_factor, 15); printf("\n");

    long i;
    timeit_start(t0);
    refine_root(root, initial, interval, convergence_factor, 20, prec);
    timeit_stop(t0);
    printf("time: %g s\n", (double) t0->cpu / 1000.0);

    FILE * fp;
    fp = fopen("zetazero.txt", "w");

    fmpz_fprint(fp, fmpr_manref(fmprb_midref(root))); fprintf(fp, "\n");
    fmpz_fprint(fp, fmpr_expref(fmprb_midref(root))); fprintf(fp, "\n");
    fmpz_fprint(fp, fmpr_manref(fmprb_radref(root))); fprintf(fp, "\n");
    fmpz_fprint(fp, fmpr_expref(fmprb_radref(root))); fprintf(fp, "\n");
    fclose(fp);

/*
    timeit_start(t0);
    refine_root(root2, initial, interval, convergence_factor, 20, prec + 10);
    timeit_stop(t0);
    printf("again time: %g s\n", (double) t0->cpu / 1000.0);
*/

    fmpcb_t s, z;
    fmpcb_init(s);
    fmpcb_init(z);
    fmprb_set(fmpcb_imagref(s), root);
    fmprb_one(fmpcb_realref(s));
    fmprb_mul_2exp_si(fmpcb_realref(s), fmpcb_realref(s), -1);

    timeit_start(t0);
    fmpcb_zeta(z, s, prec + 10);
    timeit_stop(t0);
    printf("zeta time: %g s\n", (double) t0->cpu / 1000.0);

    

    fmprb_printd(root, 100); printf("\n\n");
    fmpcb_printd(z, 10); printf("\n\n");

    printf("ok: %d, %d\n", fmprb_contains(root, root2), fmpcb_contains_zero(z));

    printf("\n");

    fmprb_clear(interval);
    fmprb_clear(root);
    fmprb_clear(root2);
    fmprb_clear(initial);
    fmpr_clear(convergence_factor);
}