1 /*
2 Rigorous numerical integration (with fast convergence for piecewise
3 holomorphic functions) using Gauss-Legendre quadrature and adaptive
4 subdivision.
5
6 Author: Fredrik Johansson.
7 This file is in the public domain.
8 */
9
10 #include <string.h>
11 #include "flint/profiler.h"
12 #include "arb_hypgeom.h"
13 #include "acb_hypgeom.h"
14 #include "acb_dirichlet.h"
15 #include "acb_modular.h"
16 #include "acb_calc.h"
17
18 /* ------------------------------------------------------------------------- */
19 /* Example integrands */
20 /* ------------------------------------------------------------------------- */
21
22 /* f(z) = sin(z) */
23 int
f_sin(acb_ptr res,const acb_t z,void * param,slong order,slong prec)24 f_sin(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
25 {
26 if (order > 1)
27 flint_abort(); /* Would be needed for Taylor method. */
28
29 acb_sin(res, z, prec);
30
31 return 0;
32 }
33
34 /* f(z) = floor(z) */
35 int
f_floor(acb_ptr res,const acb_t z,void * param,slong order,slong prec)36 f_floor(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
37 {
38 if (order > 1)
39 flint_abort(); /* Would be needed for Taylor method. */
40
41 acb_real_floor(res, z, order != 0, prec);
42
43 return 0;
44 }
45
46 /* f(z) = sqrt(1-z^2) */
47 int
f_circle(acb_ptr res,const acb_t z,void * param,slong order,slong prec)48 f_circle(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
49 {
50 if (order > 1)
51 flint_abort(); /* Would be needed for Taylor method. */
52
53 acb_one(res);
54 acb_submul(res, z, z, prec);
55 acb_real_sqrtpos(res, res, order != 0, prec);
56
57 return 0;
58 }
59
60 /* f(z) = 1/(1+z^2) */
61 int
f_atanderiv(acb_ptr res,const acb_t z,void * param,slong order,slong prec)62 f_atanderiv(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
63 {
64 if (order > 1)
65 flint_abort(); /* Would be needed for Taylor method. */
66
67 acb_mul(res, z, z, prec);
68 acb_add_ui(res, res, 1, prec);
69 acb_inv(res, res, prec);
70
71 return 0;
72 }
73
74 /* f(z) = sin(z + exp(z)) -- Rump's oscillatory example */
75 int
f_rump(acb_ptr res,const acb_t z,void * param,slong order,slong prec)76 f_rump(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
77 {
78 if (order > 1)
79 flint_abort(); /* Would be needed for Taylor method. */
80
81 acb_exp(res, z, prec);
82 acb_add(res, res, z, prec);
83 acb_sin(res, res, prec);
84
85 return 0;
86 }
87
88 /* f(z) = |z^4 + 10z^3 + 19z^2 - 6z - 6| exp(z) (for real z) --
89 Helfgott's integral on MathOverflow */
90 int
f_helfgott(acb_ptr res,const acb_t z,void * param,slong order,slong prec)91 f_helfgott(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
92 {
93 if (order > 1)
94 flint_abort(); /* Would be needed for Taylor method. */
95
96 acb_add_si(res, z, 10, prec);
97 acb_mul(res, res, z, prec);
98 acb_add_si(res, res, 19, prec);
99 acb_mul(res, res, z, prec);
100 acb_add_si(res, res, -6, prec);
101 acb_mul(res, res, z, prec);
102 acb_add_si(res, res, -6, prec);
103
104 acb_real_abs(res, res, order != 0, prec);
105
106 if (acb_is_finite(res))
107 {
108 acb_t t;
109 acb_init(t);
110 acb_exp(t, z, prec);
111 acb_mul(res, res, t, prec);
112 acb_clear(t);
113 }
114
115 return 0;
116 }
117
118 /* f(z) = zeta(z) */
119 int
f_zeta(acb_ptr res,const acb_t z,void * param,slong order,slong prec)120 f_zeta(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
121 {
122 if (order > 1)
123 flint_abort(); /* Would be needed for Taylor method. */
124
125 acb_zeta(res, z, prec);
126
127 return 0;
128 }
129
130 /* f(z) = z sin(1/z), assume on real interval */
131 int
f_essing2(acb_ptr res,const acb_t z,void * param,slong order,slong prec)132 f_essing2(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
133 {
134 if (order > 1)
135 flint_abort(); /* Would be needed for Taylor method. */
136
137 if ((order == 0) && acb_is_real(z) && arb_contains_zero(acb_realref(z)))
138 {
139 /* todo: arb_zero_pm_one, arb_unit_interval? */
140 acb_zero(res);
141 mag_one(arb_radref(acb_realref(res)));
142 }
143 else
144 {
145 acb_inv(res, z, prec);
146 acb_sin(res, res, prec);
147 }
148
149 acb_mul(res, res, z, prec);
150
151 return 0;
152 }
153
154 /* f(z) = sin(1/z), assume on real interval */
155 int
f_essing(acb_ptr res,const acb_t z,void * param,slong order,slong prec)156 f_essing(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
157 {
158 if (order > 1)
159 flint_abort(); /* Would be needed for Taylor method. */
160
161 if ((order == 0) && acb_is_real(z) && arb_contains_zero(acb_realref(z)))
162 {
163 /* todo: arb_zero_pm_one, arb_unit_interval? */
164 acb_zero(res);
165 mag_one(arb_radref(acb_realref(res)));
166 }
167 else
168 {
169 acb_inv(res, z, prec);
170 acb_sin(res, res, prec);
171 }
172
173 return 0;
174 }
175
176 /* f(z) = exp(-z) z^1000 */
177 int
f_factorial1000(acb_ptr res,const acb_t z,void * param,slong order,slong prec)178 f_factorial1000(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
179 {
180 acb_t t;
181
182 if (order > 1)
183 flint_abort(); /* Would be needed for Taylor method. */
184
185 acb_init(t);
186 acb_pow_ui(t, z, 1000, prec);
187 acb_neg(res, z);
188 acb_exp(res, res, prec);
189 acb_mul(res, res, t, prec);
190 acb_clear(t);
191
192 return 0;
193 }
194
195 /* f(z) = gamma(z) */
196 int
f_gamma(acb_ptr res,const acb_t z,void * param,slong order,slong prec)197 f_gamma(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
198 {
199 if (order > 1)
200 flint_abort(); /* Would be needed for Taylor method. */
201
202 acb_gamma(res, z, prec);
203
204 return 0;
205 }
206
207 /* f(z) = sin(z) + exp(-200-z^2) */
208 int
f_sin_plus_small(acb_ptr res,const acb_t z,void * param,slong order,slong prec)209 f_sin_plus_small(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
210 {
211 acb_t t;
212
213 if (order > 1)
214 flint_abort(); /* Would be needed for Taylor method. */
215
216 acb_init(t);
217 acb_mul(t, z, z, prec);
218 acb_add_ui(t, t, 200, prec);
219 acb_neg(t, t);
220 acb_exp(t, t, prec);
221 acb_sin(res, z, prec);
222 acb_add(res, res, t, prec);
223 acb_clear(t);
224
225 return 0;
226 }
227
228 /* f(z) = exp(z) */
229 int
f_exp(acb_ptr res,const acb_t z,void * param,slong order,slong prec)230 f_exp(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
231 {
232 if (order > 1)
233 flint_abort(); /* Would be needed for Taylor method. */
234
235 acb_exp(res, z, prec);
236
237 return 0;
238 }
239
240 /* f(z) = exp(-z^2) */
241 int
f_gaussian(acb_ptr res,const acb_t z,void * param,slong order,slong prec)242 f_gaussian(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
243 {
244 if (order > 1)
245 flint_abort(); /* Would be needed for Taylor method. */
246
247 acb_mul(res, z, z, prec);
248 acb_neg(res, res);
249 acb_exp(res, res, prec);
250
251 return 0;
252 }
253
254 int
f_monster(acb_ptr res,const acb_t z,void * param,slong order,slong prec)255 f_monster(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
256 {
257 acb_t t;
258
259 if (order > 1)
260 flint_abort(); /* Would be needed for Taylor method. */
261
262 acb_init(t);
263
264 acb_exp(t, z, prec);
265 acb_real_floor(res, t, order != 0, prec);
266
267 if (acb_is_finite(res))
268 {
269 acb_sub(res, t, res, prec);
270 acb_add(t, t, z, prec);
271 acb_sin(t, t, prec);
272 acb_mul(res, res, t, prec);
273 }
274
275 acb_clear(t);
276
277 return 0;
278 }
279
280 /* f(z) = sech(10(x-0.2))^2 + sech(100(x-0.4))^4 + sech(1000(x-0.6))^6 */
281 int
f_spike(acb_ptr res,const acb_t z,void * param,slong order,slong prec)282 f_spike(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
283 {
284 acb_t a, b, c;
285
286 if (order > 1)
287 flint_abort(); /* Would be needed for Taylor method. */
288
289 acb_init(a);
290 acb_init(b);
291 acb_init(c);
292
293 acb_mul_ui(a, z, 10, prec);
294 acb_sub_ui(a, a, 2, prec);
295 acb_sech(a, a, prec);
296 acb_pow_ui(a, a, 2, prec);
297
298 acb_mul_ui(b, z, 100, prec);
299 acb_sub_ui(b, b, 40, prec);
300 acb_sech(b, b, prec);
301 acb_pow_ui(b, b, 4, prec);
302
303 acb_mul_ui(c, z, 1000, prec);
304 acb_sub_ui(c, c, 600, prec);
305 acb_sech(c, c, prec);
306 acb_pow_ui(c, c, 6, prec);
307
308 acb_add(res, a, b, prec);
309 acb_add(res, res, c, prec);
310
311 acb_clear(a);
312 acb_clear(b);
313 acb_clear(c);
314
315 return 0;
316 }
317
318 /* f(z) = sech(z) */
319 int
f_sech(acb_ptr res,const acb_t z,void * param,slong order,slong prec)320 f_sech(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
321 {
322 if (order > 1)
323 flint_abort(); /* Would be needed for Taylor method. */
324
325 acb_sech(res, z, prec);
326
327 return 0;
328 }
329
330 /* f(z) = sech^3(z) */
331 int
f_sech3(acb_ptr res,const acb_t z,void * param,slong order,slong prec)332 f_sech3(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
333 {
334 if (order > 1)
335 flint_abort(); /* Would be needed for Taylor method. */
336
337 acb_sech(res, z, prec);
338 acb_cube(res, res, prec);
339
340 return 0;
341 }
342
343 /* f(z) = -log(z) / (1 + z) */
344 int
f_log_div1p(acb_ptr res,const acb_t z,void * param,slong order,slong prec)345 f_log_div1p(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
346 {
347 acb_t t;
348
349 if (order > 1)
350 flint_abort(); /* Would be needed for Taylor method. */
351
352 acb_init(t);
353 acb_add_ui(t, z, 1, prec);
354 acb_log(res, z, prec);
355 acb_div(res, res, t, prec);
356 acb_neg(res, res);
357
358 acb_clear(t);
359
360 return 0;
361 }
362
363 /* f(z) = z exp(-z) / (1 + exp(-z)) */
364 int
f_log_div1p_transformed(acb_ptr res,const acb_t z,void * param,slong order,slong prec)365 f_log_div1p_transformed(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
366 {
367 acb_t t;
368
369 if (order > 1)
370 flint_abort(); /* Would be needed for Taylor method. */
371
372 acb_init(t);
373 acb_neg(t, z);
374 acb_exp(t, t, prec);
375 acb_add_ui(res, t, 1, prec);
376 acb_div(res, t, res, prec);
377 acb_mul(res, res, z, prec);
378
379 acb_clear(t);
380
381 return 0;
382 }
383
384 int
f_elliptic_p_laurent_n(acb_ptr res,const acb_t z,void * param,slong order,slong prec)385 f_elliptic_p_laurent_n(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
386 {
387 slong n;
388 acb_t tau;
389
390 if (order > 1)
391 flint_abort(); /* Would be needed for Taylor method. */
392
393 n = ((slong *)(param))[0];
394
395 acb_init(tau);
396 acb_onei(tau);
397
398 acb_modular_elliptic_p(res, z, tau, prec);
399
400 acb_pow_si(tau, z, -n - 1, prec);
401 acb_mul(res, res, tau, prec);
402
403 acb_clear(tau);
404
405 return 0;
406 }
407
408 /* f(z) = zeta'(z) / zeta(z) */
409 int
f_zeta_frac(acb_ptr res,const acb_t z,void * param,slong order,slong prec)410 f_zeta_frac(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
411 {
412 acb_struct t[2];
413
414 if (order > 1)
415 flint_abort(); /* Would be needed for Taylor method. */
416
417 acb_init(t);
418 acb_init(t + 1);
419
420 acb_dirichlet_zeta_jet(t, z, 0, 2, prec);
421 acb_div(res, t + 1, t, prec);
422
423 acb_clear(t);
424 acb_clear(t + 1);
425
426 return 0;
427 }
428
429 int
f_lambertw(acb_ptr res,const acb_t z,void * param,slong order,slong prec)430 f_lambertw(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
431 {
432 acb_t t;
433
434 if (order > 1)
435 flint_abort(); /* Would be needed for Taylor method. */
436
437 acb_init(t);
438
439 prec = FLINT_MIN(prec, acb_rel_accuracy_bits(z) + 10);
440
441 if (order != 0)
442 {
443 /* check for branch cut */
444 arb_const_e(acb_realref(t), prec);
445 acb_inv(t, t, prec);
446 acb_add(t, t, z, prec);
447
448 if (arb_contains_zero(acb_imagref(t)) &&
449 arb_contains_nonpositive(acb_realref(t)))
450 {
451 acb_indeterminate(t);
452 }
453 }
454
455 if (acb_is_finite(t))
456 {
457 fmpz_t k;
458 fmpz_init(k);
459 acb_lambertw(res, z, k, 0, prec);
460 fmpz_clear(k);
461 }
462 else
463 {
464 acb_indeterminate(res);
465 }
466
467 acb_clear(t);
468
469 return 0;
470 }
471
472 /* f(z) = max(sin(z), cos(z)) */
473 int
f_max_sin_cos(acb_ptr res,const acb_t z,void * param,slong order,slong prec)474 f_max_sin_cos(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
475 {
476 acb_t s, c;
477
478 if (order > 1)
479 flint_abort(); /* Would be needed for Taylor method. */
480
481 acb_init(s);
482 acb_init(c);
483
484 acb_sin_cos(s, c, z, prec);
485 acb_real_max(res, s, c, order != 0, prec);
486
487 acb_clear(s);
488 acb_clear(c);
489
490 return 0;
491 }
492
493 /* f(z) = erf(z/sqrt(0.0002)*0.5 +1.5)*exp(-z), example provided by Silviu-Ioan Filip */
494 int
f_erf_bent(acb_ptr res,const acb_t z,void * param,slong order,slong prec)495 f_erf_bent(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
496 {
497 acb_t t;
498
499 if (order > 1)
500 flint_abort(); /* Would be needed for Taylor method. */
501
502 acb_init(t);
503
504 acb_set_ui(t, 1250);
505 acb_sqrt(t, t, prec);
506 acb_mul(t, t, z, prec);
507 acb_set_d(res, 1.5);
508 acb_add(res, res, t, prec);
509 acb_hypgeom_erf(res, res, prec);
510
511 acb_neg(t, z);
512 acb_exp(t, t, prec);
513 acb_mul(res, res, t, prec);
514
515 acb_clear(t);
516
517 return 0;
518 }
519
520 /* f(z) = Ai(z) */
521 int
f_airy_ai(acb_ptr res,const acb_t z,void * param,slong order,slong prec)522 f_airy_ai(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
523 {
524 if (order > 1)
525 flint_abort(); /* Would be needed for Taylor method. */
526
527
528 acb_hypgeom_airy(res, NULL, NULL, NULL, z, prec);
529
530 return 0;
531 }
532
533 int
f_horror(acb_ptr res,const acb_t z,void * param,slong order,slong prec)534 f_horror(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
535 {
536 acb_t s, t;
537
538 if (order > 1)
539 flint_abort(); /* Would be needed for Taylor method. */
540
541 acb_init(s);
542 acb_init(t);
543
544 acb_real_floor(res, z, order != 0, prec);
545
546 if (acb_is_finite(res))
547 {
548 acb_sub(res, z, res, prec);
549 acb_set_d(t, 0.5);
550 acb_sub(res, res, t, prec);
551 acb_sin_cos(s, t, z, prec);
552 acb_real_max(s, s, t, order != 0, prec);
553 acb_mul(res, res, s, prec);
554 }
555
556 acb_clear(s);
557 acb_clear(t);
558
559 return 0;
560 }
561
562 /* f(z) = sqrt(z) */
563 int
f_sqrt(acb_ptr res,const acb_t z,void * param,slong order,slong prec)564 f_sqrt(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
565 {
566 if (order > 1)
567 flint_abort(); /* Would be needed for Taylor method. */
568
569 acb_sqrt_analytic(res, z, order != 0, prec);
570
571 return 0;
572 }
573
574 /* f(z) = 1/sqrt(z) */
575 int
f_rsqrt(acb_ptr res,const acb_t z,void * param,slong order,slong prec)576 f_rsqrt(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
577 {
578 if (order > 1)
579 flint_abort(); /* Would be needed for Taylor method. */
580
581 acb_rsqrt_analytic(res, z, order != 0, prec);
582
583 return 0;
584 }
585
586 /* f(z) = rgamma(z) */
587 int
f_rgamma(acb_ptr res,const acb_t z,void * param,slong order,slong prec)588 f_rgamma(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
589 {
590 if (order > 1)
591 flint_abort(); /* Would be needed for Taylor method. */
592
593 acb_rgamma(res, z, prec);
594
595 return 0;
596 }
597
598 /* f(z) = exp(-z^2+iz) */
599 int
f_gaussian_twist(acb_ptr res,const acb_t z,void * param,slong order,slong prec)600 f_gaussian_twist(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
601 {
602 if (order > 1)
603 flint_abort(); /* Would be needed for Taylor method. */
604
605 acb_mul_onei(res, z);
606 acb_submul(res, z, z, prec);
607 acb_exp(res, res, prec);
608
609 return 0;
610 }
611
612 /* f(z) = exp(-z) Ai(-z) */
613 int
f_exp_airy(acb_ptr res,const acb_t z,void * param,slong order,slong prec)614 f_exp_airy(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
615 {
616 acb_t t;
617
618 if (order > 1)
619 flint_abort(); /* Would be needed for Taylor method. */
620
621 acb_init(t);
622 acb_neg(t, z);
623 acb_hypgeom_airy(res, NULL, NULL, NULL, t, prec);
624 acb_exp(t, t, prec);
625 acb_mul(res, res, t, prec);
626 acb_clear(t);
627
628 return 0;
629 }
630
631 /* f(z) = z sin(z) / (1 + cos(z)^2) */
632 int
f_sin_cos_frac(acb_ptr res,const acb_t z,void * param,slong order,slong prec)633 f_sin_cos_frac(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
634 {
635 acb_t s, c;
636
637 if (order > 1)
638 flint_abort(); /* Would be needed for Taylor method. */
639
640 acb_init(s);
641 acb_init(c);
642
643 acb_sin_cos(s, c, z, prec);
644 acb_mul(c, c, c, prec);
645 acb_add_ui(c, c, 1, prec);
646 acb_mul(s, s, z, prec);
647 acb_div(res, s, c, prec);
648
649 acb_clear(s);
650 acb_clear(c);
651
652 return 0;
653 }
654
655 /* f(z) = sin((1/1000 + (1-z)^2)^(-3/2)), example from Mioara Joldes' thesis
656 (suggested by Nicolas Brisebarre) */
657 int
f_sin_near_essing(acb_ptr res,const acb_t z,void * param,slong order,slong prec)658 f_sin_near_essing(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
659 {
660 acb_t t, u;
661
662 if (order > 1)
663 flint_abort(); /* Would be needed for Taylor method. */
664
665 acb_init(t);
666 acb_init(u);
667
668 acb_sub_ui(t, z, 1, prec);
669 acb_neg(t, t);
670 acb_mul(t, t, t, prec);
671 acb_one(u);
672 acb_div_ui(u, u, 1000, prec);
673 acb_add(t, t, u, prec);
674 acb_set_d(u, -1.5);
675 acb_pow_analytic(t, t, u, order != 0, prec);
676 acb_sin(res, t, prec);
677
678 acb_clear(t);
679 acb_clear(u);
680
681 return 0;
682 }
683
684 /* f(z) = exp(-z) (I_0(z/k))^k, from Bruno Salvy */
685 int
f_scaled_bessel(acb_ptr res,const acb_t z,void * param,slong order,slong prec)686 f_scaled_bessel(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
687 {
688 acb_t nu;
689
690 ulong k = ((ulong *) param)[0];
691
692 acb_init(nu);
693 acb_div_ui(res, z, k, prec);
694 acb_hypgeom_bessel_i_scaled(res, nu, res, prec);
695 acb_pow_ui(res, res, k, prec);
696 acb_clear(nu);
697
698 return 0;
699 }
700
701 /*
702 Bound for scaled Bessel function: 2/(2 pi x)^(1/2)
703 Bound for tail of integral: 2 N (k / (pi N))^(k / 2) / (k - 2).
704 */
705 void
scaled_bessel_tail_bound(arb_t b,ulong k,const arb_t N,slong prec)706 scaled_bessel_tail_bound(arb_t b, ulong k, const arb_t N, slong prec)
707 {
708 arb_const_pi(b, prec);
709 arb_mul(b, b, N, prec);
710 arb_ui_div(b, k, b, prec);
711 arb_sqrt(b, b, prec);
712 arb_pow_ui(b, b, k, prec);
713 arb_mul(b, b, N, prec);
714 arb_mul_ui(b, b, 2, prec);
715 arb_div_ui(b, b, k - 2, prec);
716 }
717
718 void
scaled_bessel_select_N(arb_t N,ulong k,slong prec)719 scaled_bessel_select_N(arb_t N, ulong k, slong prec)
720 {
721 slong e;
722 double f = log(k/3.14159265358979)/log(2);
723
724 e = 1;
725 while ((k / 2.0) * (e - f) - e < prec + 5)
726 e++;
727
728 arb_one(N);
729 arb_mul_2exp_si(N, N, e);
730 }
731
732 /* ------------------------------------------------------------------------- */
733 /* Main test program */
734 /* ------------------------------------------------------------------------- */
735
736 #define NUM_INTEGRALS 37
737
738 const char * descr[NUM_INTEGRALS] =
739 {
740 "int_0^100 sin(x) dx",
741 "4 int_0^1 1/(1+x^2) dx",
742 "2 int_0^{inf} 1/(1+x^2) dx (using domain truncation)",
743 "4 int_0^1 sqrt(1-x^2) dx",
744 "int_0^8 sin(x+exp(x)) dx",
745 "int_1^101 floor(x) dx",
746 "int_0^1 |x^4+10x^3+19x^2-6x-6| exp(x) dx",
747 "1/(2 pi i) int zeta(s) ds (closed path around s = 1)",
748 "int_0^1 sin(1/x) dx (slow convergence, use -heap and/or -tol)",
749 "int_0^1 x sin(1/x) dx (slow convergence, use -heap and/or -tol)",
750 "int_0^10000 x^1000 exp(-x) dx",
751 "int_1^{1+1000i} gamma(x) dx",
752 "int_{-10}^{10} sin(x) + exp(-200-x^2) dx",
753 "int_{-1020}^{-1010} exp(x) dx (use -tol 0 for relative error)",
754 "int_0^{inf} exp(-x^2) dx (using domain truncation)",
755 "int_0^1 sech(10(x-0.2))^2 + sech(100(x-0.4))^4 + sech(1000(x-0.6))^6 dx",
756 "int_0^8 (exp(x)-floor(exp(x))) sin(x+exp(x)) dx (use higher -eval)",
757 "int_0^{inf} sech(x) dx (using domain truncation)",
758 "int_0^{inf} sech^3(x) dx (using domain truncation)",
759 "int_0^1 -log(x)/(1+x) dx (using domain truncation)",
760 "int_0^{inf} x exp(-x)/(1+exp(-x)) dx (using domain truncation)",
761 "int_C wp(x)/x^(11) dx (contour for 10th Laurent coefficient of Weierstrass p-function)",
762 "N(1000) = count zeros with 0 < t <= 1000 of zeta(s) using argument principle",
763 "int_0^{1000} W_0(x) dx",
764 "int_0^pi max(sin(x), cos(x)) dx",
765 "int_{-1}^1 erf(x/sqrt(0.0002)*0.5+1.5)*exp(-x) dx",
766 "int_{-10}^10 Ai(x) dx",
767 "int_0^10 (x-floor(x)-1/2) max(sin(x),cos(x)) dx",
768 "int_{-1-i}^{-1+i} sqrt(x) dx",
769 "int_0^{inf} exp(-x^2+ix) dx (using domain truncation)",
770 "int_0^{inf} exp(-x) Ai(-x) dx (using domain truncation)",
771 "int_0^pi x sin(x) / (1 + cos(x)^2) dx",
772 "int_0^3 sin(0.001 + (1-x)^2)^(-3/2)) dx (slow convergence, use higher -eval)",
773 "int_0^{inf} exp(-x) I_0(x/3)^3 dx (using domain truncation)",
774 "int_0^{inf} exp(-x) I_0(x/15)^{15} dx (using domain truncation)",
775 "int_{-1-i}^{-1+i} 1/sqrt(x) dx",
776 "int_0^{inf} 1/gamma(x) dx (using domain truncation)",
777 };
778
main(int argc,char * argv[])779 int main(int argc, char *argv[])
780 {
781 acb_t s, t, a, b;
782 mag_t tol;
783 slong prec, goal;
784 slong N;
785 ulong k;
786 int integral, ifrom, ito;
787 int i, twice, havegoal, havetol;
788 acb_calc_integrate_opt_t options;
789
790 ifrom = ito = -1;
791
792 for (i = 1; i < argc; i++)
793 {
794 if (!strcmp(argv[i], "-i"))
795 {
796 if (!strcmp(argv[i+1], "all"))
797 {
798 ifrom = 0;
799 ito = NUM_INTEGRALS - 1;
800 }
801 else
802 {
803 ifrom = ito = atol(argv[i+1]);
804 if (ito < 0 || ito >= NUM_INTEGRALS)
805 flint_abort();
806 }
807 }
808 }
809
810 if (ifrom == -1)
811 {
812 flint_printf("Compute integrals using acb_calc_integrate.\n");
813 flint_printf("Usage: integrals -i n [-prec p] [-tol eps] [-twice] [...]\n\n");
814 flint_printf("-i n - compute integral n (0 <= n <= %d), or \"-i all\"\n", NUM_INTEGRALS - 1);
815 flint_printf("-prec p - precision in bits (default p = 64)\n");
816 flint_printf("-goal p - approximate relative accuracy goal (default p)\n");
817 flint_printf("-tol eps - approximate absolute error goal (default 2^-p)\n");
818 flint_printf("-twice - run twice (to see overhead of computing nodes)\n");
819 flint_printf("-heap - use heap for subinterval queue\n");
820 flint_printf("-verbose - show information\n");
821 flint_printf("-verbose2 - show more information\n");
822 flint_printf("-deg n - use quadrature degree up to n\n");
823 flint_printf("-eval n - limit number of function evaluations to n\n");
824 flint_printf("-depth n - limit subinterval queue size to n\n\n");
825 flint_printf("Implemented integrals:\n");
826 for (integral = 0; integral < NUM_INTEGRALS; integral++)
827 flint_printf("I%d = %s\n", integral, descr[integral]);
828 flint_printf("\n");
829 return 1;
830 }
831
832 acb_calc_integrate_opt_init(options);
833
834 prec = 64;
835 twice = 0;
836 goal = 0;
837 havetol = havegoal = 0;
838
839 acb_init(a);
840 acb_init(b);
841 acb_init(s);
842 acb_init(t);
843 mag_init(tol);
844
845 for (i = 1; i < argc; i++)
846 {
847 if (!strcmp(argv[i], "-prec"))
848 {
849 prec = atol(argv[i+1]);
850 }
851 else if (!strcmp(argv[i], "-twice"))
852 {
853 twice = 1;
854 }
855 else if (!strcmp(argv[i], "-goal"))
856 {
857 goal = atol(argv[i+1]);
858 if (goal < 0)
859 {
860 flint_printf("expected goal >= 0\n");
861 return 1;
862 }
863 havegoal = 1;
864 }
865 else if (!strcmp(argv[i], "-tol"))
866 {
867 arb_t x;
868 arb_init(x);
869 arb_set_str(x, argv[i+1], 10);
870 arb_get_mag(tol, x);
871 arb_clear(x);
872 havetol = 1;
873 }
874 else if (!strcmp(argv[i], "-deg"))
875 {
876 options->deg_limit = atol(argv[i+1]);
877 }
878 else if (!strcmp(argv[i], "-eval"))
879 {
880 options->eval_limit = atol(argv[i+1]);
881 }
882 else if (!strcmp(argv[i], "-depth"))
883 {
884 options->depth_limit = atol(argv[i+1]);
885 }
886 else if (!strcmp(argv[i], "-verbose"))
887 {
888 options->verbose = 1;
889 }
890 else if (!strcmp(argv[i], "-verbose2"))
891 {
892 options->verbose = 2;
893 }
894 else if (!strcmp(argv[i], "-heap"))
895 {
896 options->use_heap = 1;
897 }
898 }
899
900 if (!havegoal)
901 goal = prec;
902
903 if (!havetol)
904 mag_set_ui_2exp_si(tol, 1, -prec);
905
906 for (integral = ifrom; integral <= ito; integral++)
907 {
908 flint_printf("I%d = %s ...\n", integral, descr[integral]);
909
910 for (i = 0; i < 1 + twice; i++)
911 {
912 TIMEIT_ONCE_START
913 switch (integral)
914 {
915 case 0:
916 acb_set_d(a, 0);
917 acb_set_d(b, 100);
918 acb_calc_integrate(s, f_sin, NULL, a, b, goal, tol, options, prec);
919 break;
920
921 case 1:
922 acb_set_d(a, 0);
923 acb_set_d(b, 1);
924 acb_calc_integrate(s, f_atanderiv, NULL, a, b, goal, tol, options, prec);
925 acb_mul_2exp_si(s, s, 2);
926 break;
927
928 case 2:
929 acb_set_d(a, 0);
930 acb_one(b);
931 acb_mul_2exp_si(b, b, goal);
932 acb_calc_integrate(s, f_atanderiv, NULL, a, b, goal, tol, options, prec);
933 arb_add_error_2exp_si(acb_realref(s), -goal);
934 acb_mul_2exp_si(s, s, 1);
935 break;
936
937 case 3:
938 acb_set_d(a, 0);
939 acb_set_d(b, 1);
940 acb_calc_integrate(s, f_circle, NULL, a, b, goal, tol, options, prec);
941 acb_mul_2exp_si(s, s, 2);
942 break;
943
944 case 4:
945 acb_set_d(a, 0);
946 acb_set_d(b, 8);
947 acb_calc_integrate(s, f_rump, NULL, a, b, goal, tol, options, prec);
948 break;
949
950 case 5:
951 acb_set_d(a, 1);
952 acb_set_d(b, 101);
953 acb_calc_integrate(s, f_floor, NULL, a, b, goal, tol, options, prec);
954 break;
955
956 case 6:
957 acb_set_d(a, 0);
958 acb_set_d(b, 1);
959 acb_calc_integrate(s, f_helfgott, NULL, a, b, goal, tol, options, prec);
960 break;
961
962 case 7:
963 acb_zero(s);
964
965 acb_set_d_d(a, -1.0, -1.0);
966 acb_set_d_d(b, 2.0, -1.0);
967 acb_calc_integrate(t, f_zeta, NULL, a, b, goal, tol, options, prec);
968 acb_add(s, s, t, prec);
969
970 acb_set_d_d(a, 2.0, -1.0);
971 acb_set_d_d(b, 2.0, 1.0);
972 acb_calc_integrate(t, f_zeta, NULL, a, b, goal, tol, options, prec);
973 acb_add(s, s, t, prec);
974
975 acb_set_d_d(a, 2.0, 1.0);
976 acb_set_d_d(b, -1.0, 1.0);
977 acb_calc_integrate(t, f_zeta, NULL, a, b, goal, tol, options, prec);
978 acb_add(s, s, t, prec);
979
980 acb_set_d_d(a, -1.0, 1.0);
981 acb_set_d_d(b, -1.0, -1.0);
982 acb_calc_integrate(t, f_zeta, NULL, a, b, goal, tol, options, prec);
983 acb_add(s, s, t, prec);
984
985 acb_const_pi(t, prec);
986 acb_div(s, s, t, prec);
987 acb_mul_2exp_si(s, s, -1);
988 acb_div_onei(s, s);
989 break;
990
991 case 8:
992 acb_set_d(a, 0);
993 acb_set_d(b, 1);
994 acb_calc_integrate(s, f_essing, NULL, a, b, goal, tol, options, prec);
995 break;
996
997 case 9:
998 acb_set_d(a, 0);
999 acb_set_d(b, 1);
1000 acb_calc_integrate(s, f_essing2, NULL, a, b, goal, tol, options, prec);
1001 break;
1002
1003 case 10:
1004 acb_set_d(a, 0);
1005 acb_set_d(b, 10000);
1006 acb_calc_integrate(s, f_factorial1000, NULL, a, b, goal, tol, options, prec);
1007 break;
1008
1009 case 11:
1010 acb_set_d_d(a, 1.0, 0.0);
1011 acb_set_d_d(b, 1.0, 1000.0);
1012 acb_calc_integrate(s, f_gamma, NULL, a, b, goal, tol, options, prec);
1013 break;
1014
1015 case 12:
1016 acb_set_d(a, -10.0);
1017 acb_set_d(b, 10.0);
1018 acb_calc_integrate(s, f_sin_plus_small, NULL, a, b, goal, tol, options, prec);
1019 break;
1020
1021 case 13:
1022 acb_set_d(a, -1020.0);
1023 acb_set_d(b, -1010.0);
1024 acb_calc_integrate(s, f_exp, NULL, a, b, goal, tol, options, prec);
1025 break;
1026
1027 case 14:
1028 acb_set_d(a, 0);
1029 acb_set_d(b, ceil(sqrt(goal * 0.693147181) + 1.0));
1030 acb_calc_integrate(s, f_gaussian, NULL, a, b, goal, tol, options, prec);
1031 acb_mul(b, b, b, prec);
1032 acb_neg(b, b);
1033 acb_exp(b, b, prec);
1034 arb_add_error(acb_realref(s), acb_realref(b));
1035 break;
1036
1037 case 15:
1038 acb_set_d(a, 0.0);
1039 acb_set_d(b, 1.0);
1040 acb_calc_integrate(s, f_spike, NULL, a, b, goal, tol, options, prec);
1041 break;
1042
1043 case 16:
1044 acb_set_d(a, 0.0);
1045 acb_set_d(b, 8.0);
1046 acb_calc_integrate(s, f_monster, NULL, a, b, goal, tol, options, prec);
1047 break;
1048
1049 case 17:
1050 acb_set_d(a, 0);
1051 acb_set_d(b, ceil(goal * 0.693147181 + 1.0));
1052 acb_calc_integrate(s, f_sech, NULL, a, b, goal, tol, options, prec);
1053 acb_neg(b, b);
1054 acb_exp(b, b, prec);
1055 acb_mul_2exp_si(b, b, 1);
1056 arb_add_error(acb_realref(s), acb_realref(b));
1057 break;
1058
1059 case 18:
1060 acb_set_d(a, 0);
1061 acb_set_d(b, ceil(goal * 0.693147181 / 3.0 + 2.0));
1062 acb_calc_integrate(s, f_sech3, NULL, a, b, goal, tol, options, prec);
1063 acb_neg(b, b);
1064 acb_mul_ui(b, b, 3, prec);
1065 acb_exp(b, b, prec);
1066 acb_mul_2exp_si(b, b, 3);
1067 acb_div_ui(b, b, 3, prec);
1068 arb_add_error(acb_realref(s), acb_realref(b));
1069 break;
1070
1071 case 19:
1072 if (goal < 0)
1073 abort();
1074 /* error bound 2^-N (1+N) when truncated at 2^-N */
1075 N = goal + FLINT_BIT_COUNT(goal);
1076 acb_one(a);
1077 acb_mul_2exp_si(a, a, -N);
1078 acb_one(b);
1079 acb_calc_integrate(s, f_log_div1p, NULL, a, b, goal, tol, options, prec);
1080 acb_set_ui(b, N + 1);
1081 acb_mul_2exp_si(b, b, -N);
1082 arb_add_error(acb_realref(s), acb_realref(b));
1083 break;
1084
1085 case 20:
1086 if (goal < 0)
1087 abort();
1088 /* error bound (N+1) exp(-N) when truncated at N */
1089 N = goal + FLINT_BIT_COUNT(goal);
1090 acb_zero(a);
1091 acb_set_ui(b, N);
1092 acb_calc_integrate(s, f_log_div1p_transformed, NULL, a, b, goal, tol, options, prec);
1093 acb_neg(b, b);
1094 acb_exp(b, b, prec);
1095 acb_mul_ui(b, b, N + 1, prec);
1096 arb_add_error(acb_realref(s), acb_realref(b));
1097 break;
1098
1099 case 21:
1100
1101 acb_zero(s);
1102
1103 N = 10;
1104
1105 acb_set_d_d(a, 0.5, -0.5);
1106 acb_set_d_d(b, 0.5, 0.5);
1107 acb_calc_integrate(t, f_elliptic_p_laurent_n, &N, a, b, goal, tol, options, prec);
1108 acb_add(s, s, t, prec);
1109
1110 acb_set_d_d(a, 0.5, 0.5);
1111 acb_set_d_d(b, -0.5, 0.5);
1112 acb_calc_integrate(t, f_elliptic_p_laurent_n, &N, a, b, goal, tol, options, prec);
1113 acb_add(s, s, t, prec);
1114
1115 acb_set_d_d(a, -0.5, 0.5);
1116 acb_set_d_d(b, -0.5, -0.5);
1117 acb_calc_integrate(t, f_elliptic_p_laurent_n, &N, a, b, goal, tol, options, prec);
1118 acb_add(s, s, t, prec);
1119
1120 acb_set_d_d(a, -0.5, -0.5);
1121 acb_set_d_d(b, 0.5, -0.5);
1122 acb_calc_integrate(t, f_elliptic_p_laurent_n, &N, a, b, goal, tol, options, prec);
1123 acb_add(s, s, t, prec);
1124
1125 acb_const_pi(t, prec);
1126 acb_div(s, s, t, prec);
1127 acb_mul_2exp_si(s, s, -1);
1128 acb_div_onei(s, s);
1129 break;
1130
1131 case 22:
1132
1133 acb_zero(s);
1134
1135 N = 1000;
1136
1137 acb_set_d_d(a, 100.0, 0.0);
1138 acb_set_d_d(b, 100.0, N);
1139 acb_calc_integrate(t, f_zeta_frac, NULL, a, b, goal, tol, options, prec);
1140 acb_add(s, s, t, prec);
1141
1142 acb_set_d_d(a, 100, N);
1143 acb_set_d_d(b, 0.5, N);
1144 acb_calc_integrate(t, f_zeta_frac, NULL, a, b, goal, tol, options, prec);
1145 acb_add(s, s, t, prec);
1146
1147 acb_div_onei(s, s);
1148 arb_zero(acb_imagref(s));
1149
1150 acb_set_ui(t, N);
1151 acb_dirichlet_hardy_theta(t, t, NULL, NULL, 1, prec);
1152 acb_add(s, s, t, prec);
1153
1154 acb_const_pi(t, prec);
1155 acb_div(s, s, t, prec);
1156 acb_add_ui(s, s, 1, prec);
1157 break;
1158
1159 case 23:
1160 acb_set_d(a, 0.0);
1161 acb_set_d(b, 1000.0);
1162 acb_calc_integrate(s, f_lambertw, NULL, a, b, goal, tol, options, prec);
1163 break;
1164
1165 case 24:
1166 acb_set_d(a, 0.0);
1167 acb_const_pi(b, prec);
1168 acb_calc_integrate(s, f_max_sin_cos, NULL, a, b, goal, tol, options, prec);
1169 break;
1170
1171 case 25:
1172 acb_set_si(a, -1);
1173 acb_set_si(b, 1);
1174 acb_calc_integrate(s, f_erf_bent, NULL, a, b, goal, tol, options, prec);
1175 break;
1176
1177 case 26:
1178 acb_set_si(a, -10);
1179 acb_set_si(b, 10);
1180 acb_calc_integrate(s, f_airy_ai, NULL, a, b, goal, tol, options, prec);
1181 break;
1182
1183 case 27:
1184 acb_set_si(a, 0);
1185 acb_set_si(b, 10);
1186 acb_calc_integrate(s, f_horror, NULL, a, b, goal, tol, options, prec);
1187 break;
1188
1189 case 28:
1190 acb_set_d_d(a, -1, -1);
1191 acb_set_d_d(b, -1, 1);
1192 acb_calc_integrate(s, f_sqrt, NULL, a, b, goal, tol, options, prec);
1193 break;
1194
1195 case 29:
1196 acb_set_d(a, 0);
1197 acb_set_d(b, ceil(sqrt(goal * 0.693147181) + 1.0));
1198 acb_calc_integrate(s, f_gaussian_twist, NULL, a, b, goal, tol, options, prec);
1199 acb_mul(b, b, b, prec);
1200 acb_neg(b, b);
1201 acb_exp(b, b, prec);
1202 arb_add_error(acb_realref(s), acb_realref(b));
1203 arb_add_error(acb_imagref(s), acb_realref(b));
1204 break;
1205
1206 case 30:
1207 acb_set_d(a, 0);
1208 acb_set_d(b, ceil(goal * 0.693147181 + 1.0));
1209 acb_calc_integrate(s, f_exp_airy, NULL, a, b, goal, tol, options, prec);
1210 acb_neg(b, b);
1211 acb_exp(b, b, prec);
1212 acb_mul_2exp_si(b, b, 1);
1213 arb_add_error(acb_realref(s), acb_realref(b));
1214 break;
1215
1216 case 31:
1217 acb_zero(a);
1218 acb_const_pi(b, prec);
1219 acb_calc_integrate(s, f_sin_cos_frac, NULL, a, b, goal, tol, options, prec);
1220 break;
1221
1222 case 32:
1223 acb_zero(a);
1224 acb_set_ui(b, 3);
1225 acb_calc_integrate(s, f_sin_near_essing, NULL, a, b, goal, tol, options, prec);
1226 break;
1227
1228 case 33:
1229 acb_zero(a);
1230 acb_zero(b);
1231 k = 3;
1232 scaled_bessel_select_N(acb_realref(b), k, prec);
1233 acb_calc_integrate(s, f_scaled_bessel, &k, a, b, goal, tol, options, prec);
1234 scaled_bessel_tail_bound(acb_realref(a), k, acb_realref(b), prec);
1235 arb_add_error(acb_realref(s), acb_realref(a));
1236 break;
1237
1238 case 34:
1239 acb_zero(a);
1240 acb_zero(b);
1241 k = 15;
1242 scaled_bessel_select_N(acb_realref(b), k, prec);
1243 acb_calc_integrate(s, f_scaled_bessel, &k, a, b, goal, tol, options, prec);
1244 scaled_bessel_tail_bound(acb_realref(a), k, acb_realref(b), prec);
1245 arb_add_error(acb_realref(s), acb_realref(a));
1246 break;
1247
1248 case 35:
1249 acb_set_d_d(a, -1, -1);
1250 acb_set_d_d(b, -1, 1);
1251 acb_calc_integrate(s, f_rsqrt, NULL, a, b, goal, tol, options, prec);
1252
1253 case 36:
1254 if (goal < 0)
1255 abort();
1256 acb_zero(a);
1257 acb_set_ui(b, 4 + (goal + 1) / 2);
1258 acb_calc_integrate(s, f_rgamma, NULL, a, b, goal, tol, options, prec);
1259 arb_add_error_2exp_si(acb_realref(s), -goal);
1260 break;
1261
1262 default:
1263 abort();
1264 }
1265
1266 TIMEIT_ONCE_STOP
1267 }
1268 flint_printf("I%d = ", integral);
1269 acb_printn(s, 3.333 * prec, 0);
1270 flint_printf("\n\n");
1271 }
1272
1273 acb_clear(a);
1274 acb_clear(b);
1275 acb_clear(s);
1276 acb_clear(t);
1277 mag_clear(tol);
1278
1279 flint_cleanup();
1280 return 0;
1281 }
1282
1283