Guest User

Untitled

a guest
Mar 18th, 2018
75
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. void H_approx(arb_t res, const acb_t z, slong prec)
  13. {
  14.     arb_t logx, pixd8;
  15.  
  16.     arb_init(logx);
  17.     arb_log(logx, acb_realref(z), prec);
  18.  
  19.     arb_init(pixd8);
  20.     arb_const_pi(pixd8, prec);
  21.     arb_mul(pixd8, pixd8, acb_realref(z), prec);
  22.     arb_mul_2exp_si(pixd8, pixd8, -3);
  23.  
  24.     arb_add_si(res, acb_imagref(z), 7, prec);
  25.     arb_mul(res, res, logx, prec);
  26.     arb_mul_2exp_si(res, res, -2);
  27.     arb_sub(res, res, pixd8, prec);
  28.     arb_exp(res, res, prec);
  29.  
  30.     arb_clear(logx);
  31.     arb_clear(pixd8);
  32. }
  33.  
  34. void integral_remainder(mag_t m, const arb_t b)
  35. {
  36.     arb_get_mag(m, b);
  37.     mag_mul_2exp_si(m, m, 2);
  38.     mag_exp_lower(m, m);
  39.     mag_expinv(m, m);
  40.     mag_mul_ui(m, m, 25);
  41. }
  42.  
  43. void integration_limit(arb_t b, const arb_t eps, slong prec)
  44. {
  45.     arb_zero(b);
  46.     arb_get_lbound_arf(arb_midref(b), eps, prec);
  47.     arb_div_ui(b, b, 25, prec);
  48.     arb_log(b, b, prec);
  49.     arb_neg(b, b);
  50.     arb_log(b, b, prec);
  51.     arb_mul_2exp_si(b, b, -2);
  52.     arb_get_ubound_arf(arb_midref(b), b, prec);
  53.     arb_get_mid_arb(b, b);
  54. }
  55.  
  56. void phi_domain(acb_t z)
  57. {
  58.     arb_ptr x, y;
  59.     slong prec;
  60.     arb_t t;
  61.  
  62.     prec = 64;
  63.  
  64.     x = acb_realref(z);
  65.     arb_zero_pm_inf(x);
  66.  
  67.     y = acb_imagref(z);
  68.     arb_zero(y);
  69.  
  70.     arb_init(t);
  71.     arb_const_pi(t, prec);
  72.     arb_div_ui(t, t, 12, prec);
  73.     arb_get_mag_lower(arb_radref(y), t);
  74.     arb_clear(t);
  75. }
  76.  
  77. int phi_domain_contains(const acb_t u)
  78. {
  79.     int result;
  80.     acb_t d;
  81.     acb_init(d);
  82.     phi_domain(d);
  83.     result = acb_contains(d, u);
  84.     acb_clear(d);
  85.     return result;
  86. }
  87.  
  88. void phi_remainder(arb_t z, const acb_t u, ulong n0, slong prec)
  89. {
  90.     if (n0 < 1)
  91.         flint_abort();
  92.  
  93.     if (phi_domain_contains(u))
  94.     {
  95.         acb_t exp4u;
  96.         arb_t b, d;
  97.  
  98.         acb_init(exp4u);
  99.         acb_mul_2exp_si(exp4u, u, 2);
  100.         acb_exp(exp4u, exp4u, prec);
  101.  
  102.         arb_init(b);
  103.         arb_const_pi(b, prec);
  104.         arb_mul_ui(b, b, n0, prec);
  105.         arb_mul_2exp_si(b, b, -1);
  106.         arb_mul(b, b, acb_realref(exp4u), prec);
  107.  
  108.         arb_init(d);
  109.         arb_neg(d, b);
  110.         arb_exp(d, d, prec);
  111.         arb_sub_ui(d, d, 1, prec);
  112.  
  113.         arb_mul_ui(z, b, n0, prec);
  114.         arb_neg(z, z);
  115.         arb_exp(z, z, prec);
  116.         arb_div(z, z, d, prec);
  117.         arb_mul_si(z, z, -50, prec);
  118.  
  119.         acb_clear(exp4u);
  120.         arb_clear(b);
  121.         arb_clear(d);
  122.     }
  123.     else
  124.     {
  125.         arb_indeterminate(z);
  126.     }
  127. }
  128.  
  129. void phi_summand(acb_t z, const acb_t u, ulong n, slong prec)
  130. {
  131.     if (phi_domain_contains(u))
  132.     {
  133.         arb_t pin2;
  134.         acb_t a, b;
  135.  
  136.         arb_init(pin2);
  137.         arb_const_pi(pin2, prec);
  138.         arb_mul_ui(pin2, pin2, n*n, prec);
  139.  
  140.         acb_init(a);
  141.         acb_mul_2exp_si(a, u, 2);
  142.         acb_exp(a, a, prec);
  143.         acb_mul_arb(a, a, pin2, prec);
  144.  
  145.         acb_init(b);
  146.         acb_mul_2exp_si(b, a, 1);
  147.         acb_sub_ui(b, b, 3, prec);
  148.  
  149.         acb_mul_ui(z, u, 5, prec);
  150.         acb_sub(z, z, a, prec);
  151.         acb_exp(z, z, prec);
  152.         acb_mul(z, z, b, prec);
  153.         acb_mul_arb(z, z, pin2, prec);
  154.  
  155.         arb_clear(pin2);
  156.         acb_clear(a);
  157.         acb_clear(b);
  158.     }
  159.     else
  160.     {
  161.         acb_indeterminate(z);
  162.     }
  163. }
  164.  
  165. void phi(acb_t z, const acb_t u, slong prec)
  166. {
  167.     if (phi_domain_contains(u))
  168.     {
  169.         ulong n;
  170.         acb_t s;
  171.         arb_t tail;
  172.         arf_t err;
  173.         mag_t m, tail_mag_lb, tail_mag_ub;
  174.  
  175.         acb_init(s);
  176.         arb_init(tail);
  177.         arf_init(err);
  178.         mag_init(tail_mag_lb);
  179.         mag_init(tail_mag_ub);
  180.         mag_init(m);
  181.  
  182.         n = 1;
  183.         acb_zero(z);
  184.         while (1)
  185.         {
  186.             phi_summand(s, u, n, prec);
  187.             acb_add(z, z, s, prec);
  188.  
  189.             n++;
  190.  
  191.             phi_remainder(tail, u, n, prec);
  192.             acb_get_rad_ubound_arf(err, z, prec);
  193.             arf_get_mag(m, err);
  194.  
  195.             arb_get_mag_lower(tail_mag_lb, tail);
  196.             if (mag_cmp(tail_mag_lb, m) <= 0)
  197.             {
  198.                 arb_get_mag(tail_mag_ub, tail);
  199.                 acb_add_error_mag(z, tail_mag_ub);
  200.                 break;
  201.             }
  202.         }
  203.  
  204.         acb_clear(s);
  205.         arb_clear(tail);
  206.         arf_clear(err);
  207.         mag_clear(tail_mag_lb);
  208.         mag_clear(tail_mag_ub);
  209.         mag_clear(m);
  210.     }
  211.     else
  212.     {
  213.         acb_indeterminate(z);
  214.     }
  215. }
  216.  
  217. int integrand(acb_ptr res, const acb_t u, void *param, slong order, slong prec)
  218. {
  219.     if (phi_domain_contains(u))
  220.     {
  221.         acb_t a, b, c;
  222.         arb_ptr t;
  223.         acb_ptr z;
  224.  
  225.         t = ((arb_ptr *)(param))[0];
  226.         z = ((acb_ptr *)(param))[1];
  227.  
  228.         acb_init(a);
  229.         acb_sqr(a, u, prec);
  230.         acb_mul_arb(a, a, t, prec);
  231.         acb_exp(a, a, prec);
  232.  
  233.         acb_init(b);
  234.         phi(b, u, prec);
  235.  
  236.         acb_init(c);
  237.         acb_mul(c, z, u, prec);
  238.         acb_cos(c, c, prec);
  239.  
  240.         acb_mul(res, a, b, prec);
  241.         acb_mul(res, res, c, prec);
  242.  
  243.         acb_clear(a);
  244.         acb_clear(b);
  245.         acb_clear(c);
  246.     }
  247.     else
  248.     {
  249.         acb_indeterminate(res);
  250.     }
  251.     return 0;
  252. }
  253.  
  254.  
  255. void H(acb_t h, const arb_t t, const acb_t z, slong prec)
  256. {
  257.     arb_t width;
  258.     acb_t a, b;
  259.     void *param[2];
  260.     slong rel_goal;
  261.     int calc_result;
  262.     mag_t m, abs_tol;
  263.  
  264.     rel_goal = prec;
  265.     param[0] = (void *) t;
  266.     param[1] = (void *) z;
  267.  
  268.     arb_init(width);
  269.     H_approx(width, z, prec);
  270.     arb_mul_2exp_si(width, width, -rel_goal);
  271.  
  272.     acb_init(a);
  273.     acb_init(b);
  274.     integration_limit(acb_realref(b), width, prec);
  275.  
  276.     mag_init(abs_tol);
  277.     arb_get_mag(abs_tol, width);
  278.  
  279.     calc_result = acb_calc_integrate(
  280.             h, integrand, param, a, b,
  281.             rel_goal, abs_tol, NULL, prec);
  282.     if (calc_result == ARB_CALC_NO_CONVERGENCE)
  283.     {
  284.         acb_indeterminate(h);
  285.         goto finish;
  286.     }
  287.  
  288.     mag_init(m);
  289.     integral_remainder(m, acb_realref(b));
  290.     acb_add_error_mag(h, m);
  291.  
  292. finish:
  293.  
  294.     acb_clear(a);
  295.     acb_clear(b);
  296.     arb_clear(width);
  297.     mag_clear(abs_tol);
  298.     mag_clear(m);
  299. }
  300.  
  301. int main(int argc, char *argv[])
  302. {
  303.     int result;
  304.     arb_t t;
  305.     acb_t z, h;
  306.     slong digits, prec, high_prec;
  307.     arb_t half, yabs;
  308.     double dx;
  309.  
  310.     arb_init(t);
  311.     acb_init(z);
  312.     acb_init(h);
  313.     arb_init(half);
  314.     arb_init(yabs);
  315.    
  316.     result = EXIT_SUCCESS;
  317.  
  318.     if (argc != 5)
  319.     {
  320.         result = EXIT_FAILURE;
  321.         goto finish;
  322.     }
  323.  
  324.     digits = atol(argv[4]);
  325.     high_prec = digits * 3.32192809488736 + 10;
  326.  
  327.     arb_one(half);
  328.     arb_mul_2exp_si(half, half, -1);
  329.  
  330.     dx = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
  331.     prec = FLINT_MAX(32, (int) (dx / 4));
  332.     acb_indeterminate(h);
  333.     while (acb_rel_accuracy_bits(h) < high_prec)
  334.     {
  335.         prec *= 2;
  336.         arb_set_str(t, argv[1], prec);
  337.         arb_set_str(acb_realref(z), argv[2], prec);
  338.         arb_set_str(acb_imagref(z), argv[3], prec);
  339.         arb_abs(yabs, acb_imagref(z));
  340.         if (!arb_lt(t, half) || !arb_lt(yabs, half))
  341.         {
  342.             result = EXIT_FAILURE;
  343.             goto finish;
  344.         }
  345.         H(h, t, z, prec);
  346.     }
  347.  
  348.     flint_printf("Re[H]: ");
  349.     arf_printd(arb_midref(acb_realref(h)), digits);
  350.     flint_printf("\n");
  351.     flint_printf("Im[H]: ");
  352.     arf_printd(arb_midref(acb_imagref(h)), digits);
  353.     flint_printf("\n");
  354.  
  355. finish:
  356.  
  357.     if (result == EXIT_FAILURE)
  358.     {
  359.         flint_printf(
  360.             "%s t x y d\n\n"
  361.             "Evaluate H_t(x + yi) to d significant digits "
  362.             "where H is the function involved "
  363.             "in the definition of the De Bruijn-Newman constant.\n"
  364.             "This implementation requires t < 0.5 and |y| < 0.5.\n",
  365.             argv[0]);
  366.     }
  367.  
  368.     arb_clear(t);
  369.     acb_clear(z);
  370.     acb_clear(h);
  371.     arb_clear(half);
  372.     arb_clear(yabs);
  373.  
  374.     flint_cleanup();
  375.     return result;
  376. }
RAW Paste Data