mirror of
https://github.com/python/cpython.git
synced 2025-08-04 17:08:35 +00:00
gh-69639: Add mixed-mode rules for complex arithmetic (C-like) (GH-124829)
"Generally, mixed-mode arithmetic combining real and complex variables should be performed directly, not by first coercing the real to complex, lest the sign of zero be rendered uninformative; the same goes for combinations of pure imaginary quantities with complex variables." (c) Kahan, W: Branch cuts for complex elementary functions. This patch implements mixed-mode arithmetic rules, combining real and complex variables as specified by C standards since C99 (in particular, there is no special version for the true division with real lhs operand). Most C compilers implementing C99+ Annex G have only these special rules (without support for imaginary type, which is going to be deprecated in C2y).
This commit is contained in:
Sergey B Kirpichev
GitHub
parent
dcf629213b
commit
987311d42e
15 changed files with 449 additions and 98 deletions
|
|
@ -8,6 +8,7 @@
|
|||
#include "Python.h"
|
||||
#include "pycore_call.h" // _PyObject_CallNoArgs()
|
||||
#include "pycore_complexobject.h" // _PyComplex_FormatAdvancedWriter()
|
||||
#include "pycore_floatobject.h" // _Py_convert_int_to_double()
|
||||
#include "pycore_long.h" // _PyLong_GetZero()
|
||||
#include "pycore_object.h" // _PyObject_Init()
|
||||
#include "pycore_pymath.h" // _Py_ADJUST_ERANGE2()
|
||||
|
|
@ -34,6 +35,20 @@ _Py_c_sum(Py_complex a, Py_complex b)
|
|||
return r;
|
||||
}
|
||||
|
||||
Py_complex
|
||||
_Py_cr_sum(Py_complex a, double b)
|
||||
{
|
||||
Py_complex r = a;
|
||||
r.real += b;
|
||||
return r;
|
||||
}
|
||||
|
||||
static inline Py_complex
|
||||
_Py_rc_sum(double a, Py_complex b)
|
||||
{
|
||||
return _Py_cr_sum(b, a);
|
||||
}
|
||||
|
||||
Py_complex
|
||||
_Py_c_diff(Py_complex a, Py_complex b)
|
||||
{
|
||||
|
|
@ -43,6 +58,23 @@ _Py_c_diff(Py_complex a, Py_complex b)
|
|||
return r;
|
||||
}
|
||||
|
||||
Py_complex
|
||||
_Py_cr_diff(Py_complex a, double b)
|
||||
{
|
||||
Py_complex r = a;
|
||||
r.real -= b;
|
||||
return r;
|
||||
}
|
||||
|
||||
Py_complex
|
||||
_Py_rc_diff(double a, Py_complex b)
|
||||
{
|
||||
Py_complex r;
|
||||
r.real = a - b.real;
|
||||
r.imag = -b.imag;
|
||||
return r;
|
||||
}
|
||||
|
||||
Py_complex
|
||||
_Py_c_neg(Py_complex a)
|
||||
{
|
||||
|
|
@ -61,6 +93,21 @@ _Py_c_prod(Py_complex a, Py_complex b)
|
|||
return r;
|
||||
}
|
||||
|
||||
Py_complex
|
||||
_Py_cr_prod(Py_complex a, double b)
|
||||
{
|
||||
Py_complex r = a;
|
||||
r.real *= b;
|
||||
r.imag *= b;
|
||||
return r;
|
||||
}
|
||||
|
||||
static inline Py_complex
|
||||
_Py_rc_prod(double a, Py_complex b)
|
||||
{
|
||||
return _Py_cr_prod(b, a);
|
||||
}
|
||||
|
||||
/* Avoid bad optimization on Windows ARM64 until the compiler is fixed */
|
||||
#ifdef _M_ARM64
|
||||
#pragma optimize("", off)
|
||||
|
|
@ -143,6 +190,64 @@ _Py_c_quot(Py_complex a, Py_complex b)
|
|||
|
||||
return r;
|
||||
}
|
||||
|
||||
Py_complex
|
||||
_Py_cr_quot(Py_complex a, double b)
|
||||
{
|
||||
Py_complex r = a;
|
||||
if (b) {
|
||||
r.real /= b;
|
||||
r.imag /= b;
|
||||
}
|
||||
else {
|
||||
errno = EDOM;
|
||||
r.real = r.imag = 0.0;
|
||||
}
|
||||
return r;
|
||||
}
|
||||
|
||||
/* an equivalent of _Py_c_quot() function, when 1st argument is real */
|
||||
Py_complex
|
||||
_Py_rc_quot(double a, Py_complex b)
|
||||
{
|
||||
Py_complex r;
|
||||
const double abs_breal = b.real < 0 ? -b.real : b.real;
|
||||
const double abs_bimag = b.imag < 0 ? -b.imag : b.imag;
|
||||
|
||||
if (abs_breal >= abs_bimag) {
|
||||
if (abs_breal == 0.0) {
|
||||
errno = EDOM;
|
||||
r.real = r.imag = 0.0;
|
||||
}
|
||||
else {
|
||||
const double ratio = b.imag / b.real;
|
||||
const double denom = b.real + b.imag * ratio;
|
||||
r.real = a / denom;
|
||||
r.imag = (-a * ratio) / denom;
|
||||
}
|
||||
}
|
||||
else if (abs_bimag >= abs_breal) {
|
||||
const double ratio = b.real / b.imag;
|
||||
const double denom = b.real * ratio + b.imag;
|
||||
assert(b.imag != 0.0);
|
||||
r.real = (a * ratio) / denom;
|
||||
r.imag = (-a) / denom;
|
||||
}
|
||||
else {
|
||||
r.real = r.imag = Py_NAN;
|
||||
}
|
||||
|
||||
if (isnan(r.real) && isnan(r.imag) && isfinite(a)
|
||||
&& (isinf(abs_breal) || isinf(abs_bimag)))
|
||||
{
|
||||
const double x = copysign(isinf(b.real) ? 1.0 : 0.0, b.real);
|
||||
const double y = copysign(isinf(b.imag) ? 1.0 : 0.0, b.imag);
|
||||
r.real = 0.0 * (a*x);
|
||||
r.imag = 0.0 * (-a*y);
|
||||
}
|
||||
|
||||
return r;
|
||||
}
|
||||
#ifdef _M_ARM64
|
||||
#pragma optimize("", on)
|
||||
#endif
|
||||
|
|
@ -474,84 +579,91 @@ complex_hash(PyComplexObject *v)
|
|||
}
|
||||
|
||||
/* This macro may return! */
|
||||
#define TO_COMPLEX(obj, c) \
|
||||
if (PyComplex_Check(obj)) \
|
||||
c = ((PyComplexObject *)(obj))->cval; \
|
||||
else if (to_complex(&(obj), &(c)) < 0) \
|
||||
#define TO_COMPLEX(obj, c) \
|
||||
if (PyComplex_Check(obj)) \
|
||||
c = ((PyComplexObject *)(obj))->cval; \
|
||||
else if (real_to_complex(&(obj), &(c)) < 0) \
|
||||
return (obj)
|
||||
|
||||
static int
|
||||
to_complex(PyObject **pobj, Py_complex *pc)
|
||||
real_to_double(PyObject **pobj, double *dbl)
|
||||
{
|
||||
PyObject *obj = *pobj;
|
||||
|
||||
pc->real = pc->imag = 0.0;
|
||||
if (PyLong_Check(obj)) {
|
||||
pc->real = PyLong_AsDouble(obj);
|
||||
if (pc->real == -1.0 && PyErr_Occurred()) {
|
||||
*pobj = NULL;
|
||||
return -1;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
if (PyFloat_Check(obj)) {
|
||||
pc->real = PyFloat_AsDouble(obj);
|
||||
return 0;
|
||||
*dbl = PyFloat_AS_DOUBLE(obj);
|
||||
}
|
||||
*pobj = Py_NewRef(Py_NotImplemented);
|
||||
return -1;
|
||||
}
|
||||
|
||||
|
||||
static PyObject *
|
||||
complex_add(PyObject *v, PyObject *w)
|
||||
{
|
||||
Py_complex result;
|
||||
Py_complex a, b;
|
||||
TO_COMPLEX(v, a);
|
||||
TO_COMPLEX(w, b);
|
||||
result = _Py_c_sum(a, b);
|
||||
return PyComplex_FromCComplex(result);
|
||||
}
|
||||
|
||||
static PyObject *
|
||||
complex_sub(PyObject *v, PyObject *w)
|
||||
{
|
||||
Py_complex result;
|
||||
Py_complex a, b;
|
||||
TO_COMPLEX(v, a);
|
||||
TO_COMPLEX(w, b);
|
||||
result = _Py_c_diff(a, b);
|
||||
return PyComplex_FromCComplex(result);
|
||||
}
|
||||
|
||||
static PyObject *
|
||||
complex_mul(PyObject *v, PyObject *w)
|
||||
{
|
||||
Py_complex result;
|
||||
Py_complex a, b;
|
||||
TO_COMPLEX(v, a);
|
||||
TO_COMPLEX(w, b);
|
||||
result = _Py_c_prod(a, b);
|
||||
return PyComplex_FromCComplex(result);
|
||||
}
|
||||
|
||||
static PyObject *
|
||||
complex_div(PyObject *v, PyObject *w)
|
||||
{
|
||||
Py_complex quot;
|
||||
Py_complex a, b;
|
||||
TO_COMPLEX(v, a);
|
||||
TO_COMPLEX(w, b);
|
||||
errno = 0;
|
||||
quot = _Py_c_quot(a, b);
|
||||
if (errno == EDOM) {
|
||||
PyErr_SetString(PyExc_ZeroDivisionError, "division by zero");
|
||||
return NULL;
|
||||
else if (_Py_convert_int_to_double(pobj, dbl) < 0) {
|
||||
return -1;
|
||||
}
|
||||
return PyComplex_FromCComplex(quot);
|
||||
return 0;
|
||||
}
|
||||
|
||||
static int
|
||||
real_to_complex(PyObject **pobj, Py_complex *pc)
|
||||
{
|
||||
pc->imag = 0.0;
|
||||
return real_to_double(pobj, &(pc->real));
|
||||
}
|
||||
|
||||
/* Complex arithmetic rules implement special mixed-mode case where combining
|
||||
a pure-real (float or int) value and a complex value is performed directly
|
||||
without first coercing the real value to a complex value.
|
||||
|
||||
Let us consider the addition as an example, assuming that ints are implicitly
|
||||
converted to floats. We have the following rules (up to variants with changed
|
||||
order of operands):
|
||||
|
||||
complex(a, b) + complex(c, d) = complex(a + c, b + d)
|
||||
float(a) + complex(b, c) = complex(a + b, c)
|
||||
|
||||
Similar rules are implemented for subtraction, multiplication and division.
|
||||
See C11's Annex G, sections G.5.1 and G.5.2.
|
||||
*/
|
||||
|
||||
#define COMPLEX_BINOP(NAME, FUNC) \
|
||||
static PyObject * \
|
||||
complex_##NAME(PyObject *v, PyObject *w) \
|
||||
{ \
|
||||
Py_complex a; \
|
||||
errno = 0; \
|
||||
if (PyComplex_Check(w)) { \
|
||||
Py_complex b = ((PyComplexObject *)w)->cval; \
|
||||
if (PyComplex_Check(v)) { \
|
||||
a = ((PyComplexObject *)v)->cval; \
|
||||
a = _Py_c_##FUNC(a, b); \
|
||||
} \
|
||||
else if (real_to_double(&v, &a.real) < 0) { \
|
||||
return v; \
|
||||
} \
|
||||
else { \
|
||||
a = _Py_rc_##FUNC(a.real, b); \
|
||||
} \
|
||||
} \
|
||||
else if (!PyComplex_Check(v)) { \
|
||||
Py_RETURN_NOTIMPLEMENTED; \
|
||||
} \
|
||||
else { \
|
||||
a = ((PyComplexObject *)v)->cval; \
|
||||
double b; \
|
||||
if (real_to_double(&w, &b) < 0) { \
|
||||
return w; \
|
||||
} \
|
||||
a = _Py_cr_##FUNC(a, b); \
|
||||
} \
|
||||
if (errno == EDOM) { \
|
||||
PyErr_SetString(PyExc_ZeroDivisionError, \
|
||||
"division by zero"); \
|
||||
return NULL; \
|
||||
} \
|
||||
return PyComplex_FromCComplex(a); \
|
||||
}
|
||||
|
||||
COMPLEX_BINOP(add, sum)
|
||||
COMPLEX_BINOP(mul, prod)
|
||||
COMPLEX_BINOP(sub, diff)
|
||||
COMPLEX_BINOP(div, quot)
|
||||
|
||||
static PyObject *
|
||||
complex_pow(PyObject *v, PyObject *w, PyObject *z)
|
||||
{
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue