MAINT: Add frexp and ldexp to npymath and generate their ufuncs

7 files changed, 91 insertions, 207 deletions

diff --git a/numpy/core/code_generators/generate_umath.py b/numpy/core/code_generators/generate_umath.py
index e3c9cf28b..9f8d4c688 100644
--- a/numpy/core/code_generators/generate_umath.py
+++ b/numpy/core/code_generators/generate_umath.py
@@ -20,6 +20,12 @@ ReorderableNone = "PyUFunc_ReorderableNone"
 class FullTypeDescr(object):
     pass
 
+class FuncNameSuffix(object):
+    """Stores the suffix to append when generating functions names.
+    """
+    def __init__(self, suffix):
+        self.suffix = suffix
+
 class TypeDescription(object):
     """Type signature for a ufunc.
 
@@ -795,6 +801,30 @@ defdict = {
           None,
           TD(flts),
           ),
+'ldexp' :
+    Ufunc(2, 1, None,
+          docstrings.get('numpy.core.umath.ldexp'),
+          None,
+          [TypeDescription('e', None, 'ei', 'e'),
+          TypeDescription('f', None, 'fi', 'f'),
+          TypeDescription('e', FuncNameSuffix('long'), 'el', 'e'),
+          TypeDescription('f', FuncNameSuffix('long'), 'fl', 'f'),
+          TypeDescription('d', None, 'di', 'd'),
+          TypeDescription('d', FuncNameSuffix('long'), 'dl', 'd'),
+          TypeDescription('g', None, 'gi', 'g'),
+          TypeDescription('g', FuncNameSuffix('long'), 'gl', 'g'),
+          ],
+          ),
+'frexp' :
+    Ufunc(1, 2, None,
+          docstrings.get('numpy.core.umath.frexp'),
+          None,
+          [TypeDescription('e', None, 'e', 'ei'),
+          TypeDescription('f', None, 'f', 'fi'),
+          TypeDescription('d', None, 'd', 'di'),
+          TypeDescription('g', None, 'g', 'gi'),
+          ],
+          )
 }
 
 if sys.version_info[0] >= 3:
@@ -854,7 +884,7 @@ def make_arrays(funcdict):
             thedict = chartotype1  # one input and one output
 
         for t in uf.type_descriptions:
-            if t.func_data not in (None, FullTypeDescr):
+            if t.func_data not in (None, FullTypeDescr) and not isinstance(t.func_data, FuncNameSuffix):
                 funclist.append('NULL')
                 astype = ''
                 if not t.astype is None:
@@ -880,6 +910,10 @@ def make_arrays(funcdict):
                 tname = english_upper(chartoname[t.type])
                 datalist.append('(void *)NULL')
                 funclist.append('%s_%s_%s_%s' % (tname, t.in_, t.out, name))
+            elif isinstance(t.func_data, FuncNameSuffix):
+                datalist.append('(void *)NULL')
+                tname = english_upper(chartoname[t.type])
+                funclist.append('%s_%s_%s' % (tname, name, t.func_data.suffix))
             else:
                 datalist.append('(void *)NULL')
                 tname = english_upper(chartoname[t.type])
diff --git a/numpy/core/code_generators/ufunc_docstrings.py b/numpy/core/code_generators/ufunc_docstrings.py
index 4d302969e..804108397 100644
--- a/numpy/core/code_generators/ufunc_docstrings.py
+++ b/numpy/core/code_generators/ufunc_docstrings.py
@@ -3324,9 +3324,6 @@ add_newdoc('numpy.core.umath', 'true_divide',
 
     """)
 
-# This doc is not currently used, but has been converted to a C string
-# that can be found in numpy/core/src/umath/umathmodule.c where the
-# frexp ufunc is constructed.
 add_newdoc('numpy.core.umath', 'frexp',
     """
     Decompose the elements of x into mantissa and twos exponent.
@@ -3372,9 +3369,6 @@ add_newdoc('numpy.core.umath', 'frexp',
 
     """)
 
-# This doc is not currently used, but has been converted to a C string
-# that can be found in numpy/core/src/umath/umathmodule.c where the
-# ldexp ufunc is constructed.
 add_newdoc('numpy.core.umath', 'ldexp',
     """
     Returns x1 * 2**x2, element-wise.
diff --git a/numpy/core/include/numpy/npy_math.h b/numpy/core/include/numpy/npy_math.h
index 6b4f622cf..855ddf7fa 100644
--- a/numpy/core/include/numpy/npy_math.h
+++ b/numpy/core/include/numpy/npy_math.h
@@ -143,6 +143,8 @@ double npy_log2(double x);
 double npy_atan2(double x, double y);
 double npy_pow(double x, double y);
 double npy_modf(double x, double* y);
+double npy_frexp(double x, int* y);
+double npy_ldexp(double n, int y);
 
 double npy_copysign(double x, double y);
 double npy_nextafter(double x, double y);
@@ -247,6 +249,8 @@ float npy_powf(float x, float y);
 float npy_fmodf(float x, float y);
 
 float npy_modff(float x, float* y);
+float npy_frexpf(float x, int* y);
+float npy_ldexpf(float x, int y);
 
 float npy_copysignf(float x, float y);
 float npy_nextafterf(float x, float y);
@@ -288,6 +292,8 @@ npy_longdouble npy_powl(npy_longdouble x, npy_longdouble y);
 npy_longdouble npy_fmodl(npy_longdouble x, npy_longdouble y);
 
 npy_longdouble npy_modfl(npy_longdouble x, npy_longdouble* y);
+npy_longdouble npy_frexpl(npy_longdouble x, int* y);
+npy_longdouble npy_ldexpl(npy_longdouble x, int y);
 
 npy_longdouble npy_copysignl(npy_longdouble x, npy_longdouble y);
 npy_longdouble npy_nextafterl(npy_longdouble x, npy_longdouble y);
diff --git a/numpy/core/src/npymath/npy_math.c.src b/numpy/core/src/npymath/npy_math.c.src
index 05af0b132..3a1be3745 100644
--- a/numpy/core/src/npymath/npy_math.c.src
+++ b/numpy/core/src/npymath/npy_math.c.src
@@ -343,6 +343,7 @@ double npy_log2(double x)
  * asinh, acosh, atanh
  *
  * hypot, atan2, pow, fmod, modf
+ * ldexp, frexp
  *
  * We assume the above are always available in their double versions.
  *
@@ -405,6 +406,26 @@ double npy_log2(double x)
 }
 #endif
 
+#ifdef ldexp@c@
+#undef ldexp@c@
+#endif
+#ifndef HAVE_LDEXP@C@
+@type@ npy_ldexp@c@(@type@ x, int exp)
+{
+    return (@type@) npy_ldexp((double)x, exp);
+}
+#endif 
+
+#ifdef frexp@c@
+#undef frexp@c@
+#endif
+#ifndef HAVE_FREXP@C@
+@type@ npy_frexp@c@(@type@ x, int* exp)
+{
+    return (@type@) npy_frexp(x, exp);
+}
+#endif
+
 /**end repeat**/
 
 
@@ -451,6 +472,20 @@ double npy_log2(double x)
 }
 #endif
 
+#ifdef HAVE_LDEXP@C@
+@type@ npy_ldexp@c@(@type@ x, int exp)
+{
+    return ldexp@c@(x, exp);
+}
+#endif 
+
+#ifdef HAVE_FREXP@C@
+@type@ npy_frexp@c@(@type@ x, int* exp)
+{
+    return frexp@c@(x, exp);
+}
+#endif
+
 /**end repeat**/
 
 
diff --git a/numpy/core/src/umath/loops.c.src b/numpy/core/src/umath/loops.c.src
index 035a27fd2..d747864f8 100644
--- a/numpy/core/src/umath/loops.c.src
+++ b/numpy/core/src/umath/loops.c.src
@@ -1743,25 +1743,22 @@ NPY_NO_EXPORT void
     }
 }
 
-#ifdef HAVE_FREXP@C@
 NPY_NO_EXPORT void
 @TYPE@_frexp(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
 {
     UNARY_LOOP_TWO_OUT {
         const @type@ in1 = *(@type@ *)ip1;
-        *((@type@ *)op1) = frexp@c@(in1, (int *)op2);
+        *((@type@ *)op1) = npy_frexp@c@(in1, (int *)op2);
     }
 }
-#endif
 
-#ifdef HAVE_LDEXP@C@
 NPY_NO_EXPORT void
 @TYPE@_ldexp(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
 {
     BINARY_LOOP {
         const @type@ in1 = *(@type@ *)ip1;
         const int in2 = *(int *)ip2;
-        *((@type@ *)op1) = ldexp@c@(in1, in2);
+        *((@type@ *)op1) = npy_ldexp@c@(in1, in2);
     }
 }
 
@@ -1778,7 +1775,7 @@ NPY_NO_EXPORT void
         const long in2 = *(long *)ip2;
         if (((int)in2) == in2) {
             /* Range OK */
-            *((@type@ *)op1) = ldexp@c@(in1, ((int)in2));
+            *((@type@ *)op1) = npy_ldexp@c@(in1, ((int)in2));
         }
         else {
             /*
@@ -1786,15 +1783,14 @@ NPY_NO_EXPORT void
              * given that exponent has less bits than npy_int.
              */
             if (in2 > 0) {
-                *((@type@ *)op1) = ldexp@c@(in1, NPY_MAX_INT);
+                *((@type@ *)op1) = npy_ldexp@c@(in1, NPY_MAX_INT);
             }
             else {
-                *((@type@ *)op1) = ldexp@c@(in1, NPY_MIN_INT);
+                *((@type@ *)op1) = npy_ldexp@c@(in1, NPY_MIN_INT);
             }
         }
     }
 }
-#endif
 
 #define @TYPE@_true_divide @TYPE@_divide
 
@@ -2059,25 +2055,22 @@ HALF_modf(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(f
     }
 }
 
-#ifdef HAVE_FREXPF
 NPY_NO_EXPORT void
 HALF_frexp(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
 {
     UNARY_LOOP_TWO_OUT {
         const float in1 = npy_half_to_float(*(npy_half *)ip1);
-        *((npy_half *)op1) = npy_float_to_half(frexpf(in1, (int *)op2));
+        *((npy_half *)op1) = npy_float_to_half(npy_frexpf(in1, (int *)op2));
     }
 }
-#endif
 
-#ifdef HAVE_LDEXPF
 NPY_NO_EXPORT void
 HALF_ldexp(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func))
 {
     BINARY_LOOP {
         const float in1 = npy_half_to_float(*(npy_half *)ip1);
         const int in2 = *(int *)ip2;
-        *((npy_half *)op1) = npy_float_to_half(ldexpf(in1, in2));
+        *((npy_half *)op1) = npy_float_to_half(npy_ldexpf(in1, in2));
     }
 }
 
@@ -2094,7 +2087,7 @@ HALF_ldexp_long(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UN
         const long in2 = *(long *)ip2;
         if (((int)in2) == in2) {
             /* Range OK */
-            *((npy_half *)op1) = npy_float_to_half(ldexpf(in1, ((int)in2)));
+            *((npy_half *)op1) = npy_float_to_half(npy_ldexpf(in1, ((int)in2)));
         }
         else {
             /*
@@ -2102,15 +2095,14 @@ HALF_ldexp_long(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UN
              * given that exponent has less bits than npy_int.
              */
             if (in2 > 0) {
-                *((npy_half *)op1) = npy_float_to_half(ldexpf(in1, NPY_MAX_INT));
+                *((npy_half *)op1) = npy_float_to_half(npy_ldexpf(in1, NPY_MAX_INT));
             }
             else {
-                *((npy_half *)op1) = npy_float_to_half(ldexpf(in1, NPY_MIN_INT));
+                *((npy_half *)op1) = npy_float_to_half(npy_ldexpf(in1, NPY_MIN_INT));
             }
         }
     }
 }
-#endif
 
 #define HALF_true_divide HALF_divide
 
diff --git a/numpy/core/src/umath/loops.h.src b/numpy/core/src/umath/loops.h.src
index fdc9230de..a6e775a3a 100644
--- a/numpy/core/src/umath/loops.h.src
+++ b/numpy/core/src/umath/loops.h.src
@@ -248,17 +248,13 @@ NPY_NO_EXPORT void
 NPY_NO_EXPORT void
 @TYPE@_modf(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func));
 
-#ifdef HAVE_FREXP@C@
 NPY_NO_EXPORT void
 @TYPE@_frexp(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func));
-#endif
 
-#ifdef HAVE_LDEXP@C@
 NPY_NO_EXPORT void
 @TYPE@_ldexp(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func));
 NPY_NO_EXPORT void
 @TYPE@_ldexp_long(char **args, npy_intp *dimensions, npy_intp *steps, void *NPY_UNUSED(func));
-#endif
 
 #define @TYPE@_true_divide @TYPE@_divide
 
diff --git a/numpy/core/src/umath/umathmodule.c b/numpy/core/src/umath/umathmodule.c
index 52aa2e48d..57b2bb239 100644
--- a/numpy/core/src/umath/umathmodule.c
+++ b/numpy/core/src/umath/umathmodule.c
@@ -201,182 +201,6 @@ ufunc_frompyfunc(PyObject *NPY_UNUSED(dummy), PyObject *args, PyObject *NPY_UNUS
  *****************************************************************************
  */
 
-/* Less automated additions to the ufuncs */
-
-static PyUFuncGenericFunction frexp_functions[] = {
-#ifdef HAVE_FREXPF
-    HALF_frexp,
-    FLOAT_frexp,
-#endif
-    DOUBLE_frexp
-#ifdef HAVE_FREXPL
-    ,LONGDOUBLE_frexp
-#endif
-};
-
-static char frexp_signatures[] = {
-#ifdef HAVE_FREXPF
-    NPY_HALF, NPY_HALF, NPY_INT,
-    NPY_FLOAT, NPY_FLOAT, NPY_INT,
-#endif
-    NPY_DOUBLE, NPY_DOUBLE, NPY_INT
-#ifdef HAVE_FREXPL
-    ,NPY_LONGDOUBLE, NPY_LONGDOUBLE, NPY_INT
-#endif
-};
-static void * blank_data[12];
-
-#if NPY_SIZEOF_LONG == NPY_SIZEOF_INT
-#define LDEXP_LONG(typ) typ##_ldexp
-#else
-#define LDEXP_LONG(typ) typ##_ldexp_long
-#endif
-
-static PyUFuncGenericFunction ldexp_functions[] = {
-#ifdef HAVE_LDEXPF
-    HALF_ldexp,
-    FLOAT_ldexp,
-    LDEXP_LONG(HALF),
-    LDEXP_LONG(FLOAT),
-#endif
-    DOUBLE_ldexp,
-    LDEXP_LONG(DOUBLE)
-#ifdef HAVE_LDEXPL
-    ,
-    LONGDOUBLE_ldexp,
-    LDEXP_LONG(LONGDOUBLE)
-#endif
-};
-
-static const char frdoc[] =
-    "    Decompose the elements of x into mantissa and twos exponent.\n"
-    "\n"
-    "    Returns (`mantissa`, `exponent`), where `x = mantissa * 2**exponent``.\n"
-    "    The mantissa is lies in the open interval(-1, 1), while the twos\n"
-    "    exponent is a signed integer.\n"
-    "\n"
-    "    Parameters\n"
-    "    ----------\n"
-    "    x : array_like\n"
-    "        Array of numbers to be decomposed.\n"
-    "    out1: ndarray, optional\n"
-    "        Output array for the mantissa. Must have the same shape as `x`.\n"
-    "    out2: ndarray, optional\n"
-    "        Output array for the exponent. Must have the same shape as `x`.\n"
-    "\n"
-    "    Returns\n"
-    "    -------\n"
-    "    (mantissa, exponent) : tuple of ndarrays, (float, int)\n"
-    "        `mantissa` is a float array with values between -1 and 1.\n"
-    "        `exponent` is an int array which represents the exponent of 2.\n"
-    "\n"
-    "    See Also\n"
-    "    --------\n"
-    "    ldexp : Compute ``y = x1 * 2**x2``, the inverse of `frexp`.\n"
-    "\n"
-    "    Notes\n"
-    "    -----\n"
-    "    Complex dtypes are not supported, they will raise a TypeError.\n"
-    "\n"
-    "    Examples\n"
-    "    --------\n"
-    "    >>> x = np.arange(9)\n"
-    "    >>> y1, y2 = np.frexp(x)\n"
-    "    >>> y1\n"
-    "    array([ 0.   ,  0.5  ,  0.5  ,  0.75 ,  0.5  ,  0.625,  0.75 ,  0.875,\n"
-    "            0.5  ])\n"
-    "    >>> y2\n"
-    "    array([0, 1, 2, 2, 3, 3, 3, 3, 4])\n"
-    "    >>> y1 * 2**y2\n"
-    "    array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.])\n"
-    "\n";
-
-
-static char ldexp_signatures[] = {
-#ifdef HAVE_LDEXPF
-    NPY_HALF, NPY_INT, NPY_HALF,
-    NPY_FLOAT, NPY_INT, NPY_FLOAT,
-    NPY_HALF, NPY_LONG, NPY_HALF,
-    NPY_FLOAT, NPY_LONG, NPY_FLOAT,
-#endif
-    NPY_DOUBLE, NPY_INT, NPY_DOUBLE,
-    NPY_DOUBLE, NPY_LONG, NPY_DOUBLE
-#ifdef HAVE_LDEXPL
-    ,NPY_LONGDOUBLE, NPY_INT, NPY_LONGDOUBLE
-    ,NPY_LONGDOUBLE, NPY_LONG, NPY_LONGDOUBLE
-#endif
-};
-
-static const char lddoc[] =
-    "    Returns x1 * 2**x2, element-wise.\n"
-    "\n"
-    "    The mantissas `x1` and twos exponents `x2` are used to construct\n"
-    "    floating point numbers ``x1 * 2**x2``.\n"
-    "\n"
-    "    Parameters\n"
-    "    ----------\n"
-    "    x1 : array_like\n"
-    "        Array of multipliers.\n"
-    "    x2 : array_like, int\n"
-    "        Array of twos exponents.\n"
-    "    out : ndarray, optional\n"
-    "        Output array for the result.\n"
-    "\n"
-    "    Returns\n"
-    "    -------\n"
-    "    y : ndarray or scalar\n"
-    "        The result of ``x1 * 2**x2``.\n"
-    "\n"
-    "    See Also\n"
-    "    --------\n"
-    "    frexp : Return (y1, y2) from ``x = y1 * 2**y2``, inverse to `ldexp`.\n"
-    "\n"
-    "    Notes\n"
-    "    -----\n"
-    "    Complex dtypes are not supported, they will raise a TypeError.\n"
-    "\n"
-    "    `ldexp` is useful as the inverse of `frexp`, if used by itself it is\n"
-    "    more clear to simply use the expression ``x1 * 2**x2``.\n"
-    "\n"
-    "    Examples\n"
-    "    --------\n"
-    "    >>> np.ldexp(5, np.arange(4))\n"
-    "    array([  5.,  10.,  20.,  40.], dtype=float32)\n"
-    "\n"
-    "    >>> x = np.arange(6)\n"
-    "    >>> np.ldexp(*np.frexp(x))\n"
-    "    array([ 0.,  1.,  2.,  3.,  4.,  5.])\n"
-    "\n";
-
-
-static void
-InitOtherOperators(PyObject *dictionary) {
-    PyObject *f;
-    int num;
-
-    num = sizeof(frexp_functions) / sizeof(frexp_functions[0]);
-    assert(sizeof(blank_data) / sizeof(blank_data[0]) >= num);
-    f = PyUFunc_FromFuncAndData(frexp_functions, blank_data,
-                                frexp_signatures, num,
-                                1, 2, PyUFunc_None, "frexp", frdoc, 0);
-    PyDict_SetItemString(dictionary, "frexp", f);
-    Py_DECREF(f);
-
-    num = sizeof(ldexp_functions) / sizeof(ldexp_functions[0]);
-    assert(sizeof(blank_data) / sizeof(blank_data[0]) >= num);
-    f = PyUFunc_FromFuncAndData(ldexp_functions, blank_data,
-                                ldexp_signatures, num,
-                                2, 1, PyUFunc_None, "ldexp", lddoc, 0);
-    PyDict_SetItemString(dictionary, "ldexp", f);
-    Py_DECREF(f);
-
-#if defined(NPY_PY3K)
-    f = PyDict_GetItemString(dictionary, "true_divide");
-    PyDict_SetItemString(dictionary, "divide", f);
-#endif
-    return;
-}
-
 NPY_VISIBILITY_HIDDEN PyObject * npy_um_str_out = NULL;
 NPY_VISIBILITY_HIDDEN PyObject * npy_um_str_subok = NULL;
 NPY_VISIBILITY_HIDDEN PyObject * npy_um_str_array_prepare = NULL;
@@ -492,8 +316,6 @@ PyMODINIT_FUNC initumath(void)
     /* Load the ufunc operators into the array module's namespace */
     InitOperators(d);
 
-    InitOtherOperators(d);
-
     PyDict_SetItemString(d, "pi", s = PyFloat_FromDouble(NPY_PI));
     Py_DECREF(s);
     PyDict_SetItemString(d, "e", s = PyFloat_FromDouble(NPY_E));
@@ -536,6 +358,11 @@ PyMODINIT_FUNC initumath(void)
     PyModule_AddObject(m, "NZERO", PyFloat_FromDouble(NPY_NZERO));
     PyModule_AddObject(m, "NAN", PyFloat_FromDouble(NPY_NAN));
 
+#if defined(NPY_PY3K)
+    s = PyDict_GetItemString(d, "true_divide");
+    PyDict_SetItemString(d, "divide", s);
+#endif
+
     s = PyDict_GetItemString(d, "conjugate");
     s2 = PyDict_GetItemString(d, "remainder");
     /* Setup the array object's numerical structures with appropriate