/*
    Copyright (C) 2018 Association des collaborateurs de D.H.J Polymath

    This is free software: you can redistribute it and/or modify it under
    the terms of the GNU Lesser General Public License (LGPL) as published
    by the Free Software Foundation; either version 2.1 of the License, or
    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
*/

#include "acb_calc.h"

void H_approx(arb_t res, const acb_t z, slong prec)
{
    arb_t logx, pixd8;

    arb_init(logx);
    arb_log(logx, acb_realref(z), prec);

    arb_init(pixd8);
    arb_const_pi(pixd8, prec);
    arb_mul(pixd8, pixd8, acb_realref(z), prec);
    arb_mul_2exp_si(pixd8, pixd8, -3);

    arb_add_si(res, acb_imagref(z), 7, prec);
    arb_mul(res, res, logx, prec);
    arb_mul_2exp_si(res, res, -2);
    arb_sub(res, res, pixd8, prec);
    arb_exp(res, res, prec);

    arb_clear(logx);
    arb_clear(pixd8);
}

void integral_remainder(mag_t m, const arb_t b)
{
    arb_get_mag(m, b);
    mag_mul_2exp_si(m, m, 2);
    mag_exp_lower(m, m);
    mag_expinv(m, m);
    mag_mul_ui(m, m, 25);
}

void integration_limit(arb_t b, const arb_t eps, slong prec)
{
    arb_zero(b);
    arb_get_lbound_arf(arb_midref(b), eps, prec);
    arb_div_ui(b, b, 25, prec);
    arb_log(b, b, prec);
    arb_neg(b, b);
    arb_log(b, b, prec);
    arb_mul_2exp_si(b, b, -2);
    arb_get_ubound_arf(arb_midref(b), b, prec);
    arb_get_mid_arb(b, b);
}

void phi_domain(acb_t z)
{
    arb_ptr x, y;
    slong prec;
    arb_t t;

    prec = 64;

    x = acb_realref(z);
    arb_zero_pm_inf(x);

    y = acb_imagref(z);
    arb_zero(y);

    arb_init(t);
    arb_const_pi(t, prec);
    arb_div_ui(t, t, 12, prec);
    arb_get_mag_lower(arb_radref(y), t);
    arb_clear(t);
}

int phi_domain_contains(const acb_t u)
{
    int result;
    acb_t d;
    acb_init(d);
    phi_domain(d);
    result = acb_contains(d, u);
    acb_clear(d);
    return result;
}

void phi_remainder(arb_t z, const acb_t u, ulong n0, slong prec)
{
    if (n0 < 1)
        flint_abort();

    if (phi_domain_contains(u))
    {
        acb_t exp4u;
        arb_t b, d;

        acb_init(exp4u);
        acb_mul_2exp_si(exp4u, u, 2);
        acb_exp(exp4u, exp4u, prec);

        arb_init(b);
        arb_const_pi(b, prec);
        arb_mul_ui(b, b, n0, prec);
        arb_mul_2exp_si(b, b, -1);
        arb_mul(b, b, acb_realref(exp4u), prec);

        arb_init(d);
        arb_neg(d, b);
        arb_exp(d, d, prec);
        arb_sub_ui(d, d, 1, prec);

        arb_mul_ui(z, b, n0, prec);
        arb_neg(z, z);
        arb_exp(z, z, prec);
        arb_div(z, z, d, prec);
        arb_mul_si(z, z, -50, prec);

        acb_clear(exp4u);
        arb_clear(b);
        arb_clear(d);
    }
    else
    {
        arb_indeterminate(z);
    }
}

void phi_summand(acb_t z, const acb_t u, ulong n, slong prec)
{
    if (phi_domain_contains(u))
    {
        arb_t pin2;
        acb_t a, b;

        arb_init(pin2);
        arb_const_pi(pin2, prec);
        arb_mul_ui(pin2, pin2, n*n, prec);

        acb_init(a);
        acb_mul_2exp_si(a, u, 2);
        acb_exp(a, a, prec);
        acb_mul_arb(a, a, pin2, prec);

        acb_init(b);
        acb_mul_2exp_si(b, a, 1);
        acb_sub_ui(b, b, 3, prec);

        acb_mul_ui(z, u, 5, prec);
        acb_sub(z, z, a, prec);
        acb_exp(z, z, prec);
        acb_mul(z, z, b, prec);
        acb_mul_arb(z, z, pin2, prec);

        arb_clear(pin2);
        acb_clear(a);
        acb_clear(b);
    }
    else
    {
        acb_indeterminate(z);
    }
}

void phi(acb_t z, const acb_t u, slong prec)
{
    if (phi_domain_contains(u))
    {
        ulong n;
        acb_t s;
        arb_t tail;
        arf_t err;
        mag_t m, tail_mag_lb, tail_mag_ub;

        acb_init(s);
        arb_init(tail);
        arf_init(err);
        mag_init(tail_mag_lb);
        mag_init(tail_mag_ub);
        mag_init(m);

        n = 1;
        acb_zero(z);
        while (1)
        {
            phi_summand(s, u, n, prec);
            acb_add(z, z, s, prec);

            n++;

            phi_remainder(tail, u, n, prec);
            acb_get_rad_ubound_arf(err, z, prec);
            arf_get_mag(m, err);

            arb_get_mag_lower(tail_mag_lb, tail);
            if (mag_cmp(tail_mag_lb, m) <= 0)
            {
                arb_get_mag(tail_mag_ub, tail);
                acb_add_error_mag(z, tail_mag_ub);
                break;
            }
        }

        acb_clear(s);
        arb_clear(tail);
        arf_clear(err);
        mag_clear(tail_mag_lb);
        mag_clear(tail_mag_ub);
        mag_clear(m);
    }
    else
    {
        acb_indeterminate(z);
    }
}

int integrand(acb_ptr res, const acb_t u, void *param, slong order, slong prec)
{
    if (phi_domain_contains(u))
    {
        acb_t a, b, c;
        arb_ptr t;
        acb_ptr z;

        t = ((arb_ptr *)(param))[0];
        z = ((acb_ptr *)(param))[1];

        acb_init(a);
        acb_sqr(a, u, prec);
        acb_mul_arb(a, a, t, prec);
        acb_exp(a, a, prec);

        acb_init(b);
        phi(b, u, prec);

        acb_init(c);
        acb_mul(c, z, u, prec);
        acb_cos(c, c, prec);

        acb_mul(res, a, b, prec);
        acb_mul(res, res, c, prec);

        acb_clear(a);
        acb_clear(b);
        acb_clear(c);
    }
    else
    {
        acb_indeterminate(res);
    }
    return 0;
}


void H(acb_t h, const arb_t t, const acb_t z, slong prec)
{
    arb_t width;
    acb_t a, b;
    void *param[2];
    slong rel_goal;
    int calc_result;
    mag_t m, abs_tol;

    rel_goal = prec;
    param[0] = (void *) t;
    param[1] = (void *) z;

    arb_init(width);
    H_approx(width, z, prec);
    arb_mul_2exp_si(width, width, -rel_goal);

    acb_init(a);
    acb_init(b);
    integration_limit(acb_realref(b), width, prec);

    mag_init(abs_tol);
    arb_get_mag(abs_tol, width);

    calc_result = acb_calc_integrate(
            h, integrand, param, a, b,
            rel_goal, abs_tol, NULL, prec);
    if (calc_result == ARB_CALC_NO_CONVERGENCE)
    {
        acb_indeterminate(h);
        goto finish;
    }

    mag_init(m);
    integral_remainder(m, acb_realref(b));
    acb_add_error_mag(h, m);

finish:

    acb_clear(a);
    acb_clear(b);
    arb_clear(width);
    mag_clear(abs_tol);
    mag_clear(m);
}

int main(int argc, char *argv[])
{
    int result;
    arb_t t;
    acb_t z, h;
    slong digits, prec, high_prec;
    arb_t half, yabs;
    double dx;

    arb_init(t);
    acb_init(z);
    acb_init(h);
    arb_init(half);
    arb_init(yabs);

    result = EXIT_SUCCESS;

    if (argc != 5)
    {
        result = EXIT_FAILURE;
        goto finish;
    }

    digits = atol(argv[4]);
    high_prec = digits * 3.32192809488736 + 10;

    arb_one(half);
    arb_mul_2exp_si(half, half, -1);

    dx = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
    prec = FLINT_MAX(32, (int) (dx / 4));
    acb_indeterminate(h);
    while (acb_rel_accuracy_bits(h) < high_prec)
    {
        prec *= 2;
        arb_set_str(t, argv[1], prec);
        arb_set_str(acb_realref(z), argv[2], prec);
        arb_set_str(acb_imagref(z), argv[3], prec);
        arb_abs(yabs, acb_imagref(z));
        if (!arb_lt(t, half) || !arb_lt(yabs, half))
        {
            result = EXIT_FAILURE;
            goto finish;
        }
        H(h, t, z, prec);
    }

    flint_printf("Re[H]: ");
    arf_printd(arb_midref(acb_realref(h)), digits);
    flint_printf("\n");
    flint_printf("Im[H]: ");
    arf_printd(arb_midref(acb_imagref(h)), digits);
    flint_printf("\n");

finish:

    if (result == EXIT_FAILURE)
    {
        flint_printf(
            "%s t x y d\n\n"
            "Evaluate H_t(x + yi) to d significant digits "
            "where H is the function involved "
            "in the definition of the De Bruijn-Newman constant.\n"
            "This implementation requires t < 0.5 and |y| < 0.5.\n",
            argv[0]);
    }

    arb_clear(t);
    acb_clear(z);
    acb_clear(h);
    arb_clear(half);
    arb_clear(yabs);

    flint_cleanup();
    return result;
}