gh-133895: correct values of cmath.cosh/sinh in case of overflows

[skirpichev]

This is a split-off gh-134995.

The C17 standard says (cis(y) is defined as cos(y) + i sin(y)):

• ccosh(+∞ + i0) returns +∞ + i0.
• ccosh(+∞ + iy) returns +∞ cis(y), for finite nonzero y. and
• csinh(+∞ + i0) returns +∞ + i0.
• csinh(+∞ + iy) returns +∞ cis(y), for positive finite y.

So far values, computed for exceptions, aren't accessible from the pure-Python world, yet we are trying to be correct in other places. The Lib/test/mathdata/cmath_testcases.txt has data points with correct numbers (see cosh0032 and sinh0032).

• Issue: Pass return value on ValueError exceptions in the cmath/math modules #133895

[skirpichev]

CC @serhiy-storchaka

[github-actions]

This PR is stale because it has been open for 30 days with no activity.

[read-the-docs-community]

Documentation build overview

📚 cpython-previews | 🛠️ Build #32820215 | 📁 Comparing d5bd0a0 against main (a7d5a6c)

  🔍 Preview build  

479 files changed · ± 448 modified · - 31 deleted

± Modified

• about.html
• bugs.html
• glossary.html
• index.html
• license.html
• py-modindex.html
• c-api/allocation.html
• c-api/apiabiversion.html
• c-api/arg.html
• c-api/buffer.html
• and 438 more...

- Deleted

• improve-page-nojs.html
• improve-page.html
• _static/tachyon-example-flamegraph.html
• _static/tachyon-example-heatmap.html
• c-api/curses.html
• c-api/interp-lifecycle.html
• c-api/picklebuffer.html
• c-api/profiling.html
• c-api/sentinel.html
• c-api/slots.html
• and 21 more...

[skirpichev]

CC @picnixz

[zware]

The cmath.rect change appears unrelated to this issue. I would approve and merge that change on its own (no issue required), but I do not have enough background to approve the rest of it (addressing the linked issue :)).

[skirpichev]

The cmath.rect change appears unrelated to this issue. I would approve and merge that change on its own (no issue required)

FYI: #154034

[skirpichev]

@vstinner, does it make sense for you?

Below attached pure-Python version of the main code and patched version, with tests against libm. Results for unpatched case:

$ ./python a.py 
csinh (720+0j) (inf+0j) (inf+nanj)
ccosh (720+0j) (inf+0j) (inf+nanj)
csinh (720-0j) (inf-0j) (inf+nanj)
ccosh (720-0j) (inf-0j) (inf+nanj)
csinh (-720+0j) (-inf+0j) (-inf+nanj)
ccosh (-720+0j) (inf-0j) (inf+nanj)
csinh (-720-0j) (-inf-0j) (-inf+nanj)
ccosh (-720-0j) (inf+0j) (inf+nanj)

Details

import cmath
import sys
from math import (isinf, isfinite, isnan, copysign, inf,
                  e, log, nan, isclose)
import ctypes
import ctypes.util
libm = ctypes.CDLL('libm.so.6')

@ctypes.util.wrap_dll_function(libm)
def csinh(z: ctypes.c_double_complex) -> ctypes.c_double_complex:
    pass
@ctypes.util.wrap_dll_function(libm)
def ccosh(z: ctypes.c_double_complex) -> ctypes.c_double_complex:
    pass

@ctypes.util.wrap_dll_function(libm)
def sinh(z: ctypes.c_double) -> ctypes.c_double:
    pass
@ctypes.util.wrap_dll_function(libm)
def sin(z: ctypes.c_double) -> ctypes.c_double:
    pass
@ctypes.util.wrap_dll_function(libm)
def cosh(z: ctypes.c_double) -> ctypes.c_double:
    pass
@ctypes.util.wrap_dll_function(libm)
def cos(z: ctypes.c_double) -> ctypes.c_double:
    pass


# Enumeration-like indices mapping to the 7x7 matrix taxonomy
ST_NINF  = 0  # -inf
ST_NEG   = 1  # finite number < 0
ST_NZERO = 2  # -0.0
ST_PZERO = 3  # +0.0
ST_POS   = 4  # finite number > 0
ST_PINF  = 5  # +inf
ST_NAN   = 6  # nan

U = -9.5426319407711027e33  # CPython placeholder

# Matrix rows = special_type(z.real), columns = special_type(z.imag)
SINH_SPECIAL_VALUES = [
    # ST_NINF (real == -inf)
    [(inf,nan), (U,U), (-inf,-0.0), (-inf,0.0), (U,U), (inf,nan), (inf,nan)],
    # ST_NEG (real < 0)
    [(nan,nan), (U,U), (U,U), (U,U), (U,U), (nan,nan),  (nan,nan)],
    # ST_NZERO (real == -0.0)
    [(0.0,nan), (U,U), (-0.0,-0.0), (-0.0,0.0), (U,U), (0.0,nan),  (0.0,nan)],
    # ST_PZERO (real == +0.0)
    [(0.0,nan), (U,U), (0.0,-0.0), (0.0,0.0), (U,U), (0.0,nan), (0.0,nan)],
    # ST_POS (real > 0)
    [(nan,nan), (U,U), (U,U), (U, U), (U,U), (nan,nan), (nan,nan)],
    # ST_PINF (real == +inf)
    [(inf,nan), (U,U), (inf,-0.0), (inf,0.0), (U,U), (inf,nan), (inf,nan)],
    # ST_NAN (real == nan)
    [(nan,nan), (nan,nan), (nan,-0.0), (nan,0.0), (nan,nan), (nan,nan), (nan,nan)]
]

COSH_SPECIAL_VALUES = [
    # ST_NINF (real == -inf)
    [(inf,nan), (U,U), (inf,0.0), (inf,-0.0), (U,U), (inf,nan), (inf,nan)],
    # ST_NEG (real < 0)
    [(nan,nan), (U,U), (U, U), (U, U), (U, U), (nan,nan), (nan,nan)],
    # ST_NZERO (real == -0.0)
    [(nan,0.0), (U, U), (1.0,0.0), (1.0,-0.0), (U, U), (nan,0.0), (nan,0.0)],
    # ST_PZERO (real == +0.0)
    [(nan,0.0), (U, U), (1.0,-0.0), (1.0,0.0), (U, U), (nan,0.0), (nan,0.0)],
    # ST_POS (real > 0)
    [(nan,nan), (U, U), (U, U), (U, U), (U, U), (nan,nan), (nan,nan)],
    # ST_PINF (real == +inf)
    [(inf,nan), (U,U), (inf,-0.0), (inf,0.0), (U,U), (inf,nan), (inf,nan)],
    # ST_NAN (real == nan)
    [(nan,nan), (nan,nan), (nan,0.0), (nan,0.0), (nan,nan), (nan,nan), (nan,nan)]
]


def special_type(x):
    if isfinite(x):
        if x:
            return ST_POS if copysign(1.0, x) == 1.0 else ST_NEG
        else:
            return ST_PZERO if copysign(1.0, x) == 1.0 else ST_NZERO
    elif isnan(x):
        return ST_NAN
    else:
        return ST_PINF if copysign(1.0, x) == 1.0 else ST_NINF


CM_LARGE_DOUBLE = sys.float_info.max/4
CM_LOG_LARGE_DOUBLE = log(CM_LARGE_DOUBLE)


def cis(x):
    return complex(cos(x), sin(x))


def py_sinh(z, patch=0):
    if not cmath.isfinite(z):
        if isinf(z.real) and isfinite(z.imag) and z.imag:
            if z.real > 0:
                return inf*cis(z.imag)
            else:
                return -inf*cis(z.imag).conjugate()

        row = special_type(x)
        col = special_type(y)
        return complex(*SINH_SPECIAL_VALUES[row][col])

    if abs(z.real) > CM_LOG_LARGE_DOUBLE:
        x_minus_one = z.real - copysign(1, z.real)
        r = complex(cos(z.imag) * sinh(x_minus_one),
                    sin(z.imag) * cosh(x_minus_one)) * e
        if patch and isnan(r.imag):
            r = complex(r.real, copysign(0.0, z.imag))
        return r
    else:
        return complex(cos(z.imag) * sinh(z.real),
                       sin(z.imag) * cosh(z.real))


def py_cosh(z, patch=0):
    if not cmath.isfinite(z):
        if isinf(z.real) and isfinite(z.imag) and z.imag:
            if z.real > 0:
                return inf*cis(z.imag)
            else:
                return inf*cis(z.imag).conjugate()

        row = special_type(x)
        col = special_type(y)
        return complex(*COSH_SPECIAL_VALUES[row][col])

    if abs(z.real) > CM_LOG_LARGE_DOUBLE:
        x_minus_one = z.real - copysign(1, z.real)
        r = complex(cos(z.imag) * cosh(x_minus_one),
                    sin(z.imag) * sinh(x_minus_one)) * e
        if patch and isnan(r.imag):
            r = complex(r.real, copysign(0.0, z.real*z.imag))
        return r
    else:
        return complex(cos(z.imag) * cosh(z.real),
                       sin(z.imag) * sinh(z.real))


def test(stmt, f, z, cr, r):
    if not stmt:
        print(f.__name__, z, cr, r)


if __name__ == "__main__":
    cases = [-inf, inf, nan, 0.0, -0.0, 0.5, 0.5,
             720, -720]  # overflows
    patch = 0
    if len(sys.argv) > 1:
        patch = int(sys.argv[1])

    for x in cases:
        for y in cases:
            z = complex(x, y)

            for pf, cf in [(py_sinh, csinh), (py_cosh, ccosh)]:
                r = pf(z, patch)
                cr = cf(z)

                if not cmath.isnan(r):
                    test(isclose(r.real, cr.real), cf, z, cr, r)
                    test(isclose(r.imag, cr.imag), cf, z, cr, r)
                else:
                    if isnan(r.real):
                        test(isnan(cr.real), cf, z, cr, r)
                        if isnan(r.imag):
                            test(isnan(cr.imag), cf, z, cr, r)
                        else:
                            test(isclose(r.imag, cr.imag), cf, z, cr, r)
                    else:
                        test(isnan(cr.imag), cf, z, cr, r)
                        test(isclose(r.real, cr.real), cf, z, cr, r)

[vstinner]

You didn't update test_cmath. Does it mean that modified code is not currently tested? If yes, would it be possible to add tests?

[skirpichev]

Does it mean that modified code is not currently tested?

Yes, currently these values aren't accessible from the pure-Python world. See PR description and the issue thread.