#include <stdio.h>
#include <stdlib.h>
#include <flint/acb.h>
#include <flint/acb_calc.h>
#include <flint/acb_hypgeom.h>
#include <flint/fmpz.h>

// Open file for writing
//    FILE *f = fopen("LP_LL_comparison_2.txt", "w");
    // if (!f) {
    //     perror("File open failed");
    //     return 1;
    // }

// Define integrand for volume integral bound by curve parameterized by f
int f_func_vol(acb_t res, const acb_t v, void *param, slong order, slong prec)
{
   // Define f(v) = sqrt(r^2 - (v + h)^2)
    acb_t h, r, fval, fpow, vph, vph2, k, l, K, L;
    acb_init(h); acb_init(r); acb_init(fval); acb_init(vph); acb_init(vph2); acb_init(l); acb_init(k); acb_init(fpow);
    acb_init(K); acb_init(L);
    acb_set(h, ((acb_t*)param)[0]);
    acb_set(r, ((acb_t*)param)[1]);
    acb_set(k,((acb_t*)param)[2]);
    acb_set(l, ((acb_t*)param)[3]);
    acb_add_ui(K, k, 1, prec);         // K = k+1
    acb_add_ui(L, l, 1, prec);
    acb_add(vph, v, h, prec);             // v + h
    acb_mul(vph2, vph, vph, prec);        // (v + h)^2
    acb_mul(fval, r, r, prec);       // r^2
    acb_sub(fval, fval, vph2, prec);      // r^2 - (v + h)^2
    acb_sqrt(fval, fval, prec);          // sqrt(...)
    acb_pow(fpow, fval, K, prec);         // f^K
    acb_pow(res, v, l, prec);           // v^l
    acb_mul(res, res, fpow, prec);      // v^l * f^K
    acb_mul(res, res, L, prec);         // v^l * f^K * (l+1)
    acb_clear(h); acb_clear(r); acb_clear(vph); acb_clear(fval); acb_clear(vph2); acb_clear(fpow); acb_clear(K); acb_clear(L);
    return 0;
}

// Define integrand for volume integral bound by curve parameterized by g
int g_func_vol(acb_t res, const acb_t u, void *param, slong order, slong prec)
{
   // Define g(u) = sqrt(rho^2 - (u + d)^2)
    acb_t d, rho, gval, gpow, upd, upd2, k, l, K, L;
    acb_init(d); acb_init(rho); acb_init(gval); acb_init(upd); acb_init(upd2); acb_init(l); acb_init(k); acb_init(gpow);
    acb_init(K); acb_init(L);
    acb_set(d, ((acb_t*)param)[4]);
    acb_set(rho, ((acb_t*)param)[5]);
    acb_set(k,((acb_t*)param)[2]);
    acb_set(l, ((acb_t*)param)[3]);
    acb_add_ui(K, k, 1, prec);         // K = k+1
    acb_add_ui(L, l, 1, prec);
    acb_add(upd, u, d, prec);             // u + d
    acb_mul(upd2, upd, upd, prec);        // (u + d)^2
    acb_mul(gval, rho, rho, prec);       // rho^2
    acb_sub(gval, gval, upd2, prec);      // rho^2 - (u + d)^2
    acb_sqrt(gval, gval, prec);          // sqrt(...)
    acb_pow(gpow, gval, L, prec);         // g^L
    acb_pow(res, u, k, prec);           // u^k
    acb_mul(res, res, gpow, prec);      // u^k * g^L
    acb_mul(res, res, K, prec);         // u^k * g^L * (k+1)
    acb_clear(d); acb_clear(rho); acb_clear(upd); acb_clear(gval); acb_clear(upd2); acb_clear(gpow); acb_clear(K); acb_clear(L);
    return 0;
}

// Define integrand for perimeter integral associated to the curve parameterized by f
int f_func(acb_t res, const acb_t v, void *param, slong order, slong prec)
{
   // Define f(v) = sqrt(r^2 - (v + h)^2)
    acb_t h, r, fval, fpow, vph, vph2, k, l, K;
    acb_init(h); acb_init(r); acb_init(fval); acb_init(vph); acb_init(vph2); acb_init(l); acb_init(k); acb_init(fpow);
    acb_init(K);
    acb_set(h, ((acb_t*)param)[0]);
    acb_set(r, ((acb_t*)param)[1]);
    acb_set(k,((acb_t*)param)[2]);
    acb_set(l, ((acb_t*)param)[3]);
    acb_sub_ui(K, k, 1, prec);         // K = k-1
    acb_add(vph, v, h, prec);             // v + h
    acb_mul(vph2, vph, vph, prec);        // (v + h)^2
    acb_mul(fval, r, r, prec);       // r^2
    acb_sub(fval, fval, vph2, prec);      // r^2 - (v + h)^2
    acb_sqrt(fval, fval, prec);          // sqrt(...)
    acb_pow(fpow, fval, K, prec);         // f^K
    acb_pow(res, v, l, prec);           // v^l
    acb_mul(res, res, fpow, prec);      // v^l * f^K
    acb_mul(res, res, r, prec);       // v^l * f^K * r
 //   acb_mul(res, res, dlength_g, prec); // u^k * g^L * sqrt[...]
    acb_clear(h); acb_clear(r); acb_clear(vph); acb_clear(fval); acb_clear(vph2); acb_clear(fpow); acb_clear(K);
    return 0;
}

// Define integrand for perimeter integral associated to the curve parameterized by g
int g_func(acb_t res, const acb_t u, void *param, slong order, slong prec)
{
    // Define g(u) = sqrt(rho^2 - (u + d)^2)
    acb_t d, rho, gval, gpow, gprime, dlength_g, upd, upd2, k, l, L;
    acb_init(d); acb_init(rho); acb_init(gval); acb_init(upd); acb_init(upd2); acb_init(l); acb_init(k); acb_init(gpow);
    acb_init(gprime); acb_init(dlength_g); acb_init(L);
    acb_set(d, ((acb_t*)param)[4]);
    acb_set(rho, ((acb_t*)param)[5]);
    acb_set(k,((acb_t*)param)[2]);
    acb_set(l, ((acb_t*)param)[3]);
    acb_sub_ui(L, l, 1, prec);         // L = l-1
    acb_add(upd, u, d, prec);             // u + d
    acb_mul(upd2, upd, upd, prec);        // (u + d)^2
    acb_mul(gval, rho, rho, prec);       // rho^2
    acb_sub(gval, gval, upd2, prec);      // rho^2 - (u + d)^2
    acb_sqrt(gval, gval, prec);          // sqrt(...)
    acb_pow(gpow, gval, L, prec);         // g^l
    // Define g'(u) = -(u + d) / sqrt(rho^2 - (u + d)^2)
    // acb_neg(upd, upd);         // - (u + d)
    // acb_div(gprime, upd, gval, prec); // g'(u) = - (u + d) / sqrt(...)
    // // Define length element = Sqrt[1+(g')^2]
    // acb_mul(dlength_g, gprime, gprime, prec);   // (g')^2
    // acb_add_si(dlength_g, dlength_g, 1, prec);  // 1 + (g')^2
    // acb_sqrt(dlength_g, dlength_g, prec);       // sqrt(...)
    // Define final integrand
    acb_pow(res, u, k, prec);           // u^k
    acb_mul(res, res, gpow, prec);      // u^k * g^L
    acb_mul(res, res, rho, prec);       // u^k * g^L * rho
    //   acb_mul(res, res, dlength_g, prec); // u^k * g^L * sqrt[...]
    acb_clear(d); acb_clear(rho); acb_clear(upd); acb_clear(gval); acb_clear(upd2); acb_clear(gpow); acb_clear(gprime); acb_clear(dlength_g);
    acb_clear(L);
    return 0;
}


// Lambda_plane

int main(int argc, char *argv[]) {
    slong prec = 1024;
slong N = (argc > 1) ? atoi(argv[1]) : 4;
slong M = (argc > 2) ? atoi(argv[2]) : 2700;

    acb_t pi, gamma1, gamma2;
    acb_t omega_n1, omega_pow, num, den_base, hyp1, hyp2;
    acb_t term1, term2, sqrt3, sqrt3_div2_pow;
    acb_t temp1, temp2, temp3, log_pi, log_den, den_pow, result, ratio;
    fmpz_t exponent;
    acb_t a, b, c, z;

    acb_init(pi); acb_const_pi(pi, prec);
    acb_init(gamma1); acb_init(gamma2);
    acb_init(omega_n1); acb_init(omega_pow);
    acb_init(num); acb_init(den_base); acb_init(hyp1); acb_init(hyp2);
    acb_init(term1); acb_init(term2); acb_init(sqrt3); acb_init(sqrt3_div2_pow);
    acb_init(temp1); acb_init(temp2); acb_init(temp3); acb_init(log_pi); acb_init(log_den);
    acb_init(den_pow); acb_init(result); acb_init(ratio);

for (slong n = N; n <= M; n++) {

    // omega[n-1] = Pi^((n-1)/2) / Gamma((n-1)/2 + 1)
    acb_set_si(temp1, n - 1);
    acb_div_ui(temp1, temp1, 2, prec);           // (n - 1)/2
    acb_add_ui(temp2, temp1, 1, prec);           // (n - 1)/2 + 1
    acb_gamma(gamma1, temp2, prec);              // Gamma((n-1)/2 + 1)
    acb_log(log_pi, pi, prec);
    acb_mul(temp2, log_pi, temp1, prec);         // log(Pi) * ((n-1)/2)
    acb_exp(omega_n1, temp2, prec);              // Pi^((n-1)/2)
    acb_div(omega_n1, omega_n1, gamma1, prec);   // omega[n-1]
    acb_set_si(temp3, 1);
    acb_div_ui(temp3, temp3, n, prec);
    acb_pow(omega_pow, omega_n1, temp3, prec);      // omega[n-1]^(1/n)

    // numerator
  acb_set_si(temp1, n - 1);
    acb_div_ui(temp1, temp1, 2, prec);           // (n-1)/2
    acb_gamma(gamma1, temp1, prec);              // Gamma((n-1)/2)

    acb_set_si(temp2, n);
    acb_div_ui(temp2, temp2, 2, prec);           // n/2
    acb_gamma(gamma2, temp2, prec);              // Gamma(n/2)

    acb_sqrt(temp1, pi, prec);
    acb_div(temp2, gamma1, gamma2, prec);
    acb_mul(term1, temp1, temp2, prec);          // Sqrt[Pi] * Gamma((n-1)/2)/Gamma(n/2)

    // 2F1[1/2, (3 - n)/2, 3/2, 1/4]
    acb_init(a); acb_init(b); acb_init(c); acb_init(z);
    acb_set_d(a, 0.5);
    acb_set_si(temp1, 3 - n);
    acb_div_ui(b, temp1, 2, prec);
    acb_set_d(c, 1.5);
    acb_set_d(z, 0.25);
    acb_hypgeom_2f1(hyp1, a, b, c, z, 0, prec);

    acb_sub(term1, term1, hyp1, prec);
    acb_set_si(temp1, n - 1);
    acb_mul(term1, term1, temp1, prec);          // full numerator

    acb_set_si(sqrt3, 3);
    acb_sqrt(sqrt3, sqrt3, prec);
    acb_div_ui(sqrt3, sqrt3, 2, prec);
    acb_pow_ui(sqrt3_div2_pow, sqrt3, n - 1, prec);
    acb_sub(num, term1, sqrt3_div2_pow, prec);

    // denominator base
    acb_set_si(temp1, n + 1);
    acb_div_ui(temp1, temp1, 2, prec);
    acb_gamma(gamma1, temp1, prec);              // Gamma((n+1)/2)

    acb_set_si(temp1, n);
    acb_div_ui(temp1, temp1, 2, prec);
    acb_add_ui(temp1, temp1, 1, prec);           // 1 + n/2
    acb_gamma(gamma2, temp1, prec);              // Gamma(1 + n/2)

    acb_sqrt(temp1, pi, prec);
    acb_div(temp2, gamma1, gamma2, prec);
    acb_mul(term2, temp1, temp2, prec);          // Sqrt[Pi] * Gamma((n+1)/2)/Gamma(1 + n/2)

    acb_set_si(temp1, 1 - n);
    acb_div_ui(b, temp1, 2, prec);
    acb_hypgeom_2f1(hyp2, a, b, c, z, 0, prec);
    acb_sub(den_base, term2, hyp2, prec);        // base of denominator

    // Raise den_base to the power (n-1)/n
    acb_log(log_den, den_base, prec);
    acb_set_d(ratio, (double)(n - 1) / n);
    acb_mul(temp1, log_den, ratio, prec);
    acb_exp(den_pow, temp1, prec);

    // Final result
    acb_mul(temp1, omega_pow, num, prec);
    acb_div(result, temp1, den_pow, prec);

    // printf("Lambdaplane[%ld] = ", n);
    // acb_printn(result, 50, 0);
    // printf("\n");

    int k_raw, l_raw;
// Lambda_Lawson
    if (n % 2 == 0) {
        k_raw = (int) n / 2 - 1;
        l_raw = (int) n / 2 - 1;
    }
    else {
        k_raw = (n - 1) / 2 - 1;
        l_raw = (n - 1) / 2;
    }
    // int k_raw = atoi(argv[1]);
    // int l_raw = atoi(argv[2]);
    slong prec = 100;

    acb_t lambda, lambdapow, theta, pi, temp, two_pi_over_3, pi_over_6, OmegaKPlusOne, OmegaLPlusOne, SquareVol;
    acb_t r, h, rho, d, rmd, rmh, k, l;
    acb_t half_exp, gamma_arg, pi_pow, gamma_val;
    acb_init(lambda); acb_init(theta); acb_init(pi); acb_init(temp); acb_init(lambdapow);
    acb_init(two_pi_over_3); acb_init(pi_over_6);
    acb_init(r); acb_init(h); acb_init(rho); acb_init(d); acb_init(rmd); acb_init(k); acb_init(l); acb_init(rmh);
    acb_init(OmegaKPlusOne); acb_init(OmegaLPlusOne); acb_init(SquareVol);
    acb_init(half_exp); acb_init(gamma_arg); acb_init(pi_pow); acb_init(gamma_val);

    // Define k and l from their raw int values
    acb_set_si(k, k_raw);
    acb_set_si(l, l_raw);


    // pi and fixed angles
    acb_const_pi(pi, prec);
    acb_mul_ui(two_pi_over_3, pi, 2, prec);
    acb_div_ui(two_pi_over_3, two_pi_over_3, 3, prec);
    acb_div_ui(pi_over_6, pi, 6, prec);

    // lambda = sqrt(k / l)
    acb_set_si(temp, k_raw);
    acb_div_ui(lambda, temp, l_raw, prec);
    acb_sqrt(lambda, lambda, prec);

    // lambda^k
    acb_pow(lambdapow, lambda, k, prec);

    // theta = arctan(lambda)
    acb_atan(theta, lambda, prec);

    // r = lambda * sec(2pi/3 - theta)
    acb_sub(temp, two_pi_over_3, theta, prec);
    acb_sec(temp, temp, prec);
    acb_mul(r, lambda, temp, prec);

    // h = lambda * tan(2pi/3 - theta) - 1
    acb_sub(temp, two_pi_over_3, theta, prec);
    acb_tan(temp, temp, prec);
    acb_mul(h, lambda, temp, prec);
    acb_sub_si(h, h, 1, prec);

    // rho = sec(pi/6 + theta)
    acb_add(temp, pi_over_6, theta, prec);
    acb_sec(rho, temp, prec);

    // d = tan(pi/6 + theta) - lambda
    acb_add(temp, pi_over_6, theta, prec);
    acb_tan(d, temp, prec);
    acb_sub(d, d, lambda, prec);

    // rmd = rho - d
    acb_sub(rmd, rho, d, prec);

    // rmh = r - h
    acb_sub(rmh, r, h, prec);


    // Compute OmegaKPlusOne
    acb_set_si(half_exp, k_raw + 1);
    acb_div_ui(half_exp, half_exp, 2, prec);              // (k+1)/2
    acb_add_ui(gamma_arg, half_exp, 1, prec);             // gamma_arg = (k+1)/2 + 1
    acb_const_pi(pi_pow, prec);                           // pi
    acb_pow(pi_pow, pi_pow, half_exp, prec);              // pi^{(k+1)/2}
    acb_gamma(gamma_val, gamma_arg, prec);                // Gamma(...)
    acb_div(OmegaKPlusOne, pi_pow, gamma_val, prec);      // Omega_k+1

    // Compute OmegaLPlusOne
    acb_set_si(half_exp, l_raw + 1);
    acb_div_ui(half_exp, half_exp, 2, prec);              // (l+1)/2
    acb_add_ui(gamma_arg, half_exp, 1, prec);             // gamma_arg = (l+1)/2 + 1
    acb_const_pi(pi_pow, prec);                           // pi
    acb_pow(pi_pow, pi_pow, half_exp, prec);              // pi^{(l+1)/2}
    acb_gamma(gamma_val, gamma_arg, prec);                // Gamma(...)
    acb_div(OmegaLPlusOne, pi_pow, gamma_val, prec);      // Omega_l+1

    acb_calc_integrate_opt_t options;
    acb_calc_integrate_opt_init(options);
    mag_t abs_tol;
    mag_init(abs_tol);
    mag_set_ui_2exp_si(abs_tol, 1, -prec);

    // Integrate (l+1) v^l * f^{k+1}
    acb_t f_res2, one;
    acb_init(f_res2); acb_init(one);
    acb_t param[6];
    acb_init(param[0]); acb_init(param[1]); acb_init(param[2]); 
    acb_init(param[3]); acb_init(param[4]); acb_init(param[5]);
    acb_set_ui(one, 1);
    acb_set(param[0],h);
    acb_set(param[1],r);
    acb_set(param[2],k);
    acb_set(param[3],l);
    acb_set(param[4],d);
    acb_set(param[5],rho);
    acb_calc_integrate(f_res2, f_func_vol, param, one, rmh, prec, abs_tol, options, prec);

    // Integrate (k+1) u^k * g^{l+1}
    acb_t g_res2;
    acb_init(g_res2);
    acb_calc_integrate(g_res2, g_func_vol, param, lambda, rmd, prec, abs_tol, options, prec);

    // Integrate v^l f^k sqrt[1+(f')^2] = v^l * f^{k-1} * r
    acb_t f_res;
    acb_init(f_res);
    acb_calc_integrate(f_res, f_func, param, one, rmh, prec, abs_tol, options, prec);  

    // Integrate u^k g^l sqrt[1+(g')^2] = u^k * g^{l-1} * rho
    acb_t g_res;
    acb_init(g_res);
    acb_calc_integrate(g_res, g_func, param, lambda, rmd, prec, abs_tol, options, prec);

    // Compute Vol^{(n-1)/n}
    acb_t vol_pow, dim;
    acb_init(vol_pow); acb_init(dim);
    acb_set_si(dim, n);                     // dim = n
 //   acb_add(dim, k, l, prec);               // k + l
 //   acb_add_ui(dim, dim, 2, prec);          // n = k + l + 2
    acb_sub_ui(vol_pow, dim, 1, prec);      // k + l + 1
    acb_div(vol_pow, vol_pow, dim, prec);   // (n-1)/n = (k + l + 1)/(k + l + 2)
    acb_mul(SquareVol, lambdapow, lambda, prec);    // lambda^{k+1}
    acb_add(SquareVol, SquareVol, g_res2, prec);    // \lambda^{k+1} * g_res2
    acb_add(SquareVol, SquareVol, f_res2, prec);    // \lambda^{k+1} * g_res2 * f_res2
    acb_mul(SquareVol, SquareVol, OmegaKPlusOne, prec);
    acb_mul(SquareVol, SquareVol, OmegaLPlusOne, prec);
    acb_pow(SquareVol, SquareVol, vol_pow, prec);

    // Compute factor towards perimeter of portion of cone that we subtract
    acb_t cone_per, lambda2, temp2;
    acb_init(cone_per); acb_init(lambda2); acb_init(temp2);
    acb_mul(lambda2, lambda, lambda, prec);     // lambda^2
    acb_add_ui(cone_per, lambda2, 1, prec);     // lambda^2 + 1
    acb_sqrt(cone_per, cone_per, prec);         // sqrt(...)
    acb_mul(cone_per, cone_per, lambdapow, prec);     // lambda^k * sqrt(...)
    acb_sub_ui(temp2, dim, 1, prec);            // n - 1
    acb_div(cone_per, cone_per, temp2, prec);   // [lambda^2 * sqrt(...)]/(n-1)

    // Compute Lambda_Lawson
    acb_t kplusone, lplusone, LambdaLawson, diff;
    acb_init(kplusone); acb_init(lplusone); acb_init(LambdaLawson); acb_init(diff);
    acb_add_ui(kplusone, k, 1, prec);           // k + 1
    acb_add_ui(lplusone, l, 1, prec);           // l + 1
    acb_add(LambdaLawson, f_res, g_res, prec);  // adds perimeter integrals from f and g
    acb_sub(LambdaLawson, LambdaLawson, cone_per, prec);        // subtracts perimeter from cone
    acb_mul(LambdaLawson, LambdaLawson, kplusone, prec);        // (...) * (k + 1)
    acb_mul(LambdaLawson, LambdaLawson, lplusone, prec);        // (...) * (k + 1)(l + 1)
    acb_mul(LambdaLawson, LambdaLawson, OmegaKPlusOne, prec);   // (...) * (k + 1)(l + 1) * omega_{k+1}
    acb_mul(LambdaLawson, LambdaLawson, OmegaLPlusOne, prec);   // (...) * (k + 1)(l + 1) * omega_{k+1} * omega_{l+1}
    acb_div(LambdaLawson, LambdaLawson, SquareVol, prec);       // (...) * (k + 1)(l + 1) * omega_{k+1} * omega_{l+1} / Vol^{(n-1)/n}
    acb_sub(diff, result, LambdaLawson, prec);                  // Lambda_plane - Lambda_Lawson
    printf("Lambdaplane(%ld) - M(%d,%d) = ", n, k_raw, l_raw); acb_printn(diff, 25, 0); printf("\n");  

    // Clear memory
    acb_clear(lambda); acb_clear(theta); acb_clear(pi); acb_clear(temp);
    acb_clear(two_pi_over_3); acb_clear(pi_over_6);
    acb_clear(r); acb_clear(h); acb_clear(rho); acb_clear(d);
    acb_clear(half_exp); acb_clear(gamma_arg);
    acb_clear(pi_pow); acb_clear(gamma_val);
    acb_clear(rmh); acb_clear(rmd);
    acb_clear(k); acb_clear(l);
    acb_clear(vol_pow); acb_clear(dim); acb_clear(SquareVol);
    acb_clear(LambdaLawson); acb_clear(kplusone); acb_clear(lplusone); acb_clear(cone_per); acb_clear(temp2);
    }
 //   fclose(f);
    //    Cleanup
    acb_clear(pi); acb_clear(gamma1); acb_clear(gamma2);
    acb_clear(omega_n1); acb_clear(omega_pow);
    acb_clear(num); acb_clear(den_base); acb_clear(hyp1); acb_clear(hyp2);
    acb_clear(term1); acb_clear(term2); acb_clear(sqrt3); acb_clear(sqrt3_div2_pow);
    acb_clear(temp1); acb_clear(temp2); acb_clear(log_pi); acb_clear(log_den);
    acb_clear(den_pow); acb_clear(result); acb_clear(ratio);
    acb_clear(a); acb_clear(b); acb_clear(c); acb_clear(z); acb_clear(temp3);
    return 0;
}