Ledger Nano X - The secure hardware wallet
SHARE
TWEET

Untitled

a guest Mar 27th, 2018 67 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  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. }
RAW Paste Data
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
Top