• API
• FAQ
• Tools
• Archive
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
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);
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);
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.
267.     if (!integral_bound_is_permissible(theta, t, a, n, prec))
268.     {
269.         result = 0;
270.         goto finish;
271.     }
272.
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);
298.     arb_mul(w1, w1, X, prec);
299.
300.     arb_init(w2);
301.     arb_mul(w2, t, theta, 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);
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);
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);
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);
377.     arb_mul(z2, z2, theta, prec);
378.
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.
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);
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);
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.
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);
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);
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.
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);
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);
798.
799.     arb_init(r);
800.     acb_abs(r, rc, prec);
801.
802.     acb_init(z1);
803.     acb_mul_arb(z1, u, t, 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);
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.
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);
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.
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.
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);
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.
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);
1064.
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);
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);
1145.
1147.     I_integral_bound(s, theta, t, x, a, n0, prec);
1148.     arb_mul(s, s, c1, prec);
1150.
1152.     I_integral_bound(s, theta, t, x, a, n0, prec);
1153.     arb_mul(s, s, c2, 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);
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.
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);
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);
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);
1276.
1278.     D_integral_bound(s, theta, t, a, n0, prec);
1279.     arb_mul(s, s, c1, prec);
1281.
1283.     D_integral_bound(s, theta, t, a, n0, prec);
1284.     arb_mul(s, s, c2, 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);
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);
1336.
1338.     D_integral(s, theta, m_estimate, t, a, n, prec);
1339.     arb_mul(s, s, c1, prec);
1341.
1343.     D_integral(s, theta, m_estimate, t, a, n, prec);
1344.     arb_mul(s, s, c2, 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);
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.
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);
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);
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);
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);
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);
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.
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);
1839.
1840.         arb_get_lbound_arf(arb_midref(x_next), x_next, prec);
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);
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);
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.
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.
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);
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.
Top