C library for arbitrary-precision ball arithmetic

Overview

Arb

Arb is a C library for arbitrary-precision interval arithmetic. It has full support for both real and complex numbers. The library is thread-safe, portable, and extensively tested. Arb is free software distributed under the GNU Lesser General Public License (LGPL), version 2.1 or later.

arb logo

Documentation: http://arblib.org

Development updates: http://fredrikj.net/blog/

Author: Fredrik Johansson [email protected]

Bug reports, feature requests and other comments are welcome in private communication, on the GitHub issue tracker, or on the FLINT mailing list [email protected].

Build Status

Code example

The following program evaluates sin(pi + exp(-10000)). Since the input to the sine function matches a root to within 4343 digits, at least 4343-digit (14427-bit) precision is needed to get an accurate result. The program repeats the evaluation at 64-bit, 128-bit, ... precision, stopping only when the result is accurate to at least 53 bits.

#include "arb.h"

int main()
{
    slong prec;
    arb_t x, y;
    arb_init(x); arb_init(y);

    for (prec = 64; ; prec *= 2)
    {
        arb_const_pi(x, prec);
        arb_set_si(y, -10000);
        arb_exp(y, y, prec);
        arb_add(x, x, y, prec);
        arb_sin(y, x, prec);
        arb_printn(y, 15, 0); printf("\n");
        if (arb_rel_accuracy_bits(y) >= 53)
            break;
    }

    arb_clear(x); arb_clear(y);
    flint_cleanup();
}

The output is:

[+/- 6.01e-19]
[+/- 2.55e-38]
[+/- 8.01e-77]
[+/- 8.64e-154]
[+/- 5.37e-308]
[+/- 3.63e-616]
[+/- 1.07e-1232]
[+/- 9.27e-2466]
[-1.13548386531474e-4343 +/- 3.91e-4358]

Each line shows a rigorous enclosure of the exact value of the expression. The program demonstrates how the user can rely on Arb's automatic error bound tracking to get an output that is guaranteed to be accurate -- no error analysis needs to be done by the user.

For more example programs, see: http://arblib.org/examples.html

Features

Besides basic arithmetic, Arb allows working with univariate polynomials, truncated power series, and matrices over both real and complex numbers.

Basic linear algebra is supported, including matrix multiplication, determinant, inverse, nonsingular solving, matrix exponential, and computation of eigenvalues and eigenvectors.

Support for polynomials and power series is quite extensive, including methods for composition, reversion, product trees, multipoint evaluation and interpolation, complex root isolation, and transcendental functions of power series.

Other features include root isolation for real functions, rigorous numerical integration of complex functions, and discrete Fourier transforms (DFTs).

Special functions

Arb can compute a wide range of transcendental and special functions, including the gamma function, polygamma functions, Riemann zeta and Hurwitz zeta function, Dirichlet L-functions, polylogarithm, error function, Gauss hypergeometric function 2F1, confluent hypergeometric functions, Bessel functions, Airy functions, Legendre functions and other orthogonal polynomials, exponential and trigonometric integrals, incomplete gamma and beta functions, Jacobi theta functions, modular functions, Weierstrass elliptic functions, complete and incomplete elliptic integrals, arithmetic-geometric mean, Bernoulli numbers, partition function, Barnes G-function, Lambert W function.

Speed

Arb uses a midpoint-radius (ball) representation of real numbers. At high precision, this allows doing interval arithmetic without significant overhead compared to plain floating-point arithmetic. Various low-level optimizations have also been implemented to reduce overhead at precisions of just a few machine words. Most operations on polynomials and power series use asymptotically fast FFT multiplication based on FLINT. Similarly, most operations on large matrices take advantage of the fast integer matrix multiplication in FLINT.

For basic arithmetic, Arb should generally be around as fast as MPFR (http://mpfr.org), though it can be a bit slower at low precision, and around twice as fast as MPFI (https://perso.ens-lyon.fr/nathalie.revol/software.html).

Transcendental functions in Arb are quite well optimized and should generally be faster than any other arbitrary-precision software currently available. The following table compares the time in seconds to evaluate the Gauss hypergeometric function 2F1(1/2, 1/4, 1, z) at the complex number z = 5^(1/2) + 7^(1/2)i, to a given number of decimal digits (Arb 2.8-git and mpmath 0.19 on an 1.90 GHz Intel i5-4300U, Mathematica 9.0 on a 3.07 GHz Intel Xeon X5675).

Digits Mathematica mpmath Arb
10 0.00066 0.00065 0.000071
100 0.0039 0.0012 0.00048
1000 0.23 1.2 0.0093
10000 42.6 84 0.56

Dependencies, installation, and interfaces

Arb depends on FLINT (http://flintlib.org/), either GMP (http://gmplib.org) or MPIR (http://mpir.org), and MPFR (http://mpfr.org).

See http://arblib.org/setup.html for instructions on building and installing Arb directly from the source code. Arb might also be available (or coming soon) as a package for your Linux distribution.

SageMath (http://sagemath.org/) includes Arb as a standard package and contains a high-level Python interface. See the SageMath documentation for RealBallField (http://doc.sagemath.org/html/en/reference/rings_numerical/sage/rings/real_arb.html) and ComplexBallField (http://doc.sagemath.org/html/en/reference/rings_numerical/sage/rings/complex_arb.html).

Nemo (http://nemocas.org/) is a computer algebra package for the Julia programming language which includes a high-level Julia interface to Arb. The Nemo installation script will create a local installation of Arb along with other dependencies.

A standalone Python interface to FLINT and Arb is also available (https://github.com/fredrik-johansson/python-flint).

A separate wrapper of transcendental functions for use with the C99 complex double type is available (https://github.com/fredrik-johansson/arbcmath).

Other third-party wrappers include:

Comments
  • flint.h or flint/flint.h?

    flint.h or flint/flint.h?

    Currently, some files in Arb (e.g., fmpr.h) include the FLINT headers as #include "flint.h". I am wondering if this should be #include "flint/flint.h" instead, as it seems that FLINT headers are commonly installed in a flint/ subdirectory (e.g., /usr/include/flint/flint.h on my setup).

    The current way of including flint.h essentially prevents to compile files including the Arb headers without adding extra -I switches to the compiler command line. Or am I missing something?

    opened by bluescarni 43
  • testsuite failure on arm*

    testsuite failure on arm*

    Hello, the new releases are having segfaults on arm platforms during testsuite...

    make[3]: Leaving directory '/<<PKGBUILDDIR>>/acb_poly'
    make[3]: Entering directory '/<<PKGBUILDDIR>>/acb_dft'
    gcc -Wdate-time -D_FORTIFY_SOURCE=2 -g -O2 -fdebug-prefix-map=/<<PKGBUILDDIR>>=. -fstack-protector-strong -Wformat -Werror=format-security -I/<<PKGBUILDDIR>> -I/usr/local/include -I/usr/local/include -I/usr/include test/t-convol.c -o ../build/acb_dft/test/t-convol -L/<<PKGBUILDDIR>> -L/usr/local/lib -L/usr/local/lib -L/usr/lib -lflint-arb -lflint -lmpfr -lgmp -lm -lpthread  -MMD -MP -MF ../build/acb_dft/test/t-convol.d -MT "../build/acb_dft/test/t-convol" -MT "../build/acb_dft/test/t-convol.d"
    gcc -Wdate-time -D_FORTIFY_SOURCE=2 -g -O2 -fdebug-prefix-map=/<<PKGBUILDDIR>>=. -fstack-protector-strong -Wformat -Werror=format-security -I/<<PKGBUILDDIR>> -I/usr/local/include -I/usr/local/include -I/usr/include test/t-dft.c -o ../build/acb_dft/test/t-dft -L/<<PKGBUILDDIR>> -L/usr/local/lib -L/usr/local/lib -L/usr/lib -lflint-arb -lflint -lmpfr -lgmp -lm -lpthread  -MMD -MP -MF ../build/acb_dft/test/t-dft.d -MT "../build/acb_dft/test/t-dft" -MT "../build/acb_dft/test/t-dft.d"
    convol....dft....free(): invalid pointer
    make[3]: *** [../Makefile.subdirs:84: ../build/acb_dft/test/t-dft_RUN] Aborted (core dumped)
    make[3]: *** Waiting for unfinished jobs....
    

    can you please have a look? see e.g. https://launchpad.net/ubuntu/+source/flint-arb/1:2.12.0-3/+build/15056296

    or the new version https://launchpadlibrarian.net/380256280/buildlog_ubuntu-cosmic-ppc64el.flint-arb_1%3A2.14.0-1~build1_BUILDING.txt.gz

    opened by LocutusOfBorg 40
  • sin_cos test suite failure

    sin_cos test suite failure

    I thought this issue was due to interference from MLK using the intel compiler, but when I installed using gcc on my univerisity linux server (icc is not installed), running the test suite gives errors for the sin_cos function.

    I started out by locally installing gmp-6.0.0a, mpfr-3.1.3, and flint-2.5.2. To install gmp, I used the following commands:

    (let MYPATH be the desired install directory, it doesn't change)

    install gmp

    ./configure --prefix=MYPATH --disable-static make -j8 make check cd make tune make speed make tune make tune | tee tune.out

    make a copy of the original gmp-mparam.h

    cp ../mpn/x86_64/core2/gmp-mparam.h ../mpn/x86_64/core2/gmp-mparam.h~

    I remove the first several lines of tune.out, so that it can be used as a new gmp-mparam.h file

    afterwards, I can then replace with the contents of tune.out

    cp tune.out ../mpn/x86_64/core2/gmp-mparam.h

    gmp needs to be rebuilt since we tuned

    cd .. && make -j8 make check

    all tests pass

    make install


    Now, I install mpfr. I navigate to it's directory, and the path for installation will be the same as MYPATH above.

    ./configure --prefix=MYPATH --with-gmp-build=MYPATH/local/gmp-6.0.0 --disable-static make -j8 cd tune make check

    ALL tests pass here

    cd tune make tune cd .. make check

    test one more time

    make install


    Now, I install flint similar to above using the commands:

    ./configure --prefix=MYPATH --with-gmp=MYPATH --with-mpfr=MYPATH --disable-static make -j8 make check

    All tests pass

    make install


    Now, I make a local repository of my fork of arb: git clone https://github.com/rickyefarr/arb.git cd arb git remote add upstream https://github.com/fredrik-johansson/arb.git git fetch upstream git checkout master git merge upstream/master ./configure --prefix=MYPATH --with-gmp=MYPATH --with-mpfr=MYPATH --with-flint=MYPATH --disable-static make -j8 make check

    The failure says:

    sin_cos....FAIL: containment (sin) a = (-448128001983 * 2^14) +/- (536870912 * 2^-323) b = (157968445807803538695138812109622709659517717 * 2^-147) +/- (541065218 * 2^-177) make[1]: *** [../build/arb/test/t-sin_cos_RUN] Aborted (core dumped)

    Have I installed something incorrectly?

    opened by rickyefarr 40
  • Experiment with function acb_dirichlet_platt_local_hardy_z_zeros

    Experiment with function acb_dirichlet_platt_local_hardy_z_zeros

    For a project, I like to compute the first 10,000 zeros of Z(t), starting at the n-th zero with n=10^k, k=5..20.

    The function acb_dirichlet_hardy_z_zeros is feasible to use up to 10^15, however for higher k Platt's method is required and this has been accommodated in acb_dirichlet_platt_local_hardy_z_zeros. This function works impressively fast and I used it to successfully confirm the results up 10^15 and also managed to compute the required 10K zeros at n=10^16, 10^17 and 10^19.

    However, for 10^18 and 10^20, the function frequently fails to produce the zeros or aborts the code with error message: "expected the first node in the list, Abort trap: 6". Example code that generates the error message is included below (n=10^18 + 200 and 100 zeros are at most computed).

    I have tried varying the count size, the prec and experimented with different starting heights around 10^18, however nothing seems to fix the issues. I do realise the code has been developed by the 'the mysterious D.H.J. Polymath' and also that the many parameters required for Platt's method are quite 'delicate' to properly establish automatically, however I'd still be grateful for any thoughts on this.

    #include "acb_dirichlet.h"
    #include "flint/profiler.h"
    
    int main()
    {
        slong i, count, prec;
     
        count = 100;
        prec = 512;
    
        arb_ptr pa;
        pa = _arb_vec_init(count);
        
        fmpz_t n;
        fmpz_init(n);
        fmpz_set_str(n, "1000000000000000200", 10);
    
        printf("Platt_local_hardy_z_zeros....\n");
    
        acb_dirichlet_platt_local_hardy_z_zeros(pa, n, count, prec);
    
        for (i = 0; i < count; i++)
        {
            arb_printd(pa+i, 30); printf("\n");
        }
    
        fmpz_clear(n);
    	
        _arb_vec_clear(pa, count);
    
        flint_cleanup();
    
        return 0;
    }
    
    opened by rudolph-git-acc 37
  • Linear solving and matrix inverse should be more clever

    Linear solving and matrix inverse should be more clever

    Gaussian elimination done directly in interval arithmetic gives poor enclosures for large matrices. An old paper that describes a few different methods that are easy to implement is Hansen & Smith, Interval Arithmetic in Matrix Computations, Part II, SIAM J. Numer. Anal., 4(1), 1–9. There are new papers that describe fancier techniques, but even the simple preconditioning trick (method 4) is a vast improvement. Example:

    from flint import *
    from mpmath import *
    
    n = 40
    A = randmatrix(n,n)
    B = randmatrix(n,1)
    AA = arb_mat(n, n, [arb(A[i//n,i%n], 1e-14) for i in range(n*n)])
    BB = arb_mat(n, 1, [arb(B[i,0], 1e-14) for i in range(n*1)])
    AI = inverse(A)   # approximate inverse with mpmath
    AAI = arb_mat(n, n, [AI[i//n,i%n] for i in range(n*n)])
    AA.solve(BB)   # naive interval solution
    (AAI*AA).solve(AAI*BB)   # preconditioned interval solution
    

    The naive solution in one instance:

    [    [+/- 1.32e+5]]
    [    [+/- 6.04e+4]]
    [    [+/- 3.00e+4]]
    [    [+/- 4.58e+4]]
    [    [+/- 2.93e+4]]
    [    [+/- 2.09e+4]]
    [    [+/- 9.04e+3]]
    [    [+/- 2.93e+3]]
    [    [+/- 2.62e+3]]
    [    [+/- 1.93e+3]]
    [    [+/- 1.50e+3]]
    [    [+/- 8.99e+2]]
    [    [+/- 1.33e+3]]
    [    [+/- 6.14e+2]]
    [    [+/- 3.23e+2]]
    [    [+/- 3.24e+2]]
    [    [+/- 2.37e+2]]
    [    [+/- 1.49e+2]]
    [    [+/- 1.34e+2]]
    [       [+/- 91.3]]
    [       [+/- 70.9]]
    [       [+/- 11.8]]
    [       [+/- 12.3]]
    [       [+/- 24.4]]
    [       [+/- 9.99]]
    [       [+/- 7.60]]
    [       [+/- 11.4]]
    [       [+/- 3.73]]
    [       [+/- 1.89]]
    [       [+/- 3.64]]
    [       [+/- 1.50]]
    [       [+/- 2.56]]
    [       [+/- 1.24]]
    [      [+/- 0.732]]
    [      [+/- 0.447]]
    [      [+/- 0.896]]
    [ [0.5 +/- 0.0759]]
    [      [+/- 0.521]]
    [      [+/- 0.212]]
    [[-0.7 +/- 0.0472]]
    

    The preconditioned solution:

    [ [0.02839853147 +/- 7.90e-12]]
    [ [-1.0500787426 +/- 4.48e-11]]
    [[-0.80417013789 +/- 2.33e-12]]
    [ [0.10241135449 +/- 6.26e-12]]
    [  [0.2737351746 +/- 2.11e-11]]
    [ [0.07527146134 +/- 5.56e-12]]
    [[-0.36775106771 +/- 6.79e-12]]
    [ [0.24625560798 +/- 2.87e-12]]
    [ [-0.4328648101 +/- 1.51e-11]]
    [  [0.1132582518 +/- 3.43e-11]]
    [  [0.6887292062 +/- 4.07e-11]]
    [ [0.33669069978 +/- 5.46e-12]]
    [ [0.38419150904 +/- 5.89e-12]]
    [[-0.19803441003 +/- 3.70e-12]]
    [[-0.06555745595 +/- 7.70e-12]]
    [  [0.2793057899 +/- 3.17e-11]]
    [[-0.22216454080 +/- 6.48e-12]]
    [[-0.94981852995 +/- 6.97e-12]]
    [[-0.44283327345 +/- 8.45e-12]]
    [ [-0.0066854861 +/- 6.38e-11]]
    [  [0.9304996076 +/- 4.09e-11]]
    [ [1.42913754750 +/- 8.34e-12]]
    [ [0.44546042722 +/- 6.36e-12]]
    [  [0.3450231038 +/- 4.13e-11]]
    [ [0.18352732449 +/- 5.73e-12]]
    [[-0.63541868806 +/- 4.49e-12]]
    [ [0.08467962923 +/- 3.24e-12]]
    [  [0.1062654616 +/- 5.91e-11]]
    [ [0.08812065325 +/- 5.94e-12]]
    [[-1.05566529972 +/- 9.66e-12]]
    [  [0.1417693049 +/- 4.45e-11]]
    [ [0.74041981272 +/- 8.82e-12]]
    [  [0.4540493110 +/- 9.64e-12]]
    [  [0.3964147231 +/- 4.49e-11]]
    [ [-0.1807095930 +/- 3.49e-11]]
    [[-0.67377800420 +/- 6.70e-12]]
    [ [0.50512672666 +/- 7.59e-12]]
    [  [0.4187979360 +/- 4.44e-11]]
    [  [0.1328337609 +/- 2.10e-11]]
    [[-0.72121349106 +/- 2.93e-12]]
    

    The approximate solving could use doubles and even BLAS where available (with careful checking whether overflow or underflow is a problem).

    Also, the preconditioned interval solution does not need to use interval Gaussian elimination; since one has almost the inverse of the identity matrix, it should be sufficient to compute some perturbation bounds in O(n^2) instead of O(n^3).

    opened by fredrik-johansson 27
  • Implementation of elliptic functions/integrals

    Implementation of elliptic functions/integrals

    Not really an issue, more like a hopefully polite request :)

    I am an enthusiastic user of the elliptic functions/integrals module in mpmath. Any plans for implementing them in arb as well?

    opened by bluescarni 25
  • Arb fails to link under vcpkg.

    Arb fails to link under vcpkg.

    The following example fails to link. Arb version is 2.21.1. Windows.

    src/main.cpp:

    #include <malloc.h>
    #include <arb.h>
    
    int main()
    {
        arb_t x;
        arb_init(x);
        arb_const_pi(x, 50 * 3.33);
        arb_printn(x, 50, 0); flint_printf("\n");
        flint_printf("Computed with arb-%s\n", arb_version);
        arb_clear(x);
    }
    

    vcpkg.json:

    {
      "name": "testarb",
      "version-string": "0.1.0",
      "dependencies": [
        "arb"
      ]
    }
    

    CMakeLists.txt:

    cmake_minimum_required(VERSION 3.15)
    
    project(testarb CXX)
    
    include_directories("./vcpkg/packages/arb_x64-windows/include")
    include_directories("./vcpkg/packages/flint_x64-windows/include")
    include_directories("./vcpkg/packages/gmp_x64-windows/include")
    include_directories("./vcpkg/packages/mpfr_x64-windows/include")
    include_directories("./vcpkg/packages/pthreads_x64-windows/include")
    
    link_directories("./vcpkg/packages/arb_x64-windows/lib")
    link_directories("./vcpkg/packages/flint_x64-windows/lib")
    link_directories("./vcpkg/packages/gmp_x64-windows/lib")
    link_directories("./vcpkg/packages/mpfr_x64-windows/lib")
    link_directories("./vcpkg/packages/pthreads_x64-windows/lib")
    
    add_executable(testarb src/main.cpp)
    target_compile_features(testarb PRIVATE cxx_std_20)
    
    target_link_libraries(testarb PRIVATE gmp flint mpfr arb)
    

    Command line:

    git clone https://github.com/Microsoft/vcpkg.git
    .\vcpkg\bootstrap-vcpkg.bat
    cmake -B build -S . -DCMAKE_TOOLCHAIN_FILE=./vcpkg/scripts/buildsystems/vcpkg.cmake
    cmake --build build
    

    Output:

    main.obj : error LNK2019: unresolved external symbol __imp_arb_version referenced in function main [C:\work\testarb\bui
    ld\testarb.vcxproj]
    arb.lib(div.c.obj) : error LNK2019: unresolved external symbol __gmpn_div_q referenced in function arf_div [C:\work\tes
    tarb\build\testarb.vcxproj]
    

    If I open gmp.lib in 7-zip, I can find __gmpn_div_qr_1 and __gmpn_div_qr_2, but not __gmpn_div_q. How can I make this example link?

    opened by kaba2 17
  • Dirichlet

    Dirichlet

    base PR for Dirichlet module, thanks for extending l_hurwitz.

    For L-functions, I am working on l_incgam (simple) and then I will add the Fourier method.

    opened by pascalmolin 17
  • enh: linear algebra functions for symmetric positive definite (spd) matrices

    enh: linear algebra functions for symmetric positive definite (spd) matrices

    This is based on Cholesky factorization, and includes docs and tests.

    As an example of the accuracy improvement, here are error magnitudes for inverses and determinants of various sizes of Pascal matrices with prec=96, using the generic inv and det vs. the spd_inv and spd_det (check mark indicates the solution is exact, and ✗ indicates utter failure).

    | n | inv | spd_inv | det | spd_det | | --- | --- | --- | --- | --- | | 2 | ✓ | ✓ | ✓ | ✓ | | 3 | ✓ | ✓ | ✓ | ✓ | | 4 | 2.198e-25 | ✓ | 3.33e-27 | ✓ | | 5 | ✓ | ✓ | ✓ | ✓ | | 6 | 1.02e-19 | ✓ | 4.597e-24 | ✓ | | 7 | 1.267e-16 | ✓ | 1.888e-22 | ✓ | | 8 | 2.271e-13 | ✓ | 1.235e-20 | ✓ | | 9 | 6.285e-14 | ✓ | 9.331e-23 | ✓ | | 10 | 2.288e-06 | ✓ | 7.317e-17 | ✓ | | 11 | 0.007739 | ✓ | 4.992e-15 | ✓ | | 12 | 39.61 | ✓ | 4.873e-13 | ✓ | | 13 | 5.275e+04 | ✓ | 1.22e-11 | ✓ | | 14 | 2.064e+09 | ✓ | 6.186e-09 | ✓ | | 15 | 1.849e+13 | ✓ | 6.991e-07 | ✓ | | 16 | 1.757e+17 | ✓ | 9.328e-05 | ✓ | | 17 | 3.521e+18 | ✓ | 2.561e-05 | ✓ | | 18 | 1.24e+27 | ✓ | 0.9967 | ✓ | | 19 | ✗ | 1.706e+30 | 1.883e+04 | 1.486e-07 | | 20 | ✗ | 4.526e+34 | 1.073e+10 | 5.248e-06 | | 21 | ✗ | 2.31e+40 | 3.194e+16 | 0.003213 | | 22 | ✗ | ✗ | 7.323e+23 | ✗ | | 23 | ✗ | ✗ | 1.751e+35 | ✗ |

    opened by argriffing 16
  • ENH: use sparsity structure more intelligently for matrix power series

    ENH: use sparsity structure more intelligently for matrix power series

    Avoid adding unnecessary matrix power truncation terms to functions of nilpotent or "partially nilpotent" matrices. An example (maybe unnecessarily contrived) is the nilpotent matrix

    A =
    [0 2^-100 0]
    [0 0 2^100]
    [0 0 0]
    

    for which the matrix exponential with prec=32 used to have large error bounds

    [1 +/- 0, 7.88860905221012e-31 +/- 2.0842e+13, 0.5 +/- 1.321e+43]
    [0 +/- 0, 1 +/- 0, 1.26765060022823e+30 +/- 2.0842e+13]
    [0 +/- 0, 0 +/- 0, 1 +/- 0]
    

    but which is now recognized as having no error:

    [1 +/- 0, 7.88860905221012e-31 +/- 0, 0.5 +/- 0]
    [0 +/- 0, 1 +/- 0, 1.26765060022823e+30 +/- 0]
    [0 +/- 0, 0 +/- 0, 1 +/- 0]
    

    The PR includes a test, but no new docs yet. It also includes some new utility functions in fmpz_mat_extras.

    opened by argriffing 15
  • improve upper incomplete gamma algorithm selection for nonnegative real s and z

    improve upper incomplete gamma algorithm selection for nonnegative real s and z

    When s and z are nonnegative real, the asymptotic expansion is used when (z-c)^8 > (a*s)^8 + (b*prec)^8 where a=1.029287542, b=0.3319411658, c=2.391097143. For reference the acb_hypgeom_u_use_asymp condition has a=0, b=0.69314718055994530942, c=0 in this notation.

    Closes https://github.com/fredrik-johansson/arb/issues/276. Closes https://github.com/fredrik-johansson/arb/issues/166.

    The magic constants were found by locating the best transition point for each point in the grid prec=16..400 s=0.5..200.5 then contour plotting the results, guessing a nonlinear model, and fitting the constants. The mean squared error for the predicted x transition point as a function of s and prec is about 1. The model looks solid enough to me that extrapolating beyond the training grid seems OK.

    Machine floats are used for speed, and I'm not sure if overflow would be a problem. This PR only deals with nonnegative real s and z, so there are still algorithm selection problems when s and z have imaginary parts. I'm not totally sure what's going on with integer values of s. Surprisingly I haven't noticed any actual transition region issues, in other words I haven't noticed a point where neither the 1f1 nor the asymp algorithms were good enough, but I haven't looked for this specifically.

    opened by p15-git-acc 14
  • Root finding for degenerate polynomials

    Root finding for degenerate polynomials

    Let me first say that I was very happy to find the capabilities of Arb within Nemo and Oscar!

    I have one application in mind where I need to track roots of a univariate polynomial f(u) for varying coefficients a_i of f. It lies in the nature of my problem that the limit polynomial f -> f_1 has at least one multiple root. I learned from the Arb documentation that internally, the Durand-Kerner method is used. According to this article, it is possible to obtain reasonable convergence also in this degenerate case.

    Playing around with the Arb methods, I found that the error bounds for the returned list of roots of f becomes significantly worse, once f approaches f_1 which has a multiple root. Usually, one would circumvent this problem by making f square free. But in my setting, it will be most convenient to have the input in inexact form, i.e. f_1 is not an exact polynomial with a multiple root, but rather an Arb-polynomial with an 'approximate multiple root'. Therefore, doing exact division beforehand is not a feasible way to go.

    Question: Can we modify the implementation of the root-finding algorithm in such a way that it returns a list of roots with higher precision in the degenerate cases? It already returns a list of roots, making no claims about their multiplicities. And this is fine for me; it's only that the precision really seems to go down.

    The following is a sample code in julia using Oscar, derived from @jhanselman's work.

    using Oscar
    
    function recursive_continuation(f::acb_poly, g::acb_poly, r::Vector{acb})
      P = parent(f)
      P === parent(g) || error("polynomials must belong to the same parent ring")
      z = gens(P)[1]
      CC = coefficient_ring(P)
      prec = precision(CC)
      n = degree(f)
      degree(g) == n || error("number of coefficients must be equal to the degree of the polynomial")
      length(r) == n || error("number of approximate roots must be equal to the degree of the polynomial")
      temp_vec = Nemo.acb_vec(n)
      temp_vec_res = Nemo.acb_vec(n)
    
      # The following is an error estimation from Neurohr's thesis.
      d = reduce(min, [abs(r[i] - r[j]) for (i, j) in filter(t-> t[1] != t[2], [a for a in Iterators.product((1:n), (1:n))])]) 
      W = [ f(r[i]) // prod([r[i] - r[j] for j in vcat((1:i - 1), i+1:n)];init = one(CC)) for i in (1:n)]
      w = reduce(max, map(t -> real(t)^2 + imag(t)^2, W))
      
      if w < d // (2*n)
        temp_vec = Nemo.acb_vec(copy(r)) # TODO: Why does fill! not work?
        dd = ccall((:acb_poly_find_roots, Nemo.libarb), Cint, (Ptr{Nemo.acb_struct}, Ref{acb_poly}, Ptr{Nemo.acb_struct}, Int, Int), temp_vec_res, g, temp_vec, 0, prec)
    
        dd == n || println("number of roots has dropped to $dd")
        result = similar(r)
        result .= Nemo.array(CC, temp_vec_res, n)
        
        Nemo.acb_vec_clear(temp_vec, n)
        Nemo.acb_vec_clear(temp_vec_res, n)
        
        return result
      else
        h = CC(1//2).*(f+g)
        return recursive_continuation(h, g, recursive_continuation(f, h, r))
      end
    end
    
    prec = 128
    CC = ComplexField(prec)
    
    d = 5
    P, a = PolynomialRing(QQ, push!(["a$i" for i in 0:d-1], "z"))
    z = pop!(a)
    F = sum([a[i]*z^(i-1) for i in 1:d]) + z^d
    
    v = [CC(i) for i in 1:d]
    
    CCu, u = CC["u"]
    
    f = evaluate(F, push!(CCu.(v), u))
    
    r = roots(f)
    
    w = CC.([20, -10, 5, 25, -3])
    g = u^d + sum([w[i]*u^(i-1) for i in 1:d])
    @show r2 = recursive_continuation(f, g, r) # This produces good results!
    g = u^d - u^(d-1) 
    @show recursive_continuation(f, g, r) # This produces worse results.
    
    opened by HechtiDerLachs 1
  • arb_pow could do much better for <a +- a> ^ <b +- c> with b >= c > 0 and a > 0

    arb_pow could do much better for ^ with b >= c > 0 and a > 0

    Consider arb_pow(z, x, y, prec) with x = <a +- a> and y nonnegative. This works as I would expect for x=0, and also for exact half integers y, including y = 0. But if a > 0 and y is not an exact half integer, then arb_pow returns a non-finite answer (because it takes the logarithm of x, then multiplies by y, then takes the exponential; but we can't represent the half-infinite interval that the logarithm requires), whereas there is a finite answer. I propose to do roughly this, after all existing special cases have been done:

        if(arf_cmp_si(arb_midref(x), 0) > 0
            && arf_cmpabs_mag(arb_midref(x), arb_radref(x)) == 0
            && arb_is_nonnegative(y)) {
            /* x is an interval of the form <a +- a>, y is an interval of the form <b +- c> with b >=
              c. 
    
              (x, y) -> x^y is nondecreasing in x for y >= 0, and it is 0 for x=0 and y>0.  The lower
              bound of the interval for x is 0, and we have points with y>0 (because the case
              arb_is_zero(y) is excluded above). So the lowerbound of the result is 0, and the
              upperbound can be found somewhere at (2*a)^y.
    
              We compute this upperbound by computing z^y with z = < 3*a/2 +- a/2 >, which has the same
              upperbound (2*a) as x, then keeping the upper bound of this result and setting the lower
              bound to 0. This necessarily drops the precision to MAG_BITS, so we might as well do that
              from the start. */
            prec = (prec > MAG_BITS) ? MAG_BITS : prec;
            arf_mul_ui(arb_midref(z), arb_midref(x), 3, prec, ARF_RND_UP);
            arf_mul_2exp_si(arb_midref(z), arb_midref(z), -1);
            mag_mul_2exp_si(arb_radref(z), arb_radref(x), -1);
    
            arb_pow(z, z, y);
    
            /* Now we keep the upper bound of z (we may need to round up) and set the lower bound to
               zero. That is, if currently, z = < zc +/- zr >, we set it to < (zc + zr)/2 +/- (zc +
               zr)/2 >. */
            arb_get_ubound_arf(arb_midref(z), z, prec);
            arf_mul_2exp_si(arb_midref(z), arb_midref(z), -1);
            arf_get_mag(arb_radref(z), arb_midref(z));
            /* Note the above is inexact (rounding up), so need to update arb_midref(z) to match again */
            arf_set_mag(arb_midref(z), arb_radref(z));
        }
    
    opened by postmath 5
  • Spheroidal wave functions

    Spheroidal wave functions

    Arb ought to have implementations of the spheroidal wave functions, which have their own chapter in DLMF: https://dlmf.nist.gov/30

    An interesting application would be to check Alain Connes's recent work on the Riemann hypothesis which uses prolate spheroidal wave functions.

    opened by fredrik-johansson 0
  • Feature request: `arf_set_str`

    Feature request: `arf_set_str`

    Since there is a function char * arf_get_str(const arf_t x, slong prec), it will be good to also define a function with the effect of

    inline int arf_set_str(arf_t res, const char * inp, slong prec)
    {
        arb_t y; arb_init(y);
        if (arb_set_str(y, inp, prec))
            return 1;
        arf_set(res, &y->mid);
        return 0;
    }
    
    opened by MacroUniverse 0
  • Entire (generalised) Bernoulli function

    Entire (generalised) Bernoulli function

    This PR is like sympy/sympy#23984, but for Arb. It implements the generalised Bernoulli function from @PeterLuschny's "An introduction to the Bernoulli function":

    $$B(s,a)=\begin{cases}1&s=0\\-s\zeta(1-s,a)&s\ne0\end{cases}$$

    as well as the ordinary version $B(s)=B(s,1)$. These functions, like the Riemann/Hurwitz zeta functions, come in arb and acb versions too.

    Note that the integer-only Bernoulli number functions (arb_bernoulli_ui, etc.) are not touched, and that this PR is inapplicable to FLINT since it has no Hurwitz zeta function.

    It might be possible to improve accuracy near the zeta pole (i.e. the origin, $s=0$) by directly using the Maclaurin series expansion of the Bernoulli function (which is entire for fixed $a$), but this is not implemented here.


    Example of usage:

    #include "arb.h"
    #include "acb.h"
    #include <stdio.h>
    
    int main()
    {
        arb_t s, a, x;
        arb_init(s);
        arb_init(a);
        arb_init(x);
        arb_set_d(s, 0.9);
        arb_set_d(a, 0.7);
        arb_bernoulli(x, s, 100);
        arb_printn(x, 50, 0);
        printf("\n");
        arb_bernoulli_gen(x, s, a, 100);
        arb_printn(x, 50, 0);
        printf("\n");
        arb_clear(s);
        arb_clear(a);
        arb_clear(x);
    
        acb_t s2, a2, x2;
        acb_init(s2);
        acb_init(a2);
        acb_init(x2);
        arb_set_d(acb_realref(s2), 0.4);
        arb_set_d(acb_imagref(s2), -1.2);
        arb_set_d(acb_realref(a2), -0.7);
        arb_set_d(acb_imagref(a2), 0.3);
        acb_bernoulli(x2, s2, 100);
        acb_printn(x2, 50, 0);
        printf("\n");
        acb_bernoulli_gen(x2, s2, a2, 100);
        acb_printn(x2, 50, 0);
        printf("\n");
        acb_clear(s2);
        acb_clear(a2);
        acb_clear(x2);
    }
    

    Output:

    [0.54273376787061753404885841995 +/- 4.12e-30]
    [0.248587684454566258768404567803 +/- 5.42e-31]
    [0.66727675711496067968005244886 +/- 2.95e-30] + [0.627508782028046109043402027631 +/- 8.85e-31]*I
    [30.583734710361609686614823596 +/- 2.22e-28] + [25.559309022539927788601797491 +/- 5.14e-28]*I
    
    opened by Parcly-Taxel 0
Owner
Fredrik Johansson
Computer algebra, high-precision arithmetic
Fredrik Johansson
IntX is a C++11 port of IntX arbitrary precision Integer library with speed, about O(N * log N) multiplication/division algorithms implementation.

IntX IntX is a C++11 port of IntX arbitrary precision Integer library with speed, about O(N * log N) multiplication/division algorithms implementation

Telepati 9 Mar 9, 2022
Arbitrary Precision provides C++ long integer types that behave as basic integer types. This library aims to be intuitive and versatile in usage, rather than fast.

Arbitrary Precision (AP) Cross-platform and cross-standard header-only arbitrary precision arithmetic library. Currently it offers integer types that

null 17 Sep 28, 2022
The ball is the trackball

The ball is the trackball This is Arduino code and 3D-printable models for a Bluetooth trackball in which all the electronics are inside the ball.

Jacek Fedoryński 129 Nov 11, 2022
keyball is split keyboard has 100% track ball

keyball Keyball is split keyboard has 100% track ball Firmware build guide Keyball46 have separate firmwares for each of PCBs w/ trackball and w/o tra

null 139 Dec 2, 2022
Simple Ball-Based Game, made using SFML.

Ball Ball is a basic windows-only, ball-based game, made using C++ SFML. The goal of the game is to claim 5 balls every 10 seconds. Ball can never be

orlando 2 Dec 20, 2021
A fun game where you don't press the red ball!

DarkBall DarkBall is a fun to play game where you can press little balls/button, but never press the red ball(or any of its friends) You can find/play

Catermelon 3 Jun 25, 2022
A remake of the classic DX-Ball built with iGraphics

DX-Ball Reincarnation This is a remake of the classic DX-Ball games made with C using the iGraphics library and SDL2. To compile from source, run the

null 2 Feb 24, 2022
✖🌱 A DirectX 12 starter repo that you could use to get the ball rolling.

DirectX 12 Seed A DirectX 12 repo you can use to get started with your own renderer. Setup First install: Git CMake Visual Studio Then type the follow

Alain Galvan 71 Nov 16, 2022
MIRACL Cryptographic SDK: Multiprecision Integer and Rational Arithmetic Cryptographic Library is a C software library that is widely regarded by developers as the gold standard open source SDK for elliptic curve cryptography (ECC).

MIRACL What is MIRACL? Multiprecision Integer and Rational Arithmetic Cryptographic Library – the MIRACL Crypto SDK – is a C software library that is

MIRACL 517 Nov 22, 2022
MIRACL Cryptographic SDK: Multiprecision Integer and Rational Arithmetic Cryptographic Library is a C software library that is widely regarded by developers as the gold standard open source SDK for elliptic curve cryptography (ECC).

MIRACL What is MIRACL? Multiprecision Integer and Rational Arithmetic Cryptographic Library – the MIRACL Crypto SDK – is a C software library that is

MIRACL 517 Nov 22, 2022
LibTomMath is a free open source portable number theoretic multiple-precision integer library written entirely in C.

libtommath This is the git repository for LibTomMath, a free open source portable number theoretic multiple-precision integer (MPI) library written en

libtom 540 Nov 23, 2022
Intel:registered: Homomorphic Encryption Acceleration Library accelerates modular arithmetic operations used in homomorphic encryption

Intel Homomorphic Encryption Acceleration Library (HEXL) Intel ®️ HEXL is an open-source library which provides efficient implementations of integer a

Intel Corporation 161 Nov 14, 2022
A single-header C/C++ library for parsing and evaluation of arithmetic expressions

ceval A C/C++ header for parsing and evaluation of arithmetic expressions. [README file is almost identical to that of the ceval library] Functions ac

e_t 9 Oct 10, 2022
A single-header C/C++ library for parsing and evaluation of arithmetic expressions

ceval A C/C++ header for parsing and evaluation of arithmetic expressions. [README file is almost identical to that of the ceval library] Functions ac

e_t 9 Oct 10, 2022
C++ Matrix -- High performance and accurate (e.g. edge cases) matrix math library with expression template arithmetic operators

Matrix This is a math and arithmetic matrix library. It has stood many years of performing in mission critical production for financial systems. It ha

Hossein Moein 71 Oct 29, 2022
A C/C++ library for parsing and evaluation of arithmetic expressions.

ceval A C/C++ header for parsing and evaluation of arithmetic expressions. Functions accessibe from main() Function Argument(s) Return Value ceval_res

e_t 6 Nov 8, 2022
A simple C++ 03/11/etc timer class for ~microsecond-precision cross-platform benchmarking. The implementation is as limited and as simple as possible to create the lowest amount of overhead.

plf_nanotimer A simple C++ 03/11/etc timer class for ~microsecond-precision cross-platform benchmarking. The implementation is as limited and simple a

Matt Bentley 100 Nov 16, 2022
An updated, 1.2.1 version of Equalizer APO running with internal double precision processing (64 bit)

EqualizerAPO - 64bit port This repo contains an updated, 1.2.1 64-bit port of EqualizerAPO - system wide equalizer for Windows. The port here is inspi

FireKahuna 42 Nov 28, 2022
Extracts high-precision mouse/pointer motion data on Windows. Good for drawing software!

window_mouse_queue This is a wrapper for GetMouseMovePointsEx function that allows to extract high-precision mouse/pointer motion data on Windows. Goo

YellowAfterlife's GameMaker Things 6 Feb 21, 2022
Calculator that suffers from floating point precision

calc A calculator that suffers from floating precision. This calculator suffers from floating point precision and isn't much more than a fun project.

atiedebee 1 Nov 28, 2021