API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Apr 5th, 2018
124
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 23.84 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_poly.h"
  11.  
  12. void acb_poly_onei(acb_poly_t res)
  13. {
  14.     acb_poly_one(res);
  15.     acb_onei(acb_poly_get_coeff_ptr(res, 0));
  16. }
  17.  
  18. void
  19. acb_poly_add_acb(acb_poly_t res,
  20.         const acb_poly_t x, const acb_t y, slong prec)
  21. {
  22.     slong len = x->length;
  23.  
  24.     if (len == 0)
  25.     {
  26.         acb_poly_set_acb(res, y);
  27.     }
  28.     else
  29.     {
  30.         acb_poly_fit_length(res, len);
  31.  
  32.         acb_add(res->coeffs, x->coeffs, y, prec);
  33.  
  34.         if (res != x)
  35.             _acb_vec_set(res->coeffs + 1, x->coeffs + 1, len - 1);
  36.  
  37.         _acb_poly_set_length(res, len);
  38.         _acb_poly_normalise(res);
  39.     }
  40. }
  41.  
  42. void
  43. acb_poly_add_arb(acb_poly_t res,
  44.         const acb_poly_t x, const arb_t y, slong prec)
  45. {
  46.     acb_t z;
  47.     acb_init(z);
  48.     acb_set_arb(z, y);
  49.     acb_poly_add_acb(res, x, z, prec);
  50.     acb_clear(z);
  51. }
  52.  
  53.  
  54. void _acb_poly_alpha1_series(acb_poly_t z, const acb_poly_t s,
  55.         slong n, slong prec)
  56. {
  57.     acb_poly_t a;
  58.     acb_poly_init(a);
  59.  
  60.     acb_poly_zero(z);
  61.  
  62.     acb_poly_inv_series(a, s, n, prec);
  63.     acb_poly_scalar_mul_2exp_si(a, a, -1);
  64.     acb_poly_add_series(z, z, a, n, prec);
  65.  
  66.     acb_poly_add_si(a, s, -1, prec);
  67.     acb_poly_inv_series(a, a, n, prec);
  68.     acb_poly_add_series(z, z, a, n, prec);
  69.  
  70.     acb_poly_log_series(a, s, n, prec);
  71.     acb_poly_scalar_mul_2exp_si(a, a, -1);
  72.     acb_poly_add_series(z, z, a, n, prec);
  73.  
  74.     {
  75.         arb_t c;
  76.         arb_ptr u;
  77.         arb_init(c);
  78.         arb_const_log_sqrt2pi(c, prec);
  79.         u = acb_realref(acb_poly_get_coeff_ptr(z, 0));
  80.         arb_sub(u, u, c, prec);
  81.         arb_clear(c);
  82.     }
  83.  
  84.     acb_poly_clear(a);
  85. }
  86.  
  87. void acb_poly_alpha1_series(acb_poly_t z, const acb_poly_t s,
  88.         slong n, slong prec)
  89. {
  90.     if (z == s)
  91.     {
  92.         acb_poly_t u;
  93.         acb_poly_init(u);
  94.         _acb_poly_alpha1_series(u, s, n, prec);
  95.         acb_poly_swap(z, u);
  96.         acb_poly_clear(u);
  97.     }
  98.     else
  99.     {
  100.         _acb_poly_alpha1_series(z, s, n, prec);
  101.     }
  102. }
  103.  
  104. void _acb_poly_H01_series(acb_poly_t z, const acb_poly_t s,
  105.         slong n, slong prec)
  106. {
  107.     arb_t c;
  108.     acb_poly_t a, b;
  109.     acb_t u;
  110.     arb_ptr x;
  111.  
  112.     acb_poly_init(a);
  113.     acb_poly_init(b);
  114.  
  115.     acb_init(u);
  116.  
  117.     arb_init(c);
  118.     arb_const_log_sqrt2pi(c, prec);
  119.  
  120.     acb_poly_zero(z);
  121.  
  122.     acb_poly_add_si(a, s, -1, prec);
  123.     acb_poly_log_series(a, a, n, prec);
  124.     acb_poly_add_series(z, z, a, n, prec);
  125.  
  126.     acb_zero(u);
  127.     x = acb_realref(u);
  128.     arb_log_ui(x, 2, prec);
  129.     arb_neg(x, x);
  130.     arb_add_ui(x, x, 1, prec);
  131.     arb_mul_2exp_si(x, x, -1);
  132.     arb_add(x, x, c, prec);
  133.     acb_poly_scalar_mul(a, s, u, prec);
  134.     acb_poly_sub_series(z, z, a, n, prec);
  135.  
  136.     x = acb_realref(acb_poly_get_coeff_ptr(z, 0));
  137.     arb_add(x, x, c, prec);
  138.  
  139.     acb_poly_add_si(b, s, 1, prec);
  140.     acb_poly_scalar_mul_2exp_si(a, s, -1);
  141.     acb_poly_log_series(a, a, n, prec);
  142.     acb_poly_scalar_mul_2exp_si(a, a, -1);
  143.     acb_poly_mullow(a, a, b, n, prec);
  144.     acb_poly_add_series(z, z, a, n, prec);
  145.  
  146.     acb_poly_exp_series(z, z, n, prec);
  147.  
  148.     acb_poly_clear(a);
  149.     acb_poly_clear(b);
  150.     acb_clear(u);
  151.     arb_clear(c);
  152. }
  153.  
  154. void acb_poly_H01_series(acb_poly_t z, const acb_poly_t s,
  155.         slong n, slong prec)
  156. {
  157.     if (z == s)
  158.     {
  159.         acb_poly_t u;
  160.         acb_poly_init(u);
  161.         _acb_poly_H01_series(u, s, n, prec);
  162.         acb_poly_swap(z, u);
  163.         acb_poly_clear(u);
  164.     }
  165.     else
  166.     {
  167.         _acb_poly_H01_series(z, s, n, prec);
  168.     }
  169. }
  170.  
  171. static void
  172. _acb_poly_s_series(acb_poly_t res,
  173.         const acb_poly_t x, const acb_poly_t y,
  174.         slong n, slong prec)
  175. {
  176.     acb_poly_t ix;
  177.  
  178.     acb_poly_init(ix);
  179.     acb_poly_onei(ix);
  180.  
  181.     acb_poly_mullow(ix, ix, x, n, prec);
  182.  
  183.     acb_poly_sub_series(res, ix, y, n, prec);
  184.     acb_poly_add_si(res, res, 1, prec);
  185.     acb_poly_scalar_mul_2exp_si(res, res, -1);
  186.  
  187.     acb_poly_clear(ix);
  188. }
  189.  
  190. static void
  191. _acb_poly_abbeff_summand_series(acb_poly_t res,
  192.         const acb_t logk, const acb_poly_t a, const acb_poly_t b,
  193.         slong n, slong prec)
  194. {
  195.     acb_poly_scalar_mul(res, a, logk, prec);
  196.     acb_poly_sub_series(res, res, b, n, prec);
  197.     acb_poly_scalar_mul(res, res, logk, prec);
  198.     acb_poly_exp_series(res, res, n, prec);
  199. }
  200.  
  201.  
  202. void
  203. abbeff_series(acb_poly_t res,
  204.         const acb_poly_t t, const acb_poly_t x, const acb_poly_t y,
  205.         slong N, slong n, slong prec)
  206. {
  207.     acb_t logk;
  208.     acb_poly_t s, sc;
  209.     acb_poly_t alph1, alph2;
  210.     acb_poly_t hs, hsc;
  211.     acb_poly_t ta1, ta2;
  212.     acb_poly_t A0, B0;
  213.     acb_poly_t a, b1, b2;
  214.     acb_poly_t summand;
  215.     acb_poly_t Asum, Bsum;
  216.     acb_poly_t A, B;
  217.  
  218.     acb_poly_init(s);
  219.     _acb_poly_s_series(s, x, y, n, prec);
  220.  
  221.     acb_poly_init(sc);
  222.     acb_poly_add_si(sc, s, -1, prec);
  223.     acb_poly_neg(sc, sc);
  224.  
  225.     acb_poly_init(alph1);
  226.     acb_poly_alpha1_series(alph1, s, n, prec);
  227.  
  228.     acb_poly_init(alph2);
  229.     acb_poly_alpha1_series(alph2, sc, n, prec);
  230.  
  231.     acb_poly_init(hs);
  232.     acb_poly_H01_series(hs, s, n, prec);
  233.  
  234.     acb_poly_init(hsc);
  235.     acb_poly_H01_series(hsc, sc, n, prec);
  236.  
  237.     acb_poly_init(ta1);
  238.     acb_poly_mullow(ta1, t, alph1, n, prec);
  239.  
  240.     acb_poly_init(ta2);
  241.     acb_poly_mullow(ta2, t, alph2, n, prec);
  242.  
  243.     acb_poly_init(A0);
  244.     acb_poly_mullow(A0, ta1, alph1, n, prec);
  245.     acb_poly_scalar_mul_2exp_si(A0, A0, -2);
  246.     acb_poly_exp_series(A0, A0, n, prec);
  247.     acb_poly_mullow(A0, A0, hs, n, prec);
  248.     acb_poly_scalar_mul_2exp_si(A0, A0, -3);
  249.  
  250.     acb_poly_init(B0);
  251.     acb_poly_mullow(B0, ta2, alph2, n, prec);
  252.     acb_poly_scalar_mul_2exp_si(B0, B0, -2);
  253.     acb_poly_exp_series(B0, B0, n, prec);
  254.     acb_poly_mullow(B0, B0, hsc, n, prec);
  255.     acb_poly_scalar_mul_2exp_si(B0, B0, -3);
  256.  
  257.     acb_poly_init(a);
  258.     acb_poly_scalar_mul_2exp_si(a, t, -2);
  259.  
  260.     acb_poly_init(b1);
  261.     acb_poly_scalar_mul_2exp_si(b1, ta1, -1);
  262.     acb_poly_add_series(b1, b1, s, n, prec);
  263.  
  264.     acb_poly_init(b2);
  265.     acb_poly_scalar_mul_2exp_si(b2, ta2, -1);
  266.     acb_poly_add_series(b2, b2, sc, n, prec);
  267.  
  268.     acb_init(logk);
  269.     acb_poly_init(summand);
  270.     acb_poly_init(Asum);
  271.     acb_poly_init(Bsum);
  272.     {
  273.         slong k;
  274.         for (k = N; k >= 2; k--)
  275.         {
  276.             acb_set_si(logk, k);
  277.             acb_log(logk, logk, prec);
  278.  
  279.             _acb_poly_abbeff_summand_series(summand, logk, a, b1, n, prec);
  280.             acb_poly_add_series(Asum, Asum, summand, n, prec);
  281.  
  282.             _acb_poly_abbeff_summand_series(summand, logk, a, b2, n, prec);
  283.             acb_poly_add_series(Bsum, Bsum, summand, n, prec);
  284.         }
  285.     }
  286.     acb_poly_add_si(Asum, Asum, 1, prec);
  287.     acb_poly_add_si(Bsum, Bsum, 1, prec);
  288.  
  289.     acb_poly_init(A);
  290.     acb_poly_mullow(A, A0, Asum, n, prec);
  291.  
  292.     acb_poly_init(B);
  293.     acb_poly_mullow(B, B0, Bsum, n, prec);
  294.  
  295.     acb_poly_div_series(res, A, B0, n, prec);
  296.     acb_poly_add_series(res, res, Bsum, n, prec);
  297.  
  298.     acb_clear(logk);
  299.     acb_poly_clear(s);
  300.     acb_poly_clear(sc);
  301.     acb_poly_clear(alph1);
  302.     acb_poly_clear(alph2);
  303.     acb_poly_clear(hs);
  304.     acb_poly_clear(hsc);
  305.     acb_poly_clear(ta1);
  306.     acb_poly_clear(ta2);
  307.     acb_poly_clear(A0);
  308.     acb_poly_clear(B0);
  309.     acb_poly_clear(a);
  310.     acb_poly_clear(b1);
  311.     acb_poly_clear(b2);
  312.     acb_poly_clear(summand);
  313.     acb_poly_clear(Asum);
  314.     acb_poly_clear(Bsum);
  315.     acb_poly_clear(A);
  316.     acb_poly_clear(B);
  317. }
  318.  
  319.  
  320.  
  321.  
  322. void arg2(arb_t res, const acb_t v, const acb_t u, slong prec)
  323. {
  324.     if (acb_is_zero(v))
  325.     {
  326.         acb_arg(res, u, prec);
  327.     }
  328.     else
  329.     {
  330.         acb_t p;
  331.  
  332.         acb_init(p);
  333.  
  334.         acb_conj(p, v);
  335.         acb_mul(p, p, u, prec);
  336.         acb_arg(res, p, prec);
  337.  
  338.         acb_clear(p);
  339.     }
  340. }
  341.  
  342. static void
  343. _abbeff(acb_t res,
  344.         const arb_t t, const arb_t x, const arb_t y,
  345.         slong N, slong prec)
  346. {
  347.     acb_poly_t pt, px, py, pr;
  348.     slong n = 1;
  349.  
  350.     acb_poly_init(pt);
  351.     acb_poly_init(px);
  352.     acb_poly_init(py);
  353.     acb_poly_init(pr);
  354.  
  355.     acb_poly_one(pt);
  356.     arb_set(acb_realref(acb_poly_get_coeff_ptr(pt, 0)), t);
  357.  
  358.     acb_poly_one(px);
  359.     arb_set(acb_realref(acb_poly_get_coeff_ptr(px, 0)), x);
  360.  
  361.     acb_poly_one(py);
  362.     arb_set(acb_realref(acb_poly_get_coeff_ptr(py, 0)), y);
  363.  
  364.     abbeff_series(pr, pt, px, py, N, n, prec);
  365.     acb_set(res, acb_poly_get_coeff_ptr(pr, 0));
  366.  
  367.     acb_poly_clear(pt);
  368.     acb_poly_clear(px);
  369.     acb_poly_clear(py);
  370.     acb_poly_clear(pr);
  371. }
  372.  
  373.  
  374. typedef struct
  375. {
  376.     acb_struct z;
  377.     acb_struct h;
  378. } meshpt_struct;
  379. typedef meshpt_struct meshpt_t[1];
  380. typedef meshpt_struct *meshpt_ptr;
  381. typedef const meshpt_struct *meshpt_srcptr;
  382.  
  383. #define meshpt_zref(p) (&(p)->z)
  384. #define meshpt_href(p) (&(p)->h)
  385.  
  386. void meshpt_init(meshpt_t p)
  387. {
  388.     acb_init(meshpt_zref(p));
  389.     acb_init(meshpt_href(p));
  390. }
  391.  
  392. void meshpt_clear(meshpt_t p)
  393. {
  394.     acb_clear(meshpt_zref(p));
  395.     acb_clear(meshpt_href(p));
  396. }
  397.  
  398. void meshpt_swap(meshpt_t a, meshpt_t b)
  399. {
  400.     acb_swap(meshpt_zref(a), meshpt_zref(b));
  401.     acb_swap(meshpt_href(a), meshpt_href(b));
  402. }
  403.  
  404.  
  405.  
  406. typedef struct
  407. {
  408.     meshpt_struct left;
  409.     meshpt_srcptr right;
  410. } meshival_struct;
  411. typedef meshival_struct meshival_t[1];
  412. typedef meshival_struct *meshival_ptr;
  413. typedef const meshival_struct *meshival_srcptr;
  414.  
  415. #define meshival_leftref(m) (&(m)->left)
  416.  
  417. void meshival_init(meshival_t m)
  418. {
  419.     meshpt_init(meshival_leftref(m));
  420.     m->right = NULL;
  421. }
  422.  
  423. void meshival_clear(meshival_t m)
  424. {
  425.     meshpt_clear(meshival_leftref(m));
  426.     m->right = NULL;
  427. }
  428.  
  429. void meshival_center(acb_t res, const meshival_t m, slong prec)
  430. {
  431.     acb_srcptr z_left = meshpt_zref(meshival_leftref(m));
  432.     acb_srcptr z_right = meshpt_zref(m->right);
  433.     acb_add(res, z_left, z_right, prec);
  434.     acb_mul_2exp_si(res, res, -1);
  435. }
  436.  
  437.  
  438. void compute_lbound(arb_t res, const acb_t fa, const acb_t fb,
  439.         const acb_t fprime, const acb_t delta, slong prec)
  440. {
  441.     arb_t afa, afb, afprime, adelta;
  442.  
  443.     arb_init(afa);
  444.     arb_init(afb);
  445.     arb_init(afprime);
  446.     arb_init(adelta);
  447.  
  448.     acb_abs(afa, fa, prec);
  449.     acb_abs(afb, fb, prec);
  450.     acb_abs(afprime, fprime, prec);
  451.     acb_abs(adelta, delta, prec);
  452.  
  453.     arb_add(res, afa, afb, prec);
  454.     arb_submul(res, afprime, adelta, prec);
  455.     arb_mul_2exp_si(res, res, -1);
  456.  
  457.     arb_clear(afa);
  458.     arb_clear(afb);
  459.     arb_clear(afprime);
  460.     arb_clear(adelta);
  461. }
  462.  
  463. int meshival_validate(const meshival_t m, const arb_t t,
  464.         const arb_t c, slong N, slong prec)
  465. {
  466.     acb_poly_t hpoly, tpoly, xpoly, ypoly;
  467.     acb_t z, f, fprime, step_delta;
  468.     arb_t lbound;
  469.     acb_srcptr z_left, z_right;
  470.     acb_srcptr h_left, h_right;
  471.     arb_srcptr x, y;
  472.     slong n = 2;
  473.     int result = 0;
  474.  
  475.     z_left = meshpt_zref(meshival_leftref(m));
  476.     z_right = meshpt_zref(m->right);
  477.  
  478.     h_left = meshpt_href(meshival_leftref(m));
  479.     h_right = meshpt_href(m->right);
  480.  
  481.     acb_init(step_delta);
  482.     acb_sub(step_delta, z_right, z_left, prec);
  483.  
  484.     acb_init(z);
  485.     acb_union(z, z_left, z_right, prec);
  486.     x = acb_realref(z);
  487.     y = acb_imagref(z);
  488.  
  489.     acb_poly_init(tpoly);
  490.     acb_poly_add_arb(tpoly, tpoly, t, prec);
  491.  
  492.     acb_poly_init(xpoly);
  493.     acb_poly_one(xpoly);
  494.     acb_poly_shift_left(xpoly, xpoly, 1);
  495.     acb_poly_add_arb(xpoly, xpoly, x, prec);
  496.  
  497.     acb_poly_init(ypoly);
  498.     acb_poly_add_arb(ypoly, ypoly, y, prec);
  499.  
  500.     acb_poly_init(hpoly);
  501.     abbeff_series(hpoly, tpoly, xpoly, ypoly, N, n, prec);
  502.  
  503.     acb_init(f);
  504.     acb_poly_get_coeff_acb(f, hpoly, 0);
  505.  
  506.     acb_init(fprime);
  507.     acb_poly_get_coeff_acb(fprime, hpoly, 1);
  508.  
  509.     arb_init(lbound);
  510.     acb_abs(lbound, f, prec);
  511.  
  512.     if (arb_gt(lbound, c))
  513.     {
  514.         result = 1;
  515.     }
  516.     else
  517.     {
  518.         compute_lbound(lbound, h_left, h_right, fprime, step_delta, prec);
  519.         if (arb_gt(lbound, c))
  520.         {
  521.             result = 1;
  522.         }
  523.     }
  524.  
  525. finish:
  526.  
  527.     acb_poly_clear(hpoly);
  528.     acb_poly_clear(tpoly);
  529.     acb_poly_clear(xpoly);
  530.     acb_poly_clear(ypoly);
  531.  
  532.     acb_clear(z);
  533.     acb_clear(f);
  534.     acb_clear(fprime);
  535.     acb_clear(step_delta);
  536.  
  537.     arb_clear(lbound);
  538.  
  539.     return result;
  540. }
  541.  
  542.  
  543. void
  544. evaluate_meshpt(meshpt_t m,
  545.         const arb_t t, const arb_t x, const arb_t y,
  546.         slong N, slong prec)
  547. {
  548.     acb_ptr z = meshpt_zref(m);
  549.     acb_ptr h = meshpt_href(m);
  550.     acb_set_arb_arb(z, x, y);
  551.     _abbeff(h, t, x, y, N, prec);
  552. }
  553.  
  554.  
  555. int validate_f(const acb_t f, const arb_t c, slong prec)
  556. {
  557.     arb_t af;
  558.  
  559.     arb_init(af);
  560.     acb_abs(af, f, prec);
  561.     int result = 1;
  562.  
  563.     if (arb_lt(af, c))
  564.     {
  565.         fprintf(stderr, "Failed to compute the winding number because "
  566.                 "|(Aeff + Beff)/Beff0| < c\n");
  567.         result = -1;
  568.         goto finish;
  569.     }
  570.  
  571.     if (!arb_gt(af, c))
  572.     {
  573.         result = 0;
  574.         goto finish;
  575.     }
  576.  
  577.     if (arb_contains_zero(acb_imagref(f)) && !arb_is_positive(acb_realref(f)))
  578.     {
  579.         result = 0;
  580.         goto finish;
  581.     }
  582.  
  583. finish:
  584.  
  585.     arb_clear(af);
  586.     return result;
  587. }
  588.  
  589. int
  590. evaluate_mesh(arb_t res,
  591.         const arb_t t, const acb_t za, const acb_t zb,
  592.         const arb_t c, slong maxdepth, slong N, slong prec)
  593. {
  594.     acb_t center;
  595.     arb_t accum, arg;
  596.     meshpt_t mb;
  597.     meshival_struct *stack;
  598.     slong i, depth;
  599.     meshival_ptr ival;
  600.     int result = 1;
  601.  
  602.     stack = flint_malloc(maxdepth * sizeof(*stack));
  603.     for (i = 0; i < maxdepth; i++)
  604.     {
  605.         meshival_init(stack + i);
  606.     }
  607.  
  608.     acb_init(center);
  609.     arb_init(accum);
  610.     arb_init(arg);
  611.     meshpt_init(mb);
  612.  
  613.     evaluate_meshpt(mb,
  614.             t, acb_realref(zb), acb_imagref(zb), N, prec);
  615.  
  616.     result = validate_f(meshpt_href(mb), c, prec);
  617.     if (result < 1)
  618.     {
  619.         arb_indeterminate(res);
  620.         goto finish;
  621.     }
  622.  
  623.     depth = 0;
  624.     ival = stack + depth;
  625.     evaluate_meshpt(meshival_leftref(ival),
  626.             t, acb_realref(za), acb_imagref(za), N, prec);
  627.     ival->right = mb;
  628.     depth++;
  629.  
  630.     result = validate_f(meshpt_href(meshival_leftref(ival)), c, prec);
  631.     if (result < 1)
  632.     {
  633.         arb_indeterminate(res);
  634.         goto finish;
  635.     }
  636.  
  637.     while (1)
  638.     {
  639.         meshival_ptr curr = stack + depth - 1;
  640.         meshival_ptr next = stack + depth;
  641.  
  642.         if (depth + 5 > maxdepth)
  643.         {
  644.             result = 0;
  645.             arb_indeterminate(res);
  646.             goto finish;
  647.         }
  648.  
  649.         if (!depth)
  650.         {
  651.             break;
  652.         }
  653.  
  654.         if (meshival_validate(curr, t, c, N, prec))
  655.         {
  656.             acb_srcptr h_left = meshpt_href(meshival_leftref(curr));
  657.             acb_srcptr h_right = meshpt_href(curr->right);
  658.  
  659.             arg2(arg, h_left, h_right, prec);
  660.             arb_add(accum, accum, arg, prec);
  661.  
  662.             if (!arb_is_finite(accum))
  663.             {
  664.                 result = 0;
  665.                 arb_indeterminate(res);
  666.                 goto finish;
  667.             }
  668.  
  669.             depth--;
  670.         }
  671.         else
  672.         {
  673.             meshival_center(center, curr, prec);
  674.  
  675.             evaluate_meshpt(meshival_leftref(next),
  676.                     t, acb_realref(center), acb_imagref(center), N, prec);
  677.  
  678.             result = validate_f(meshpt_href(meshival_leftref(next)), c, prec);
  679.             if (result < 1)
  680.             {
  681.                 arb_indeterminate(res);
  682.                 goto finish;
  683.             }
  684.  
  685.             meshpt_swap(meshival_leftref(curr), meshival_leftref(next));
  686.  
  687.             next->right = meshival_leftref(curr);
  688.  
  689.             depth++;
  690.         }
  691.     }
  692.  
  693.     arb_set(res, accum);
  694.  
  695. finish:
  696.  
  697.     acb_clear(center);
  698.     arb_clear(accum);
  699.     arb_clear(arg);
  700.     meshpt_clear(mb);
  701.  
  702.     for (i = 0; i < maxdepth; i++)
  703.     {
  704.         meshival_clear(stack + i);
  705.     }
  706.     flint_free(stack);
  707.  
  708.     return result;
  709. }
  710.  
  711. slong get_N(const arb_t t, const arb_t x, slong prec)
  712. {
  713.     arb_t pi, u;
  714.     slong N;
  715.     slong result;
  716.  
  717.     arb_init(pi);
  718.     arb_init(u);
  719.  
  720.     arb_const_pi(pi, prec);
  721.     arb_mul(u, pi, t, prec);
  722.     arb_mul_2exp_si(u, u, -2);
  723.     arb_add(u, u, x, prec);
  724.     arb_div(u, u, pi, prec);
  725.     arb_mul_2exp_si(u, u, -2);
  726.     arb_sqrt(u, u, prec);
  727.     arb_floor(u, u, prec);
  728.  
  729.     N = (slong) arf_get_d(arb_midref(u), ARF_RND_DOWN);
  730.  
  731.     if (arb_contains_si(u, N) &&
  732.         !arb_contains_si(u, N-1) &&
  733.         !arb_contains_si(u, N+1))
  734.     {
  735.         result = N;
  736.     }
  737.     else
  738.     {
  739.         fprintf(stderr, "Unexpected error: could not compute N\n");
  740.         flint_abort();
  741.     }
  742.    
  743.     arb_clear(pi);
  744.     arb_clear(u);
  745.  
  746.     return result;
  747. }
  748.  
  749.  
  750.  
  751. int main(int argc, char *argv[])
  752. {
  753.     arb_t c, t, xa, xb, ya, yb;
  754.     arb_t arg, arg_total;
  755.     const char *N_str, *c_str, *t_str, *xa_str, *xb_str, *ya_str, *yb_str;
  756.     acb_t z1, z2;
  757.     slong N_user, N, prec;
  758.     slong prec_base, prec_incr;
  759.     slong maxdepth = 50;
  760.     int mesh_result;
  761.     int result = EXIT_SUCCESS;
  762.     int usage_error = 0;
  763.  
  764.     arb_init(c);
  765.     arb_init(t);
  766.     arb_init(xa);
  767.     arb_init(xb);
  768.     arb_init(ya);
  769.     arb_init(yb);
  770.     acb_init(z1);
  771.     acb_init(z2);
  772.  
  773.     arb_init(arg);
  774.     arb_init(arg_total);
  775.  
  776.     if (argc != 8)
  777.     {
  778.         usage_error = 1;
  779.         result = EXIT_FAILURE;
  780.         goto finish;
  781.     }
  782.  
  783.     N_str = argv[1];
  784.     c_str = argv[2];
  785.     t_str = argv[3];
  786.     xa_str = argv[4];
  787.     xb_str = argv[5];
  788.     ya_str = argv[6];
  789.     yb_str = argv[7];
  790.  
  791.     prec_base = 32;
  792.     prec_incr = 32;
  793.  
  794.     N_user = atol(N_str);
  795.     if (N_user > 0)
  796.     {
  797.         N = N_user;
  798.     }
  799.     else
  800.     {
  801.         prec = 128;
  802.         slong N_low, N_high;
  803.         arb_set_str(t, t_str, prec);
  804.         arb_set_str(xa, xa_str, prec);
  805.         arb_set_str(xb, xb_str, prec);
  806.         if (arb_lt(xb, xa))
  807.         {
  808.             usage_error = 1;
  809.             result = EXIT_FAILURE;
  810.             goto finish;
  811.         }
  812.         N_low = get_N(t, xa, prec);
  813.         N_high = get_N(t, xb, prec);
  814.         if (N_low == N_high)
  815.         {
  816.             N = N_low;
  817.         }
  818.         else
  819.         {
  820.             fprintf(stderr, "Imputed N is ambiguous: [%ld, %ld].\n",
  821.                     N_low, N_high);
  822.             fprintf(stderr, "Either specify a positive integer N explicitly, "
  823.                     "or choose a smaller range of x.\n");
  824.             usage_error = 1;
  825.             result = EXIT_FAILURE;
  826.             goto finish;
  827.         }
  828.     }
  829.  
  830.     prec = prec_base;
  831.     while (1)
  832.     {
  833.         arb_set_str(c, c_str, prec);
  834.         arb_set_str(t, t_str, prec);
  835.         arb_set_str(xa, xa_str, prec);
  836.         arb_set_str(xb, xb_str, prec);
  837.         arb_set_str(ya, ya_str, prec);
  838.         arb_set_str(yb, yb_str, prec);
  839.         if (arb_lt(xb, xa) || arb_lt(yb, ya))
  840.         {
  841.             usage_error = 1;
  842.             result = EXIT_FAILURE;
  843.             goto finish;
  844.         }
  845.         acb_set_arb_arb(z1, xa, ya);
  846.         acb_set_arb_arb(z2, xb, ya);
  847.         mesh_result = evaluate_mesh(arg, t, z1, z2, c, maxdepth, N, prec);
  848.         if (mesh_result == -1)
  849.         {
  850.             result = EXIT_FAILURE;
  851.             goto finish;
  852.         }
  853.         else if (mesh_result == 1)
  854.         {
  855.             break;
  856.         }
  857.         prec += prec_incr;
  858.     }
  859.     arb_add(arg_total, arg_total, arg, prec);
  860.  
  861.     prec = prec_base;
  862.     while (1)
  863.     {
  864.         arb_set_str(c, c_str, prec);
  865.         arb_set_str(t, t_str, prec);
  866.         arb_set_str(xa, xa_str, prec);
  867.         arb_set_str(xb, xb_str, prec);
  868.         arb_set_str(ya, ya_str, prec);
  869.         arb_set_str(yb, yb_str, prec);
  870.         if (arb_lt(xb, xa) || arb_lt(yb, ya))
  871.         {
  872.             usage_error = 1;
  873.             result = EXIT_FAILURE;
  874.             goto finish;
  875.         }
  876.         acb_set_arb_arb(z1, xb, ya);
  877.         acb_set_arb_arb(z2, xb, yb);
  878.         mesh_result = evaluate_mesh(arg, t, z1, z2, c, maxdepth, N, prec);
  879.         if (mesh_result == -1)
  880.         {
  881.             result = EXIT_FAILURE;
  882.             goto finish;
  883.         }
  884.         else if (mesh_result == 1)
  885.         {
  886.             break;
  887.         }
  888.         prec += prec_incr;
  889.     }
  890.     arb_add(arg_total, arg_total, arg, prec);
  891.  
  892.     prec = prec_base;
  893.     while (1)
  894.     {
  895.         arb_set_str(c, c_str, prec);
  896.         arb_set_str(t, t_str, prec);
  897.         arb_set_str(xa, xa_str, prec);
  898.         arb_set_str(xb, xb_str, prec);
  899.         arb_set_str(ya, ya_str, prec);
  900.         arb_set_str(yb, yb_str, prec);
  901.         if (arb_lt(xb, xa) || arb_lt(yb, ya) || arb_is_negative(c))
  902.         {
  903.             usage_error = 1;
  904.             result = EXIT_FAILURE;
  905.             goto finish;
  906.         }
  907.         acb_set_arb_arb(z1, xb, yb);
  908.         acb_set_arb_arb(z2, xa, yb);
  909.         mesh_result = evaluate_mesh(arg, t, z1, z2, c, maxdepth, N, prec);
  910.         if (mesh_result == -1)
  911.         {
  912.             result = EXIT_FAILURE;
  913.             goto finish;
  914.         }
  915.         else if (mesh_result == 1)
  916.         {
  917.             break;
  918.         }
  919.         prec += prec_incr;
  920.     }
  921.     arb_add(arg_total, arg_total, arg, prec);
  922.  
  923.     prec = prec_base;
  924.     while (1)
  925.     {
  926.         arb_set_str(c, c_str, prec);
  927.         arb_set_str(t, t_str, prec);
  928.         arb_set_str(xa, xa_str, prec);
  929.         arb_set_str(xb, xb_str, prec);
  930.         arb_set_str(ya, ya_str, prec);
  931.         arb_set_str(yb, yb_str, prec);
  932.         if (arb_lt(xb, xa) || arb_lt(yb, ya))
  933.         {
  934.             usage_error = 1;
  935.             result = EXIT_FAILURE;
  936.             goto finish;
  937.         }
  938.         acb_set_arb_arb(z1, xa, yb);
  939.         acb_set_arb_arb(z2, xa, ya);
  940.         mesh_result = evaluate_mesh(arg, t, z1, z2, c, maxdepth, N, prec);
  941.         if (mesh_result == -1)
  942.         {
  943.             result = EXIT_FAILURE;
  944.             goto finish;
  945.         }
  946.         else if (mesh_result == 1)
  947.         {
  948.             break;
  949.         }
  950.         prec += prec_incr;
  951.     }
  952.     arb_add(arg_total, arg_total, arg, prec);
  953.  
  954.     {
  955.         arb_t pi;
  956.         arb_init(pi);
  957.         arb_const_pi(pi, prec);
  958.         arb_div(arg_total, arg_total, pi, prec);
  959.         arb_mul_2exp_si(arg_total, arg_total, -1);
  960.         arb_clear(pi);
  961.     }
  962.  
  963.     flint_printf("(Aeff+Beff)/Beff0 winding number script\n");
  964.  
  965.     if (N_user <= 0)
  966.     {
  967.         flint_printf("imputed N: %ld\n", N);
  968.     }
  969.  
  970.     {
  971.         double d = arf_get_d(arb_midref(arg_total), ARF_RND_NEAR);
  972.         slong nroots = (int) (d + 0.5);
  973.         flint_printf("winding number: ");
  974.         if (arb_contains_si(arg_total, nroots) &&
  975.             !arb_contains_si(arg_total, nroots - 1) &&
  976.             !arb_contains_si(arg_total, nroots + 1))
  977.         {
  978.             flint_printf("%ld\n", nroots);
  979.         }
  980.         else
  981.         {
  982.             arb_printd(arg_total, 10);
  983.             flint_printf("\n");
  984.         }
  985.     }
  986.  
  987. finish:
  988.  
  989.     if (usage_error && result == EXIT_FAILURE)
  990.     {
  991.         flint_printf("Usage:\n");
  992.         flint_printf("%s N c t xa xb ya yb\n\n", argv[0]);
  993.         flint_printf(
  994.     "This script computes the winding number about 0 of f(p).\n"
  995.     "Let r be a rectangle in the complex plane\n"
  996.     "r : [xa <= x <= xb] + i[ya <= y <= yb]\n"
  997.     "and let p be the counterclockwise path along the boundary of r\n"
  998.     "starting and ending at (xa + i*ya).\n"
  999.     "Let f(z) = (Aeff_t(z) + Beff_t(z)) / Beff0_t(z).\n"
  1000.     "Additionally, f(p) is required to stay away from 0 so\n"
  1001.     "the script will either abort or will fail to terminate\n"
  1002.     "if |f(z)| <= c at any point z on the path p.\n"
  1003.     "The Nth partial sum is used in Aeff and Beff.\n");
  1004.     }
  1005.  
  1006.     arb_clear(c);
  1007.     arb_clear(t);
  1008.     arb_clear(xa);
  1009.     arb_clear(xb);
  1010.     arb_clear(ya);
  1011.     arb_clear(yb);
  1012.  
  1013.     arb_clear(arg);
  1014.     arb_clear(arg_total);
  1015.  
  1016.     acb_clear(z1);
  1017.     acb_clear(z2);
  1018.  
  1019.     flint_cleanup();
  1020.     return result;
  1021. }
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!