API tools faq
paste
Guest User

Untitled

a guest
Mar 27th, 2018
163
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 45.29 KB | None | 0 0
raw download clone embed print report
  1. /*
  2.     Copyright (C) 2018 Association des collaborateurs de D.H.J Polymath
  3.  
  4.     This is free software: you can redistribute it and/or modify it under
  5.     the terms of the GNU Lesser General Public License (LGPL) as published
  6.     by the Free Software Foundation; either version 2.1 of the License, or
  7.     (at your option) any later version.  See <http://www.gnu.org/licenses/>.
  8. */
  9.  
  10. #include "acb_calc.h"
  11.  
  12.  
  13. void eulerian_polynomial(fmpz_poly_t res, slong len)
  14. {
  15.     if (len)
  16.     {
  17.         slong n, m;
  18.         fmpz *curr, *prev, *tmp;
  19.  
  20.         prev = _fmpz_vec_init(len);
  21.         curr = _fmpz_vec_init(len);
  22.  
  23.         fmpz_one(curr + 0);
  24.         for (n = 1; n < len+1; n++)
  25.         {
  26.             tmp = curr; curr = prev; prev = tmp;
  27.             fmpz_one(curr + 0);
  28.             fmpz_one(curr + n - 1);
  29.             for (m = 1; m < n-1; m++)
  30.             {
  31.                 fmpz_mul_si(curr + m, prev + m - 1, n - m);
  32.                 fmpz_addmul_ui(curr + m, prev + m, m + 1);
  33.             }
  34.         }
  35.  
  36.         fmpz_poly_zero(res);
  37.         for (m = 0; m < len; m++)
  38.         {
  39.             fmpz_poly_set_coeff_fmpz(res, m, curr + m);
  40.         }
  41.  
  42.         _fmpz_vec_clear(prev, len);
  43.         _fmpz_vec_clear(curr, len);
  44.     }
  45.     else
  46.     {
  47.         fmpz_poly_zero(res);
  48.     }
  49. }
  50.  
  51. static void
  52. generalized_geometric_remainder(arb_t res,
  53.         const arb_t c, ulong n0, ulong k, slong prec)
  54. {
  55.     fmpz_poly_t zpoly;
  56.     arb_poly_t poly;
  57.     arb_t cc, s, d;
  58.     fmpz_t bin, n0p;
  59.     slong i, j;
  60.  
  61.     fmpz_poly_init(zpoly);
  62.     arb_poly_init(poly);
  63.  
  64.     arb_init(cc);
  65.     arb_one(cc);
  66.     arb_sub(cc, cc, c, prec);
  67.     arb_inv(cc, cc, prec);
  68.  
  69.     fmpz_init(bin);
  70.     fmpz_init(n0p);
  71.     arb_init(s);
  72.  
  73.     arb_init(d);
  74.     arb_one(d);
  75.  
  76.     arb_zero(res);
  77.     for (i = 0; i < k+1; i++)
  78.     {
  79.         arb_mul(d, d, cc, prec);
  80.  
  81.         fmpz_bin_uiui(bin, k, i);
  82.         fmpz_set_ui(n0p, n0);
  83.         fmpz_pow_ui(n0p, n0p, k - i);
  84.         fmpz_mul(n0p, n0p, bin);
  85.  
  86.         if (i)
  87.         {
  88.             eulerian_polynomial(zpoly, i);
  89.             fmpz_poly_shift_left(zpoly, zpoly, 1);
  90.             arb_poly_set_fmpz_poly(poly, zpoly, prec);
  91.             arb_poly_evaluate(s, poly, c, prec);
  92.         }
  93.         else
  94.         {
  95.             arb_one(s);
  96.         }
  97.         arb_mul_fmpz(s, s, n0p, prec);
  98.         arb_mul(s, s, d, prec);
  99.         arb_add(res, res, s, prec);
  100.     }
  101.  
  102.     arb_poly_clear(poly);
  103.     fmpz_poly_clear(zpoly);
  104.     arb_clear(cc);
  105.     arb_clear(s);
  106.     arb_clear(d);
  107.     fmpz_clear(bin);
  108.     fmpz_clear(n0p);
  109. }
  110.  
  111.  
  112.  
  113.  
  114. static void
  115. _gamma_helper(arb_t res, const arb_t theta, slong n, slong prec)
  116. {
  117.     arb_t beta;
  118.  
  119.     arb_init(beta);
  120.     arb_const_pi(beta, prec);
  121.     arb_mul_si(beta, beta, n*n, prec);
  122.  
  123.     arb_mul_2exp_si(res, theta, 2);
  124.     arb_cos(res, res, prec);
  125.     arb_mul(res, res, beta, prec);
  126.  
  127.     arb_clear(beta);
  128. }
  129.  
  130. static void
  131. _gamma(arb_t res, const arb_t theta, const arb_t X, slong n, slong prec)
  132. {
  133.     arb_t g;
  134.  
  135.     arb_init(g);
  136.     _gamma_helper(g, theta, n, prec);
  137.  
  138.     arb_mul_2exp_si(res, X, 2);
  139.     arb_exp(res, res, prec);
  140.     arb_mul(res, res, g, prec);
  141.  
  142.     arb_clear(g);
  143. }
  144.  
  145. void
  146. integral_remainder_bound_permissibility_lhs(arb_t res,
  147.         const arb_t theta, const arb_t X, slong n, slong prec)
  148. {
  149.     _gamma(res, theta, X, n, prec);
  150. }
  151.  
  152. void
  153. integral_bound_permissibility_lhs(arb_t res,
  154.         const arb_t theta, slong n, slong prec)
  155. {
  156.     _gamma_helper(res, theta, n, prec);
  157. }
  158.  
  159.  
  160. void
  161. integral_remainder_bound_permissibility_rhs(arb_t res,
  162.         const arb_t t, const arb_t a, const arb_t X, slong prec)
  163. {
  164.     arb_t r1, r2;
  165.  
  166.     arb_init(r1);
  167.     arb_init(r2);
  168.  
  169.     arb_mul_2exp_si(r1, t, -1);
  170.  
  171.     arb_mul(r2, t, X, prec);
  172.     arb_mul_2exp_si(r2, r2, 1);
  173.     arb_add(r2, r2, a, prec);
  174.     arb_mul_2exp_si(r2, r2, -2);
  175.  
  176.     arb_max(res, r1, r2, prec);
  177.  
  178.     arb_clear(r1);
  179.     arb_clear(r2);
  180. }
  181.  
  182. void
  183. integral_bound_permissibility_rhs(arb_t res,
  184.         const arb_t t, const arb_t a, slong prec)
  185. {
  186.     arb_t r1, r2;
  187.  
  188.     arb_init(r1);
  189.     arb_mul_2exp_si(r1, t, -1);
  190.  
  191.     arb_init(r2);
  192.     arb_mul_2exp_si(r2, a, -2);
  193.  
  194.     arb_max(res, r1, r2, prec);
  195.  
  196.     arb_clear(r1);
  197.     arb_clear(r2);
  198. }
  199.  
  200. int integral_remainder_bound_is_permissible(const arb_t theta, const arb_t t,
  201.         const arb_t a, slong n, const arb_t X, slong prec)
  202. {
  203.     arb_t lhs, rhs;
  204.     int result;
  205.  
  206.     arb_init(lhs);
  207.     arb_init(rhs);
  208.  
  209.     integral_remainder_bound_permissibility_lhs(lhs, theta, X, n, prec);
  210.     integral_remainder_bound_permissibility_rhs(rhs, t, a, X, prec);
  211.  
  212.     result = arb_gt(lhs, rhs);
  213.  
  214.     arb_clear(lhs);
  215.     arb_clear(rhs);
  216.  
  217.     return result;
  218. }
  219.  
  220.  
  221. int integral_bound_is_permissible(const arb_t theta, const arb_t t,
  222.         const arb_t a, slong n, slong prec)
  223. {
  224.     arb_t lhs, rhs;
  225.     int result;
  226.  
  227.     arb_init(lhs);
  228.     arb_init(rhs);
  229.  
  230.     integral_bound_permissibility_lhs(lhs, theta, n, prec);
  231.     integral_bound_permissibility_rhs(rhs, t, a, prec);
  232.  
  233.     result = arb_gt(lhs, rhs);
  234.  
  235.     arb_clear(lhs);
  236.     arb_clear(rhs);
  237.  
  238.     return result;
  239. }
  240.  
  241. int integral_bounds_are_permissible(
  242.         const arb_t theta, const arb_t t,
  243.         const arb_t y, slong n, slong prec)
  244. {
  245.     arb_t a;
  246.     int result = 1;
  247.  
  248.     arb_init(a);
  249.  
  250.     arb_sub_ui(a, y, 9, prec);
  251.     arb_neg(a, a);
  252.     if (!integral_bound_is_permissible(theta, t, a, n, prec))
  253.     {
  254.         result = 0;
  255.         goto finish;
  256.     }
  257.  
  258.     arb_sub_ui(a, y, 5, prec);
  259.     arb_neg(a, a);
  260.     if (!integral_bound_is_permissible(theta, t, a, n, prec))
  261.     {
  262.         result = 0;
  263.         goto finish;
  264.     }
  265.  
  266.     arb_add_ui(a, y, 9, prec);
  267.     if (!integral_bound_is_permissible(theta, t, a, n, prec))
  268.     {
  269.         result = 0;
  270.         goto finish;
  271.     }
  272.  
  273.     arb_add_ui(a, y, 5, prec);
  274.     if (!integral_bound_is_permissible(theta, t, a, n, prec))
  275.     {
  276.         result = 0;
  277.         goto finish;
  278.     }
  279.  
  280. finish:
  281.  
  282.     arb_clear(a);
  283.     return result;
  284. }
  285.  
  286.  
  287. static void
  288. _integral_remainder_helper(arb_t res1, arb_t res2,
  289.         const arb_t theta, const arb_t t, const arb_t x,
  290.         const arb_t a, slong n, const arb_t X, slong prec)
  291. {
  292.     arb_t z1, z2;
  293.     arb_t w1, w2;
  294.  
  295.     arb_init(w1);
  296.     arb_mul(w1, t, X, prec);
  297.     arb_add(w1, w1, a, prec);
  298.     arb_mul(w1, w1, X, prec);
  299.  
  300.     arb_init(w2);
  301.     arb_mul(w2, t, theta, prec);
  302.     arb_add(w2, w2, x, prec);
  303.     arb_mul(w2, w2, theta, prec);
  304.  
  305.     arb_init(z1);
  306.     _gamma(z1, theta, X, n, prec);
  307.  
  308.     arb_init(z2);
  309.     arb_mul(z2, t, X, prec);
  310.     arb_mul_2exp_si(z2, z2, 1);
  311.     arb_add(z2, z2, a, prec);
  312.  
  313.     arb_sub(res1, w1, z1, prec);
  314.     arb_sub(res1, res1, w2, prec);
  315.     arb_exp(res1, res1, prec);
  316.  
  317.     arb_mul_2exp_si(res2, z1, 2);
  318.     arb_sub(res2, res2, z2, prec);
  319.  
  320.     arb_clear(z1);
  321.     arb_clear(z2);
  322.     arb_clear(w1);
  323.     arb_clear(w2);
  324. }
  325.  
  326. static void
  327. _D_integral_remainder_helper(arb_t res1, arb_t res2,
  328.         const arb_t theta, const arb_t t,
  329.         const arb_t a, slong n, const arb_t X, slong prec)
  330. {
  331.     arb_t z1, z2;
  332.     arb_t w1;
  333.  
  334.     arb_init(w1);
  335.     arb_mul(w1, t, X, prec);
  336.     arb_add(w1, w1, a, prec);
  337.     arb_mul(w1, w1, X, prec);
  338.  
  339.     arb_init(z1);
  340.     _gamma(z1, theta, X, n, prec);
  341.  
  342.     arb_init(z2);
  343.     arb_mul(z2, t, X, prec);
  344.     arb_mul_2exp_si(z2, z2, 1);
  345.     arb_add(z2, z2, a, prec);
  346.  
  347.     arb_sub(res1, w1, z1, prec);
  348.     arb_exp(res1, res1, prec);
  349.  
  350.     arb_mul_2exp_si(res2, z1, 2);
  351.     arb_sub(res2, res2, z2, prec);
  352.  
  353.     arb_clear(z1);
  354.     arb_clear(z2);
  355.     arb_clear(w1);
  356. }
  357.  
  358. static void
  359. _integral_helper(arb_t res1, arb_t res2,
  360.         const arb_t theta, const arb_t t, const arb_t x,
  361.         const arb_t a, slong n, slong prec)
  362. {
  363.     arb_t beta, z1, z2;
  364.  
  365.     arb_init(beta);
  366.     arb_const_pi(beta, prec);
  367.     arb_mul_si(beta, beta, n*n, prec);
  368.  
  369.     arb_init(z1);
  370.     arb_mul_2exp_si(z1, theta, 2);
  371.     arb_cos(z1, z1, prec);
  372.     arb_mul(z1, z1, beta, prec);
  373.  
  374.     arb_init(z2);
  375.     arb_mul(z2, t, theta, prec);
  376.     arb_add(z2, z2, x, prec);
  377.     arb_mul(z2, z2, theta, prec);
  378.  
  379.     arb_add(res1, z1, z2, prec);
  380.     arb_neg(res1, res1);
  381.     arb_exp(res1, res1, prec);
  382.  
  383.     arb_mul_2exp_si(res2, z1, 2);
  384.     arb_sub(res2, res2, a, prec);
  385.  
  386.     arb_clear(beta);
  387.     arb_clear(z1);
  388.     arb_clear(z2);
  389. }
  390.  
  391. static void
  392. _D_integral_helper(arb_t res1, arb_t res2,
  393.         const arb_t theta, const arb_t t,
  394.         const arb_t a, slong n, slong prec)
  395. {
  396.     arb_t beta, z1;
  397.  
  398.     arb_init(beta);
  399.     arb_const_pi(beta, prec);
  400.     arb_mul_si(beta, beta, n*n, prec);
  401.  
  402.     arb_init(z1);
  403.     arb_mul_2exp_si(z1, theta, 2);
  404.     arb_cos(z1, z1, prec);
  405.     arb_mul(z1, z1, beta, prec);
  406.  
  407.     arb_neg(res1, z1);
  408.     arb_exp(res1, res1, prec);
  409.  
  410.     arb_mul_2exp_si(res2, z1, 2);
  411.     arb_sub(res2, res2, a, prec);
  412.  
  413.     arb_clear(beta);
  414.     arb_clear(z1);
  415. }
  416.  
  417.  
  418. void
  419. I_integral_remainder_bound(arb_t res,
  420.         const arb_t theta, const arb_t t, const arb_t x,
  421.         const arb_t a, slong n, const arb_t X, slong prec)
  422. {
  423.     arb_t den;
  424.     arb_init(den);
  425.     _integral_remainder_helper(res, den, theta, t, x, a, n, X, prec);
  426.     arb_div(res, res, den, prec);
  427.     arb_clear(den);
  428. }
  429.  
  430. void I_integral_bound(arb_t res,
  431.         const arb_t theta, const arb_t t, const arb_t x,
  432.         const arb_t a, slong n, slong prec)
  433. {
  434.     arb_t den;
  435.     arb_init(den);
  436.     _integral_helper(res, den, theta, t, x, a, n, prec);
  437.     arb_div(res, res, den, prec);
  438.     arb_clear(den);
  439. }
  440.  
  441. void
  442. D_integral_remainder_bound(arb_t res,
  443.         const arb_t theta, const arb_t t,
  444.         const arb_t a, slong n, const arb_t X, slong prec)
  445. {
  446.     arb_t u, v;
  447.     arb_t q;
  448.  
  449.     arb_init(u);
  450.     arb_init(v);
  451.     _D_integral_remainder_helper(u, v, theta, t, a, n, X, prec);
  452.     arb_inv(v, v, prec);
  453.  
  454.     arb_init(q);
  455.     arb_hypot(q, X, theta, prec);
  456.  
  457.     arb_add(res, q, v, prec);
  458.     arb_mul(res, res, u, prec);
  459.     arb_mul(res, res, v, prec);
  460.  
  461.     arb_clear(u);
  462.     arb_clear(v);
  463.     arb_clear(q);
  464. }
  465.  
  466. void D_integral_bound(arb_t res,
  467.         const arb_t theta, const arb_t t,
  468.         const arb_t a, slong n, slong prec)
  469. {
  470.     arb_t u, v;
  471.  
  472.     arb_init(u);
  473.     arb_init(v);
  474.  
  475.     _D_integral_helper(u, v, theta, t, a, n, prec);
  476.     arb_inv(v, v, prec);
  477.     arb_add(res, theta, v, prec);
  478.     arb_mul(res, res, u, prec);
  479.     arb_mul(res, res, v, prec);
  480.  
  481.     arb_clear(u);
  482.     arb_clear(v);
  483. }
  484.  
  485.  
  486.  
  487.  
  488.  
  489.  
  490.  
  491.  
  492.  
  493.  
  494.  
  495.  
  496. void
  497. _acb_get_rad_ubound_mag(mag_t res, const acb_t z, slong prec)
  498. {
  499.     arf_t err;
  500.     arf_init(err);
  501.     acb_get_rad_ubound_arf(err, z, prec);
  502.     arf_get_mag(res, err);
  503.     arf_clear(err);
  504. }
  505.  
  506.  
  507. typedef struct
  508. {
  509.     arb_t t;
  510.     arb_t x;
  511.     arb_t a;
  512.     slong n;
  513. } IJ_integrand_param_struct;
  514. typedef IJ_integrand_param_struct IJ_integrand_param_t[1];
  515. typedef IJ_integrand_param_struct *IJ_integrand_param_ptr;
  516. typedef const IJ_integrand_param_struct *IJ_integrand_param_srcptr;
  517.  
  518. static void
  519. IJ_integrand_param_init(IJ_integrand_param_t p,
  520.         const arb_t t, const arb_t x, const arb_t a, slong n)
  521. {
  522.     arb_init(p->t);
  523.     arb_init(p->x);
  524.     arb_init(p->a);
  525.     arb_set(p->t, t);
  526.     arb_set(p->x, x);
  527.     arb_set(p->a, a);
  528.     p->n = n;
  529. }
  530.  
  531. static void
  532. IJ_integrand_param_clear(IJ_integrand_param_t p)
  533. {
  534.     arb_clear(p->t);
  535.     arb_clear(p->x);
  536.     arb_clear(p->a);
  537. }
  538.  
  539.  
  540.  
  541. typedef struct
  542. {
  543.     arb_t t;
  544.     arb_t theta;
  545.     arb_t a;
  546.     slong n;
  547. } D_integrand_param_struct;
  548. typedef D_integrand_param_struct D_integrand_param_t[1];
  549. typedef D_integrand_param_struct *D_integrand_param_ptr;
  550. typedef const D_integrand_param_struct *D_integrand_param_srcptr;
  551.  
  552. static void
  553. D_integrand_param_init(D_integrand_param_t p,
  554.         const arb_t t, const arb_t theta, const arb_t a, slong n)
  555. {
  556.     arb_init(p->t);
  557.     arb_init(p->theta);
  558.     arb_init(p->a);
  559.     arb_set(p->t, t);
  560.     arb_set(p->theta, theta);
  561.     arb_set(p->a, a);
  562.     p->n = n;
  563. }
  564.  
  565. static void
  566. D_integrand_param_clear(D_integrand_param_t p)
  567. {
  568.     arb_clear(p->t);
  569.     arb_clear(p->theta);
  570.     arb_clear(p->a);
  571. }
  572.  
  573.  
  574.  
  575.  
  576.  
  577. int IJ_integrand_domain_contains(const acb_t u)
  578. {
  579.     int result;
  580.     acb_t domain;
  581.     arb_ptr x, y;
  582.  
  583.     acb_init(domain);
  584.     x = acb_realref(domain);
  585.     y = acb_imagref(domain);
  586.  
  587.     arb_zero_pm_inf(x);
  588.  
  589.     mag_const_pi_lower(arb_radref(y));
  590.     mag_mul_2exp_si(arb_radref(y), arb_radref(y), -3);
  591.  
  592.     result = acb_contains(domain, u);
  593.  
  594.     acb_clear(domain);
  595.  
  596.     return result;
  597. }
  598.  
  599. static void
  600. I_integrand(acb_t res, const acb_t u,
  601.         const arb_t t, const arb_t x, const arb_t a, slong n, slong prec)
  602. {
  603.     arb_t beta;
  604.     acb_t b, ibu, tuu;
  605.  
  606.     arb_init(beta);
  607.     arb_const_pi(beta, prec);
  608.     arb_mul_si(beta, beta, n*n, prec);
  609.  
  610.     acb_init(b);
  611.     arb_set(acb_realref(b), x);
  612.     arb_neg(acb_imagref(b), a);
  613.  
  614.     acb_init(ibu);
  615.     acb_onei(ibu);
  616.     acb_mul(ibu, ibu, b, prec);
  617.     acb_mul(ibu, ibu, u, prec);
  618.  
  619.     acb_init(tuu);
  620.     acb_sqr(tuu, u, prec);
  621.     acb_mul_arb(tuu, tuu, t, prec);
  622.  
  623.     acb_mul_2exp_si(res, u, 2);
  624.     acb_exp(res, res, prec);
  625.     acb_mul_arb(res, res, beta, prec);
  626.     acb_neg(res, res);
  627.     acb_add(res, res, tuu, prec);
  628.     acb_add(res, res, ibu, prec);
  629.     acb_exp(res, res, prec);
  630.  
  631.     arb_clear(beta);
  632.     acb_clear(b);
  633.     acb_clear(ibu);
  634.     acb_clear(tuu);
  635. }
  636.  
  637.  
  638. static int
  639. f_I_integrand(acb_ptr res, const acb_t u, void * param,
  640.         slong order, slong prec)
  641. {
  642.     IJ_integrand_param_srcptr p;
  643.  
  644.     if (order > 1)
  645.         flint_abort();
  646.  
  647.     if (!IJ_integrand_domain_contains(u))
  648.     {
  649.         acb_indeterminate(res);
  650.         return 0;
  651.     }
  652.  
  653.     p = param;
  654.     I_integrand(res, u, p->t, p->x, p->a, p->n, prec);
  655.  
  656.     return 0;
  657. }
  658.  
  659.  
  660. void
  661. I_proper_integral(acb_t res,
  662.         const arb_t lower_limit,
  663.         const arb_t upper_limit,
  664.         const arb_t theta,
  665.         const mag_t abs_estimate,
  666.         const arb_t t, const arb_t x, const arb_t a, slong n, slong prec)
  667. {
  668.     IJ_integrand_param_t p;
  669.     acb_t a_lim, b_lim;
  670.     acb_calc_integrate_opt_t options;
  671.     int calc_result;
  672.     mag_t abs_tol;
  673.     slong rel_goal = prec;
  674.  
  675.     IJ_integrand_param_init(p, t, x, a, n);
  676.  
  677.     acb_init(a_lim);
  678.     acb_init(b_lim);
  679.     acb_set_arb_arb(a_lim, lower_limit, theta);
  680.     acb_set_arb_arb(b_lim, upper_limit, theta);
  681.  
  682.     mag_init(abs_tol);
  683.     mag_mul_2exp_si(abs_tol, abs_estimate, -prec);
  684.  
  685.     acb_calc_integrate_opt_init(options);
  686.     options->depth_limit = 100*prec;
  687.     options->use_heap = 1;
  688.     options->verbose = 0;
  689.     calc_result = acb_calc_integrate(
  690.             res, f_I_integrand, p, a_lim, b_lim,
  691.             rel_goal, abs_tol, options, prec);
  692.     if (calc_result == ARB_CALC_NO_CONVERGENCE)
  693.     {
  694.         acb_indeterminate(res);
  695.     }
  696.  
  697.     IJ_integrand_param_clear(p);
  698.     acb_clear(a_lim);
  699.     acb_clear(b_lim);
  700.     mag_clear(abs_tol);
  701. }
  702.  
  703.  
  704. void
  705. I_integral(acb_t res,
  706.         const arb_t theta, const mag_t m_estimate,
  707.         const arb_t t, const arb_t x, const arb_t a, slong n, slong prec)
  708. {
  709.     arb_t lower, upper, remainder;
  710.     acb_t value, partial;
  711.     mag_t m_partial_sum, tail_mag_lb, tail_mag_ub;
  712.     slong r, p, q, tmp;
  713.  
  714.     arb_init(lower);
  715.     arb_init(upper);
  716.     arb_init(remainder);
  717.  
  718.     acb_init(value);
  719.     acb_init(partial);
  720.  
  721.     arb_init(remainder);
  722.  
  723.     mag_init(m_partial_sum);
  724.     mag_init(tail_mag_lb);
  725.     mag_init(tail_mag_ub);
  726.  
  727.     r = 0;
  728.     p = 0;
  729.     q = 1;
  730.     while (1)
  731.     {
  732.         tmp = p + q;
  733.         arb_set_ui(lower, r);
  734.         arb_set_ui(upper, r + q);
  735.         arb_mul_2exp_si(lower, lower, -4);
  736.         arb_mul_2exp_si(upper, upper, -4);
  737.        
  738.         r = r + q;
  739.         p = q;
  740.         q = tmp;
  741.  
  742.         I_proper_integral(value, lower, upper,
  743.                 theta, m_estimate, t, x, a, n, prec);
  744.         acb_add(partial, partial, value, prec);
  745.  
  746.         if (integral_remainder_bound_is_permissible(
  747.                     theta, t, a, n, upper, prec))
  748.         {
  749.             I_integral_remainder_bound(
  750.                     remainder, theta, t, x, a, n, upper, prec);
  751.             arb_get_mag(tail_mag_ub, remainder);
  752.  
  753.             _acb_get_rad_ubound_mag(m_partial_sum, partial, prec);
  754.             arb_get_mag_lower(tail_mag_lb, remainder);
  755.  
  756.             if (mag_cmp(tail_mag_lb, m_partial_sum) <= 0)
  757.             {
  758.                 acb_set(res, partial);
  759.                 acb_add_error_mag(res, tail_mag_ub);
  760.                 break;
  761.             }
  762.         }
  763.     }
  764.  
  765.     arb_clear(lower);
  766.     arb_clear(upper);
  767.  
  768.     acb_clear(value);
  769.     acb_clear(partial);
  770.  
  771.     arb_clear(remainder);
  772.  
  773.     mag_clear(m_partial_sum);
  774.     mag_clear(tail_mag_lb);
  775.     mag_clear(tail_mag_ub);
  776. }
  777.  
  778.  
  779.  
  780.  
  781.  
  782.  
  783. static void
  784. D_integrand(acb_t res, const acb_t u,
  785.         const arb_t t, const arb_t theta, const arb_t a, slong n, slong prec)
  786. {
  787.     arb_t beta, c4, r;
  788.     acb_t z1, z2, rc;
  789.  
  790.     arb_init(beta);
  791.     arb_const_pi(beta, prec);
  792.     arb_mul_si(beta, beta, n*n, prec);
  793.  
  794.     acb_init(rc);
  795.     acb_onei(rc);
  796.     acb_mul(rc, rc, u, prec);
  797.     acb_add_arb(rc, rc, theta, prec);
  798.  
  799.     arb_init(r);
  800.     acb_abs(r, rc, prec);
  801.  
  802.     acb_init(z1);
  803.     acb_mul_arb(z1, u, t, prec);
  804.     acb_add_arb(z1, z1, a, prec);
  805.     acb_mul(z1, z1, u, prec);
  806.  
  807.     arb_init(c4);
  808.     arb_mul_2exp_si(c4, theta, 2);
  809.     arb_cos(c4, c4, prec);
  810.  
  811.     acb_init(z2);
  812.     acb_mul_2exp_si(z2, u, 2);
  813.     acb_exp(z2, z2, prec);
  814.     acb_mul_arb(z2, z2, beta, prec);
  815.     acb_mul_arb(z2, z2, c4, prec);
  816.  
  817.     acb_sub(res, z1, z2, prec);
  818.     acb_exp(res, res, prec);
  819.     acb_mul_arb(res, res, r, prec);
  820.  
  821.     arb_clear(beta);
  822.     arb_clear(c4);
  823.     arb_clear(r);
  824.  
  825.     acb_clear(z1);
  826.     acb_clear(z2);
  827.     acb_clear(rc);
  828. }
  829.  
  830. static int
  831. f_D_integrand(acb_ptr res, const acb_t u, void * param,
  832.         slong order, slong prec)
  833. {
  834.     D_integrand_param_srcptr p;
  835.  
  836.     if (order > 1)
  837.         flint_abort();
  838.  
  839.     p = param;
  840.     D_integrand(res, u, p->t, p->theta, p->a, p->n, prec);
  841.  
  842.     return 0;
  843. }
  844.  
  845. void
  846. D_proper_integral(arb_t res,
  847.         const arb_t lower_limit,
  848.         const arb_t upper_limit,
  849.         const arb_t theta,
  850.         const mag_t abs_estimate,
  851.         const arb_t t, const arb_t a, slong n, slong prec)
  852. {
  853.     D_integrand_param_t p;
  854.     acb_t a_lim, b_lim, value;
  855.     acb_calc_integrate_opt_t options;
  856.     int calc_result;
  857.     mag_t abs_tol;
  858.     slong rel_goal = prec;
  859.  
  860.     D_integrand_param_init(p, t, theta, a, n);
  861.  
  862.     acb_init(a_lim);
  863.     acb_init(b_lim);
  864.     acb_set_arb(a_lim, lower_limit);
  865.     acb_set_arb(b_lim, upper_limit);
  866.  
  867.     mag_init(abs_tol);
  868.     mag_mul_2exp_si(abs_tol, abs_estimate, -prec);
  869.  
  870.     acb_calc_integrate_opt_init(options);
  871.     options->depth_limit = 100*prec;
  872.     options->use_heap = 1;
  873.     options->verbose = 0;
  874.  
  875.     acb_init(value);
  876.     calc_result = acb_calc_integrate(
  877.             value, f_D_integrand, p, a_lim, b_lim,
  878.             rel_goal, abs_tol, options, prec);
  879.     if (calc_result == ARB_CALC_NO_CONVERGENCE)
  880.     {
  881.         acb_indeterminate(value);
  882.     }
  883.  
  884.     arb_set(res, acb_realref(value));
  885.  
  886.     D_integrand_param_clear(p);
  887.     acb_clear(a_lim);
  888.     acb_clear(b_lim);
  889.     acb_clear(value);
  890.     mag_clear(abs_tol);
  891. }
  892.  
  893. void
  894. D_integral(arb_t res,
  895.         const arb_t theta, const mag_t m_estimate,
  896.         const arb_t t, const arb_t a, slong n, slong prec)
  897. {
  898.     arb_t lower, upper, remainder;
  899.     arb_t value, partial;
  900.     mag_t m_partial_sum, tail_mag_lb, tail_mag_ub;
  901.     slong r, p, q, tmp;
  902.  
  903.     arb_init(lower);
  904.     arb_init(upper);
  905.     arb_init(remainder);
  906.  
  907.     arb_init(value);
  908.     arb_init(partial);
  909.  
  910.     arb_init(remainder);
  911.  
  912.     mag_init(m_partial_sum);
  913.     mag_init(tail_mag_lb);
  914.     mag_init(tail_mag_ub);
  915.  
  916.     r = 0;
  917.     p = 0;
  918.     q = 1;
  919.     while (1)
  920.     {
  921.         tmp = p + q;
  922.         arb_set_ui(lower, r);
  923.         arb_set_ui(upper, r + q);
  924.         arb_mul_2exp_si(lower, lower, -4);
  925.         arb_mul_2exp_si(upper, upper, -4);
  926.        
  927.         r = r + q;
  928.         p = q;
  929.         q = tmp;
  930.  
  931.         D_proper_integral(value, lower, upper,
  932.                 theta, m_estimate, t, a, n, prec);
  933.         arb_add(partial, partial, value, prec);
  934.  
  935.         if (integral_remainder_bound_is_permissible(
  936.                     theta, t, a, n, upper, prec))
  937.         {
  938.             D_integral_remainder_bound(
  939.                     remainder, theta, t, a, n, upper, prec);
  940.             arb_get_mag(tail_mag_ub, remainder);
  941.  
  942.             mag_set(m_partial_sum, arb_radref(partial));
  943.             arb_get_mag_lower(tail_mag_lb, remainder);
  944.  
  945.             if (mag_cmp(tail_mag_lb, m_partial_sum) <= 0)
  946.             {
  947.                 arb_set(res, partial);
  948.                 arb_add_error_mag(res, tail_mag_ub);
  949.                 break;
  950.             }
  951.         }
  952.     }
  953.  
  954.     arb_clear(lower);
  955.     arb_clear(upper);
  956.  
  957.     arb_clear(value);
  958.     arb_clear(partial);
  959.  
  960.     arb_clear(remainder);
  961.  
  962.     mag_clear(m_partial_sum);
  963.     mag_clear(tail_mag_lb);
  964.     mag_clear(tail_mag_ub);
  965. }
  966.  
  967.  
  968.  
  969.  
  970.  
  971.  
  972.  
  973. void
  974. H_approximation(arb_t res, const arb_t x, const arb_t y, slong prec)
  975. {
  976.     if (arb_contains_nonpositive(x))
  977.     {
  978.         arb_one(res);
  979.     }
  980.     else
  981.     {
  982.         arb_t logx, pixd8;
  983.  
  984.         arb_init(logx);
  985.         arb_log(logx, x, prec);
  986.  
  987.         arb_init(pixd8);
  988.         arb_const_pi(pixd8, prec);
  989.         arb_mul(pixd8, pixd8, x, prec);
  990.         arb_mul_2exp_si(pixd8, pixd8, -3);
  991.  
  992.         arb_add_si(res, y, 7, prec);
  993.         arb_mul(res, res, logx, prec);
  994.         arb_mul_2exp_si(res, res, -2);
  995.         arb_sub(res, res, pixd8, prec);
  996.         arb_exp(res, res, prec);
  997.  
  998.         arb_clear(logx);
  999.         arb_clear(pixd8);
  1000.     }
  1001. }
  1002.  
  1003.  
  1004. void
  1005. get_theta(arb_t res, const acb_t z, slong prec)
  1006. {
  1007.     arb_t pid8;
  1008.  
  1009.     arb_init(pid8);
  1010.     arb_const_pi(pid8, prec);
  1011.     arb_mul_2exp_si(pid8, pid8, -3);
  1012.  
  1013.     arb_add_si(res, acb_imagref(z), 9, prec);
  1014.     arb_atan2(res, res, acb_realref(z), prec);
  1015.     arb_mul_2exp_si(res, res, -2);
  1016.     arb_sub(res, pid8, res, prec);
  1017.  
  1018.     arb_clear(pid8);
  1019. }
  1020.  
  1021.  
  1022.  
  1023. static void
  1024. H_summand(acb_t res,
  1025.         const arb_t theta, mag_t m_estimate,
  1026.         const arb_t t, const arb_t x, const arb_t y, slong n, slong prec)
  1027. {
  1028.     arb_t a;
  1029.     acb_t s, total;
  1030.     arb_t beta, c1, c2;
  1031.  
  1032.     arb_init(a);
  1033.     acb_init(s);
  1034.     acb_init(total);
  1035.  
  1036.     arb_init(beta);
  1037.     arb_const_pi(beta, prec);
  1038.     arb_mul_si(beta, beta, n*n, prec);
  1039.  
  1040.     arb_init(c1);
  1041.     arb_sqr(c1, beta, prec);
  1042.  
  1043.     arb_init(c2);
  1044.     arb_mul_si(c2, beta, 3, prec);
  1045.     arb_mul_2exp_si(c2, c2, -1);
  1046.  
  1047.     arb_sub_ui(a, y, 9, prec);
  1048.     arb_neg(a, a);
  1049.     I_integral(s, theta, m_estimate, t, x, a, n, prec);
  1050.     acb_mul_arb(s, s, c1, prec);
  1051.     acb_add(total, total, s, prec);
  1052.  
  1053.     arb_sub_ui(a, y, 5, prec);
  1054.     arb_neg(a, a);
  1055.     I_integral(s, theta, m_estimate, t, x, a, n, prec);
  1056.     acb_mul_arb(s, s, c2, prec);
  1057.     acb_sub(total, total, s, prec);
  1058.  
  1059.     arb_add_ui(a, y, 9, prec);
  1060.     I_integral(s, theta, m_estimate, t, x, a, n, prec);
  1061.     acb_conj(s, s);
  1062.     acb_mul_arb(s, s, c1, prec);
  1063.     acb_add(total, total, s, prec);
  1064.  
  1065.     arb_add_ui(a, y, 5, prec);
  1066.     I_integral(s, theta, m_estimate, t, x, a, n, prec);
  1067.     acb_conj(s, s);
  1068.     acb_mul_arb(s, s, c2, prec);
  1069.     acb_sub(total, total, s, prec);
  1070.  
  1071.     acb_set(res, total);
  1072.  
  1073.     arb_clear(a);
  1074.     acb_clear(s);
  1075.     acb_clear(total);
  1076.  
  1077.     arb_clear(beta);
  1078.     arb_clear(c1);
  1079.     arb_clear(c2);
  1080. }
  1081.  
  1082. static void
  1083. _summation_alpha(arb_t res, const arb_t theta, slong n0, slong prec)
  1084. {
  1085.     arb_t c4;
  1086.  
  1087.     arb_init(c4);
  1088.     arb_mul_2exp_si(c4, theta, 2);
  1089.     arb_cos(c4, c4, prec);
  1090.  
  1091.     arb_const_pi(res, prec);
  1092.     arb_mul(res, res, c4, prec);
  1093.     arb_mul_si(res, res, n0, prec);
  1094.     arb_neg(res, res);
  1095.     arb_exp(res, res, prec);
  1096.  
  1097.     arb_clear(c4);
  1098. }
  1099.  
  1100. static void
  1101. I_summation_remainder_bound(arb_t res,
  1102.         const arb_t theta, const arb_t t, const arb_t x, const arb_t y,
  1103.         slong n0, slong prec)
  1104. {
  1105.     arb_t pi, c1, c2, g2, g4;
  1106.     arb_t alpha, a;
  1107.     arb_t s, total;
  1108.  
  1109.     arb_init(alpha);
  1110.     _summation_alpha(alpha, theta, n0, prec);
  1111.  
  1112.     arb_init(g2);
  1113.     generalized_geometric_remainder(g2, alpha, n0, 2, prec);
  1114.  
  1115.     arb_init(g4);
  1116.     generalized_geometric_remainder(g4, alpha, n0, 4, prec);
  1117.  
  1118.     arb_init(pi);
  1119.     arb_const_pi(pi, prec);
  1120.  
  1121.     arb_init(c1);
  1122.     arb_sqr(c1, pi, prec);
  1123.     arb_mul(c1, c1, g4, prec);
  1124.  
  1125.     arb_init(c2);
  1126.     arb_mul_si(c2, pi, 3, prec);
  1127.     arb_mul_2exp_si(c2, c2, -1);
  1128.     arb_mul(c2, c2, g2, prec);
  1129.  
  1130.     arb_init(a);
  1131.     arb_init(total);
  1132.     arb_init(s);
  1133.  
  1134.     arb_sub_ui(a, y, 9, prec);
  1135.     arb_neg(a, a);
  1136.     I_integral_bound(s, theta, t, x, a, n0, prec);
  1137.     arb_mul(s, s, c1, prec);
  1138.     arb_add(total, total, s, prec);
  1139.  
  1140.     arb_sub_ui(a, y, 5, prec);
  1141.     arb_neg(a, a);
  1142.     I_integral_bound(s, theta, t, x, a, n0, prec);
  1143.     arb_mul(s, s, c2, prec);
  1144.     arb_add(total, total, s, prec);
  1145.  
  1146.     arb_add_ui(a, y, 9, prec);
  1147.     I_integral_bound(s, theta, t, x, a, n0, prec);
  1148.     arb_mul(s, s, c1, prec);
  1149.     arb_add(total, total, s, prec);
  1150.  
  1151.     arb_add_ui(a, y, 5, prec);
  1152.     I_integral_bound(s, theta, t, x, a, n0, prec);
  1153.     arb_mul(s, s, c2, prec);
  1154.     arb_add(total, total, s, prec);
  1155.  
  1156.     arb_set(res, total);
  1157.  
  1158.     arb_clear(pi);
  1159.     arb_clear(c1);
  1160.     arb_clear(c2);
  1161.     arb_clear(g2);
  1162.     arb_clear(g4);
  1163.  
  1164.     arb_clear(alpha);
  1165.     arb_clear(a);
  1166.     arb_clear(s);
  1167.     arb_clear(total);
  1168. }
  1169.  
  1170.  
  1171. void H(acb_t res,
  1172.         const arb_t theta, mag_t m_estimate,
  1173.         const arb_t t, const arb_t x, const arb_t y, slong prec)
  1174. {
  1175.     acb_t h, htotal;
  1176.     arb_t remainder;
  1177.     mag_t m_partial_sum, tail_mag_lb, tail_mag_ub;
  1178.     slong n;
  1179.  
  1180.     acb_init(h);
  1181.     acb_init(htotal);
  1182.  
  1183.     arb_init(remainder);
  1184.  
  1185.     mag_init(m_partial_sum);
  1186.     mag_init(tail_mag_lb);
  1187.     mag_init(tail_mag_ub);
  1188.  
  1189.     n = 1;
  1190.     while (1)
  1191.     {
  1192.         H_summand(h, theta, m_estimate, t, x, y, n, prec);
  1193.         acb_add(htotal, htotal, h, prec);
  1194.  
  1195.         n++;
  1196.  
  1197.         if (integral_bounds_are_permissible(theta, t, y, n, prec))
  1198.         {
  1199.             I_summation_remainder_bound(remainder, theta, t, x, y, n, prec);
  1200.             arb_get_mag(tail_mag_ub, remainder);
  1201.  
  1202.             _acb_get_rad_ubound_mag(m_partial_sum, htotal, prec);
  1203.             arb_get_mag_lower(tail_mag_lb, remainder);
  1204.  
  1205.             if (mag_cmp(tail_mag_lb, m_partial_sum) <= 0)
  1206.             {
  1207.                 acb_set(res, htotal);
  1208.                 acb_add_error_mag(res, tail_mag_ub);
  1209.                 break;
  1210.             }
  1211.         }
  1212.     }
  1213.  
  1214.     if (arb_is_zero(x) || arb_is_zero(y))
  1215.     {
  1216.         arb_zero(acb_imagref(res));
  1217.     }
  1218.  
  1219.     acb_clear(h);
  1220.     acb_clear(htotal);
  1221.  
  1222.     arb_clear(remainder);
  1223.  
  1224.     mag_clear(m_partial_sum);
  1225.     mag_clear(tail_mag_lb);
  1226.     mag_clear(tail_mag_ub);
  1227. }
  1228.  
  1229.  
  1230.  
  1231. static void
  1232. D_summation_remainder_bound(arb_t res,
  1233.         const arb_t theta, const arb_t t, const arb_t y,
  1234.         slong n0, slong prec)
  1235. {
  1236.     arb_t pi, c1, c2, g2, g4;
  1237.     arb_t alpha, a;
  1238.     arb_t s, total;
  1239.  
  1240.     arb_init(alpha);
  1241.     _summation_alpha(alpha, theta, n0, prec);
  1242.  
  1243.     arb_init(g2);
  1244.     generalized_geometric_remainder(g2, alpha, n0, 2, prec);
  1245.  
  1246.     arb_init(g4);
  1247.     generalized_geometric_remainder(g4, alpha, n0, 4, prec);
  1248.  
  1249.     arb_init(pi);
  1250.     arb_const_pi(pi, prec);
  1251.  
  1252.     arb_init(c1);
  1253.     arb_sqr(c1, pi, prec);
  1254.     arb_mul(c1, c1, g4, prec);
  1255.  
  1256.     arb_init(c2);
  1257.     arb_mul_si(c2, pi, 3, prec);
  1258.     arb_mul_2exp_si(c2, c2, -1);
  1259.     arb_mul(c2, c2, g2, prec);
  1260.  
  1261.     arb_init(a);
  1262.     arb_init(total);
  1263.     arb_init(s);
  1264.  
  1265.     arb_sub_ui(a, y, 9, prec);
  1266.     arb_neg(a, a);
  1267.     D_integral_bound(s, theta, t, a, n0, prec);
  1268.     arb_mul(s, s, c1, prec);
  1269.     arb_add(total, total, s, prec);
  1270.  
  1271.     arb_sub_ui(a, y, 5, prec);
  1272.     arb_neg(a, a);
  1273.     D_integral_bound(s, theta, t, a, n0, prec);
  1274.     arb_mul(s, s, c2, prec);
  1275.     arb_add(total, total, s, prec);
  1276.  
  1277.     arb_add_ui(a, y, 9, prec);
  1278.     D_integral_bound(s, theta, t, a, n0, prec);
  1279.     arb_mul(s, s, c1, prec);
  1280.     arb_add(total, total, s, prec);
  1281.  
  1282.     arb_add_ui(a, y, 5, prec);
  1283.     D_integral_bound(s, theta, t, a, n0, prec);
  1284.     arb_mul(s, s, c2, prec);
  1285.     arb_add(total, total, s, prec);
  1286.  
  1287.     arb_set(res, total);
  1288.  
  1289.     arb_clear(pi);
  1290.     arb_clear(c1);
  1291.     arb_clear(c2);
  1292.     arb_clear(g2);
  1293.     arb_clear(g4);
  1294.  
  1295.     arb_clear(alpha);
  1296.     arb_clear(a);
  1297.     arb_clear(s);
  1298.     arb_clear(total);
  1299. }
  1300.  
  1301. static void
  1302. D_summand(arb_t res,
  1303.         const arb_t theta, mag_t m_estimate,
  1304.         const arb_t t, const arb_t y, slong n, slong prec)
  1305. {
  1306.     arb_t a;
  1307.     arb_t s, total;
  1308.     arb_t beta, c1, c2;
  1309.  
  1310.     arb_init(a);
  1311.     arb_init(s);
  1312.     arb_init(total);
  1313.  
  1314.     arb_init(beta);
  1315.     arb_const_pi(beta, prec);
  1316.     arb_mul_si(beta, beta, n*n, prec);
  1317.  
  1318.     arb_init(c1);
  1319.     arb_sqr(c1, beta, prec);
  1320.  
  1321.     arb_init(c2);
  1322.     arb_mul_si(c2, beta, 3, prec);
  1323.     arb_mul_2exp_si(c2, c2, -1);
  1324.  
  1325.     arb_sub_ui(a, y, 9, prec);
  1326.     arb_neg(a, a);
  1327.     D_integral(s, theta, m_estimate, t, a, n, prec);
  1328.     arb_mul(s, s, c1, prec);
  1329.     arb_add(total, total, s, prec);
  1330.  
  1331.     arb_sub_ui(a, y, 5, prec);
  1332.     arb_neg(a, a);
  1333.     D_integral(s, theta, m_estimate, t, a, n, prec);
  1334.     arb_mul(s, s, c2, prec);
  1335.     arb_add(total, total, s, prec);
  1336.  
  1337.     arb_add_ui(a, y, 9, prec);
  1338.     D_integral(s, theta, m_estimate, t, a, n, prec);
  1339.     arb_mul(s, s, c1, prec);
  1340.     arb_add(total, total, s, prec);
  1341.  
  1342.     arb_add_ui(a, y, 5, prec);
  1343.     D_integral(s, theta, m_estimate, t, a, n, prec);
  1344.     arb_mul(s, s, c2, prec);
  1345.     arb_add(total, total, s, prec);
  1346.  
  1347.     arb_set(res, total);
  1348.  
  1349.     arb_clear(a);
  1350.     arb_clear(s);
  1351.     arb_clear(total);
  1352.  
  1353.     arb_clear(beta);
  1354.     arb_clear(c1);
  1355.     arb_clear(c2);
  1356. }
  1357.  
  1358.  
  1359. void
  1360. D(arb_t res,
  1361.         const arb_t theta, mag_t m_estimate,
  1362.         const arb_t t, const arb_t y, slong prec)
  1363. {
  1364.     arb_t d, dtotal;
  1365.     arb_t remainder;
  1366.     mag_t m_partial_sum, tail_mag_lb, tail_mag_ub;
  1367.     slong n;
  1368.  
  1369.     mag_init(m_partial_sum);
  1370.     mag_init(tail_mag_lb);
  1371.     mag_init(tail_mag_ub);
  1372.  
  1373.     arb_init(d);
  1374.     arb_init(dtotal);
  1375.     arb_init(remainder);
  1376.  
  1377.     n = 1;
  1378.     while (1)
  1379.     {
  1380.         D_summand(d, theta, m_estimate, t, y, n, prec);
  1381.         arb_add(dtotal, dtotal, d, prec);
  1382.  
  1383.         n++;
  1384.  
  1385.         if (integral_bounds_are_permissible(theta, t, y, n, prec))
  1386.         {
  1387.             D_summation_remainder_bound(remainder, theta, t, y, n, prec);
  1388.             arb_get_mag(tail_mag_ub, remainder);
  1389.  
  1390.             mag_set(m_partial_sum, arb_radref(dtotal));
  1391.             arb_get_mag_lower(tail_mag_lb, remainder);
  1392.  
  1393.             if (mag_cmp(tail_mag_lb, m_partial_sum) <= 0)
  1394.             {
  1395.                 arb_set(res, dtotal);
  1396.                 arb_add_error_mag(res, tail_mag_ub);
  1397.                 break;
  1398.             }
  1399.         }
  1400.     }
  1401.  
  1402.     arb_clear(d);
  1403.     arb_clear(dtotal);
  1404.     arb_clear(remainder);
  1405.  
  1406.     mag_clear(m_partial_sum);
  1407.     mag_clear(tail_mag_lb);
  1408.     mag_clear(tail_mag_ub);
  1409. }
  1410.  
  1411.  
  1412.  
  1413.  
  1414.  
  1415.  
  1416.  
  1417.  
  1418.  
  1419. typedef struct
  1420. {
  1421.     acb_struct z;
  1422.     acb_struct h;
  1423.     arb_struct d;
  1424.     arb_struct theta;
  1425. } xmeshpt_struct;
  1426. typedef xmeshpt_struct xmeshpt_t[1];
  1427. typedef xmeshpt_struct *xmeshpt_ptr;
  1428. typedef const xmeshpt_struct *xmeshpt_srcptr;
  1429.  
  1430. #define xmeshpt_zref(p) (&(p)->z)
  1431. #define xmeshpt_href(p) (&(p)->h)
  1432. #define xmeshpt_dref(p) (&(p)->d)
  1433. #define xmeshpt_thetaref(p) (&(p)->theta)
  1434.  
  1435. void xmeshpt_init(xmeshpt_t p)
  1436. {
  1437.     acb_init(xmeshpt_zref(p));
  1438.     acb_init(xmeshpt_href(p));
  1439.     arb_init(xmeshpt_dref(p));
  1440.     arb_init(xmeshpt_thetaref(p));
  1441. }
  1442.  
  1443. void xmeshpt_clear(xmeshpt_t p)
  1444. {
  1445.     acb_clear(xmeshpt_zref(p));
  1446.     acb_clear(xmeshpt_href(p));
  1447.     arb_clear(xmeshpt_dref(p));
  1448.     arb_clear(xmeshpt_thetaref(p));
  1449. }
  1450.  
  1451. void xmeshpt_swap(xmeshpt_t a, xmeshpt_t b)
  1452. {
  1453.     acb_swap(xmeshpt_zref(a), xmeshpt_zref(b));
  1454.     acb_swap(xmeshpt_href(a), xmeshpt_href(b));
  1455.     arb_swap(xmeshpt_dref(a), xmeshpt_dref(b));
  1456.     arb_swap(xmeshpt_thetaref(a), xmeshpt_thetaref(b));
  1457. }
  1458.  
  1459.  
  1460.  
  1461. typedef struct
  1462. {
  1463.     acb_struct z;
  1464.     acb_struct h;
  1465. } ymeshpt_struct;
  1466. typedef ymeshpt_struct ymeshpt_t[1];
  1467. typedef ymeshpt_struct *ymeshpt_ptr;
  1468. typedef const ymeshpt_struct *ymeshpt_srcptr;
  1469.  
  1470. #define ymeshpt_zref(p) (&(p)->z)
  1471. #define ymeshpt_href(p) (&(p)->h)
  1472.  
  1473. void ymeshpt_init(ymeshpt_t p)
  1474. {
  1475.     acb_init(ymeshpt_zref(p));
  1476.     acb_init(ymeshpt_href(p));
  1477. }
  1478.  
  1479. void ymeshpt_clear(ymeshpt_t p)
  1480. {
  1481.     acb_clear(ymeshpt_zref(p));
  1482.     acb_clear(ymeshpt_href(p));
  1483. }
  1484.  
  1485. void ymeshpt_swap(ymeshpt_t a, ymeshpt_t b)
  1486. {
  1487.     acb_swap(ymeshpt_zref(a), ymeshpt_zref(b));
  1488.     acb_swap(ymeshpt_href(a), ymeshpt_href(b));
  1489. }
  1490.  
  1491.  
  1492.  
  1493. typedef struct
  1494. {
  1495.     ymeshpt_struct left;
  1496.     ymeshpt_srcptr right;
  1497. } ymeshival_struct;
  1498. typedef ymeshival_struct ymeshival_t[1];
  1499. typedef ymeshival_struct *ymeshival_ptr;
  1500. typedef const ymeshival_struct *ymeshival_srcptr;
  1501.  
  1502. #define ymeshival_leftref(m) (&(m)->left)
  1503.  
  1504. void ymeshival_init(ymeshival_t m)
  1505. {
  1506.     ymeshpt_init(ymeshival_leftref(m));
  1507.     m->right = NULL;
  1508. }
  1509.  
  1510. void ymeshival_clear(ymeshival_t m)
  1511. {
  1512.     ymeshpt_clear(ymeshival_leftref(m));
  1513.     m->right = NULL;
  1514. }
  1515.  
  1516. void ymeshival_center(acb_t res, const ymeshival_t m, slong prec)
  1517. {
  1518.     acb_srcptr z_left = ymeshpt_zref(ymeshival_leftref(m));
  1519.     acb_srcptr z_right = ymeshpt_zref(m->right);
  1520.     acb_add(res, z_left, z_right, prec);
  1521.     acb_mul_2exp_si(res, res, -1);
  1522. }
  1523.  
  1524.  
  1525. int ymeshival_validate(const ymeshival_t m, const arb_t t, slong prec)
  1526. {
  1527.     acb_t z, center;
  1528.     arb_t theta;
  1529.     arb_t estimate;
  1530.     mag_t m_estimate;
  1531.     acb_t step_delta;
  1532.     arb_t step_width;
  1533.     arb_t u, d;
  1534.     arb_t lhs, rhs;
  1535.     acb_srcptr z_left, z_right;
  1536.     acb_srcptr h_left, h_right;
  1537.     arb_t abs_h_left, abs_h_right;
  1538.     arb_srcptr x, y;
  1539.     int result;
  1540.  
  1541.     z_left = ymeshpt_zref(ymeshival_leftref(m));
  1542.     z_right = ymeshpt_zref(m->right);
  1543.  
  1544.     h_left = ymeshpt_href(ymeshival_leftref(m));
  1545.     h_right = ymeshpt_href(m->right);
  1546.  
  1547.     acb_init(center);
  1548.     ymeshival_center(center, m, prec);
  1549.  
  1550.     arb_init(theta);
  1551.     get_theta(theta, center, prec);
  1552.     mag_zero(arb_radref(theta));
  1553.  
  1554.     acb_init(step_delta);
  1555.     acb_sub(step_delta, z_right, z_left, prec);
  1556.  
  1557.     arb_init(step_width);
  1558.     acb_abs(step_width, step_delta, prec);
  1559.  
  1560.     arb_init(estimate);
  1561.     mag_init(m_estimate);
  1562.     H_approximation(estimate, acb_realref(center), acb_imagref(center), prec);
  1563.     arb_get_mag(m_estimate, estimate);
  1564.  
  1565.     acb_init(z);
  1566.     acb_union(z, z_left, z_right, prec);
  1567.     x = acb_realref(z);
  1568.     y = acb_imagref(z);
  1569.  
  1570.     arb_init(u);
  1571.     arb_init(d);
  1572.     D(d, theta, m_estimate, t, y, prec);
  1573.     arb_sqr(u, theta, prec);
  1574.     arb_mul(u, u, t, prec);
  1575.     arb_neg(u, u);
  1576.     arb_exp(u, u, prec);
  1577.     arb_mul(d, d, u, prec);
  1578.  
  1579.     arb_init(lhs);
  1580.     arb_mul(lhs, x, theta, prec);
  1581.     arb_neg(lhs, lhs);
  1582.     arb_exp(lhs, lhs, prec);
  1583.     arb_mul(lhs, lhs, d, prec);
  1584.     arb_mul(lhs, lhs, step_width, prec);
  1585.  
  1586.     arb_init(rhs);
  1587.     arb_init(abs_h_left);
  1588.     arb_init(abs_h_right);
  1589.     acb_abs(abs_h_left, h_left, prec);
  1590.     acb_abs(abs_h_right, h_right, prec);
  1591.     arb_max(rhs, abs_h_left, abs_h_right, prec);
  1592.  
  1593.     result = arb_lt(lhs, rhs);
  1594.  
  1595.     acb_clear(z);
  1596.     acb_clear(center);
  1597.     arb_clear(theta);
  1598.  
  1599.     acb_clear(step_delta);
  1600.     arb_clear(step_width);
  1601.  
  1602.     arb_clear(u);
  1603.     arb_clear(d);
  1604.  
  1605.     arb_clear(estimate);
  1606.     mag_clear(m_estimate);
  1607.  
  1608.     arb_clear(lhs);
  1609.     arb_clear(rhs);
  1610.  
  1611.     arb_clear(abs_h_left);
  1612.     arb_clear(abs_h_right);
  1613.  
  1614.     return result;
  1615. }
  1616.  
  1617.  
  1618. void
  1619. evaluate_ymeshpt(ymeshpt_t m,
  1620.         const arb_t t, const arb_t x, const arb_t y,
  1621.         const arb_t htheta_max, slong prec)
  1622. {
  1623.     arb_t htheta;
  1624.     arb_t estimate;
  1625.     mag_t m_estimate;
  1626.  
  1627.     acb_ptr z = ymeshpt_zref(m);
  1628.     acb_ptr h = ymeshpt_href(m);
  1629.  
  1630.     acb_set_arb_arb(z, x, y);
  1631.  
  1632.     arb_init(htheta);
  1633.     get_theta(htheta, z, prec);
  1634.     arb_min(htheta, htheta, htheta_max, prec);
  1635.     mag_zero(arb_radref(htheta));
  1636.  
  1637.     arb_init(estimate);
  1638.     mag_init(m_estimate);
  1639.     H_approximation(estimate, x, y, prec);
  1640.     arb_get_mag(m_estimate, estimate);
  1641.  
  1642.     H(h, htheta, m_estimate, t, x, y, prec);
  1643.  
  1644.     arb_clear(htheta);
  1645.     arb_clear(estimate);
  1646.     mag_clear(m_estimate);
  1647. }
  1648.  
  1649.  
  1650.  
  1651. void
  1652. evaluate_xmeshpt(xmeshpt_t m,
  1653.         const arb_t t, const arb_t x, const arb_t y,
  1654.         const arb_t htheta_max, slong prec)
  1655. {
  1656.     arb_t htheta, u;
  1657.     arb_t estimate;
  1658.     mag_t m_estimate;
  1659.  
  1660.     acb_ptr z = xmeshpt_zref(m);
  1661.     acb_ptr h = xmeshpt_href(m);
  1662.     arb_ptr d = xmeshpt_dref(m);
  1663.     arb_ptr theta = xmeshpt_thetaref(m);
  1664.  
  1665.     acb_set_arb_arb(z, x, y);
  1666.  
  1667.     get_theta(theta, z, prec);
  1668.  
  1669.     arb_init(htheta);
  1670.     arb_min(htheta, theta, htheta_max, prec);
  1671.     mag_zero(arb_radref(htheta));
  1672.  
  1673.     arb_init(estimate);
  1674.     mag_init(m_estimate);
  1675.     H_approximation(estimate, x, y, prec);
  1676.     arb_get_mag(m_estimate, estimate);
  1677.  
  1678.     arb_init(u);
  1679.     D(d, theta, m_estimate, t, y, prec);
  1680.     arb_sqr(u, theta, prec);
  1681.     arb_mul(u, u, t, prec);
  1682.     arb_neg(u, u);
  1683.     arb_exp(u, u, prec);
  1684.     arb_mul(d, d, u, prec);
  1685.  
  1686.     H(h, htheta, m_estimate, t, x, y, prec);
  1687.  
  1688.  
  1689.     arb_clear(htheta);
  1690.     arb_clear(u);
  1691.  
  1692.     arb_clear(estimate);
  1693.     mag_clear(m_estimate);
  1694. }
  1695.  
  1696. int
  1697. evaluate_ymesh(arb_t res,
  1698.         const arb_t t, const arb_t x, const arb_t ya, const arb_t yb,
  1699.         const arb_t htheta_max, slong maxdepth, slong prec)
  1700. {
  1701.     acb_t center;
  1702.     arb_t accum, arg;
  1703.     acb_t h_ratio;
  1704.     ymeshpt_t mb;
  1705.     ymeshival_struct *stack;
  1706.     slong i, depth;
  1707.     ymeshival_ptr ival;
  1708.     int result = 1;
  1709.  
  1710.     stack = flint_malloc(maxdepth * sizeof(*stack));
  1711.     for (i = 0; i < maxdepth; i++)
  1712.     {
  1713.         ymeshival_init(stack + i);
  1714.     }
  1715.  
  1716.     acb_init(center);
  1717.     arb_init(accum);
  1718.     arb_init(arg);
  1719.     acb_init(h_ratio);
  1720.     ymeshpt_init(mb);
  1721.  
  1722.     evaluate_ymeshpt(mb, t, x, yb, htheta_max, prec);
  1723.  
  1724.     depth = 0;
  1725.     ival = stack + depth;
  1726.     evaluate_ymeshpt(ymeshival_leftref(ival), t, x, ya, htheta_max, prec);
  1727.     ival->right = mb;
  1728.     depth++;
  1729.  
  1730.     while (1)
  1731.     {
  1732.         ymeshival_ptr curr = stack + depth - 1;
  1733.         ymeshival_ptr next = stack + depth;
  1734.  
  1735.         if (depth + 5 > maxdepth)
  1736.         {
  1737.             result = 0;
  1738.             arb_indeterminate(res);
  1739.             goto finish;
  1740.         }
  1741.  
  1742.         if (!depth)
  1743.         {
  1744.             break;
  1745.         }
  1746.  
  1747.         if (ymeshival_validate(curr, t, prec))
  1748.         {
  1749.             acb_srcptr h_left = ymeshpt_href(ymeshival_leftref(curr));
  1750.             acb_srcptr h_right = ymeshpt_href(curr->right);
  1751.  
  1752.             acb_div(h_ratio, h_right, h_left, prec);
  1753.  
  1754.             acb_arg(arg, h_ratio, prec);
  1755.  
  1756.             arb_add(accum, accum, arg, prec);
  1757.  
  1758.             if (!arb_is_finite(accum))
  1759.             {
  1760.                 result = 0;
  1761.                 arb_indeterminate(res);
  1762.                 goto finish;
  1763.             }
  1764.  
  1765.             depth--;
  1766.         }
  1767.         else
  1768.         {
  1769.             ymeshival_center(center, curr, prec);
  1770.  
  1771.             evaluate_ymeshpt(ymeshival_leftref(next),
  1772.                     t, acb_realref(center), acb_imagref(center),
  1773.                     htheta_max, prec);
  1774.  
  1775.             ymeshpt_swap(
  1776.                     ymeshival_leftref(curr),
  1777.                     ymeshival_leftref(next));
  1778.  
  1779.             next->right = ymeshival_leftref(curr);
  1780.  
  1781.             depth++;
  1782.         }
  1783.     }
  1784.  
  1785.     arb_set(res, accum);
  1786.  
  1787. finish:
  1788.  
  1789.     acb_clear(center);
  1790.     arb_clear(accum);
  1791.     arb_clear(arg);
  1792.     acb_clear(h_ratio);
  1793.     ymeshpt_clear(mb);
  1794.  
  1795.     for (i = 0; i < maxdepth; i++)
  1796.     {
  1797.         ymeshival_clear(stack + i);
  1798.     }
  1799.     flint_free(stack);
  1800.  
  1801.     return result;
  1802. }
  1803.  
  1804.  
  1805. int
  1806. evaluate_xmesh(arb_t res,
  1807.         const arb_t t, const arb_t xa, const arb_t xb,
  1808.         const arb_t y, const arb_t htheta_max, slong prec)
  1809. {
  1810.     arb_t step, x_next;
  1811.     arb_t accum, arg;
  1812.     acb_t h_ratio;
  1813.     arb_t c1, c2;
  1814.     xmeshpt_t ma, mb;
  1815.     int result = 1;
  1816.     int finished_flag = 0;
  1817.  
  1818.     arb_init(step);
  1819.     arb_init(x_next);
  1820.     arb_init(accum);
  1821.     arb_init(arg);
  1822.     acb_init(h_ratio);
  1823.     arb_init(c1);
  1824.     arb_init(c2);
  1825.     xmeshpt_init(ma);
  1826.     xmeshpt_init(mb);
  1827.  
  1828.     evaluate_xmeshpt(ma, t, xa, y, htheta_max, prec);
  1829.  
  1830.     while (1)
  1831.     {
  1832.         arb_mul(c1, xmeshpt_thetaref(ma), acb_realref(xmeshpt_zref(ma)), prec);
  1833.         arb_exp(c1, c1, prec);
  1834.         acb_abs(c2, xmeshpt_href(ma), prec);
  1835.         arb_div(step, c1, xmeshpt_dref(ma), prec);
  1836.         arb_mul(step, step, c2, prec);
  1837.         arb_mul_2exp_si(step, step, -1);
  1838.         arb_add(x_next, acb_realref(xmeshpt_zref(ma)), step, prec);
  1839.  
  1840.         arb_get_lbound_arf(arb_midref(x_next), x_next, prec);
  1841.         mag_zero(arb_radref(x_next));
  1842.  
  1843.         if (arb_lt(x_next, xb))
  1844.         {
  1845.             evaluate_xmeshpt(mb, t, x_next, y, htheta_max, prec);
  1846.         }
  1847.         else
  1848.         {
  1849.             evaluate_xmeshpt(mb, t, xb, y, htheta_max, prec);
  1850.             finished_flag = 1;
  1851.         }
  1852.  
  1853.         acb_div(h_ratio, xmeshpt_href(mb), xmeshpt_href(ma), prec);
  1854.         acb_arg(arg, h_ratio, prec);
  1855.         arb_add(accum, accum, arg, prec);
  1856.  
  1857.  
  1858.         if (!arb_is_finite(accum))
  1859.         {
  1860.             result = 0;
  1861.             break;
  1862.         }
  1863.  
  1864.         if (finished_flag)
  1865.         {
  1866.             break;
  1867.         }
  1868.         else
  1869.         {
  1870.             xmeshpt_swap(ma, mb);
  1871.         }
  1872.     }
  1873.  
  1874.     arb_set(res, accum);
  1875.  
  1876.     arb_clear(step);
  1877.     arb_clear(x_next);
  1878.     arb_clear(accum);
  1879.     arb_clear(arg);
  1880.     acb_clear(h_ratio);
  1881.     arb_clear(c1);
  1882.     arb_clear(c2);
  1883.     xmeshpt_clear(ma);
  1884.     xmeshpt_clear(mb);
  1885.  
  1886.     return result;
  1887. }
  1888.  
  1889.  
  1890.  
  1891. int main(int argc, char *argv[])
  1892. {
  1893.     arb_t t, xa, xb, ya, yb;
  1894.     arb_t arg, arg_total, htheta_max;
  1895.     const char *t_str, *xa_str, *xb_str, *ya_str, *yb_str;
  1896.     slong prec;
  1897.     slong prec_base = 64;
  1898.     slong prec_incr = 32;
  1899.     slong maxdepth = 50;
  1900.     int result = EXIT_SUCCESS;
  1901.  
  1902.     arb_init(t);
  1903.     arb_init(xa);
  1904.     arb_init(xb);
  1905.     arb_init(ya);
  1906.     arb_init(yb);
  1907.  
  1908.     arb_init(arg);
  1909.     arb_init(arg_total);
  1910.     arb_init(htheta_max);
  1911.  
  1912.     if (argc != 6)
  1913.     {
  1914.         result = EXIT_FAILURE;
  1915.         goto finish;
  1916.     }
  1917.  
  1918.     t_str = argv[1];
  1919.     xa_str = argv[2];
  1920.     xb_str = argv[3];
  1921.     ya_str = argv[4];
  1922.     yb_str = argv[5];
  1923.  
  1924.     prec_base = 64;
  1925.     prec_incr = 32;
  1926.  
  1927.     arb_set_str(htheta_max, "0.369", 32);
  1928.  
  1929.     prec = prec_base / 2;
  1930.     while (1)
  1931.     {
  1932.         prec *= 2;
  1933.  
  1934.         arb_set_str(t, t_str, prec);
  1935.         arb_set_str(xa, xa_str, prec);
  1936.         arb_set_str(xb, xb_str, prec);
  1937.         arb_set_str(ya, ya_str, prec);
  1938.         arb_set_str(yb, yb_str, prec);
  1939.  
  1940.         if (arb_lt(xb, xa) || arb_lt(yb, ya))
  1941.         {
  1942.             result = EXIT_FAILURE;
  1943.             goto finish;
  1944.         }
  1945.  
  1946.         if (evaluate_ymesh(arg, t, xa, ya, yb, htheta_max, maxdepth, prec))
  1947.         {
  1948.             break;
  1949.         }
  1950.     }
  1951.  
  1952.     arb_neg(arg, arg);
  1953.     arb_add(arg_total, arg_total, arg, prec);
  1954.  
  1955.     prec = prec_base / 2;
  1956.     while (1)
  1957.     {
  1958.         prec *= 2;
  1959.  
  1960.         arb_set_str(t, t_str, prec);
  1961.         arb_set_str(xa, xa_str, prec);
  1962.         arb_set_str(xb, xb_str, prec);
  1963.         arb_set_str(ya, ya_str, prec);
  1964.         arb_set_str(yb, yb_str, prec);
  1965.  
  1966.         if (arb_lt(xb, xa) || arb_lt(yb, ya))
  1967.         {
  1968.             result = EXIT_FAILURE;
  1969.             goto finish;
  1970.         }
  1971.  
  1972.         if (evaluate_ymesh(arg, t, xb, ya, yb, htheta_max, maxdepth, prec))
  1973.         {
  1974.             break;
  1975.         }
  1976.     }
  1977.  
  1978.     arb_add(arg_total, arg_total, arg, prec);
  1979.  
  1980.     prec = prec_base / 2;
  1981.     while (1)
  1982.     {
  1983.         prec *= 2;
  1984.  
  1985.         arb_set_str(t, t_str, prec);
  1986.         arb_set_str(xa, xa_str, prec);
  1987.         arb_set_str(xb, xb_str, prec);
  1988.         arb_set_str(ya, ya_str, prec);
  1989.         arb_set_str(yb, yb_str, prec);
  1990.  
  1991.         if (arb_lt(xb, xa) || arb_lt(yb, ya))
  1992.         {
  1993.             result = EXIT_FAILURE;
  1994.             goto finish;
  1995.         }
  1996.  
  1997.         if (evaluate_xmesh(arg, t, xa, xb, ya, htheta_max, prec))
  1998.         {
  1999.             break;
  2000.         }
  2001.     }
  2002.  
  2003.     arb_add(arg_total, arg_total, arg, prec);
  2004.  
  2005.     prec = prec_base / 2;
  2006.     while (1)
  2007.     {
  2008.         prec *= 2;
  2009.  
  2010.         arb_set_str(t, t_str, prec);
  2011.         arb_set_str(xa, xa_str, prec);
  2012.         arb_set_str(xb, xb_str, prec);
  2013.         arb_set_str(ya, ya_str, prec);
  2014.         arb_set_str(yb, yb_str, prec);
  2015.  
  2016.         if (arb_lt(xb, xa) || arb_lt(yb, ya))
  2017.         {
  2018.             result = EXIT_FAILURE;
  2019.             goto finish;
  2020.         }
  2021.  
  2022.         if (evaluate_xmesh(arg, t, xa, xb, yb, htheta_max, prec))
  2023.         {
  2024.             break;
  2025.         }
  2026.     }
  2027.  
  2028.     arb_neg(arg, arg);
  2029.     arb_add(arg_total, arg_total, arg, prec);
  2030.  
  2031.     {
  2032.         arb_t pi;
  2033.         arb_init(pi);
  2034.         arb_const_pi(pi, prec);
  2035.         arb_div(arg_total, arg_total, pi, prec);
  2036.         arb_mul_2exp_si(arg_total, arg_total, -1);
  2037.         arb_clear(pi);
  2038.     }
  2039.  
  2040.     {
  2041.         double d = arf_get_d(arb_midref(arg_total), ARF_RND_NEAR);
  2042.         slong nroots = (int) (d + 0.5);
  2043.         if (arb_contains_si(arg_total, nroots) &&
  2044.             !arb_contains_si(arg_total, nroots - 1) &&
  2045.             !arb_contains_si(arg_total, nroots + 1))
  2046.         {
  2047.             flint_printf("%ld\n", nroots);
  2048.         }
  2049.         else
  2050.         {
  2051.             arb_printd(arg_total, 10);
  2052.             flint_printf("\n");
  2053.         }
  2054.     }
  2055.  
  2056. finish:
  2057.  
  2058.     if (result == EXIT_FAILURE)
  2059.     {
  2060.         flint_printf("%s t xa xb ya yb\n\n", argv[0]);
  2061.         flint_printf(
  2062.                 "Count the roots of H_t(x + yi).\n"
  2063.                 "xa < x < xb\n"
  2064.                 "ya < y < yb\n");
  2065.     }
  2066.  
  2067.     arb_clear(t);
  2068.     arb_clear(xa);
  2069.     arb_clear(xb);
  2070.     arb_clear(ya);
  2071.     arb_clear(yb);
  2072.  
  2073.     arb_clear(arg);
  2074.     arb_clear(arg_total);
  2075.     arb_clear(htheta_max);
  2076.  
  2077.     flint_cleanup();
  2078.     return result;
  2079. }
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!