Giter Club home page Giter Club logo

arb's Introduction

See FLINT

2023: Arb has been merged into FLINT. The present repository is archived and will no longer be updated. See

https://github.com/flintlib/flint/

for new developments!

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 (https://github.com/Nemocas/Nemo.jl/) 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:

arb's People

Contributors

ahrvoje avatar akobel avatar albinahlback avatar alexbrc avatar arbguest avatar david-berghaus avatar davidlowryduda avatar fandango96 avatar fredrik-johansson avatar hannorein avatar isuruf avatar jdemeyer avatar joel-dahne avatar jpflori avatar jwbober avatar kiwifb avatar mezzarobba avatar mkoeppe avatar orlitzky avatar oscardssmith avatar p15-git-acc avatar pascalmolin avatar postmath avatar rickyefarr avatar rwst avatar saraedum avatar texnokrates avatar thofma avatar tthsqe12 avatar wbhart avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

arb's Issues

Polygamma functions

Extend Stirling series code to compute polygamma functions (i.e. starting the series from some derivative), and optimize for a small number of derivatives by using a direct recurrence instead of binary splitting.

hilbert_class_poly spam in tests

Running make check gives many lines that look like test failures, but which I think are more like debugging printfs in non-test code.

hilbert_class_poly failed at 8837953 bits of precision

These hilbert_class_poly tests are taking a while, but I haven't tried to time them. I'm running the development branch of arb with the lib{flint, gmp, mpfr}-dev that are packaged for Ubuntu 15.04.

log1mexp -- adding yet another function to the wish list

log1mexp(x) = log(1 - exp(-x))
According to this CRAN note it should be evaluated using either log1p or expm1 depending on some threshold value of x. Maybe I will add this. Just an editorial comment, the name is a little funky because you might think it would compute log(1 - exp(x)).

Output radius smaller than radius of input

Consider the following code snippet:


include "acb.h"

void main()
{
acb_t s,z;
acb_init(s);
acb_init(z);

arb_const_pi(acb_realref(s), 1000);
arb_const_pi(acb_imagref(s), 1000);

acb_zeta(z,s,1000);
acb_printd(z,20);
printf("\n");

acb_clear(s);
acb_clear(z);   

}


The output above will give:
(0.9096902551448360325 - 0.065857401879019999596j) +/- (9.94e-302, 9.82e-302j)

Shouldn't the radius of output always be larger (or equal to) the radius of the approximate machine precision midpoint?

Setting the value of z to be 3+5*i, will result in
(0.91252658899897131011 + 0.050842871074571362072j) +/- (7.88e-305, 1.81e-305j)
which is more drastic.

Maybe I'm wrong, but it appears as though the exponent is moving in the opposite direction.

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?

Further improve computation of roots of unity

Chebyshev polynomials (with fast evaluation) may be faster than complex exponentiation. Also maybe use a lookup table of 1/q values for small q and do binary exponentiation to get to p/q.

Implement `mag_get_mpfr`

Currently, there is a function arf_get_mpfr, but in order to obtain an mpfr from a mag, it seems to be necessary to use mag_get_fmpr and fmpr_get_mpfr.

It would be convenient to be able to do this in one step by a function mag_get_mpfr, say.

Compute polylogarithms efficiently in special cases

Rational functions should be used when s is a negative integer. The series expansion in log(z) should be used when z is close to 1. Functional equations should be used for the dilogarithm. There are other special cases. See Crandall's paper.

Faster sine and cosine of power series

For the sine and cosine of power series, a single exponential (with inverse) should be faster when the order gets large enough, at least when the expansion point actually is complex (a case distinction might be needed). Requires benchmarking.

Comparison methods

Add methods arb_eq, arb_ne, arb_lt, arb_le, arb_gt, arb_ge implementing mathematically meaningful comparisons (interpreted on the extended real line)

(X op Y)

is true if

(x op y) holds for all x in X, for all y in Y

where presence of NaN is handled correctly (causing the result to be false).

Make Check Fails Because libarb.so.0 is not found

With the newest changes to the configure file, the shared object file libarb.so.0.0.0 is created. When running the tests, it results in failure due to "libarb.so.0" not found.

Also, when using the --prefix= flag with ./configure, the include files are placed in the ./include file (as expected) but the library files are now placed above the lib directory.

Improve Documentation of _acb_poly_zeta_em_choose_param

Currently the documentation for this function is (taken from website):

void _acb_poly_zeta_em_choose_param(mag_t bound, ulong * N, ulong * M, const acb_t s, const acb_t a, long d, long target, long prec)

Chooses N and M for Euler-Maclaurin summation of the Hurwitz zeta function, using a default algorithm.

Most of the parameters are self explanatory but the parameters d and target are not explained.

./configure ignores LDFLAGS

We need to pass the proper rpath options to the linker, for example by:

export LDFLAGS="-L${FLINT_DIR}/lib -Wl,-rpath=${FLINT_DIR}/lib -L${GMP_DIR}/lib -Wl,-rpath=${GMP_DIR}/lib -L${MPFR_DIR}/lib -Wl,-rpath=${MPFR_DIR}/lib"

But the LDFLAGS variable gets ignored, unlike autotools. So we then need to patch the final .so library manually with patchelf. See hashdist/hashstack#232.

The solution is to allow the user to specify LDFLAGS, just like you allow the user to specify CFLAGS.

This bug applies to Flint as well (flintlib/flint#66).

Faster argument reduction for log gamma

Write a faster logarithmic rising factorial (with correct branch cuts) for reducing the complex log gamma function. Also implement the logarithmic reflection formula.

Make root-finding more robust

In _fmpcb_poly_refine_roots_durand_kerner / _acb_poly_refine_roots_durand_kerner, it's possible for a divide by zero to occur. In this case the root-finding never converges (possibly leading to a segmentation fault when the precision is increased many times).

A polynomial that triggers this is:

0 604462842248924277768128 -133175711609 1

A workaround is to remove the zero root.

It's not clear whether this problem can occur for polynomials without trailing zero coefficients (in principle, most likely, not sure about practice). At the very least, the code should be made more robust: avoid the division by zero, and perhaps just set the root approximation to a random value in that case.

_arb_vec_norm seems to compute the squared norm

This can't be right can it? It's like this in the implementation and in the docs. Should the function name be changed to _arb_vec_squared_norm or maybe a square root should be taken? I couldn't find where this function is actually called.

poly_roots returns low accuracy

I tried the poly_roots example, reading the docs first:

If -refine d is passed, the roots are refined to an absolute
tolerance better than 10^(-d).

and using it:

$ ./poly_roots -print 30 p 50
prec=53: 0 isolated roots | cpu/wall(s): 0.08 0.085
prec=106: 50 isolated roots | cpu/wall(s): 0.06 0.059
done!
(-0.998866404420071052872997501108 + -3.09168790452712139074186521034e-18j)  +/-  (3.46e-15, 3.46e-15j)
(-0.994031969432090699513049676088 + 4.55853633110982516428872887659e-18j)  +/-  (1.06e-14, 1.06e-14j)
(-0.985354084048005917634215112514 + -8.60630010071774494920947139187e-18j)  +/-  (1.52e-14, 1.52e-14j)
(-0.972864385106692071159319738104 + 3.38386716312239390021536623006e-17j)  +/-  (1.57e-14, 1.57e-14j)
(-0.956610955242807949996900440343 + 5.45523778230809524962583938682e-18j)  +/-  (1.17e-14, 1.17e-14j)
(-0.936656618944877934533226198517 + 2.09284038499267596686478631427e-18j)  +/-  (7.86e-15, 7.86e-15j)
(-0.91307855665579189245635086212 + -2.46419863027644378825880794859e-18j)  +/-  (4.02e-15, 4.02e-15j)
(-0.885967979523613050912782675287 + 1.28655456602927425464451102937e-18j)  +/-  (1.76e-15, 1.76e-15j)
(-0.855429769429946084292451608376 + -1.37030378937902583483317169162e-18j)  +/-  (6.37e-16, 6.37e-16j)
(-0.821582070859335948219539209478 + 5.87437576856758209097429539944e-21j)  +/-  (2.08e-16, 2.08e-16j)
(-0.784555832900399263917567894563 + -6.18904593730578494178647537457e-20j)  +/-  (5.33e-17, 5.33e-17j)
(-0.744494302226068538281710950413 + -4.2589207006935182581870433993e-21j)  +/-  (1.21e-17, 1.21e-17j)
(-0.701552468706822251092289919753 + 8.22040063472779484038226475286e-22j)  +/-  (2.45e-18, 2.45e-18j)
(-0.655896465685439360780709128553 + 9.25999490785417670856128676698e-23j)  +/-  (3.89e-19, 3.89e-19j)
(-0.607702927184950239180429248324 + 6.81185141768706493633810771142e-23j)  +/-  (5.26e-20, 5.26e-20j)
(-0.557158304514650054315513239705 + 1.24569451988933598598007685593e-23j)  +/-  (6.69e-21, 6.69e-21j)
(-0.504458144907464201651459692139 + -4.1527094009297494785120667117e-26j)  +/-  (6.8e-22, 6.8e-22j)
(-0.449806334974038789147131398609 + -2.85438003288030529714473117607e-26j)  +/-  (6.04e-23, 6.04e-23j)
(-0.393414311897565127394229248133 + -3.66480380945852035019023630225e-27j)  +/-  (4.36e-24, 4.36e-24j)
(-0.335500245419437356836988257522 + 3.03389812654565888179360934898e-28j)  +/-  (3.04e-25, 3.04e-25j)
(-0.276288193779531990327645278528 + 1.02653955429918324828636981052e-29j)  +/-  (1.72e-26, 1.72e-26j)
(-0.216007236876041756847284532617 + -1.26865142760366896330150055325e-30j)  +/-  (9.86e-28, 9.86e-28j)
(-0.154890589998145902071628620941 + 3.26802839025889866667840408734e-32j)  +/-  (4.43e-29, 4.43e-29j)
(-0.0931747015600861408544503776396 + -2.78785302260664264672814381419e-34j)  +/-  (2.39e-30, 2.39e-30j)
(-0.0310983383271888761123289896659 + -3.1897750430736787953652257393e-34j)  +/-  (8.81e-32, 8.81e-32j)
(0.0310983383271888761123289896659 + 2.8009970016647681211252053531e-34j)  +/-  (8.78e-32, 8.78e-32j)
(0.0931747015600861408544503776396 + 1.22221190190830635856314273573e-33j)  +/-  (2.07e-30, 2.07e-30j)
(0.154890589998145902071628620941 + 2.98543075705758790343667925774e-32j)  +/-  (4.41e-29, 4.41e-29j)
(0.216007236876041756847284532619 + 9.45458090435984279846572806304e-31j)  +/-  (9.56e-28, 9.56e-28j)
(0.27628819377953199032764527854 + -1.35372986452557078080916816695e-29j)  +/-  (1.65e-26, 1.65e-26j)
(0.335500245419437356836988257408 + -5.60875091397265461462228546794e-28j)  +/-  (2.83e-25, 2.83e-25j)
(0.393414311897565127394229253316 + -2.31280383510649099901256794077e-27j)  +/-  (4.32e-24, 4.32e-24j)
(0.449806334974038789147131420251 + 3.04539523904888183844687913822e-27j)  +/-  (6.01e-23, 6.01e-23j)
(0.50445814490746420165145991668 + -8.26340692631503763934751187039e-25j)  +/-  (6.53e-22, 6.53e-22j)
(0.55715830451465005431551349539 + -2.20782064902516410343299466383e-25j)  +/-  (6.39e-21, 6.39e-21j)
(0.607702927184950239180375059144 + 3.52369474412386921582747697242e-23j)  +/-  (5.03e-20, 5.03e-20j)
(0.655896465685439360780829533945 + 3.12006908837958824279213297135e-23j)  +/-  (3.69e-19, 3.69e-19j)
(0.701552468706822251088434086683 + 5.27955271960739429466084120095e-21j)  +/-  (2.7e-18, 2.7e-18j)
(0.744494302226068538249512965338 + 4.12071007284606577164103053378e-20j)  +/-  (1.26e-17, 1.26e-17j)
(0.784555832900399263850757760055 + 1.13208297672283036312605253738e-19j)  +/-  (5.33e-17, 5.33e-17j)
(0.821582070859335948057064280635 + 3.33452537239634780370291429568e-20j)  +/-  (2.17e-16, 2.17e-16j)
(0.855429769429946085084070052295 + -4.95826764776914415213743709223e-19j)  +/-  (6.22e-16, 6.22e-16j)
(0.885967979523613050233022377457 + 1.59211931742218982489995745465e-18j)  +/-  (1.76e-15, 1.76e-15j)
(0.913078556655791898217434330239 + -1.78935719993836398214582758251e-18j)  +/-  (4.1e-15, 4.1e-15j)
(0.936656618944877930419888295104 + 1.20442224701355206872853432221e-17j)  +/-  (7.78e-15, 7.78e-15j)
(0.956610955242807941189670550424 + -2.208836970104941730490344241e-17j)  +/-  (1.29e-14, 1.29e-14j)
(0.972864385106692059568175324642 + -6.48091560049456287693354433304e-18j)  +/-  (1.49e-14, 1.49e-14j)
(0.985354084048005884919504220873 + 2.39614708155537040712644619796e-17j)  +/-  (1.49e-14, 1.49e-14j)
(0.994031969432090687304406847174 + 3.11204286556621781340063626733e-18j)  +/-  (9.97e-15, 9.97e-15j)
(0.99886640442007105147887935094 + 1.52915873310933633632672092345e-17j)  +/-  (3.49e-15, 3.49e-15j)
cpu/wall(s): 0.14 0.145

As you can see, I requested 30 digits, but I only got ~1e-15 accuracy for most of the roots. So the documentation is not consistent with the results, or I misunderstood something. I am using 7c6b3fb.

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?

Exact Rational Values

If one desires to set arb_t variable to be an integer, a precision is not required to do so. On the other hand, in order for a user to set an arb_t variable to be rational, one must specify a precision.

Would it be better to allow rational numbers to be stored exactly?

dev guide?

I'm wondering how to set up an arb development environment, especially with respect to building and testing? I've made a couple of PRs already but with a pretty inefficient workflow that involves rebuilding and retesting everything all the time.

Below I'm pasting the kind of script I'm starting to write to work around this. Is there a better way? If not, that is OK.

from __future__ import print_function, division

import subprocess
import argparse
import os

def main(subdirectory):

    # Find the full path to the subdirectory.
    cwd = os.getcwd()
    p = os.path.join(cwd, 'build', subdirectory, 'test')

    # Identify files in the test directory ending with ".d".
    # These are dependency files which correspond to the executable tests.
    testfiles = []
    for filename in os.listdir(p):
        base, extension = os.path.splitext(filename)
        if extension == '.d':
            testfiles.append(base)
    if testfiles:
        print('found', len(testfiles), 'test scripts')
        # Define the environment variables for the test scripts.
        scriptenv = os.environ.copy()
        scriptenv['LD_LIBRARY_PATH'] = cwd
        for testfile in testfiles:
            name = os.path.join(p, testfile)
            proc = subprocess.Popen([name], env=scriptenv)
            proc.wait()
    else:
        print('did not find any test scripts')

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('--subdirectory', '-s', required=True,
            help='a subdirectory to test (e.g. acb_mat)')
    args = parser.parse_args()
    main(args.subdirectory)

Document and clean up fmpz_extras

The fmpz_extras module is completely undocumented. There are also some functions in mag.h that should be moved to fmpz_extras. Some function renaming is probably in order. Ultimately, many of the fmpz_extras functions should be moved to flint.

Add incomplete gamma functions

Gamma(s,z)
gamma(s,z)
gamma*(s,z) (without the z^s factor)

Series expansions with respect to z.

Maybe also regularized incomplete gamma functions.

Gauss-Legendre Integration

I have some code that generates Gauss-Legendre quadrature rules. As of right now, the code returns a struct that contains the length, a node vector, and a weight vector. Is this something that would be useful to be added into the library? If so, any suggestions for changes would be appreciated.

The file is located at the following link:

https://github.com/Therickaman/arb/blob/77bacde5db14d2e9d97c77658b32a480415f071b/gaussLegendre.c

The code is pretty rough, but needs discussion.

Should the radius be set to zero for the weights and nodes, so that they are pretty much treated as exact? Should the return structure be changed? To integrate a user defined function f using arb what should the format be when passed to functions that integrate?

/usr/local/include/nmod_vec.h:37:26: fatal error: ulong_extras.h

Hi, I have Ubuntu 14.04 LTS. After installation : https://pl.wikibooks.org/wiki/Programowanie_w_systemie_UNIX/ARB#Instalacja
I have tried to compile a program from http://fredrikj.net/arb/setup.html

/*
gcc -I/usr/local/include/flint/ -L/usr/local/include/flint/ -larb c.c
gcc -larb c.c
gcc -lflint -lgmp -lmpfr -larb -I/usr/include/flint c.c
gcc -I/usr/local/include -I/usr/local/include/flint/ -L/usr/local/lib -larb c.c
gcc -lflint -lgmp -lmpfr -larb -I/usr/local/include -I/usr/local/include/flint/ -L/usr/local/lib c.c
*/

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);
}

Nothing works :

gcc -larb c.c
In file included from /usr/local/include/fmpz.h:45:0,
from /usr/local/include/fmpr.h:37,
from /usr/local/include/fmprb.h:29,
from /usr/local/include/arb.h:36,
from c.c:16:
/usr/local/include/nmod_vec.h:37:26: fatal error: ulong_extras.h: Nie ma takiego pliku ani katalogu
#include "ulong_extras.h"
^
compilation terminated.
a@zalman:/c/varia/arb/1$ gcc -I/usr/local/include -L/usr/local/lib -larb c.c
In file included from /usr/local/include/fmpz.h:45:0,
from /usr/local/include/fmpr.h:37,
from /usr/local/include/fmprb.h:29,
from /usr/local/include/arb.h:36,
from c.c:16:
/usr/local/include/nmod_vec.h:37:26: fatal error: ulong_extras.h: Nie ma takiego pliku ani katalogu
#include "ulong_extras.h"
^
compilation terminated.
a@zalman:
/c/varia/arb/1$ locate ulong_extras.h
/home/a/flint2/ulong_extras.h
/usr/local/include/flint/ulong_extras.h
a@zalman:~/c/varia/arb/1$ gcc -I/usr/local/include -I/usr/local/include/flint/ -L/usr/local/lib -larb c.c
/tmp/cc9MgIvo.o: In function _fmpz_demote': c.c:(.text+0x28): undefined reference to_fmpz_clear_mpz'
/tmp/cc9MgIvo.o: In function arf_clear': c.c:(.text+0xfc): undefined reference to_arf_demote'
/tmp/cc9MgIvo.o: In function arb_printn': c.c:(.text+0x181): undefined reference toarb_fprintn'
/tmp/cc9MgIvo.o: In function main': c.c:(.text+0x1a8): undefined reference toarb_const_pi'
c.c:(.text+0x1cd): undefined reference to flint_printf' c.c:(.text+0x1d4): undefined reference toarb_version'
c.c:(.text+0x1e6): undefined reference to flint_printf' collect2: error: ld returned 1 exit status a@zalman:~/c/varia/arb/1$ gcc -lflint -lgmp -lmpfr -larb -I/usr/local/include -I/usr/local/include/flint/ -L/usr/local/lib c.c /tmp/ccdsyzvn.o: In function_fmpz_demote':
c.c:(.text+0x28): undefined reference to _fmpz_clear_mpz' /tmp/ccdsyzvn.o: In functionarf_clear':
c.c:(.text+0xfc): undefined reference to _arf_demote' /tmp/ccdsyzvn.o: In functionarb_printn':
c.c:(.text+0x181): undefined reference to arb_fprintn' /tmp/ccdsyzvn.o: In functionmain':
c.c:(.text+0x1a8): undefined reference to arb_const_pi' c.c:(.text+0x1cd): undefined reference toflint_printf'
c.c:(.text+0x1d4): undefined reference to arb_version' c.c:(.text+0x1e6): undefined reference toflint_printf'
collect2: error: ld returned 1 exit status

But I can compile arb program :

cd arb
a@zalman:/arb$ build/examples/keiper_li 100
zeta: cpu/wall(s): 0.002 0.002
virt/peak/res/peak(MB): 25.09 25.09 5.61 5.61
log: cpu/wall(s): 0 0.001
gamma: cpu/wall(s): 0 0
binomial transform: cpu/wall(s): 0.001 0
0: -0.69314718055994530941723212145817656807550013436029 +/- 6.8425e-49
1: 0.023095708966121033814310247906495291621932127152268 +/- 3.9868e-48
2: 0.046172867614023335192864243096033943387066108314182 +/- 1.2738e-47
3: 0.069212973518108267930497348872601068994212026394777 +/- 4.371e-47
4: 0.092197619873060409647627872409439018065541673496812 +/- 1.4868e-46
5: 0.11510854289223549048622128109857276671349132305427 +/- 5.0294e-46
6: 0.13792766871372988290416713700341666356138966082885 +/- 1.7102e-45
7: 0.16063715965299421294040287257385366292282442055073 +/- 5.7197e-45
8: 0.1832194596433825790819393177472185984899809829119 +/- 1.8768e-44
9: 0.20565733870917046170289387421343304741236553444664 +/- 6.1276e-44
10: 0.22793393631931577436930340573684453380748386010511 +/- 2.0045e-43
91: 1.1359793878546436002006561104313275551533628689671 +/- 1.4607e-12
92: 1.1415669830327489959098935595739035602228202255345 +/- 3.3161e-12
93: 1.1471530138314698130651717448646005362180579682481 +/- 7.578e-12
94: 1.1527358323988920619177722932285812811855471805232 +/- 1.7414e-11
95: 1.1583134879564705721273560069201683049927946600968 +/- 4.019e-11
96: 1.1638837494169548403983835123477214984893522087425 +/- 9.3048e-11
97: 1.1694441288315003718368092993985140290769315690077 +/- 2.1584e-10
98: 1.1749919055490018710276285598764024158910189825872 +/- 5.0109e-10
99: 1.1805241509689506316631106075371800417317838090829 +/- 1.1631e-09
100: 1.1860377537679132994976271782020433469152756094002 +/- 2.6968e-09
virt/peak/res/peak(MB): 25.09 25.09 6.03 6.03
a@zalman:
/arb$

Any tips ?
TIA

Adam

Alternative algorithms for Taylor shift

The addition and convolution algorithms for polynomial Taylor shifts should be provided. At least, testing should be done to determine how they compare to divide-and-conquer in terms of speed and numerical stability.

Issues while including bernoulli.h

I've a following simple code which calls the bernoulli_fmprb_ui_zeta function.

#include "fmprb.h"
#include "bernoulli.h"

int main()
{
    fmprb_t y;
    fmprb_init(y);
    bernoulli_fmprb_ui_zeta(y, 9, 50 * 3.33);
    fmprb_printd(y, 50); printf("\n");
    fmprb_clear(y);
}

However, it raises an error

In file included from /usr/local/include/arb.h:31:0,
                 from /usr/local/include/bernoulli.h:36,
                 from check.cpp:2:
/usr/local/include/arf.h: In function ‘void arf_set(arf_struct*, const arf_struct*)’:
/usr/local/include/arf.h:392:13: error: invalid conversion from ‘void*’ to ‘mp_ptr {aka long unsigned int*}’ [-fpermissive]
             ARF_GET_MPN_WRITE(yptr, n, y);

The error occurs because of including bernoulli.h

I'm compiling it using the command:

g++ -I/usr/local/include/flint check.cpp -o check  -larb -lflint;

@fredrik-johansson can you please look into this issue or correct me if I'm doing anything wrong?

Power where the base contains zero

The power 0^c should arguably be finite when Re(c) > 0. The same goes for b^c where b is a ball containing zero, as well as the series expansion of 0^(c+x). This affects acb_poly_powsum_series, among other functions.

This should allow avoiding some indeterminate output for the Hurwitz zeta function, e.g. zeta(-2+x,-3) (but need to verify that such output actually is meaningful).

Add memory manager for arb types

This is a proposal and I'd like to get feedback about this idea. I've noticed that for operations that are repeated removing local variable initialization and clearing, can lead to quite a speed gain. On the other hand, doing this is quite impractical because it leads to unreadable code and very hard to track errors.

I think that a possible solution would be to create a stack for each arb type. Whenever a variable is initialized, if the stack is empty then initialize as normal (perhaps even add a few more initialized variables to the stack). If not, pop a variable off the stack to be used for the variable. Eventually, the stack would be large enough to guarantee local variables are not re-initialized. Especially for polynomial types, this would avoid many calls to the system functions, malloc and free...which causes the kernel to become involved (if I'm not mistaken).

Even if the speed gain is negligible for some types, another added benefit to doing this is the fact that custom data structures would need to be created for each type (stack, queue, linked list, etc.). Having these structures available to the user could really help with implementing more complex functions. Also, having these structures customized for each type avoids unnecessary casts to and from a void pointer type.

As another example, the built in function for derivatives of zeta can re-use already completed computations such as log(n+a), for a fixed and n an non-negative integer.

I just wanted to get feedback to see what everyone thought about this idea. Perhaps I'm incorrect about speed gains.

Polynomial middle product

There should be a mulmid implementation for arb_poly and acb_poly. This can speed up Newton iteration (for division and other functions).

The basecase algorithm is trivial, of course. The more interesting case is the blockwise algorithm. Can the scaling be optimized when computing middle products?

laguerre_l testsuite failure on i386

On this Linux 32-bit system:

Linux arando 3.13.0-40-generic #69-Ubuntu SMP Thu Nov 13 17:56:26 UTC 2014 i686 i686 i686 GNU/Linux

testing arb-2.8.0 gives

laguerre_l....FAIL: consistency

iter = 444

n = (4 + 0j)  +/-  (0, 0j)

m = (-4 + 0j)  +/-  (0, 0j)

z = (-4503462724239360.00001525855623 - 0.499999999999999999999999996769j)  +/-  (7e-61, 6.91e-77j)

res1 = (1.71385883572304972605522896443e+61 + 7.61129353418828371335823987317e+45j)  +/-  (4, 1.78e-15j)

res2 = (1.5222587068376557286128225538e+46 + 5.07029412715334721112563931956e+30j)  +/-  (1.62e-17, 5.65e-33j)

res3 = (nan + 0j)  +/-  (inf, 0j)

s = (nan + 5.07029412715334721112563931956e+30j)  +/-  (inf, 2.03e-31j)

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?

Testsuite fails on ppc64le

On the POWER8 ppc64le system

Linux sardonis 3.19.0-15-generic #15-Ubuntu SMP Thu Apr 16 23:32:13 UTC 2015 ppc64le ppc64le ppc64le GNU/Linux

I get the following reproducible test suite failure with arb 2.7.0:

make[1]: Entering directory '/home/jdemeyer/sage/local/var/tmp/sage/build/arb-2.7.0.p0/src/arf'
[...]
abs_bound_le_2exp_fmpz....PASS
complex_sqr....FAIL!
prec = 698, rnd = 0

a = (6339285465397756945139231368132496554128019744394747208275350942065595714223327282882004073897576134671584866440587073589025503475990026413326511053883607915299395946163061310275939887721284702577049600905500587493134063493826068090023461124398687901692686283320433763252148270486039504843224058218698609943517081716892849130905696941657348763839640087212158083454683180059633412954189769703715978964001161003786946680425738899972518273473719571986450603741721249271434462214845897250532971870021425902922261780901847636017986587369648930793121015247658270597452971859383674118378417538423342960709853716020960438322933655567693176273354358777068473650766777407727302868991 * 2^-2235)

b = (20769147822475905021182172033122559 * 2^-113)

e1 = (1315013909709381559722549835740549443786829800454131997155171828694211243086620368814692247495737405770735789621415903963936315776768232380505812958907618889969868245794264502147656964862275963676046558834458623 * 2^-698)

f1 = (1315031467311928423708944301303640969760880181222735669239969305516422962189234000124364732949570062519188408040675832296657303229505388769788357929947309747087029353794143201921833529433945488716348391784710143 * 2^-696)

e2 = (-323517714603253067973887569535021297144815213232427194808985473711617 * 2^-226)

f2 = (1315031467311928423708944301303640969760880181222735669239969305516422962189234000124364732949570062519188408040675832296657303229505388769788357929947309747087029353794143201921833529433945488716348391784710143 * 2^-696)

r1 = 3, r2 = 3
../Makefile.subdirs:84: recipe for target '../build/arf/test/t-complex_sqr_RUN' failed
make[1]: *** [../build/arf/test/t-complex_sqr_RUN] Aborted
make[1]: Leaving directory '/home/jdemeyer/sage/local/var/tmp/sage/build/arb-2.7.0.p0/src/arf'

The compiler is vanilla GCC 4.9.3. I compiled arb and its dependencies (mpir, mpfr, flint) with -O0 to rule out compiler errors and it did not make a difference.

Running Sage with this version of arb also gives errors in the Sage test suite in functions arb_digamma, arb_const_catalan and arb_polylog_si.

Docs for fmprb.h uclear (to me)

I am reading the docs (http://fredrikj.net/arb/fmprb.html):

The precision parameter passed to each function roughly indicates the precision to which calculations on the midpoint are carried out (operations on the radius are always done using a fixed, small precision.)

Here is an example of such a function that accepts prec:

void fmprb_div(fmprb_t z, const fmprb_t x, const fmprb_t y, long prec)

So the x and y ball floating point numbers each contain some precision, e.g.
x = 1 \pm 3 and y = 2 \pm 0.01. What prec should I, as a user, set to prec if I want to calculate z = x/y? It says that the result, i.e. z, will be rounded to prec bits. So if I set prec too low, then I'll lose accuracy. What if I set prec too high? If I want Arb to do its best and return me the answer as accurate as possible (given the algorithms that it implements), what should I set for prec? Will it return an error if prec is too high? Or will it just return z with lower accuracy if I set prec too high?

It would be nice to update the documentation which answers all these questions.

fprint

Arb can print but not fprint its objects. If this is added, should foo_print_bar(...) call foo_fprint_bar(stdout, ...) or should the print/fprint functions be implemented separately with code duplication?

[doc] What is the precision of the midpoint after arb_t operations?

The documentation of arb_t doesn't mention what the precision (number of bits) of the midpoint of a ball is after an operation. It says that the internal working precision is not well-defined, but I think at least the output precision should be well-defined.

I guess that after calling some function arb_foo(arb_t x, ..., long prec), the midpoint of x will be represented as prec bits. Is this guess right? If so, it should be documented.

Improve string conversion

String conversion code should be done without relying on MPFR (to support big exponents, etc.). Would also be nice to have some way to convert a number to a decimal expansion in which all digits are guaranteed correct (or a decimal interval), and pretty-printing such as 3.14159? instead of 3.1415953283490 +/- 1e-6.

Unused parameter warning

I tend to use as many diagnostic/warning flags as possible when compiling C/C++ code, and I noticed this warning while compiling code using Arb:

/home/yardbird/repos/arb/fmpr.h:410:46: warning: unused parameter ‘yexp’ [-Wunused-parameter]

Would it be possible to silence this warning? In C++ it is possible to remove the name yexp from the parameter list of the function - so that instead of

const fmpz_t yexp

there is only

const fmpz_t

But I am not sure this is legal in C. Otherwise a cast to void of yexp in the body of the function would suffice.

How to best expose Arb numbers in C++

We want to have good support for Arb in CSymPy, essentially allow users to use Arb to evaluate expressions, and it will give the result with a given accuracy.

What is the best way to expose the Arb numbers in C++? Should we create some lightweight class, similar to mpz_class in GMP (and that should be part of Arb itself, just like mpz_class is part of GMP itself)? Or should we just create a more heavyweight class like ArbFloat, which would internally have a pointer to the Arb C datastructure, and then allow some operations on that (this class would only be part of CSymPy)?

Add missing trigonometric and hyperbolic functions

The following functions should be available for arb, acb, arb_poly, acb_poly: sin, cos, tan, sinh, cosh, tanh, asin, acos, atan, asinh, acosh, atanh

Extra functions: sin_cos, sinh_cosh, sin_pi, cos_pi, sin_cos_pi

Maybe also: sec, csc, cot, sech, csch, coth, asec, acsc, acot, asech, acsch, acoth

Other functions that could be useful: sinc, sinc_pi, cosm1

Makefile problems on OS X 10.10

I've tried to compile Arb (standalone) on OS X 10.10.2,
compiler Apple LLVM version 6.0 (clang-600.0.56) (based on LLVM 3.5svn) or gcc-4.9.

Everything seems fine until the compilation ends and this line of the makefile is executed: https://github.com/fredrik-johansson/arb/blob/master/Makefile.in#L111

The problem seems to be that $(OBJS) is empty, and in fact if you remove the line
(e.g sed -i.breworig 's/$(QUIET_AR) $(AR) rcs libarb.a $(OBJS);//' Makefile)
everything goes through.

As expectable make check fails:

test/t-epsilon_arg.c: In function 'acb_modular_epsilon_arg_naive':
test/t-epsilon_arg.c:52:9: warning: implicit declaration of function 'fmpq_dedekind_sum' [-Wimplicit-function-declaration]
         fmpq_dedekind_sum(t, d, c);
         ^
Undefined symbols for architecture x86_64:
  "_fmpq_dedekind_sum", referenced from:
      _main in ccR5NFtc.o
ld: symbol(s) not found for architecture x86_64

Any idea how to fix it?

add more acb-related 'set' and 'contains' convenience functions

Currently there's a function called acb_set_fmpz_fmpz that sets the real part and the imaginary part of a complex ball to the provided integers. There are also functions like acb_contains_fmpz.

Here are some similarly named but hypothetical functions that I'd like:
acb_contains_fmpz_fmpz
acb_set_fmpq_fmpq
acb_contains_fmpq_fmpq
acb_mat_set_fmpq_mat_fmpq_mat
acb_mat_contains_fmpq_mat_fmpq_mat

Any opinions on this idea? Right now I think the acb_mat/test/t-mul.c test doesn't actually involve any numbers with nonzero imaginary parts. My motivation for posting this issue was that when I looked at improving the coverage of that test I started to want more convenience functions.

Precomputed inverses

Algorithm like vector-scalar division and basecase power series division could be sped up by conditionally precomputing inverses.

It should be conditional, because if the number is small (in bit size) and the precision is high enough, it is better to just divide than to multiply by the inverse. For a complex number a+bi, one can still speed up this division by precomputing the extended denominator a^2+b^2.

relative exponential functions (aka exprel aka phi functions)?

There is a family of relative exponential functions in gsl. I wonder if it's worth adding some of these to arb? I think the real part of exprel(i*x) is sinc(x), and sinc is listed in #35.

index name hypgeom elementary
0 exp(z) 1F1(1, 1, z) exp(z)
1 exprel(z) 1F1(1; 2; z) (exp(z) - 1)/z
2 exprel_2(z) 1F1(1; 3; z) 2*(exp(z) - z - 1) / z^2
n exprel_n(n, z) 1F1(1; n+1; z) (n!/z^n) * (exp(x) - sum_k=0^(n-1) x^k/k!)

They are also related to exponential integrators where they are called phi functions, for example in
http://core.ac.uk/download/pdf/1633681.pdf
http://www1.maths.leeds.ac.uk/~jitse/phikrylov.pdf

A recent preprint using exprel:
http://arxiv.org/pdf/1504.01804.pdf

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.