1.. _examples: 2 3Example programs 4=============================================================================== 5 6.. highlight:: text 7 8The *examples* directory 9(https://github.com/fredrik-johansson/arb/tree/master/examples) 10contains several complete C programs, which are documented below. Running:: 11 12 make examples 13 14will compile the programs and place the binaries in ``build/examples``. 15 16pi.c 17------------------------------------------------------------------------------- 18 19This program computes `\pi` to an accuracy of roughly *n* decimal digits 20by calling the :func:`arb_const_pi` function with a 21working precision of roughly `n \log_2(10)` bits. 22 23Sample output, computing `\pi` to one million digits:: 24 25 > build/examples/pi 1000000 26 computing pi with a precision of 3321933 bits... cpu/wall(s): 0.58 0.586 27 virt/peak/res/peak(MB): 28.24 36.84 8.86 15.56 28 [3.14159265358979323846{...999959 digits...}42209010610577945815 +/- 3e-1000000] 29 30The program prints an interval guaranteed to contain `\pi`, and where 31all displayed digits are correct up to an error of plus or minus 32one unit in the last place (see :func:`arb_printn`). 33By default, only the first and last few digits are printed. 34Pass 0 as a second argument to print all digits (or pass *m* to 35print *m* + 1 leading and *m* trailing digits, as above with 36the default *m* = 20). 37 38hilbert_matrix.c 39------------------------------------------------------------------------------- 40 41Given an input integer *n*, this program accurately computes the 42determinant of the *n* by *n* Hilbert matrix. 43Hilbert matrices are notoriously ill-conditioned: although the 44entries are close to unit magnitude, the determinant `h_n` 45decreases superexponentially (nearly as `1/4^{n^2}`) as 46a function of *n*. 47This program automatically doubles the working precision 48until the ball computed for `h_n` by :func:`arb_mat_det` 49does not contain zero. 50 51Sample output:: 52 53 $ build/examples/hilbert_matrix 200 54 prec=20: [+/- 1.32e-335] 55 prec=40: [+/- 1.63e-545] 56 prec=80: [+/- 1.30e-933] 57 prec=160: [+/- 3.62e-1926] 58 prec=320: [+/- 1.81e-4129] 59 prec=640: [+/- 3.84e-8838] 60 prec=1280: [2.955454297e-23924 +/- 8.29e-23935] 61 success! 62 cpu/wall(s): 8.494 8.513 63 virt/peak/res/peak(MB): 134.98 134.98 111.57 111.57 64 65Called with ``-eig n``, instead of computing the determinant, 66the program computes the smallest eigenvalue of the Hilbert matrix 67(in fact, it isolates all eigenvalues and prints the smallest eigenvalue):: 68 69 $ build/examples/hilbert_matrix -eig 50 70 prec=20: nan 71 prec=40: nan 72 prec=80: nan 73 prec=160: nan 74 prec=320: nan 75 prec=640: [1.459157797e-74 +/- 2.49e-84] 76 success! 77 cpu/wall(s): 1.84 1.841 78 virt/peak/res/peak(MB): 33.97 33.97 10.51 10.51 79 80keiper_li.c 81------------------------------------------------------------------------------- 82 83Given an input integer *n*, this program rigorously computes numerical 84values of the Keiper-Li coefficients 85`\lambda_0, \ldots, \lambda_n`. The Keiper-Li coefficients 86have the property that `\lambda_n > 0` for all `n > 0` if and only if the 87Riemann hypothesis is true. This program was used for the record 88computations described in [Joh2013]_ (the paper describes 89the algorithm in some more detail). 90 91The program takes the following parameters:: 92 93 keiper_li n [-prec prec] [-threads num_threads] [-out out_file] 94 95The program prints the first and last few coefficients. It can optionally 96write all the computed data to a file. The working precision defaults 97to a value that should give all the coefficients to a few digits of 98accuracy, but can optionally be set higher (or lower). 99On a multicore system, using several threads results in faster 100execution. 101 102Sample output:: 103 104 > build/examples/keiper_li 1000 -threads 2 105 zeta: cpu/wall(s): 0.4 0.244 106 virt/peak/res/peak(MB): 167.98 294.69 5.09 7.43 107 log: cpu/wall(s): 0.03 0.038 108 gamma: cpu/wall(s): 0.02 0.016 109 binomial transform: cpu/wall(s): 0.01 0.018 110 0: -0.69314718055994530941723212145817656807550013436026 +/- 6.5389e-347 111 1: 0.023095708966121033814310247906495291621932127152051 +/- 2.0924e-345 112 2: 0.046172867614023335192864243096033943387066108314123 +/- 1.674e-344 113 3: 0.0692129735181082679304973488726010689942120263932 +/- 5.0219e-344 114 4: 0.092197619873060409647627872409439018065541673490213 +/- 2.0089e-343 115 5: 0.11510854289223549048622128109857276671349132303596 +/- 1.0044e-342 116 6: 0.13792766871372988290416713700341666356138966078654 +/- 6.0264e-342 117 7: 0.16063715965299421294040287257385366292282442046163 +/- 2.1092e-341 118 8: 0.18321945964338257908193931774721859848998098273432 +/- 8.4368e-341 119 9: 0.20565733870917046170289387421343304741236553410044 +/- 7.5931e-340 120 10: 0.22793393631931577436930340573684453380748385942738 +/- 7.5931e-339 121 991: 2.3196617961613367928373899656994682562101430813341 +/- 2.461e-11 122 992: 2.3203766239254884035349896518332550233162909717288 +/- 9.5363e-11 123 993: 2.321092061239733282811659116333262802034375592414 +/- 1.8495e-10 124 994: 2.3218073540188462110258826121503870112747188888893 +/- 3.5907e-10 125 995: 2.3225217392815185726928702951225314023773358152533 +/- 6.978e-10 126 996: 2.3232344485814623873333223609413703912358283071281 +/- 1.3574e-09 127 997: 2.3239447114886014522889542667580382034526509232475 +/- 2.6433e-09 128 998: 2.3246517591032700808344143240352605148856869322209 +/- 5.1524e-09 129 999: 2.3253548275861382119812576052060526988544993162101 +/- 1.0053e-08 130 1000: 2.3260531616864664574065046940832238158044982041872 +/- 3.927e-08 131 virt/peak/res/peak(MB): 170.18 294.69 7.51 7.51 132 133logistic.c 134------------------------------------------------------------------------------- 135 136This program computes the *n*-th iterate of the logistic map defined 137by `x_{n+1} = r x_n (1 - x_n)` where `r` and `x_0` are given. 138It takes the following parameters:: 139 140 logistic n [x_0] [r] [digits] 141 142The inputs `x_0`, *r* and *digits* default to 0.5, 3.75 and 10 respectively. 143The computation is automatically restarted with doubled precision 144until the result is accurate to *digits* decimal digits. 145 146Sample output:: 147 148 > build/examples/logistic 10 149 Trying prec=64 bits...success! 150 cpu/wall(s): 0 0.001 151 x_10 = [0.6453672908 +/- 3.10e-11] 152 153 > build/examples/logistic 100 154 Trying prec=64 bits...ran out of accuracy at step 18 155 Trying prec=128 bits...ran out of accuracy at step 53 156 Trying prec=256 bits...success! 157 cpu/wall(s): 0 0 158 x_100 = [0.8882939923 +/- 1.60e-11] 159 160 > build/examples/logistic 10000 161 Trying prec=64 bits...ran out of accuracy at step 18 162 Trying prec=128 bits...ran out of accuracy at step 53 163 Trying prec=256 bits...ran out of accuracy at step 121 164 Trying prec=512 bits...ran out of accuracy at step 256 165 Trying prec=1024 bits...ran out of accuracy at step 525 166 Trying prec=2048 bits...ran out of accuracy at step 1063 167 Trying prec=4096 bits...ran out of accuracy at step 2139 168 Trying prec=8192 bits...ran out of accuracy at step 4288 169 Trying prec=16384 bits...ran out of accuracy at step 8584 170 Trying prec=32768 bits...success! 171 cpu/wall(s): 0.859 0.858 172 x_10000 = [0.8242048008 +/- 4.35e-11] 173 174 > build/examples/logistic 1234 0.1 3.99 30 175 Trying prec=64 bits...ran out of accuracy at step 0 176 Trying prec=128 bits...ran out of accuracy at step 10 177 Trying prec=256 bits...ran out of accuracy at step 76 178 Trying prec=512 bits...ran out of accuracy at step 205 179 Trying prec=1024 bits...ran out of accuracy at step 461 180 Trying prec=2048 bits...ran out of accuracy at step 974 181 Trying prec=4096 bits...success! 182 cpu/wall(s): 0.009 0.009 183 x_1234 = [0.256445391958651410579677945635 +/- 3.92e-31] 184 185real_roots.c 186------------------------------------------------------------------------------- 187 188This program isolates the roots of a function on the interval `(a,b)` 189(where *a* and *b* are input as double-precision literals) 190using the routines in the :ref:`arb_calc <arb-calc>` module. 191The program takes the following arguments:: 192 193 real_roots function a b [-refine d] [-verbose] [-maxdepth n] [-maxeval n] [-maxfound n] [-prec n] 194 195The following functions (specified by an integer code) are implemented: 196 197 * 0 - `Z(x)` (Riemann-Siegel Z-function) 198 * 1 - `\sin(x)` 199 * 2 - `\sin(x^2)` 200 * 3 - `\sin(1/x)` 201 * 4 - `\operatorname{Ai}(x)` (Airy function) 202 * 5 - `\operatorname{Ai}'(x)` (Airy function) 203 * 6 - `\operatorname{Bi}(x)` (Airy function) 204 * 7 - `\operatorname{Bi}'(x)` (Airy function) 205 206The following options are available: 207 208 * ``-refine d``: If provided, after isolating the roots, attempt to refine 209 the roots to *d* digits of accuracy using a few bisection steps followed 210 by Newton's method with adaptive precision, and then print them. 211 212 * ``-verbose``: Print more information. 213 214 * ``-maxdepth n``: Stop searching after *n* recursive subdivisions. 215 216 * ``-maxeval n``: Stop searching after approximately *n* function evaluations 217 (the actual number evaluations will be a small multiple of this). 218 219 * ``-maxfound n``: Stop searching after having found *n* isolated roots. 220 221 * ``-prec n``: Working precision to use for the root isolation. 222 223With *function* 0, the program isolates roots of the Riemann zeta function 224on the critical line, and guarantees that no roots are missed 225(there are more efficient ways to do this, but it is a nice example):: 226 227 > build/examples/real_roots 0 0.0 50.0 -verbose 228 interval: [0, 50] 229 maxdepth = 30, maxeval = 100000, maxfound = 100000, low_prec = 30 230 found isolated root in: [14.111328125, 14.16015625] 231 found isolated root in: [20.99609375, 21.044921875] 232 found isolated root in: [25, 25.048828125] 233 found isolated root in: [30.419921875, 30.4443359375] 234 found isolated root in: [32.91015625, 32.958984375] 235 found isolated root in: [37.548828125, 37.59765625] 236 found isolated root in: [40.91796875, 40.966796875] 237 found isolated root in: [43.310546875, 43.3349609375] 238 found isolated root in: [47.998046875, 48.0224609375] 239 found isolated root in: [49.755859375, 49.7802734375] 240 --------------------------------------------------------------- 241 Found roots: 10 242 Subintervals possibly containing undetected roots: 0 243 Function evaluations: 3058 244 cpu/wall(s): 0.202 0.202 245 virt/peak/res/peak(MB): 26.12 26.14 2.76 2.76 246 247Find just one root and refine it to approximately 75 digits:: 248 249 > build/examples/real_roots 0 0.0 50.0 -maxfound 1 -refine 75 250 interval: [0, 50] 251 maxdepth = 30, maxeval = 100000, maxfound = 1, low_prec = 30 252 refined root (0/8): 253 [14.134725141734693790457251983562470270784257115699243175685567460149963429809 +/- 2.57e-76] 254 255 --------------------------------------------------------------- 256 Found roots: 1 257 Subintervals possibly containing undetected roots: 7 258 Function evaluations: 761 259 cpu/wall(s): 0.055 0.056 260 virt/peak/res/peak(MB): 26.12 26.14 2.75 2.75 261 262Find the first few roots of an Airy function and refine them to 50 digits each:: 263 264 > build/examples/real_roots 4 -10 0 -refine 50 265 interval: [-10, 0] 266 maxdepth = 30, maxeval = 100000, maxfound = 100000, low_prec = 30 267 refined root (0/6): 268 [-9.022650853340980380158190839880089256524677535156083 +/- 4.85e-52] 269 270 refined root (1/6): 271 [-7.944133587120853123138280555798268532140674396972215 +/- 1.92e-52] 272 273 refined root (2/6): 274 [-6.786708090071758998780246384496176966053882477393494 +/- 3.84e-52] 275 276 refined root (3/6): 277 [-5.520559828095551059129855512931293573797214280617525 +/- 1.05e-52] 278 279 refined root (4/6): 280 [-4.087949444130970616636988701457391060224764699108530 +/- 2.46e-52] 281 282 refined root (5/6): 283 [-2.338107410459767038489197252446735440638540145672388 +/- 1.48e-52] 284 285 --------------------------------------------------------------- 286 Found roots: 6 287 Subintervals possibly containing undetected roots: 0 288 Function evaluations: 200 289 cpu/wall(s): 0.003 0.003 290 virt/peak/res/peak(MB): 26.12 26.14 2.24 2.24 291 292Find roots of `\sin(x^2)` on `(0,100)`. The algorithm cannot isolate 293the root at `x = 0` (it is at the endpoint of the interval, and in any 294case a root of multiplicity higher than one). The failure is reported:: 295 296 > build/examples/real_roots 2 0 100 297 interval: [0, 100] 298 maxdepth = 30, maxeval = 100000, maxfound = 100000, low_prec = 30 299 --------------------------------------------------------------- 300 Found roots: 3183 301 Subintervals possibly containing undetected roots: 1 302 Function evaluations: 34058 303 cpu/wall(s): 0.032 0.032 304 virt/peak/res/peak(MB): 26.32 26.37 2.04 2.04 305 306This does not miss any roots:: 307 308 > build/examples/real_roots 2 1 100 309 interval: [1, 100] 310 maxdepth = 30, maxeval = 100000, maxfound = 100000, low_prec = 30 311 --------------------------------------------------------------- 312 Found roots: 3183 313 Subintervals possibly containing undetected roots: 0 314 Function evaluations: 34039 315 cpu/wall(s): 0.023 0.023 316 virt/peak/res/peak(MB): 26.32 26.37 2.01 2.01 317 318Looking for roots of `\sin(1/x)` on `(0,1)`, the algorithm finds many roots, 319but will never find all of them since there are infinitely many:: 320 321 > build/examples/real_roots 3 0.0 1.0 322 interval: [0, 1] 323 maxdepth = 30, maxeval = 100000, maxfound = 100000, low_prec = 30 324 --------------------------------------------------------------- 325 Found roots: 10198 326 Subintervals possibly containing undetected roots: 24695 327 Function evaluations: 202587 328 cpu/wall(s): 0.171 0.171 329 virt/peak/res/peak(MB): 28.39 30.38 4.05 4.05 330 331Remark: the program always computes rigorous containing intervals 332for the roots, but the accuracy after refinement could be less than *d* digits. 333 334poly_roots.c 335------------------------------------------------------------------------------- 336 337This program finds the complex roots of an integer polynomial 338by calling :func:`arb_fmpz_poly_complex_roots`, which in turn calls 339:func:`acb_poly_find_roots` with increasing 340precision until the roots certainly have been isolated. 341The program takes the following arguments:: 342 343 poly_roots [-refine d] [-print d] <poly> 344 345 Isolates all the complex roots of a polynomial with integer coefficients. 346 347 If -refine d is passed, the roots are refined to a relative tolerance 348 better than 10^(-d). By default, the roots are only computed to sufficient 349 accuracy to isolate them. The refinement is not currently done efficiently. 350 351 If -print d is passed, the computed roots are printed to d decimals. 352 By default, the roots are not printed. 353 354 The polynomial can be specified by passing the following as <poly>: 355 356 a <n> Easy polynomial 1 + 2x + ... + (n+1)x^n 357 t <n> Chebyshev polynomial T_n 358 u <n> Chebyshev polynomial U_n 359 p <n> Legendre polynomial P_n 360 c <n> Cyclotomic polynomial Phi_n 361 s <n> Swinnerton-Dyer polynomial S_n 362 b <n> Bernoulli polynomial B_n 363 w <n> Wilkinson polynomial W_n 364 e <n> Taylor series of exp(x) truncated to degree n 365 m <n> <m> The Mignotte-like polynomial x^n + (100x+1)^m, n > m 366 coeffs <c0 c1 ... cn> c0 + c1 x + ... + cn x^n 367 368 Concatenate to multiply polynomials, e.g.: p 5 t 6 coeffs 1 2 3 369 for P_5(x)*T_6(x)*(1+2x+3x^2) 370 371This finds the roots of the Wilkinson polynomial with roots at the 372positive integers 1, 2, ..., 100:: 373 374 > build/examples/poly_roots -print 15 w 100 375 computing squarefree factorization... 376 cpu/wall(s): 0.001 0.001 377 roots with multiplicity 1 378 searching for 100 roots, 100 deflated 379 prec=32: 0 isolated roots | cpu/wall(s): 0.098 0.098 380 prec=64: 0 isolated roots | cpu/wall(s): 0.247 0.247 381 prec=128: 0 isolated roots | cpu/wall(s): 0.498 0.497 382 prec=256: 0 isolated roots | cpu/wall(s): 0.713 0.713 383 prec=512: 100 isolated roots | cpu/wall(s): 0.104 0.105 384 done! 385 [1.00000000000000 +/- 3e-20] 386 [2.00000000000000 +/- 3e-19] 387 [3.00000000000000 +/- 1e-19] 388 [4.00000000000000 +/- 1e-19] 389 [5.00000000000000 +/- 1e-19] 390 ... 391 [96.0000000000000 +/- 1e-17] 392 [97.0000000000000 +/- 1e-17] 393 [98.0000000000000 +/- 3e-17] 394 [99.0000000000000 +/- 3e-17] 395 [100.000000000000 +/- 3e-17] 396 cpu/wall(s): 1.664 1.664 397 398This finds the roots of a Bernoulli polynomial which has both real 399and complex roots:: 400 401 > build/examples/poly_roots -refine 100 -print 20 b 16 402 computing squarefree factorization... 403 cpu/wall(s): 0.001 0 404 roots with multiplicity 1 405 searching for 16 roots, 16 deflated 406 prec=32: 16 isolated roots | cpu/wall(s): 0.006 0.006 407 prec=64: 16 isolated roots | cpu/wall(s): 0.001 0.001 408 prec=128: 16 isolated roots | cpu/wall(s): 0.001 0.001 409 prec=256: 16 isolated roots | cpu/wall(s): 0.001 0.002 410 prec=512: 16 isolated roots | cpu/wall(s): 0.002 0.001 411 done! 412 [-0.94308706466055783383 +/- 2.02e-21] 413 [-0.75534059252067985752 +/- 2.70e-21] 414 [-0.24999757119077421009 +/- 4.27e-21] 415 [0.24999757152512726002 +/- 4.43e-21] 416 [0.75000242847487273998 +/- 4.43e-21] 417 [1.2499975711907742101 +/- 1.43e-20] 418 [1.7553405925206798575 +/- 1.74e-20] 419 [1.9430870646605578338 +/- 3.21e-20] 420 [-0.99509334829256233279 +/- 9.42e-22] + [0.44547958157103608805 +/- 3.59e-21]*I 421 [-0.99509334829256233279 +/- 9.42e-22] + [-0.44547958157103608805 +/- 3.59e-21]*I 422 [1.9950933482925623328 +/- 1.10e-20] + [0.44547958157103608805 +/- 3.59e-21]*I 423 [1.9950933482925623328 +/- 1.10e-20] + [-0.44547958157103608805 +/- 3.59e-21]*I 424 [-0.92177327714429290564 +/- 4.68e-21] + [-1.0954360955079385542 +/- 1.71e-21]*I 425 [-0.92177327714429290564 +/- 4.68e-21] + [1.0954360955079385542 +/- 1.71e-21]*I 426 [1.9217732771442929056 +/- 3.54e-20] + [1.0954360955079385542 +/- 1.71e-21]*I 427 [1.9217732771442929056 +/- 3.54e-20] + [-1.0954360955079385542 +/- 1.71e-21]*I 428 cpu/wall(s): 0.011 0.012 429 430Roots are automatically separated by multiplicity by performing an initial 431squarefree factorization:: 432 433 > build/examples/poly_roots -print 5 p 5 p 5 t 7 coeffs 1 5 10 10 5 1 434 computing squarefree factorization... 435 cpu/wall(s): 0 0 436 roots with multiplicity 1 437 searching for 6 roots, 3 deflated 438 prec=32: 3 isolated roots | cpu/wall(s): 0 0.001 439 done! 440 [-0.97493 +/- 2.10e-6] 441 [-0.78183 +/- 1.49e-6] 442 [-0.43388 +/- 3.75e-6] 443 [0.43388 +/- 3.75e-6] 444 [0.78183 +/- 1.49e-6] 445 [0.97493 +/- 2.10e-6] 446 roots with multiplicity 2 447 searching for 4 roots, 2 deflated 448 prec=32: 2 isolated roots | cpu/wall(s): 0 0 449 done! 450 [-0.90618 +/- 1.56e-7] 451 [-0.53847 +/- 6.91e-7] 452 [0.53847 +/- 6.91e-7] 453 [0.90618 +/- 1.56e-7] 454 roots with multiplicity 3 455 searching for 1 roots, 0 deflated 456 prec=32: 0 isolated roots | cpu/wall(s): 0 0 457 done! 458 0 459 roots with multiplicity 5 460 searching for 1 roots, 1 deflated 461 prec=32: 1 isolated roots | cpu/wall(s): 0 0 462 done! 463 -1.0000 464 cpu/wall(s): 0 0.001 465 466complex_plot.c 467------------------------------------------------------------------------------- 468 469This program plots one of the predefined functions over a complex 470interval `[x_a, x_b] + [y_a, y_b]i` using domain coloring, at 471a resolution of *xn* times *yn* pixels. 472 473The program takes the parameters:: 474 475 complex_plot [-range xa xb ya yb] [-size xn yn] <func> 476 477Defaults parameters are `[-10,10] + [-10,10]i` and *xn* = *yn* = 512. 478 479A color function can be selected with -color. Valid options 480are 0 (phase=hue, magnitude=brightness) and 1 (phase only, 481white-gold-black-blue-white counterclockwise). 482 483The output is written to ``arbplot.ppm``. If you have ImageMagick, 484run ``convert arbplot.ppm arbplot.png`` to get a PNG. 485 486Function codes ``<func>`` are: 487 488 * ``gamma`` - Gamma function 489 * ``digamma`` - Digamma function 490 * ``lgamma`` - Logarithmic gamma function 491 * ``zeta`` - Riemann zeta function 492 * ``erf`` - Error function 493 * ``ai`` - Airy function Ai 494 * ``bi`` - Airy function Bi 495 * ``besselj`` - Bessel function `J_0` 496 * ``bessely`` - Bessel function `Y_0` 497 * ``besseli`` - Bessel function `I_0` 498 * ``besselk`` - Bessel function `K_0` 499 * ``modj`` - Modular j-function 500 * ``modeta`` - Dedekind eta function 501 * ``barnesg`` - Barnes G-function 502 * ``agm`` - Arithmetic geometric mean 503 504The function is just sampled at point values; no attempt is made to resolve 505small features by adaptive subsampling. 506 507For example, the following plots the Riemann zeta function around 508a portion of the critical strip with imaginary part between 100 and 140:: 509 510 > build/examples/complex_plot zeta -range -10 10 100 140 -size 256 512 511 512lvalue.c 513------------------------------------------------------------------------------- 514 515This program evaluates Dirichlet L-functions. It takes the following input:: 516 517 > build/examples/lvalue 518 lvalue [-character q n] [-re a] [-im b] [-prec p] [-z] [-deflate] [-len l] 519 520 Print value of Dirichlet L-function at s = a+bi. 521 Default a = 0.5, b = 0, p = 53, (q, n) = (1, 0) (Riemann zeta) 522 [-z] - compute Z(s) instead of L(s) 523 [-deflate] - remove singular term at s = 1 524 [-len l] - compute l terms in Taylor series at s 525 526Evaluating the Riemann zeta function and 527the Dirichlet beta function at `s = 2`:: 528 529 > build/examples/lvalue -re 2 -prec 128 530 L(s) = [1.64493406684822643647241516664602518922 +/- 4.37e-39] 531 cpu/wall(s): 0.001 0.001 532 virt/peak/res/peak(MB): 26.86 26.88 2.05 2.05 533 534 > build/examples/lvalue -character 4 3 -re 2 -prec 128 535 L(s) = [0.91596559417721901505460351493238411077 +/- 7.86e-39] 536 cpu/wall(s): 0.002 0.003 537 virt/peak/res/peak(MB): 26.86 26.88 2.31 2.31 538 539Evaluating the L-function for character number 101 modulo 1009 540at `s = 1/2` and `s = 1`:: 541 542 > build/examples/lvalue -character 1009 101 543 L(s) = [-0.459256562383872 +/- 5.24e-16] + [1.346937111206009 +/- 3.03e-16]*I 544 cpu/wall(s): 0.012 0.012 545 virt/peak/res/peak(MB): 26.86 26.88 2.30 2.30 546 547 > build/examples/lvalue -character 1009 101 -re 1 548 L(s) = [0.657952586112728 +/- 6.02e-16] + [1.004145273214022 +/- 3.10e-16]*I 549 cpu/wall(s): 0.017 0.018 550 virt/peak/res/peak(MB): 26.86 26.88 2.30 2.30 551 552Computing the first few coefficients in the Laurent series of the 553Riemann zeta function at `s = 1`:: 554 555 > build/examples/lvalue -re 1 -deflate -len 8 556 L(s) = [0.577215664901532861 +/- 5.29e-19] 557 L'(s) = [0.072815845483676725 +/- 2.68e-19] 558 [x^2] L(s+x) = [-0.004845181596436159 +/- 3.87e-19] 559 [x^3] L(s+x) = [-0.000342305736717224 +/- 4.20e-19] 560 [x^4] L(s+x) = [9.6890419394471e-5 +/- 2.40e-19] 561 [x^5] L(s+x) = [-6.6110318108422e-6 +/- 4.51e-20] 562 [x^6] L(s+x) = [-3.316240908753e-7 +/- 3.85e-20] 563 [x^7] L(s+x) = [1.0462094584479e-7 +/- 7.78e-21] 564 cpu/wall(s): 0.003 0.004 565 virt/peak/res/peak(MB): 26.86 26.88 2.30 2.30 566 567Evaluating the Riemann zeta function near the first nontrivial root:: 568 569 > build/examples/lvalue -re 0.5 -im 14.134725 570 L(s) = [1.76743e-8 +/- 1.93e-14] + [-1.110203e-7 +/- 2.84e-14]*I 571 cpu/wall(s): 0.001 0.001 572 virt/peak/res/peak(MB): 26.86 26.88 2.31 2.31 573 574 > build/examples/lvalue -z -re 14.134725 -prec 200 575 Z(s) = [-1.12418349839417533300111494358128257497862927935658e-7 +/- 4.62e-58] 576 cpu/wall(s): 0.001 0.001 577 virt/peak/res/peak(MB): 26.86 26.88 2.57 2.57 578 579 > build/examples/lvalue -z -re 14.134725 -len 4 580 Z(s) = [-1.124184e-7 +/- 7.00e-14] 581 Z'(s) = [0.793160414884 +/- 4.09e-13] 582 [x^2] Z(s+x) = [0.065164586492 +/- 5.39e-13] 583 [x^3] Z(s+x) = [-0.020707762705 +/- 5.37e-13] 584 cpu/wall(s): 0.002 0.003 585 virt/peak/res/peak(MB): 26.86 26.88 2.57 2.57 586 587lcentral.c 588------------------------------------------------------------------------------- 589 590This program computes the central value `L(1/2)` for each Dirichlet L-function 591character modulo *q* for each *q* in the range *qmin* to *qmax*. Usage:: 592 593 > build/examples/lcentral 594 Computes central values (s = 0.5) of Dirichlet L-functions. 595 596 usage: build/examples/lcentral [--quiet] [--check] [--prec <bits>] qmin qmax 597 598The first few values:: 599 600 > build/examples/lcentral 1 8 601 3,2: [0.48086755769682862618122006324 +/- 7.35e-30] 602 4,3: [0.66769145718960917665869092930 +/- 1.62e-30] 603 5,2: [0.76374788011728687822451215264 +/- 2.32e-30] + [0.21696476751886069363858659310 +/- 3.06e-30]*I 604 5,4: [0.23175094750401575588338366176 +/- 2.21e-30] 605 5,3: [0.76374788011728687822451215264 +/- 2.32e-30] + [-0.21696476751886069363858659310 +/- 3.06e-30]*I 606 7,3: [0.71394334376831949285993820742 +/- 1.21e-30] + [0.47490218277139938263745243935 +/- 4.52e-30]*I 607 7,2: [0.31008936259836766059195052534 +/- 5.29e-30] + [-0.07264193137017790524562171245 +/- 5.48e-30]*I 608 7,6: [1.14658566690370833367712697646 +/- 1.95e-30] 609 7,4: [0.31008936259836766059195052534 +/- 5.29e-30] + [0.07264193137017790524562171245 +/- 5.48e-30]*I 610 7,5: [0.71394334376831949285993820742 +/- 1.21e-30] + [-0.47490218277139938263745243935 +/- 4.52e-30]*I 611 8,5: [0.37369171291254730738158695002 +/- 4.01e-30] 612 8,3: [1.10042140952554837756713576997 +/- 3.37e-30] 613 cpu/wall(s): 0.002 0.003 614 virt/peak/res/peak(MB): 26.32 26.34 2.35 2.35 615 616Testing a large *q*:: 617 618 > build/examples/lcentral --quiet --check --prec 256 100000 100000 619 cpu/wall(s): 1.668 1.667 620 virt/peak/res/peak(MB): 35.67 46.66 11.67 22.61 621 622It is conjectured that the central value never vanishes. Running with ``--check`` 623verifies that the interval certainly is nonzero. This can fail with 624insufficient precision:: 625 626 > build/examples/lcentral --check --prec 15 100000 100000 627 100000,71877: [0.1 +/- 0.0772] + [+/- 0.136]*I 628 100000,90629: [2e+0 +/- 0.106] + [+/- 0.920]*I 629 100000,28133: [+/- 0.811] + [-2e+0 +/- 0.501]*I 630 100000,3141: [0.8 +/- 0.0407] + [-0.1 +/- 0.0243]*I 631 100000,53189: [4.0 +/- 0.0826] + [+/- 0.107]*I 632 100000,53253: [1.9 +/- 0.0855] + [-3.9 +/- 0.0681]*I 633 Value could be zero! 634 100000,53381: [+/- 0.0329] + [+/- 0.0413]*I 635 Aborted 636 637integrals.c 638------------------------------------------------------------------------------- 639 640This program computes integrals using :func:`acb_calc_integrate`. 641Invoking the program without parameters shows usage:: 642 643 > build/examples/integrals 644 Compute integrals using acb_calc_integrate. 645 Usage: integrals -i n [-prec p] [-tol eps] [-twice] [...] 646 647 -i n - compute integral n (0 <= n <= 23), or "-i all" 648 -prec p - precision in bits (default p = 64) 649 -goal p - approximate relative accuracy goal (default p) 650 -tol eps - approximate absolute error goal (default 2^-p) 651 -twice - run twice (to see overhead of computing nodes) 652 -heap - use heap for subinterval queue 653 -verbose - show information 654 -verbose2 - show more information 655 -deg n - use quadrature degree up to n 656 -eval n - limit number of function evaluations to n 657 -depth n - limit subinterval queue size to n 658 659 Implemented integrals: 660 I0 = int_0^100 sin(x) dx 661 I1 = 4 int_0^1 1/(1+x^2) dx 662 I2 = 2 int_0^{inf} 1/(1+x^2) dx (using domain truncation) 663 I3 = 4 int_0^1 sqrt(1-x^2) dx 664 I4 = int_0^8 sin(x+exp(x)) dx 665 I5 = int_1^101 floor(x) dx 666 I6 = int_0^1 |x^4+10x^3+19x^2-6x-6| exp(x) dx 667 I7 = 1/(2 pi i) int zeta(s) ds (closed path around s = 1) 668 I8 = int_0^1 sin(1/x) dx (slow convergence, use -heap and/or -tol) 669 I9 = int_0^1 x sin(1/x) dx (slow convergence, use -heap and/or -tol) 670 I10 = int_0^10000 x^1000 exp(-x) dx 671 I11 = int_1^{1+1000i} gamma(x) dx 672 I12 = int_{-10}^{10} sin(x) + exp(-200-x^2) dx 673 I13 = int_{-1020}^{-1010} exp(x) dx (use -tol 0 for relative error) 674 I14 = int_0^{inf} exp(-x^2) dx (using domain truncation) 675 I15 = int_0^1 sech(10(x-0.2))^2 + sech(100(x-0.4))^4 + sech(1000(x-0.6))^6 dx 676 I16 = int_0^8 (exp(x)-floor(exp(x))) sin(x+exp(x)) dx (use higher -eval) 677 I17 = int_0^{inf} sech(x) dx (using domain truncation) 678 I18 = int_0^{inf} sech^3(x) dx (using domain truncation) 679 I19 = int_0^1 -log(x)/(1+x) dx (using domain truncation) 680 I20 = int_0^{inf} x exp(-x)/(1+exp(-x)) dx (using domain truncation) 681 I21 = int_C wp(x)/x^(11) dx (contour for 10th Laurent coefficient of Weierstrass p-function) 682 I22 = N(1000) = count zeros with 0 < t <= 1000 of zeta(s) using argument principle 683 I23 = int_0^{1000} W_0(x) dx 684 I24 = int_0^pi max(sin(x), cos(x)) dx 685 I25 = int_{-1}^1 erf(x/sqrt(0.0002)*0.5+1.5)*exp(-x) dx 686 I26 = int_{-10}^10 Ai(x) dx 687 I27 = int_0^10 (x-floor(x)-1/2) max(sin(x),cos(x)) dx 688 I28 = int_{-1-i}^{-1+i} sqrt(x) dx 689 I29 = int_0^{inf} exp(-x^2+ix) dx (using domain truncation) 690 I30 = int_0^{inf} exp(-x) Ai(-x) dx (using domain truncation) 691 I31 = int_0^pi x sin(x) / (1 + cos(x)^2) dx 692 693A few examples:: 694 695 build/examples/integrals -i 4 696 I4 = int_0^8 sin(x+exp(x)) dx ... 697 cpu/wall(s): 0.02 0.02 698 I4 = [0.34740017265725 +/- 3.95e-15] 699 700 > build/examples/integrals -i 3 -prec 333 -tol 1e-80 701 I3 = 4 int_0^1 sqrt(1-x^2) dx ... 702 cpu/wall(s): 0.024 0.024 703 I3 = [3.141592653589793238462643383279502884197169399375105820974944592307816406286209 +/- 4.24e-79] 704 705 > build/examples/integrals -i 9 -heap 706 I9 = int_0^1 x sin(1/x) dx (slow convergence, use -heap and/or -tol) ... 707 cpu/wall(s): 0.019 0.018 708 I9 = [0.3785300 +/- 3.17e-8] 709 710fpwrap.c 711------------------------------------------------------------------------------- 712 713This program demonstrates calling the floating-point wrapper:: 714 715 > build/examples/fpwrap 716 zeta(2) = 1.644934066848226 717 zeta(0.5 + 123i) = 0.006252861175594465 + 0.08206030514520983i 718 719 720.. highlight:: c 721 722