Merge pull request #20465 from charris/backport-20360

 BUG: Force ``npymath` ` to respect ``npy_longdouble``

numpy/core/include/numpy/npy_common.h

@@ -356,14 +356,31 @@ typedef unsigned long npy_ulonglong;
 typedef unsigned char npy_bool;
 #define NPY_FALSE 0
 #define NPY_TRUE 1
-
-
+/*
+ * `NPY_SIZEOF_LONGDOUBLE` isn't usually equal to sizeof(long double).
+ * In some certain cases, it may forced to be equal to sizeof(double)
+ * even against the compiler implementation and the same goes for
+ * `complex long double`.
+ *
+ * Therefore, avoid `long double`, use `npy_longdouble` instead,
+ * and when it comes to standard math functions make sure of using
+ * the double version when `NPY_SIZEOF_LONGDOUBLE` == `NPY_SIZEOF_DOUBLE`.
+ * For example:
+ *   npy_longdouble *ptr, x;
+ *   #if NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE
+ *       npy_longdouble r = modf(x, ptr);
+ *   #else
+ *       npy_longdouble r = modfl(x, ptr);
+ *   #endif
+ *
+ * See https://github.com/numpy/numpy/issues/20348
+ */
 #if NPY_SIZEOF_LONGDOUBLE == NPY_SIZEOF_DOUBLE
-        typedef double npy_longdouble;
-        #define NPY_LONGDOUBLE_FMT "g"
+    #define NPY_LONGDOUBLE_FMT "g"
+    typedef double npy_longdouble;
 #else
-        typedef long double npy_longdouble;
-        #define NPY_LONGDOUBLE_FMT "Lg"
+    #define NPY_LONGDOUBLE_FMT "Lg"
+    typedef long double npy_longdouble;
 #endif
 
 #ifndef Py_USING_UNICODE

numpy/core/src/multiarray/_multiarray_tests.c.src

@@ -2093,7 +2093,7 @@ PrintFloat_Printf_g(PyObject *obj, int precision)
     }
     else if (PyArray_IsScalar(obj, LongDouble)) {
         npy_longdouble x = PyArrayScalar_VAL(obj, LongDouble);
-        PyOS_snprintf(str, sizeof(str), "%.*Lg", precision, x);
+        PyOS_snprintf(str, sizeof(str), "%.*" NPY_LONGDOUBLE_FMT, precision, x);
     }
     else{
         double val = PyFloat_AsDouble(obj);

numpy/core/src/npymath/npy_math_internal.h.src

@@ -480,10 +480,16 @@ NPY_INPLACE @type@ npy_frexp@c@(@type@ x, int* exp)
 
 /**begin repeat
  * #type = npy_longdouble, npy_double, npy_float#
+ * #TYPE = LONGDOUBLE, DOUBLE, FLOAT#
  * #c = l,,f#
  * #C = L,,F#
  */
-
+#undef NPY__FP_SFX
+#if NPY_SIZEOF_@TYPE@ == NPY_SIZEOF_DOUBLE
+    #define NPY__FP_SFX(X) X
+#else
+    #define NPY__FP_SFX(X) NPY_CAT(X, @c@)
+#endif
 /*
  * On arm64 macOS, there's a bug with sin, cos, and tan where they don't
  * raise "invalid" when given INFINITY as input.
@@ -509,7 +515,7 @@ NPY_INPLACE @type@ npy_@kind@@c@(@type@ x)
         return (x - x);
     }
 #endif
-    return @kind@@c@(x);
+    return NPY__FP_SFX(@kind@)(x);
 }
 #endif
 
@@ -524,7 +530,7 @@ NPY_INPLACE @type@ npy_@kind@@c@(@type@ x)
 #ifdef HAVE_@KIND@@C@
 NPY_INPLACE @type@ npy_@kind@@c@(@type@ x, @type@ y)
 {
-    return @kind@@c@(x, y);
+    return NPY__FP_SFX(@kind@)(x, y);
 }
 #endif
 /**end repeat1**/
@@ -555,21 +561,21 @@ npy_@kind@@c@(@type@ x, @type@ y)
 #ifdef HAVE_MODF@C@
 NPY_INPLACE @type@ npy_modf@c@(@type@ x, @type@ *iptr)
 {
-    return modf@c@(x, iptr);
+    return NPY__FP_SFX(modf)(x, iptr);
 }
 #endif
 
 #ifdef HAVE_LDEXP@C@
 NPY_INPLACE @type@ npy_ldexp@c@(@type@ x, int exp)
 {
-    return ldexp@c@(x, exp);
+    return NPY__FP_SFX(ldexp)(x, exp);
 }
 #endif
 
 #ifdef HAVE_FREXP@C@
 NPY_INPLACE @type@ npy_frexp@c@(@type@ x, int* exp)
 {
-    return frexp@c@(x, exp);
+    return NPY__FP_SFX(frexp)(x, exp);
 }
 #endif
 
@@ -592,10 +598,10 @@ NPY_INPLACE @type@ npy_cbrt@c@(@type@ x)
 #else
 NPY_INPLACE @type@ npy_cbrt@c@(@type@ x)
 {
-    return cbrt@c@(x);
+    return NPY__FP_SFX(cbrt)(x);
 }
 #endif
-
+#undef NPY__FP_SFX
 /**end repeat**/
 
 
@@ -605,10 +611,16 @@ NPY_INPLACE @type@ npy_cbrt@c@(@type@ x)
 
 /**begin repeat
  * #type = npy_float, npy_double, npy_longdouble#
+ * #TYPE = FLOAT, DOUBLE, LONGDOUBLE#
  * #c = f, ,l#
  * #C = F, ,L#
  */
-
+#undef NPY__FP_SFX
+#if NPY_SIZEOF_@TYPE@ == NPY_SIZEOF_DOUBLE
+    #define NPY__FP_SFX(X) X
+#else
+    #define NPY__FP_SFX(X) NPY_CAT(X, @c@)
+#endif
 @type@ npy_heaviside@c@(@type@ x, @type@ h0)
 {
     if (npy_isnan(x)) {
@@ -625,10 +637,10 @@ NPY_INPLACE @type@ npy_cbrt@c@(@type@ x)
     }
 }
 
-#define LOGE2    NPY_LOGE2@c@
-#define LOG2E    NPY_LOG2E@c@
-#define RAD2DEG  (180.0@c@/NPY_PI@c@)
-#define DEG2RAD  (NPY_PI@c@/180.0@c@)
+#define LOGE2    NPY__FP_SFX(NPY_LOGE2)
+#define LOG2E    NPY__FP_SFX(NPY_LOG2E)
+#define RAD2DEG  (NPY__FP_SFX(180.0)/NPY__FP_SFX(NPY_PI))
+#define DEG2RAD  (NPY__FP_SFX(NPY_PI)/NPY__FP_SFX(180.0))
 
 NPY_INPLACE @type@ npy_rad2deg@c@(@type@ x)
 {
@@ -783,7 +795,7 @@ npy_divmod@c@(@type@ a, @type@ b, @type@ *modulus)
 #undef LOG2E
 #undef RAD2DEG
 #undef DEG2RAD
-
+#undef NPY__FP_SFX
 /**end repeat**/
 
 /**begin repeat