arb/acb_hypgeom/test/t-coulomb.c

/*
    Copyright (C) 2019 Fredrik Johansson

    This file is part of Arb.

    Arb 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_hypgeom.h"

int main()
{
    slong iter;
    flint_rand_t state;

    flint_printf("coulomb....");
    fflush(stdout);

    flint_randinit(state);

    for (iter = 0; iter < 2000 * arb_test_multiplier(); iter++)
    {
        acb_t z, t, u, eta, l;
        acb_t F1, G1, Hpos1, Hneg1;
        acb_t F2, G2, Hpos2, Hneg2;
        slong prec1, prec2;
        unsigned int mask;

        acb_init(z); acb_init(t); acb_init(u); acb_init(eta); acb_init(l);
        acb_init(F1); acb_init(G1); acb_init(Hpos1); acb_init(Hneg1);
        acb_init(F2); acb_init(G2); acb_init(Hpos2); acb_init(Hneg2);

        prec1 = 2 + n_randint(state, 200);
        prec2 = 2 + n_randint(state, 200);

        acb_randtest_param(eta, state, 1 + n_randint(state, 200), 1 + n_randint(state, 10));
        acb_randtest_param(l, state, 1 + n_randint(state, 200), 1 + n_randint(state, 10));
        acb_randtest_param(z, state, 1 + n_randint(state, 200), 1 + n_randint(state, 100));
        acb_randtest_param(z, state, 1 + n_randint(state, 200), 1 + n_randint(state, 100));
        acb_randtest_param(t, state, 1 + n_randint(state, 200), 1 + n_randint(state, 100));
        acb_add(z, z, t, 1000);
        acb_sub(z, z, t, 1000);

        acb_hypgeom_coulomb(F1, G1, Hpos1, Hneg1, l, eta, z, prec1);
        acb_hypgeom_coulomb(F2, G2, Hpos2, Hneg2, l, eta, z, prec2);

        if (!acb_overlaps(F1, F2) || !acb_overlaps(G1, G2) ||
            !acb_overlaps(Hpos1, Hpos2) || !acb_overlaps(Hneg1, Hneg2))
        {
            flint_printf("FAIL: consistency\n\n");
            flint_printf("l = "); acb_printn(l, 30, 0); flint_printf("\n\n");
            flint_printf("eta = "); acb_printn(eta, 30, 0); flint_printf("\n\n");
            flint_printf("z = "); acb_printn(z, 30, 0); flint_printf("\n\n");
            flint_printf("F1 = "); acb_printn(F1, 30, 0); flint_printf("\n");
            flint_printf("F2 = "); acb_printn(F2, 30, 0); flint_printf("\n\n");
            flint_printf("G1 = "); acb_printn(G1, 30, 0); flint_printf("\n");
            flint_printf("G2 = "); acb_printn(G2, 30, 0); flint_printf("\n\n");
            flint_printf("Hpos1 = "); acb_printn(Hpos1, 30, 0); flint_printf("\n");
            flint_printf("Hpos2 = "); acb_printn(Hpos2, 30, 0); flint_printf("\n\n");
            flint_printf("Hneg1 = "); acb_printn(Hneg1, 30, 0); flint_printf("\n");
            flint_printf("Hneg2 = "); acb_printn(Hneg2, 30, 0); flint_printf("\n\n");
            flint_abort();
        }

        acb_mul_onei(t, F1);
        acb_add(t, t, G1, prec1);
        acb_div_onei(u, F1);
        acb_add(u, u, G1, prec1);

        if (!acb_overlaps(t, Hpos1) || !acb_overlaps(u, Hneg1))
        {
            flint_printf("FAIL: complex\n\n");
            flint_printf("l = "); acb_printn(l, 30, 0); flint_printf("\n\n");
            flint_printf("eta = "); acb_printn(eta, 30, 0); flint_printf("\n\n");
            flint_printf("z = "); acb_printn(z, 30, 0); flint_printf("\n\n");
            flint_printf("F1 = "); acb_printn(F1, 30, 0); flint_printf("\n");
            flint_printf("G1 = "); acb_printn(G1, 30, 0); flint_printf("\n");
            flint_printf("Hpos1 = "); acb_printn(Hpos1, 30, 0); flint_printf("\n");
            flint_printf("Hneg1 = "); acb_printn(Hneg1, 30, 0); flint_printf("\n");
            flint_printf("t = "); acb_printn(t, 30, 0); flint_printf("\n");
            flint_printf("u = "); acb_printn(u, 30, 0); flint_printf("\n");
            flint_abort();
        }

        mask = n_randlimb(state);

        if (n_randint(state, 2) == 0)
        {
            acb_randtest_param(t, state, 1 + n_randint(state, 200), 1 + n_randint(state, 100));
            acb_add(z, z, t, prec2);
            acb_sub(z, z, t, prec2);
        }

        /* also test aliasing */
        acb_set(G2, z);

        acb_hypgeom_coulomb((mask & 1) ? F2 : NULL,
                            (mask & 2) ? G2 : NULL,
                            (mask & 4) ? Hpos2 : NULL,
                            (mask & 8) ? Hneg2 : NULL, l, eta, G2, prec2);

        if (((mask & 1) && !acb_overlaps(F1, F2)) ||
            ((mask & 2) && !acb_overlaps(G1, G2)) ||
            ((mask & 4) && !acb_overlaps(Hpos1, Hpos2)) ||
            ((mask & 8) && !acb_overlaps(Hneg1, Hneg2)))
        {
            flint_printf("FAIL: consistency (mask)\n\n");
            flint_printf("mask = %u\n\n", mask);
            flint_printf("l = "); acb_printn(l, 30, 0); flint_printf("\n\n");
            flint_printf("eta = "); acb_printn(eta, 30, 0); flint_printf("\n\n");
            flint_printf("z = "); acb_printn(z, 30, 0); flint_printf("\n\n");
            flint_printf("F1 = "); acb_printn(F1, 30, 0); flint_printf("\n");
            flint_printf("F2 = "); acb_printn(F2, 30, 0); flint_printf("\n\n");
            flint_printf("G1 = "); acb_printn(G1, 30, 0); flint_printf("\n");
            flint_printf("G2 = "); acb_printn(G2, 30, 0); flint_printf("\n\n");
            flint_printf("Hpos1 = "); acb_printn(Hpos1, 30, 0); flint_printf("\n");
            flint_printf("Hpos2 = "); acb_printn(Hpos2, 30, 0); flint_printf("\n\n");
            flint_printf("Hneg1 = "); acb_printn(Hneg1, 30, 0); flint_printf("\n");
            flint_printf("Hneg2 = "); acb_printn(Hneg2, 30, 0); flint_printf("\n\n");
            flint_abort();
        }

        /* Check F_{l-1} G_{l} - F_{l} G_{l-1} = l (l^2 + eta^2)^(1/2) */
        acb_sub_ui(t, l, 1, prec2);
        acb_hypgeom_coulomb(F2, G2, NULL, NULL, t, eta, z, prec2);

        acb_mul_onei(t, eta);
        acb_add(t, t, l, prec2);
        acb_rsqrt(t, t, prec2);
        acb_div_onei(u, eta);
        acb_add(u, u, l, prec2);
        acb_rsqrt(u, u, prec2);
        acb_mul(u, t, u, prec2);
        acb_mul(u, u, l, prec2);

        acb_mul(t, F2, G1, prec2);
        acb_submul(t, F1, G2, prec2);

        if (!acb_overlaps(t, u))
        {
            flint_printf("FAIL: cross-product\n\n");
            flint_printf("l = "); acb_printn(l, 30, 0); flint_printf("\n\n");
            flint_printf("eta = "); acb_printn(eta, 30, 0); flint_printf("\n\n");
            flint_printf("z = "); acb_printn(z, 30, 0); flint_printf("\n\n");
            flint_printf("F1 = "); acb_printn(F1, 30, 0); flint_printf("\n");
            flint_printf("G1 = "); acb_printn(G1, 30, 0); flint_printf("\n");
            flint_printf("F2 = "); acb_printn(F2, 30, 0); flint_printf("\n");
            flint_printf("G2 = "); acb_printn(G2, 30, 0); flint_printf("\n");
            flint_printf("t = "); acb_printn(t, 30, 0); flint_printf("\n");
            flint_printf("u = "); acb_printn(u, 30, 0); flint_printf("\n");
            flint_abort();
        }

        acb_clear(z); acb_clear(t); acb_clear(u); acb_clear(eta); acb_clear(l);
        acb_clear(F1); acb_clear(G1); acb_clear(Hpos1); acb_clear(Hneg1);
        acb_clear(F2); acb_clear(G2); acb_clear(Hpos2); acb_clear(Hneg2);
    }

    flint_randclear(state);
    flint_cleanup();
    flint_printf("PASS\n");
    return EXIT_SUCCESS;
}