path: root/numpy/core/src/umath/funcs.inc.src

/* -*- c -*- */

/*
 * This file is for the definitions of the non-c99 functions used in ufuncs.
 * All the complex ufuncs are defined here along with a smattering of real and
 * object functions.
 */

#define NPY_NO_DEPRECATED_API NPY_API_VERSION
#include "npy_pycompat.h"
#include "npy_import.h"


/*
 *****************************************************************************
 **                        PYTHON OBJECT FUNCTIONS                          **
 *****************************************************************************
 */

static PyObject *
Py_square(PyObject *o)
{
    return PyNumber_Multiply(o, o);
}

static PyObject *
Py_get_one(PyObject *NPY_UNUSED(o))
{
    return PyInt_FromLong(1);
}

static PyObject *
Py_reciprocal(PyObject *o)
{
    PyObject *one = PyInt_FromLong(1);
    PyObject *result;

    if (!one) {
        return NULL;
    }
#if defined(NPY_PY3K)
    result = PyNumber_TrueDivide(one, o);
#else
    result = PyNumber_Divide(one, o);
#endif
    Py_DECREF(one);
    return result;
}

/*
 * Define numpy version of PyNumber_Power as binary function.
 */
static PyObject *
npy_ObjectPower(PyObject *x, PyObject *y)
{
    return PyNumber_Power(x, y, Py_None);
}

/**begin repeat
 * #Kind = Max, Min#
 * #OP = Py_GE, Py_LE#
 */
static PyObject *
npy_Object@Kind@(PyObject *i1, PyObject *i2)
{
    PyObject *result;
    int cmp;

    cmp = PyObject_RichCompareBool(i1, i2, @OP@);
    if (cmp < 0) {
        return NULL;
    }
    if (cmp == 1) {
        result = i1;
    }
    else {
        result = i2;
    }
    Py_INCREF(result);
    return result;
}
/**end repeat**/

/* Emulates Python's 'a or b' behavior */
static PyObject *
npy_ObjectLogicalOr(PyObject *i1, PyObject *i2)
{
    if (i1 == NULL) {
        Py_XINCREF(i2);
        return i2;
    }
    else if (i2 == NULL) {
        Py_INCREF(i1);
        return i1;
    }
    else {
        int retcode = PyObject_IsTrue(i1);
        if (retcode == -1) {
            return NULL;
        }
        else if (retcode) {
            Py_INCREF(i1);
            return i1;
        }
        else {
            Py_INCREF(i2);
            return i2;
        }
    }
}

/* Emulates Python's 'a and b' behavior */
static PyObject *
npy_ObjectLogicalAnd(PyObject *i1, PyObject *i2)
{
    if (i1 == NULL) {
        return NULL;
    }
    else if (i2 == NULL) {
        return NULL;
    }
    else {
        int retcode = PyObject_IsTrue(i1);
        if (retcode == -1) {
            return NULL;
        }
        else if (!retcode) {
            Py_INCREF(i1);
            return i1;
        }
        else {
            Py_INCREF(i2);
            return i2;
        }
    }
}


/* Emulates Python's 'not b' behavior */
static PyObject *
npy_ObjectLogicalNot(PyObject *i1)
{
    if (i1 == NULL) {
        return NULL;
    }
    else {
        int retcode = PyObject_Not(i1);
        if (retcode == -1) {
            return NULL;
        }
        else if (retcode) {
            Py_INCREF(Py_True);
            return Py_True;
        }
        else {
            Py_INCREF(Py_False);
            return Py_False;
        }
    }
}

static PyObject *
npy_ObjectFloor(PyObject *obj) {
    PyObject *math_floor_func = NULL;

    npy_cache_import("math", "floor", &math_floor_func);
    if (math_floor_func == NULL) {
        return NULL;
    }
    return PyObject_CallFunction(math_floor_func, "O", obj);
}

static PyObject *
npy_ObjectCeil(PyObject *obj) {
    PyObject *math_ceil_func = NULL;

    npy_cache_import("math", "ceil", &math_ceil_func);
    if (math_ceil_func == NULL) {
        return NULL;
    }
    return PyObject_CallFunction(math_ceil_func, "O", obj);
}

static PyObject *
npy_ObjectTrunc(PyObject *obj) {
    PyObject *math_trunc_func = NULL;

    npy_cache_import("math", "trunc", &math_trunc_func);
    if (math_trunc_func == NULL) {
        return NULL;
    }
    return PyObject_CallFunction(math_trunc_func, "O", obj);
}

static PyObject *
npy_ObjectGCD(PyObject *i1, PyObject *i2)
{
    PyObject *gcd = NULL;

    /* use math.gcd if available, and valid on the provided types */
#if PY_VERSION_HEX >= 0x03050000
    {
        static PyObject *math_gcd_func = NULL;

        npy_cache_import("math", "gcd", &math_gcd_func);
        if (math_gcd_func == NULL) {
            return NULL;
        }
        gcd = PyObject_CallFunction(math_gcd_func, "OO", i1, i2);
        if (gcd != NULL) {
            return gcd;
        }
        /* silence errors, and fall back on pure-python gcd */
        PyErr_Clear();
    }
#endif

    /* otherwise, use our internal one, written in python */
    {
        static PyObject *internal_gcd_func = NULL;

        npy_cache_import("numpy.core._internal", "_gcd", &internal_gcd_func);
        if (internal_gcd_func == NULL) {
            return NULL;
        }
        gcd = PyObject_CallFunction(internal_gcd_func, "OO", i1, i2);
        if (gcd == NULL) {
            return NULL;
        }
        /* _gcd has some unusual behaviour regarding sign */
        return PyNumber_Absolute(gcd);
    }
}

static PyObject *
npy_ObjectLCM(PyObject *i1, PyObject *i2)
{
    /* lcm(a, b) = abs(a // gcd(a, b) * b) */

    PyObject *gcd = npy_ObjectGCD(i1, i2);
    PyObject *tmp;
    if(gcd == NULL) {
        return NULL;
    }
    /* Floor divide preserves integer types - we know the division will have
     * no remainder
     */
    tmp = PyNumber_FloorDivide(i1, gcd);
    if(tmp == NULL) {
        return NULL;
    }

    tmp = PyNumber_Multiply(tmp, i2);
    if(tmp == NULL) {
        return NULL;
    }

    /* even though we fix gcd to be positive, we need to do it again here */
    return PyNumber_Absolute(tmp);
}


static PyObject *
npy_ObjectClip(PyObject *arr, PyObject *min, PyObject *max) {
    PyObject *o = npy_ObjectMax(arr, min);
    if (o == NULL) {
        return NULL;
    }
    Py_SETREF(o, npy_ObjectMin(o, max));
    return o;
}

/*
 *****************************************************************************
 **                           COMPLEX FUNCTIONS                             **
 *****************************************************************************
 */


/*
 * Don't pass structures between functions (only pointers) because how
 * structures are passed is compiler dependent and could cause segfaults if
 * umath_ufunc_object.inc is compiled with a different compiler than an
 * extension that makes use of the UFUNC API
 */

/**begin repeat
 *
 * #ctype = npy_cfloat, npy_cdouble, npy_clongdouble#
 * #ftype = npy_float, npy_double, npy_longdouble#
 * #c = f, ,l#
 */

static void
nc_neg@c@(@ctype@ *a, @ctype@ *r)
{
    r->real = -a->real;
    r->imag = -a->imag;
    return;
}

static void
nc_pos@c@(@ctype@ *a, @ctype@ *r)
{
    r->real = +a->real;
    r->imag = +a->imag;
    return;
}

static void
nc_sqrt@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_csqrt@c@(*x);
    return;
}

static void
nc_rint@c@(@ctype@ *x, @ctype@ *r)
{
    r->real = npy_rint@c@(x->real);
    r->imag = npy_rint@c@(x->imag);
}

static void
nc_log@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_clog@c@(*x);
    return;
}

static void
nc_log1p@c@(@ctype@ *x, @ctype@ *r)
{
    @ftype@ l = npy_hypot@c@(x->real + 1,x->imag);
    r->imag = npy_atan2@c@(x->imag, x->real + 1);
    r->real = npy_log@c@(l);
    return;
}

static void
nc_exp@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_cexp@c@(*x);
    return;
}

static void
nc_exp2@c@(@ctype@ *x, @ctype@ *r)
{
    @ctype@ a;
    a.real = x->real*NPY_LOGE2@c@;
    a.imag = x->imag*NPY_LOGE2@c@;
    nc_exp@c@(&a, r);
    return;
}

static void
nc_expm1@c@(@ctype@ *x, @ctype@ *r)
{
    @ftype@ a = npy_exp@c@(x->real);
    r->real = a*npy_cos@c@(x->imag) - 1.0@c@;
    r->imag = a*npy_sin@c@(x->imag);
    return;
}

static void
nc_pow@c@(@ctype@ *a, @ctype@ *b, @ctype@ *r)
{
   *r = npy_cpow@c@(*a, *b);
    return;
}

static void
nc_acos@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_cacos@c@(*x);
    return;
}

static void
nc_acosh@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_cacosh@c@(*x);
    return;
}

static void
nc_asin@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_casin@c@(*x);
    return;
}


static void
nc_asinh@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_casinh@c@(*x);
    return;
}

static void
nc_atan@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_catan@c@(*x);
    return;
}

static void
nc_atanh@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_catanh@c@(*x);
    return;
}

static void
nc_cos@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_ccos@c@(*x);
    return;
}

static void
nc_cosh@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_ccosh@c@(*x);
    return;
}

static void
nc_log10@c@(@ctype@ *x, @ctype@ *r)
{
    nc_log@c@(x, r);
    r->real *= NPY_LOG10E@c@;
    r->imag *= NPY_LOG10E@c@;
    return;
}

static void
nc_log2@c@(@ctype@ *x, @ctype@ *r)
{
    nc_log@c@(x, r);
    r->real *= NPY_LOG2E@c@;
    r->imag *= NPY_LOG2E@c@;
    return;
}

static void
nc_sin@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_csin@c@(*x);
    return;
}

static void
nc_sinh@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_csinh@c@(*x);
    return;
}

static void
nc_tan@c@(@ctype@ *x, @ctype@ *r)
{
   *r = npy_ctan@c@(*x);
   return;
}

static void
nc_tanh@c@(@ctype@ *x, @ctype@ *r)
{
    *r = npy_ctanh@c@(*x);
    return;
}

/**end repeat**/