numerical failure in exp_integral

fchapoton commented 1 year ago

found by the patchbot in some tickets

sage -t --long --warn-long 61.1 --random-seed=242205119125380898007386286571458929899 src/sage/functions/exp_integral.py

giving

File "src/sage/functions/exp_integral.py", line 1493, in sage.functions.exp_integral.exponential_integral_1
Failed example:
    for prec in [20..128]:  # long time (15s on sage.math, 2013)
        R = RealField(prec)
        S = RealField(prec+64)
        a = R.random_element(-15,10).exp()
        n = 2^ZZ.random_element(14)
        x = exponential_integral_1(a, n)
        y = exponential_integral_1(S(a), n)
        c = RDF(4 * max(1.0, y[0]))
        for i in range(n):
            e = float(abs(S(x[i]) - y[i]) << prec)
            if e >= c:
                print("exponential_integral_1(%s, %s)[%s] with precision %s has error of %s >= %s"%(a, n, i, prec, e, c))
Expected nothing
Got:
    exponential_integral_1(49.05921476503182064969544573696770370, 1)[0] with precision 127 has error of 6.330959330559951 >= 4.0

Component: numerical

Issue created by migration from https://trac.sagemath.org/ticket/34845

fchapoton commented 1 year ago

still there in SageMath version 10.2.beta3,

collares commented 1 month ago

Still a problem with Sage 10.4:

sage -t --long --random-seed=289960408163339990509890784516587877557 /nix/store/m58x4sba0rmcbgj45x5vb4f9zhwc0jfg-sage-src-10.4/src/sage/functions/exp_integral.py
**********************************************************************
File "/nix/store/m58x4sba0rmcbgj45x5vb4f9zhwc0jfg-sage-src-10.4/src/sage/functions/exp_integral.py", line 1508, in sage.functions.exp_integral.exponential_integral_1
Failed example:
    for prec in [20..128]:            # long time (15s on sage.math, 2013), needs sage.libs.pari
        R = RealField(prec)
        S = RealField(prec+64)
        a = R.random_element(-15,10).exp()
        n = 2^ZZ.random_element(14)
        x = exponential_integral_1(a, n)
        y = exponential_integral_1(S(a), n)
        c = RDF(4 * max(1.0, y[0]))
        for i in range(n):
            e = float(abs(S(x[i]) - y[i]) << prec)
            if e >= c:
                print("exponential_integral_1(%s, %s)[%s] with precision %s has error of %s >= %s"%(a, n, i, prec, e, c))
Expected nothing
Got:
    exponential_integral_1(40.20025815043087219312615077426751077, 1)[0] with precision 127 has error of 6.194629742503663 >= 4.0

sagemath / sage

numerical failure in exp_integral #34845