API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Mar 25th, 2018
105
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 20.52 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. static void
  13. _acb_get_rad_ubound_mag(mag_t res, const acb_t z, slong prec)
  14. {
  15.     arf_t err;
  16.     arf_init(err);
  17.     acb_get_rad_ubound_arf(err, z, prec);
  18.     arf_get_mag(res, err);
  19.     arf_clear(err);
  20. }
  21.  
  22. typedef struct
  23. {
  24.     acb_t z;
  25.     arb_t t;
  26.     arb_t T;
  27. } thirdbound_integrand_param_struct;
  28. typedef thirdbound_integrand_param_struct thirdbound_integrand_param_t[1];
  29. typedef thirdbound_integrand_param_struct *thirdbound_integrand_param_ptr;
  30. typedef const thirdbound_integrand_param_struct *thirdbound_integrand_param_srcptr;
  31.  
  32. static void
  33. thirdbound_integrand_param_init(thirdbound_integrand_param_t p,
  34.         const acb_t z, const arb_t t, const arb_t T)
  35. {
  36.     acb_init(p->z);
  37.     arb_init(p->t);
  38.     arb_init(p->T);
  39.     acb_set(p->z, z);
  40.     arb_set(p->t, t);
  41.     arb_set(p->T, T);
  42. }
  43.  
  44. static void
  45. thirdbound_integrand_param_clear(thirdbound_integrand_param_t p)
  46. {
  47.     acb_clear(p->z);
  48.     arb_clear(p->t);
  49.     arb_clear(p->T);
  50. }
  51.  
  52. void thirdbound_Q(arb_t res, const arb_t t, const arb_t y,
  53.         const arb_t X, const arb_t T, slong prec)
  54. {
  55.     arb_t z1, z2, z3;
  56.  
  57.     arb_init(z1);
  58.     arb_mul_2exp_si(z1, X, 1);
  59.     arb_add(z1, z1, y, prec);
  60.     arb_sub_ui(z1, z1, 1, prec);
  61.     arb_div(z1, z1, t, prec);
  62.  
  63.     arb_init(z2);
  64.     arb_set_ui(z2, 3);
  65.     arb_div(z2, z2, T, prec);
  66.     arb_mul_2exp_si(z2, z2, -2);
  67.  
  68.     arb_init(z3);
  69.     arb_log_hypot(z3, X, T, prec);
  70.     arb_mul_2exp_si(z3, z3, -1);
  71.  
  72.     arb_const_log_sqrt2pi(res, prec);
  73.     arb_add(res, res, z1, prec);
  74.     arb_sub(res, res, z2, prec);
  75.     arb_sub(res, res, z3, prec);
  76.  
  77.     arb_clear(z1);
  78.     arb_clear(z2);
  79.     arb_clear(z3);
  80. }
  81.  
  82. void thirdbound_F(acb_t res, const acb_t s, slong prec)
  83. {
  84.     acb_t sd2, sm1;
  85.     acb_t z1, z2, z3;
  86.  
  87.     acb_init(sd2);
  88.     acb_mul_2exp_si(sd2, s, -1);
  89.  
  90.     acb_init(sm1);
  91.     acb_sub_ui(sm1, s, 1, prec);
  92.  
  93.     acb_init(z1);
  94.     acb_log(z1, sm1, prec);
  95.  
  96.     acb_init(z2);
  97.     acb_const_pi(z2, prec);
  98.     acb_log(z2, z2, prec);
  99.     acb_mul(z2, z2, sd2, prec);
  100.  
  101.     acb_init(z3);
  102.     acb_log(z3, sd2, prec);
  103.     acb_mul(z3, z3, sm1, prec);
  104.     acb_mul_2exp_si(z3, z3, -1);
  105.  
  106.     acb_log(res, s, prec);
  107.     acb_add(res, res, z1, prec);
  108.     acb_sub(res, res, z2, prec);
  109.     acb_add(res, res, z3, prec);
  110.     acb_sub(res, res, sd2, prec);
  111.  
  112.     acb_clear(sd2);
  113.     acb_clear(sm1);
  114.  
  115.     acb_clear(z1);
  116.     acb_clear(z2);
  117.     acb_clear(z3);
  118. }
  119.  
  120. void riemann_xi(acb_t res, const acb_t s, slong prec)
  121. {
  122.     acb_t pi, z1, z2, z3, z4;
  123.  
  124.     acb_init(pi);
  125.     acb_const_pi(pi, prec);
  126.  
  127.     acb_init(z1);
  128.     acb_sub_ui(z1, s, 1, prec);
  129.     acb_mul(z1, z1, s, prec);
  130.     acb_mul_2exp_si(z1, z1, -1);
  131.  
  132.     acb_init(z2);
  133.     acb_mul_2exp_si(z2, s, -1);
  134.     acb_neg(z2, z2);
  135.     acb_pow(z2, pi, z2, prec);
  136.  
  137.     acb_init(z3);
  138.     acb_mul_2exp_si(z3, s, -1);
  139.     acb_gamma(z3, z3, prec);
  140.  
  141.     acb_init(z4);
  142.     acb_zeta(z4, s, prec);
  143.  
  144.     acb_mul(res, z1, z2, prec);
  145.     acb_mul(res, res, z3, prec);
  146.     acb_mul(res, res, z4, prec);
  147.  
  148.     acb_clear(pi);
  149.     acb_clear(z1);
  150.     acb_clear(z2);
  151.     acb_clear(z3);
  152.     acb_clear(z4);
  153. }
  154.  
  155.  
  156.  
  157.  
  158. static void
  159. _thirdbound_remainder_helper(arb_t res,
  160.         const arb_t t, const acb_t z,
  161.         const arb_t X, const arb_t T, slong prec)
  162. {
  163.     arb_t z1, z2, z3, z4;
  164.     acb_t u;
  165.     arb_srcptr x, y;
  166.  
  167.     x = acb_realref(z);
  168.     y = acb_imagref(z);
  169.  
  170.     arb_init(z1);
  171.     arb_mul_2exp_si(z1, x, -1);
  172.     arb_sub(z1, T, z1, prec);
  173.     arb_sqr(z1, z1, prec);
  174.     arb_div(z1, z1, t, prec);
  175.     arb_mul_2exp_si(z1, z1, -2);
  176.  
  177.     arb_init(z2);
  178.     arb_set_str(z2, "0.33", prec);
  179.     arb_sub(z2, T, z2, prec);
  180.     arb_mul_ui(z2, z2, 6, prec);
  181.     arb_inv(z2, z2, prec);
  182.  
  183.     acb_init(u);
  184.     acb_set_arb_arb(u, X, T);
  185.     thirdbound_F(u, u, prec);
  186.  
  187.     arb_init(z3);
  188.     arb_set(z3, acb_realref(u));
  189.  
  190.     arb_init(z4);
  191.     arb_sub_ui(z4, y, 1, prec);
  192.     arb_mul_2exp_si(z4, z4, -1);
  193.     arb_add(z4, X, z4, prec);
  194.     arb_sqr(z4, z4, prec);
  195.     arb_div(z4, z4, t, prec);
  196.  
  197.     arb_add(res, z1, z2, prec);
  198.     arb_add(res, res, z3, prec);
  199.     arb_sub(res, res, z4, prec);
  200.     arb_exp(res, res, prec);
  201.  
  202.     arb_clear(z1);
  203.     arb_clear(z2);
  204.     arb_clear(z3);
  205.     arb_clear(z4);
  206.  
  207.     acb_clear(u);
  208. }
  209.  
  210.  
  211. void thirdbound_H_remainder(arb_t res,
  212.         const arb_t t, const acb_t z,
  213.         const arb_t X, const arb_t T, slong prec)
  214. {
  215.     arb_t q, u;
  216.     arb_srcptr x, y;
  217.  
  218.     x = acb_realref(z);
  219.     y = acb_imagref(z);
  220.  
  221.     arb_init(q);
  222.     thirdbound_Q(q, t, y, X, T, prec);
  223.  
  224.     arb_init(u);
  225.     arb_mul_2exp_si(u, t, 5);
  226.     arb_sqrt(u, u, prec);
  227.     arb_mul(u, u, q, prec);
  228.  
  229.     _thirdbound_remainder_helper(res, t, z, X, T, prec);
  230.  
  231.     arb_div(res, res, u, prec);
  232.  
  233.     arb_clear(q);
  234.     arb_clear(u);
  235. }
  236.  
  237.  
  238. void thirdbound_Hprime_remainder(arb_t res,
  239.         const arb_t t, const acb_t z,
  240.         const arb_t X, const arb_t T, slong prec)
  241. {
  242.     arb_t qinv, u;
  243.     arb_t z1, z2, w;
  244.     arb_srcptr x, y;
  245.  
  246.     x = acb_realref(z);
  247.     y = acb_imagref(z);
  248.  
  249.     arb_init(qinv);
  250.     thirdbound_Q(qinv, t, y, X, T, prec);
  251.     arb_inv(qinv, qinv, prec);
  252.  
  253.     arb_init(u);
  254.     arb_pow_ui(u, t, 3, prec);
  255.     arb_mul_2exp_si(u, u, 5);
  256.     arb_sqrt(u, u, prec);
  257.  
  258.     arb_init(z1);
  259.     arb_mul_2exp_si(z1, x, -1);
  260.     arb_sub(z1, T, z1, prec);
  261.     arb_abs(z1, z1);
  262.  
  263.     arb_init(z2);
  264.     arb_sub_ui(z2, y, 1, prec);
  265.     arb_mul_2exp_si(z2, z2, -1);
  266.     arb_add(z2, z2, X, prec);
  267.    
  268.     arb_init(w);
  269.     arb_add(w, z1, z2, prec);
  270.     arb_add(w, w, qinv, prec);
  271.     arb_mul(w, w, qinv, prec);
  272.  
  273.     _thirdbound_remainder_helper(res, t, z, X, T, prec);
  274.     arb_mul(res, res, w, prec);
  275.     arb_div(res, res, u, prec);
  276.  
  277.     arb_clear(qinv);
  278.     arb_clear(u);
  279.     arb_clear(z1);
  280.     arb_clear(z2);
  281.     arb_clear(w);
  282. }
  283.  
  284.  
  285.  
  286.  
  287. static void
  288. _thirdbound_integrand_helper(acb_t res1, acb_t res2,
  289.         const acb_t s, const acb_t z, const arb_t T, slong prec)
  290. {
  291.     acb_t onei, iT;
  292.     acb_t z1, z2;
  293.     arb_srcptr x, y;
  294.  
  295.     x = acb_realref(z);
  296.     y = acb_imagref(z);
  297.  
  298.     acb_init(onei);
  299.     acb_onei(onei);
  300.  
  301.     acb_init(iT);
  302.     acb_mul_arb(iT, onei, T, prec);
  303.  
  304.     acb_init(z1);
  305.     acb_add(z1, s, iT, prec);
  306.     riemann_xi(z1, z1, prec);
  307.  
  308.     acb_init(z2);
  309.     acb_mul_arb(z2, onei, x, prec);
  310.     acb_sub_arb(z2, z2, y, prec);
  311.     acb_add_ui(z2, z2, 1, prec);
  312.     acb_mul_2exp_si(z2, z2, -1);
  313.     acb_sub(z2, iT, z2, prec);
  314.  
  315.     acb_swap(res1, z1);
  316.     acb_swap(res2, z2);
  317.  
  318.     acb_clear(z1);
  319.     acb_clear(z2);
  320.  
  321.     acb_clear(onei);
  322.     acb_clear(iT);
  323. }
  324.  
  325.  
  326.  
  327. void thirdbound_H_integrand(acb_t res, const acb_t z,
  328.         const arb_t t, const acb_t s, const arb_t T, slong prec)
  329. {
  330.     acb_t h1, h2, z1, z2, w1, w2;
  331.  
  332.     acb_init(h1);
  333.     acb_init(h2);
  334.     _thirdbound_integrand_helper(h1, h2, s, z, T, prec);
  335.  
  336.     acb_init(z1);
  337.     acb_add(z1, s, h2, prec);
  338.     acb_sqr(z1, z1, prec);
  339.     acb_neg(z1, z1);
  340.     acb_div_arb(z1, z1, t, prec);
  341.     acb_exp(z1, z1, prec);
  342.  
  343.     acb_init(z2);
  344.     acb_sub(z2, h2, s, prec);
  345.     acb_add_ui(z2, z2, 1, prec);
  346.     acb_sqr(z2, z2, prec);
  347.     acb_neg(z2, z2);
  348.     acb_div_arb(z2, z2, t, prec);
  349.     acb_exp(z2, z2, prec);
  350.  
  351.     acb_init(w1);
  352.     acb_mul(w1, h1, z1, prec);
  353.  
  354.     acb_init(w2);
  355.     acb_conj(w2, h1);
  356.     acb_mul(w2, w2, z2, prec);
  357.  
  358.     acb_add(res, w1, w2, prec);
  359.  
  360.     acb_clear(h1);
  361.     acb_clear(h2);
  362.     acb_clear(z1);
  363.     acb_clear(z2);
  364.     acb_clear(w1);
  365.     acb_clear(w2);
  366. }
  367.  
  368. void thirdbound_Hprime_integrand(acb_t res, const acb_t z,
  369.         const arb_t t, const acb_t s, const arb_t T, slong prec)
  370. {
  371.     acb_t h1, h2, u1, u2, z1, z2, w1, w2;
  372.  
  373.     acb_init(h1);
  374.     acb_init(h2);
  375.     _thirdbound_integrand_helper(h1, h2, s, z, T, prec);
  376.  
  377.     acb_init(u1);
  378.     acb_add(u1, s, h2, prec);
  379.  
  380.     acb_init(u2);
  381.     acb_add_ui(u2, h2, 1, prec);
  382.     acb_sub(u2, u2, s, prec);
  383.  
  384.     acb_init(z1);
  385.     acb_sqr(z1, u1, prec);
  386.     acb_neg(z1, z1);
  387.     acb_div_arb(z1, z1, t, prec);
  388.     acb_exp(z1, z1, prec);
  389.  
  390.     acb_init(z2);
  391.     acb_sqr(z2, u2, prec);
  392.     acb_neg(z2, z2);
  393.     acb_div_arb(z2, z2, t, prec);
  394.     acb_exp(z2, z2, prec);
  395.  
  396.     acb_init(w1);
  397.     acb_mul(w1, h1, u1, prec);
  398.     acb_mul(w1, w1, z1, prec);
  399.  
  400.     acb_init(w2);
  401.     acb_conj(w2, h1);
  402.     acb_mul(w2, w2, u2, prec);
  403.     acb_mul(w2, w2, z2, prec);
  404.  
  405.     acb_add(res, w1, w2, prec);
  406.  
  407.     acb_clear(h1);
  408.     acb_clear(h2);
  409.     acb_clear(u1);
  410.     acb_clear(u2);
  411.     acb_clear(z1);
  412.     acb_clear(z2);
  413.     acb_clear(w1);
  414.     acb_clear(w2);
  415. }
  416.  
  417. static int
  418. f_thirdbound_H_integrand(acb_ptr res, const acb_t u, void * param,
  419.         slong order, slong prec)
  420. {
  421.     thirdbound_integrand_param_srcptr p;
  422.  
  423.     if (order > 1)
  424.         flint_abort();
  425.  
  426.     p = param;
  427.     thirdbound_H_integrand(res, p->z, p->t, u, p->T, prec);
  428.  
  429.     return 0;
  430. }
  431.  
  432. static int
  433. f_thirdbound_Hprime_integrand(acb_ptr res, const acb_t u, void * param,
  434.         slong order, slong prec)
  435. {
  436.     thirdbound_integrand_param_srcptr p;
  437.  
  438.     if (order > 1)
  439.         flint_abort();
  440.  
  441.     p = param;
  442.     thirdbound_Hprime_integrand(res, p->z, p->t, u, p->T, prec);
  443.  
  444.     return 0;
  445. }
  446.  
  447. void
  448. thirdbound_H_proper_integral(acb_t res,
  449.         const mag_t abs_estimate,
  450.         const arb_t lower_limit, const arb_t upper_limit,
  451.         const acb_t z, const arb_t t, const arb_t T, slong prec)
  452. {
  453.     thirdbound_integrand_param_t p;
  454.     acb_t a_lim, b_lim;
  455.     acb_calc_integrate_opt_t options;
  456.     int calc_result;
  457.     mag_t abs_tol;
  458.     slong rel_goal = prec;
  459.  
  460.     thirdbound_integrand_param_init(p, z, t, T);
  461.  
  462.     acb_init(a_lim);
  463.     acb_init(b_lim);
  464.     acb_set_arb(a_lim, lower_limit);
  465.     acb_set_arb(b_lim, upper_limit);
  466.  
  467.     mag_init(abs_tol);
  468.     mag_mul_2exp_si(abs_tol, abs_estimate, -prec);
  469.  
  470.     acb_calc_integrate_opt_init(options);
  471.     options->verbose = 0;
  472.     calc_result = acb_calc_integrate(
  473.             res, f_thirdbound_H_integrand, p, a_lim, b_lim,
  474.             rel_goal, abs_tol, options, prec);
  475.     if (calc_result == ARB_CALC_NO_CONVERGENCE)
  476.     {
  477.         acb_indeterminate(res);
  478.     }
  479.     else
  480.     {
  481.         arb_t u;
  482.         arb_init(u);
  483.         arb_const_pi(u, prec);
  484.         arb_mul(u, u, t, prec);
  485.         arb_sqrt(u, u, prec);
  486.         acb_div_arb(res, res, u, prec);
  487.         acb_mul_2exp_si(res, res, -3);
  488.         arb_clear(u);
  489.     }
  490.  
  491.     thirdbound_integrand_param_clear(p);
  492.     acb_clear(a_lim);
  493.     acb_clear(b_lim);
  494.     mag_clear(abs_tol);
  495. }
  496.  
  497.  
  498. void
  499. thirdbound_Hprime_proper_integral(acb_t res,
  500.         const mag_t abs_estimate,
  501.         const arb_t lower_limit, const arb_t upper_limit,
  502.         const acb_t z, const arb_t t, const arb_t T, slong prec)
  503. {
  504.     thirdbound_integrand_param_t p;
  505.     acb_t a_lim, b_lim;
  506.     acb_calc_integrate_opt_t options;
  507.     int calc_result;
  508.     mag_t abs_tol;
  509.     slong rel_goal = prec;
  510.  
  511.     thirdbound_integrand_param_init(p, z, t, T);
  512.  
  513.     acb_init(a_lim);
  514.     acb_init(b_lim);
  515.     acb_set_arb(a_lim, lower_limit);
  516.     acb_set_arb(b_lim, upper_limit);
  517.  
  518.     mag_init(abs_tol);
  519.     mag_mul_2exp_si(abs_tol, abs_estimate, -prec);
  520.  
  521.     acb_calc_integrate_opt_init(options);
  522.     options->verbose = 0;
  523.     calc_result = acb_calc_integrate(
  524.             res, f_thirdbound_Hprime_integrand, p, a_lim, b_lim,
  525.             rel_goal, abs_tol, options, prec);
  526.     if (calc_result == ARB_CALC_NO_CONVERGENCE)
  527.     {
  528.         acb_indeterminate(res);
  529.     }
  530.     else
  531.     {
  532.         fmpq_t q;
  533.         arb_t sqrtpi, tpow;
  534.         acb_t onei;
  535.  
  536.         fmpq_init(q);
  537.         fmpq_set_si(q, 3, 2);
  538.  
  539.         arb_init(sqrtpi);
  540.         arb_const_sqrt_pi(sqrtpi, prec);
  541.  
  542.         arb_init(tpow);
  543.         arb_pow_fmpq(tpow, t, q, prec);
  544.  
  545.         acb_init(onei);
  546.         acb_onei(onei);
  547.  
  548.         acb_mul(res, res, onei, prec);
  549.         acb_div_arb(res, res, sqrtpi, prec);
  550.         acb_div_arb(res, res, tpow, prec);
  551.         acb_mul_2exp_si(res, res, -3);
  552.  
  553.         fmpq_clear(q);
  554.         arb_clear(sqrtpi);
  555.         arb_clear(tpow);
  556.         acb_clear(onei);
  557.     }
  558.  
  559.     thirdbound_integrand_param_clear(p);
  560.     acb_clear(a_lim);
  561.     acb_clear(b_lim);
  562.     mag_clear(abs_tol);
  563. }
  564.  
  565.  
  566. int thirdbound_H_remainder_is_permissible(const arb_t t, const arb_t y,
  567.         const arb_t X, const arb_t T, slong prec)
  568. {
  569.     arb_t q, zero, one, two, half;
  570.     int result;
  571.  
  572.     arb_init(q);
  573.     arb_init(zero);
  574.     arb_init(one);
  575.     arb_init(two);
  576.     arb_init(half);
  577.  
  578.     arb_zero(zero);
  579.     arb_one(one);
  580.     arb_set_ui(two, 2);
  581.     arb_mul_2exp_si(half, one, -1);
  582.     thirdbound_Q(q, t, y, X, T, prec);
  583.  
  584.     result = (
  585.             arb_gt(X, one) &&
  586.             arb_ge(y, zero) &&
  587.             arb_ge(T, two) &&
  588.             arb_le(t, half) &&
  589.             arb_gt(q, zero));
  590.  
  591.     arb_clear(q);
  592.     arb_clear(zero);
  593.     arb_clear(one);
  594.     arb_clear(two);
  595.     arb_clear(half);
  596.  
  597.     return result;
  598. }
  599.  
  600. void
  601. thirdbound_H_integral(acb_t res,
  602.         const mag_t m_estimate,
  603.         const acb_t z, const arb_t t, const arb_t T, slong prec)
  604. {
  605.     arb_t lower, upper, remainder;
  606.     acb_t value, partial;
  607.     mag_t m_partial_sum, tail_mag_lb, tail_mag_ub;
  608.     slong r, p, q, tmp;
  609.  
  610.     arb_init(lower);
  611.     arb_init(upper);
  612.     arb_init(remainder);
  613.  
  614.     acb_init(value);
  615.     acb_init(partial);
  616.  
  617.     arb_init(remainder);
  618.  
  619.     mag_init(m_partial_sum);
  620.     mag_init(tail_mag_lb);
  621.     mag_init(tail_mag_ub);
  622.  
  623.     r = 0;
  624.     p = 0;
  625.     q = 1;
  626.     while (1)
  627.     {
  628.         tmp = p + q;
  629.         arb_set_ui(lower, 1 + r);
  630.         arb_set_ui(upper, 1 + r + q);
  631.         arb_mul_2exp_si(lower, lower, -1);
  632.         arb_mul_2exp_si(upper, upper, -1);
  633.        
  634.         r = r + q;
  635.         p = q;
  636.         q = tmp;
  637.  
  638.         thirdbound_H_proper_integral(value,
  639.                 m_estimate, lower, upper, z, t, T, prec);
  640.         acb_add(partial, partial, value, prec);
  641.  
  642.         if (thirdbound_H_remainder_is_permissible(
  643.                     t, acb_imagref(z), upper, T, prec))
  644.         {
  645.             thirdbound_H_remainder(remainder,
  646.                     t, z, upper, T, prec);
  647.             arb_get_mag(tail_mag_ub, remainder);
  648.  
  649.             _acb_get_rad_ubound_mag(m_partial_sum, partial, prec);
  650.             arb_get_mag_lower(tail_mag_lb, remainder);
  651.  
  652.             if (mag_cmp(tail_mag_lb, m_partial_sum) <= 0)
  653.             {
  654.                 acb_set(res, partial);
  655.                 acb_add_error_mag(res, tail_mag_ub);
  656.                 break;
  657.             }
  658.         }
  659.     }
  660.  
  661.     arb_clear(lower);
  662.     arb_clear(upper);
  663.  
  664.     acb_clear(value);
  665.     acb_clear(partial);
  666.  
  667.     arb_clear(remainder);
  668.  
  669.     mag_clear(m_partial_sum);
  670.     mag_clear(tail_mag_lb);
  671.     mag_clear(tail_mag_ub);
  672. }
  673.  
  674.  
  675. void
  676. thirdbound_Hprime_integral(acb_t res,
  677.         const mag_t m_estimate,
  678.         const acb_t z, const arb_t t, const arb_t T, slong prec)
  679. {
  680.     arb_t lower, upper, remainder;
  681.     acb_t value, partial;
  682.     mag_t m_partial_sum, tail_mag_lb, tail_mag_ub;
  683.     slong r, p, q, tmp;
  684.  
  685.     arb_init(lower);
  686.     arb_init(upper);
  687.     arb_init(remainder);
  688.  
  689.     acb_init(value);
  690.     acb_init(partial);
  691.  
  692.     arb_init(remainder);
  693.  
  694.     mag_init(m_partial_sum);
  695.     mag_init(tail_mag_lb);
  696.     mag_init(tail_mag_ub);
  697.  
  698.     r = 0;
  699.     p = 0;
  700.     q = 1;
  701.     while (1)
  702.     {
  703.         tmp = p + q;
  704.         arb_set_ui(lower, 1 + r);
  705.         arb_set_ui(upper, 1 + r + q);
  706.         arb_mul_2exp_si(lower, lower, -1);
  707.         arb_mul_2exp_si(upper, upper, -1);
  708.        
  709.         r = r + q;
  710.         p = q;
  711.         q = tmp;
  712.  
  713.         thirdbound_Hprime_proper_integral(value,
  714.                 m_estimate, lower, upper, z, t, T, prec);
  715.         acb_add(partial, partial, value, prec);
  716.  
  717.         if (thirdbound_H_remainder_is_permissible(
  718.                     t, acb_imagref(z), upper, T, prec))
  719.         {
  720.             thirdbound_Hprime_remainder(remainder,
  721.                     t, z, upper, T, prec);
  722.             arb_get_mag(tail_mag_ub, remainder);
  723.  
  724.             _acb_get_rad_ubound_mag(m_partial_sum, partial, prec);
  725.             arb_get_mag_lower(tail_mag_lb, remainder);
  726.  
  727.             if (mag_cmp(tail_mag_lb, m_partial_sum) <= 0)
  728.             {
  729.                 acb_set(res, partial);
  730.                 acb_add_error_mag(res, tail_mag_ub);
  731.                 break;
  732.             }
  733.         }
  734.     }
  735.  
  736.     arb_clear(lower);
  737.     arb_clear(upper);
  738.  
  739.     acb_clear(value);
  740.     acb_clear(partial);
  741.  
  742.     arb_clear(remainder);
  743.  
  744.     mag_clear(m_partial_sum);
  745.     mag_clear(tail_mag_lb);
  746.     mag_clear(tail_mag_ub);
  747. }
  748.  
  749. static void
  750. H_approximation(arb_t res, const arb_t x, const arb_t y, slong prec)
  751. {
  752.     if (arb_contains_nonpositive(x))
  753.     {
  754.         arb_one(res);
  755.     }
  756.     else
  757.     {
  758.         arb_t logx, pixd8;
  759.  
  760.         arb_init(logx);
  761.         arb_log(logx, x, prec);
  762.  
  763.         arb_init(pixd8);
  764.         arb_const_pi(pixd8, prec);
  765.         arb_mul(pixd8, pixd8, x, prec);
  766.         arb_mul_2exp_si(pixd8, pixd8, -3);
  767.  
  768.         arb_add_si(res, y, 7, prec);
  769.         arb_mul(res, res, logx, prec);
  770.         arb_mul_2exp_si(res, res, -2);
  771.         arb_sub(res, res, pixd8, prec);
  772.         arb_exp(res, res, prec);
  773.  
  774.         arb_clear(logx);
  775.         arb_clear(pixd8);
  776.     }
  777. }
  778.  
  779. int main(int argc, char *argv[])
  780. {
  781.     int result;
  782.     arb_t t, T;
  783.     acb_t z, h;
  784.     arb_t half, four;
  785.     slong digits, prec, high_prec;
  786.     arb_t estimate;
  787.     mag_t m_estimate;
  788.     slong n;
  789.     arb_ptr x, y;
  790.  
  791.     arb_init(t);
  792.     arb_init(T);
  793.  
  794.     acb_init(z);
  795.     acb_init(h);
  796.  
  797.     arb_init(half);
  798.     arb_init(four);
  799.  
  800.     arb_init(estimate);
  801.     mag_init(m_estimate);
  802.  
  803.     x = acb_realref(z);
  804.     y = acb_imagref(z);
  805.  
  806.     arb_one(half);
  807.     arb_mul_2exp_si(half, half, -1);
  808.     arb_set_ui(four, 4);
  809.    
  810.     result = EXIT_SUCCESS;
  811.  
  812.     if (argc != 6)
  813.     {
  814.         result = EXIT_FAILURE;
  815.         goto finish;
  816.     }
  817.  
  818.     n = atol(argv[4]);
  819.     if (n != 0 && n != 1)
  820.     {
  821.         result = EXIT_FAILURE;
  822.         goto finish;
  823.     }
  824.  
  825.     digits = atol(argv[5]);
  826.     high_prec = digits * 3.32192809488736 + 10;
  827.  
  828.     prec = 16;
  829.     acb_indeterminate(h);
  830.     while (arb_rel_accuracy_bits(acb_realref(h)) < high_prec ||
  831.            arb_rel_accuracy_bits(acb_imagref(h)) < high_prec)
  832.     {
  833.         prec *= 2;
  834.         arb_set_str(t, argv[1], prec);
  835.         arb_set_str(acb_realref(z), argv[2], prec);
  836.         arb_set_str(acb_imagref(z), argv[3], prec);
  837.  
  838.         if (!arb_is_nonnegative(x) ||
  839.             !arb_is_nonnegative(y) ||
  840.             !arb_le(t, half))
  841.         {
  842.             result = EXIT_FAILURE;
  843.             goto finish;
  844.         }
  845.  
  846.         arb_const_pi(T, prec);
  847.         arb_mul_2exp_si(T, T, -1);
  848.         arb_add(T, T, x, prec);
  849.         arb_mul_2exp_si(T, T, -1);
  850.         arb_max(T, T, four, prec);
  851.  
  852.         H_approximation(estimate, x, y, prec);
  853.         arb_get_mag(m_estimate, estimate);
  854.  
  855.         if (n == 0)
  856.         {
  857.             thirdbound_H_integral(h, m_estimate, z, t, T, prec);
  858.             if (arb_is_zero(x) || arb_is_zero(y))
  859.             {
  860.                 arb_zero(acb_imagref(h));
  861.             }
  862.         }
  863.         else if (n == 1)
  864.         {
  865.             thirdbound_Hprime_integral(h, m_estimate, z, t, T, prec);
  866.             if (arb_is_zero(x))
  867.             {
  868.                 arb_zero(acb_realref(h));
  869.             }
  870.             if (arb_is_zero(y))
  871.             {
  872.                 arb_zero(acb_imagref(h));
  873.             }
  874.         }
  875.         else
  876.         {
  877.             flint_abort();
  878.         }
  879.  
  880.     }
  881.  
  882.     flint_printf("Re: ");
  883.     arf_printd(arb_midref(acb_realref(h)), digits);
  884.     flint_printf("\n");
  885.     flint_printf("Im: ");
  886.     arf_printd(arb_midref(acb_imagref(h)), digits);
  887.     flint_printf("\n");
  888.  
  889. finish:
  890.  
  891.     if (result == EXIT_FAILURE)
  892.     {
  893.         flint_printf("Usage:\n");
  894.         flint_printf("%s t x y n d\n\n", argv[0]);
  895.         flint_printf(
  896.                 "Evaluate the nth derivative of H_t(x + yi) "
  897.                 "to d significant digits "
  898.                 "where H is the function involved "
  899.                 "in the definition of the De Bruijn-Newman constant.\n"
  900.                 "Requires x >= 0 and y >= 0 and t <= 1/2 and n in {0, 1}.\n");
  901.     }
  902.  
  903.     arb_clear(t);
  904.     arb_clear(T);
  905.     arb_clear(half);
  906.     arb_clear(four);
  907.  
  908.     acb_clear(z);
  909.     acb_clear(h);
  910.  
  911.     arb_clear(estimate);
  912.     mag_clear(m_estimate);
  913.  
  914.     flint_cleanup();
  915.     return result;
  916. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
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!