/*
    Copyright (C) 2013 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_poly.h"

#define EPS 1e-13
#define NUM_DERIVS 3
#define NUM_TESTS 56

const double polylog_testdata[NUM_TESTS][10] = {
  {-2.0, 0.0, 0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0},
  {-2.0, 0.0, -0.5, 0.0,
    -0.0740740740740740741, 0.0,
    -0.158799035455864473, 0.0,
    0.00278425812292977954, 0.0},
  {-2.0, 0.0, 0.99609375, 0.0,
    33358080.0, 0.0,
    -215693527.719090182, 0.0,
    703924885.809675425, 0.0},
  {-2.0, 0.0, 1.00390625, 0.0,
    -33751296.0, 0.0,
    218367905.967311029, -106032823.56283664,
    -546519337.1290267, 686023009.262281684},
  {-2.0, 0.0, -4.25, 0.0,
    0.0954540546377281071, 0.0,
    -0.0297817299750221174, 0.0,
    -0.196648773789941696, 0.0},
  {-2.0, 0.0, 3.5, 0.0,
    -1.008, 0.0,
    0.691243661008350802, -3.19575256895409845,
    4.53340070033764442, 2.22882281653983385},
  {-2.0, 0.0, 0.0, -1.0,
    0.0, 0.5,
    -0.852556797635011582, -0.25189536315108886,
    0.345626503121598963, -0.46292146649778986},
  {-2.0, 0.0, 0.125, 0.125,
    0.071424, 0.280832,
    0.0646851488459364134, -0.118051747326526959,
    -0.0396751348737248634, 0.0452067417426917285},
  {-1.0, 0.0, 0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0},
  {-1.0, 0.0, -0.5, 0.0,
    -0.222222222222222222, 0.0,
    -0.131883201632587537, 0.0,
    0.0199225628389378777, 0.0},
  {-1.0, 0.0, 0.99609375, 0.0,
    65280.0, 0.0,
    -389461.499432291984, 0.0,
    1182814.77437627002, 0.0},
  {-1.0, 0.0, 1.00390625, 0.0,
    65792.0, 0.0,
    -392773.095775796478, 206691.925664168516,
    868954.597422132011, -1233932.55215831045},
  {-1.0, 0.0, -4.25, 0.0,
    -0.154195011337868481, 0.0,
    -0.474824989764482489, 0.0,
    -0.21428706813151403, 0.0},
  {-1.0, 0.0, 3.5, 0.0,
    0.56, 0.0,
    -0.322735710701998933, 2.00176023742981729,
    -3.09698193717149541, -0.395213225234415973},
  {-1.0, 0.0, 0.0, -1.0,
    -0.5, 0.0,
    -0.183855151549436048, -0.58312180806163756,
    0.254077251327014729, 0.0351429519611278982},
  {-1.0, 0.0, 0.125, 0.125,
    0.1088, 0.1984,
    0.0192018441183911, -0.0548368329876085499,
    -0.0115465156713130216, 0.020954643261003461},
  {0.0, 0.0, 0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0},
  {0.0, 0.0, -0.5, 0.0,
    -0.333333333333333333, 0.0,
    -0.0904116649848191551, 0.0,
    0.0199028637728058217, 0.0},
  {0.0, 0.0, 0.99609375, 0.0,
    255.0, 0.0,
    -1269.73106874463507, 0.0,
    3359.60296533314103, 0.0},
  {0.0, 0.0, 1.00390625, 0.0,
    -257.0, 0.0,
    1273.86116119868176, -805.817494984366455,
    -2113.9649224353639, 4004.84182305521042},
  {0.0, 0.0, -4.25, 0.0,
    -0.809523809523809524, 0.0,
    -0.795769271870623909, 0.0,
    -0.0965372637237445457, 0.0},
  {0.0, 0.0, 3.5, 0.0,
    -1.4, 0.0,
    -1.78816683065098841, -2.50773109725857055,
    1.83790031889197745, -2.01262260402527512},
  {0.0, 0.0, 0.0, -1.0,
    -0.5, -0.5,
    0.120782237635245222, -0.391594392706836776,
    0.0669079578849868665, 0.116416441749488305},
  {0.0, 0.0, 0.125, 0.125,
    0.12, 0.16,
    0.00583206536388767829, -0.0256940749384562616,
    -0.00344424777009813297, 0.00965046958017889366},
  {-0.5, 0.0, 0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0},
  {-0.5, 0.0, -0.5, 0.0,
    -0.283012810746506023, 0.0,
    -0.111080563146094489, 0.0,
    0.0211430885590715694, 0.0},
  {-0.5, 0.0, 0.99609375, 0.0,
    3619.14124113924805, 0.0,
    -20195.2825251289493, 0.0,
    58032.3050555790807, 0.0},
  {-0.5, 0.0, 1.00390625, 0.0,
    -0.207985517845521836, 3640.61848390552728,
    -11437.7014354502895, -20327.8206114649109,
    63861.4335397379033, 40487.2831296029499},
  {-0.5, 0.0, -4.25, 0.0,
    -0.441427843410865187, 0.0,
    -0.665497904840929966, 0.0,
    -0.162969337925229287, 0.0},
  {-0.5, 0.0, 3.5, 0.0,
    -0.232043460812781422, 0.63203565652521718,
    -2.44367789279529917, 0.119367210264842707,
    -0.763914228371074801, -2.81228485141461073},
  {-0.5, 0.0, 0.0, -1.0,
    -0.537549381115898973, -0.27517974122882025,
    0.0158915835170861709, -0.505981074388456133,
    0.147999289996985447, 0.105632732120539676},
  {-0.5, 0.0, 0.125, 0.125,
    0.116021225799463628, 0.175602658223910358,
    0.0105539353881710306, -0.0374742285737420805,
    -0.0062875342272019156, 0.0142081398192247589},
  {1.0, 0.0, 0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0},
  {1.0, 0.0, -0.5, 0.0,
    -0.405465108108164382, 0.0,
    -0.0556552111040376179, 0.0,
    0.0145754167941056437, 0.0},
  {1.0, 0.0, 0.99609375, 0.0,
    5.54517744447956248, 0.0,
    -13.0766609768056118, 0.0,
    24.0940545588450487, 0.0},
  {1.0, 0.0, 1.00390625, 0.0,
    5.54517744447956248, -3.14159265358979324,
    -8.16844417312495553, 15.6134381896777968,
    -0.38789687903338657, -36.2148484363837508},
  {1.0, 0.0, -4.25, 0.0,
    -1.65822807660353238, 0.0,
    -0.861140718490863904, 0.0,
    0.0240108114276114004, 0.0},
  {1.0, 0.0, 3.5, 0.0,
    -0.916290731874155065, -3.14159265358979324,
    2.57284574234941582, -2.52133906787958123,
    1.44241441233798233, 1.57208443290129947},
  {1.0, 0.0, 0.0, -1.0,
    -0.346573590279972655, -0.78539816339744831,
    0.160292055087885226, -0.192901316796912429,
    -0.0115066749133700615, 0.0770708622146679417},
  {1.0, 0.0, 0.125, 0.125,
    0.123430038965762899, 0.141897054604163923,
    0.00180689203856968544, -0.0122114794810442984,
    -0.00105066979284378462, 0.00449853805632375062},
  {2.0, 0.0, 0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0},
  {2.0, 0.0, -0.5, 0.0,
    -0.448414206923646202, 0.0,
    -0.0320237060779612787, 0.0,
    0.0092536571135294828, 0.0},
  {2.0, 0.0, 0.99609375, 0.0,
    1.61932072927810507, 0.0,
    -0.863010139919679163, 0.0,
    0.825515949454683347, 0.0},
  {2.0, 0.0, 1.00390625, 0.0,
    1.6704551616562229, -0.0122479400888173039,
    -0.992661633138610541, 0.0731191212424626706,
    1.04852574935488925, -0.214307793003452771},
  {2.0, 0.0, -4.25, 0.0,
    -2.46898738747224279, 0.0,
    -0.739287387997829118, 0.0,
    0.0874470644522297429, 0.0},
  {2.0, 0.0, 3.5, 0.0,
    2.1959348115757808, -3.93567093851438969,
    2.87640729249480637, 0.777030723254221356,
    -0.81939570865486996, 1.19241843763256774},
  {2.0, 0.0, 0.0, -1.0,
    -0.205616758356028305, -0.915965594177219015,
    0.117193531789480469, -0.0815807361165927951,
    -0.0255408233794369879, 0.0372076062178391175},
  {2.0, 0.0, 0.125, 0.125,
    0.124500148439393858, 0.133240721203487057,
    0.000569113907620539642, -0.00588251562210034232,
    -0.000326714635669985778, 0.00213086661076891839},
  {-2.0, 2.0, 0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0,
    0.0, 0.0},
  {-2.0, 2.0, -0.5, 0.0,
    -0.132203419506655614, -0.390689361106274866,
    -0.256343261725123592, 0.106175177622871437,
    0.0749290443999118563, 0.0319207027895210173},
  {-2.0, 2.0, 0.99609375, 0.0,
    13841369.6929212372, -8339597.62269406988,
    -87251051.2148652244, 65225350.0315463319,
    273254167.909661771, -247710863.316563213},
  {-2.0, 2.0, 1.00390625, 0.0,
    -26028.6661639721745, 15961.1559569935978,
    113847.512122984734, -206358.723605581203,
    6578.40068463757013, 908971.828153761815},
  {-2.0, 2.0, -4.25, 0.0,
    1.26205464628714387, 0.274313692011971484,
    0.442242006131113713, -1.61799252235272836,
    -0.927785295033312282, -0.172205994085466244},
  {-2.0, 2.0, 3.5, 0.0,
    0.0870998538268457367, 0.0980699787750225838,
    0.152937181587789832, -0.0211963252018165386,
    0.0431921445114569281, -0.0695183434413184998},
  {-2.0, 2.0, 0.0, -1.0,
    0.199312975007554535, -0.0775592701795263651,
    -0.0210851901426345675, -0.167326332453265841,
    -0.0632997708310078762, -0.0392019115527122943},
  {-2.0, 2.0, 0.125, 0.125,
    0.315234332693229892, 0.165559827416958489,
    -0.16530664867351669, -0.0392852937947218337,
    0.0780220225006149953, 0.0215283250955342792},
};

int main()
{
    slong iter;
    flint_rand_t state;

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

    flint_randinit(state);

    /* check particular values against table */
    {
        acb_t s, z, t;
        acb_ptr w1, w2;
        slong i, j, prec;

        acb_init(s);
        acb_init(z);
        acb_init(t);
        w1 = _acb_vec_init(NUM_DERIVS);
        w2 = _acb_vec_init(NUM_DERIVS);

        for (prec = 32; prec <= 512; prec *= 4)
        {
            for (i = 0; i < NUM_TESTS; i++)
            {
                acb_zero(s);
                arf_set_d(arb_midref(acb_realref(s)), polylog_testdata[i][0]);
                arf_set_d(arb_midref(acb_imagref(s)), polylog_testdata[i][1]);

                acb_zero(z);
                arf_set_d(arb_midref(acb_realref(z)), polylog_testdata[i][2]);
                arf_set_d(arb_midref(acb_imagref(z)), polylog_testdata[i][3]);

                _acb_poly_polylog_cpx_small(w1, s, z, NUM_DERIVS, prec);
                _acb_poly_polylog_cpx_zeta(w2, s, z, NUM_DERIVS, prec);

                for (j = 0; j < NUM_DERIVS; j++)
                {
                    arf_set_d(arb_midref(acb_realref(t)), polylog_testdata[i][4+2*j]);
                    mag_set_d(arb_radref(acb_realref(t)), fabs(polylog_testdata[i][4+2*j]) * EPS);
                    arf_set_d(arb_midref(acb_imagref(t)), polylog_testdata[i][4+2*j+1]);
                    mag_set_d(arb_radref(acb_imagref(t)), fabs(polylog_testdata[i][4+2*j+1]) * EPS);

                    if (!acb_overlaps(w1 + j, t) || !acb_overlaps(w2 + j, t)
                        || !acb_overlaps(w1 + j, w2 + j))
                    {
                        flint_printf("FAIL\n\n");
                        flint_printf("s = "); acb_printd(s, 15); flint_printf("\n\n");
                        flint_printf("z = "); acb_printd(z, 15); flint_printf("\n\n");
                        flint_printf("t = "); acb_printd(t, 15); flint_printf("\n\n");
                        flint_printf("w1 = "); acb_printd(w1 + j, 15); flint_printf("\n\n");
                        flint_printf("w2 = "); acb_printd(w2 + j, 15); flint_printf("\n\n");
                        flint_abort();
                    }
                }
            }
        }

        _acb_vec_clear(w1, NUM_DERIVS);
        _acb_vec_clear(w2, NUM_DERIVS);
        acb_clear(s);
        acb_clear(z);
        acb_clear(t);
    }

    for (iter = 0; iter < 250 * arb_test_multiplier(); iter++)
    {
        acb_t s, z;
        acb_ptr w1, w2;
        slong i, len1, len2, prec1, prec2;

        acb_init(s);
        acb_init(z);

        if (n_randint(state, 2))
            acb_randtest(s, state, 1 + n_randint(state, 300), 3);
        else
            acb_set_si(s, n_randint(state, 20) - 10);

        switch (n_randint(state, 3))
        {
            case 0:
                acb_randtest(z, state, 1 + n_randint(state, 300), 10);
                break;
            case 1:
                arb_randtest(acb_realref(z), state, 1 + n_randint(state, 300), 10);
                break;
            case 2:
                acb_one(z);
                break;
        }

        prec1 = 2 + n_randint(state, 300);
        prec2 = prec1 + 30;
        len1 = 1 + n_randint(state, 20);
        len2 = 1 + n_randint(state, 20);

        w1 = _acb_vec_init(len1);
        w2 = _acb_vec_init(len2);

        switch (n_randint(state, 5))
        {
            case 0:
                _acb_poly_polylog_cpx_small(w1, s, z, len1, prec1);
                break;
            case 1:
                _acb_poly_polylog_cpx_zeta(w1, s, z, len1, prec1);
                break;
            default:
                _acb_poly_polylog_cpx(w1, s, z, len1, prec1);
                break;
        }

        switch (n_randint(state, 5))
        {
            case 0:
                _acb_poly_polylog_cpx_small(w2, s, z, len2, prec2);
                break;
            case 1:
                _acb_poly_polylog_cpx_zeta(w2, s, z, len2, prec2);
                break;
            default:
                _acb_poly_polylog_cpx(w2, s, z, len2, prec2);
                break;
        }

        for (i = 0; i < FLINT_MIN(len1, len2); i++)
        {
            if (!acb_overlaps(w1 + i, w2 + i))
            {
                flint_printf("FAIL: overlap\n\n");
                flint_printf("iter = %wd\n", iter);
                flint_printf("len1 = %wd, len2 = %wd, i = %wd\n\n", len1, len2, i);
                flint_printf("s = "); acb_printd(s, prec1 / 3.33); flint_printf("\n\n");
                flint_printf("z = "); acb_printd(z, prec1 / 3.33); flint_printf("\n\n");
                flint_printf("w1 = "); acb_printd(w1 + i, prec1 / 3.33); flint_printf("\n\n");
                flint_printf("w2 = "); acb_printd(w2 + i, prec2 / 3.33); flint_printf("\n\n");
                flint_abort();
            }
        }

        acb_clear(s);
        acb_clear(z);
        _acb_vec_clear(w1, len1);
        _acb_vec_clear(w2, len2);
    }

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