[Numpy-svn] r6391 - branches/coremath/numpy/core/src

numpy-svn@scip... numpy-svn@scip...
Wed Feb 18 11:28:39 CST 2009


Author: cdavid
Date: 2009-02-18 11:28:30 -0600 (Wed, 18 Feb 2009)
New Revision: 6391

Added:
   branches/coremath/numpy/core/src/npy_math.c.src
   branches/coremath/numpy/core/src/npy_math.h
Log:
Start working on core math library itself.

Added: branches/coremath/numpy/core/src/npy_math.c.src
===================================================================
--- branches/coremath/numpy/core/src/npy_math.c.src	2009-02-18 17:28:03 UTC (rev 6390)
+++ branches/coremath/numpy/core/src/npy_math.c.src	2009-02-18 17:28:30 UTC (rev 6391)
@@ -0,0 +1,304 @@
+/*
+ * vim:syntax=c
+ * A small module to implement missing C99 math capabilities required by numpy
+ *
+ * Please keep this independant of python ! Only basic types (npy_longdouble)
+ * can be used, otherwise, pure C, without any use of Python facilities
+ *
+ * How to add a function to this section
+ * -------------------------------------
+ *
+ * Say you want to add `foo`, these are the steps and the reasons for them.
+ *
+ * 1) Add foo to the appropriate list in the configuration system. The
+ *    lists can be found in numpy/core/setup.py lines 63-105. Read the
+ *    comments that come with them, they are very helpful.
+ *
+ * 2) The configuration system will define a macro HAVE_FOO if your function
+ *    can be linked from the math library. The result can depend on the
+ *    optimization flags as well as the compiler, so can't be known ahead of
+ *    time. If the function can't be linked, then either it is absent, defined
+ *    as a macro, or is an intrinsic (hardware) function.
+ *
+ *    i) Undefine any possible macros:
+ *
+ *    #ifdef foo
+ *    #undef foo
+ *    #endif
+ *
+ *    ii) Avoid as much as possible to declare any function here. Declaring
+ *    functions is not portable: some platforms define some function inline
+ *    with a non standard identifier, for example, or may put another
+ *    idendifier which changes the calling convention of the function. If you
+ *    really have to, ALWAYS declare it for the one platform you are dealing
+ *    with:
+ *
+ *    Not ok:
+ *        double exp(double a);
+ *
+ *    Ok:
+ *        #ifdef SYMBOL_DEFINED_WEIRD_PLATFORM
+ *        double exp(double);
+ *        #endif
+ */
+
+#include <Python.h>
+#include <math.h>
+
+#include "config.h"
+#include "npy_math.h"
+
+/*
+ *****************************************************************************
+ **                     BASIC MATH FUNCTIONS                                **
+ *****************************************************************************
+ */
+
+/* Original code by Konrad Hinsen.  */
+#ifndef HAVE_EXPM1
+double expm1(double x)
+{
+    double u = exp(x);
+    if (u == 1.0) {
+        return x;
+    } else if (u-1.0 == -1.0) {
+        return -1;
+    } else {
+        return (u-1.0) * x/log(u);
+    }
+}
+#endif
+
+#ifndef HAVE_LOG1P
+double log1p(double x)
+{
+    double u = 1. + x;
+    if (u == 1.0) {
+        return x;
+    } else {
+        return log(u) * x / (u - 1);
+    }
+}
+#endif
+
+#ifndef HAVE_HYPOT
+double hypot(double x, double y)
+{
+    double yx;
+
+    x = fabs(x);
+    y = fabs(y);
+    if (x < y) {
+        double temp = x;
+        x = y;
+        y = temp;
+    }
+    if (x == 0.)
+        return 0.;
+    else {
+        yx = y/x;
+        return x*sqrt(1.+yx*yx);
+    }
+}
+#endif
+
+#ifndef HAVE_ACOSH
+double acosh(double x)
+{
+    return 2*log(sqrt((x+1.0)/2)+sqrt((x-1.0)/2));
+}
+#endif
+
+#ifndef HAVE_ASINH
+double asinh(double xx)
+{
+    double x, d;
+    int sign;
+    if (xx < 0.0) {
+        sign = -1;
+        x = -xx;
+    }
+    else {
+        sign = 1;
+        x = xx;
+    }
+    if (x > 1e8) {
+        d = x;
+    } else {
+        d = sqrt(x*x + 1);
+    }
+    return sign*log1p(x*(1.0 + x/(d+1)));
+}
+#endif
+
+#ifndef HAVE_ATANH
+double atanh(double x)
+{
+    if (x > 0) {
+        return -0.5*log1p(-2.0*x/(1.0 + x));
+    }
+    else {
+        return 0.5*log1p(2.0*x/(1.0 - x));
+    }
+}
+#endif
+
+#ifndef HAVE_RINT
+double rint(double x)
+{
+    double y, r;
+
+    y = floor(x);
+    r = x - y;
+
+    if (r > 0.5) goto rndup;
+
+    /* Round to nearest even */
+    if (r==0.5) {
+        r = y - 2.0*floor(0.5*y);
+        if (r==1.0) {
+        rndup:
+            y+=1.0;
+        }
+    }
+    return y;
+}
+#endif
+
+#ifndef HAVE_TRUNC
+double trunc(double x)
+{
+    return x < 0 ? ceil(x) : floor(x);
+}
+#endif
+
+#ifndef HAVE_EXP2
+#define LOG2 0.69314718055994530943
+double exp2(double x)
+{
+    return exp(LOG2*x);
+}
+#undef LOG2
+#endif
+
+#ifndef HAVE_LOG2
+#define INVLOG2 1.4426950408889634074
+double log2(double x)
+{
+    return INVLOG2*log(x);
+}
+#undef INVLOG2
+#endif
+
+/*
+ *****************************************************************************
+ **                     IEEE 754 FPU HANDLING                               **
+ *****************************************************************************
+ */
+#if !defined(HAVE_DECL_ISNAN)
+    # define isnan(x) ((x) != (x))
+#endif
+
+/* VS 2003 with /Ox optimizes (x)-(x) to 0, which is not IEEE compliant. So we
+ * force (x) + (-x), which seems to work. */
+#if !defined(HAVE_DECL_ISFINITE)
+    # define isfinite(x) !isnan((x) + (-x))
+#endif
+
+#if !defined(HAVE_DECL_ISINF)
+#define isinf(x) (!isfinite(x) && !isnan(x))
+#endif
+
+#if !defined(HAVE_DECL_SIGNBIT)
+    #include "_signbit.c"
+    # define signbit(x) \
+              (sizeof (x) == sizeof (long double) ? signbit_ld (x) \
+               : sizeof (x) == sizeof (double) ? signbit_d (x) \
+               : signbit_f (x))
+
+static int signbit_f (float x)
+{
+    return signbit_d((double)x);
+}
+
+static int signbit_ld (long double x)
+{
+    return signbit_d((double)x);
+}
+#endif
+
+/*
+ * if C99 extensions not available then define dummy functions that use the
+ * double versions for
+ *
+ * sin, cos, tan
+ * sinh, cosh, tanh,
+ * fabs, floor, ceil, rint, trunc
+ * sqrt, log10, log, exp, expm1
+ * asin, acos, atan,
+ * asinh, acosh, atanh
+ *
+ * hypot, atan2, pow, fmod, modf
+ *
+ * We assume the above are always available in their double versions.
+ *
+ * NOTE: some facilities may be available as macro only  instead of functions.
+ * For simplicity, we define our own functions and undef the macros. We could
+ * instead test for the macro, but I am lazy to do that for now.
+ */
+
+/**begin repeat
+ * #type = npy_longdouble, float#
+ * #TYPE = NPY_LONGDOUBLE, FLOAT#
+ * #c = l,f#
+ * #C = L,F#
+ */
+
+/**begin repeat1
+ * #kind = sin,cos,tan,sinh,cosh,tanh,fabs,floor,ceil,rint,trunc,sqrt,log10,
+ *         log,exp,expm1,asin,acos,atan,asinh,acosh,atanh,log1p,exp2,log2#
+ * #KIND = SIN,COS,TAN,SINH,COSH,TANH,FABS,FLOOR,CEIL,RINT,TRUNC,SQRT,LOG10,
+ *         LOG,EXP,EXPM1,ASIN,ACOS,ATAN,ASINH,ACOSH,ATANH,LOG1P,EXP2,LOG2#
+ */
+
+#ifdef @kind@@c@
+#undef @kind@@c@
+#endif
+#ifndef HAVE_@KIND@@C@
+@type@ @kind@@c@(@type@ x)
+{
+    return (@type@) @kind@((double)x);
+}
+#endif
+
+/**end repeat1**/
+
+/**begin repeat1
+ * #kind = atan2,hypot,pow,fmod#
+ * #KIND = ATAN2,HYPOT,POW,FMOD#
+ */
+#ifdef @kind@@c@
+#undef @kind@@c@
+#endif
+#ifndef HAVE_@KIND@@C@
+@type@ @kind@@c@(@type@ x, @type@ y)
+{
+    return (@type@) @kind@((double)x, (double) y);
+}
+#endif
+/**end repeat1**/
+
+#ifdef modf@c@
+#undef modf@c@
+#endif
+#ifndef HAVE_MODF@C@
+@type@ modf@c@(@type@ x, @type@ *iptr)
+{
+    double niptr;
+    double y = modf((double)x, &niptr);
+    *iptr = (@type@) niptr;
+    return (@type@) y;
+}
+#endif
+
+/**end repeat**/

Added: branches/coremath/numpy/core/src/npy_math.h
===================================================================
--- branches/coremath/numpy/core/src/npy_math.h	2009-02-18 17:28:03 UTC (rev 6390)
+++ branches/coremath/numpy/core/src/npy_math.h	2009-02-18 17:28:30 UTC (rev 6391)
@@ -0,0 +1,121 @@
+#ifndef __NPY_MATH_C99_H_
+#define __NPY_MATH_C99_H_
+
+#include <math.h>
+#include <numpy/npy_common.h>
+/*
+ * C99 double math funcs
+ */
+double npy_expm1(double x);
+double npy_log1p(double x);
+double npy_hypot(double x, double y);
+double npy_acosh(double x);
+double npy_asinh(double xx);
+double npy_atanh(double x);
+double npy_rint(double x);
+double npy_trunc(double x);
+double npy_exp2(double x);
+double npy_log2(double x);
+
+/*
+ * IEEE 754 fpu handling. Those are guaranteed to be macros
+ */
+#ifndef NPY_HAVE_DECL_ISNAN
+        #define npy_isnan(x) _npy_isnan((x))
+#else
+        #define npy_isnan(x) isnan((x))
+#endif
+
+#ifndef NPY_HAVE_DECL_ISFINITE
+        #define npy_isfinite(x) _npy_isfinite((x))
+#else
+        #define npy_isfinite(x) isfinite((x))
+#endif
+
+#ifndef NPY_HAVE_DECL_ISFINITE
+        #define npy_isinf(x) _npy_isinf((x))
+#else
+        #define npy_isinf(x) isinf((x))
+#endif
+
+#ifndef NPY_HAVE_DECL_SIGNBIT
+        #define npy_signbit(x) _npy_signbit((x))
+#else
+        #define npy_signbit(x) signbit((x))
+#endif
+
+/*
+ * float C99 math functions
+ */
+
+float npy_sinf(float x);
+float npy_cosf(float x);
+float npy_tanf(float x);
+float npy_sinhf(float x);
+float npy_coshf(float x);
+float npy_tanhf(float x);
+float npy_fabsf(float x);
+float npy_floorf(float x);
+float npy_ceilf(float x);
+float npy_rintf(float x);
+float npy_truncf(float x);
+float npy_sqrtf(float x);
+float npy_log10f(float x);
+float npy_logf(float x);
+float npy_expf(float x);
+float npy_expm1f(float x);
+float npy_asinf(float x);
+float npy_acosf(float x);
+float npy_atanf(float x);
+float npy_asinhf(float x);
+float npy_acoshf(float x);
+float npy_atanhf(float x);
+float npy_log1pf(float x);
+float npy_exp2f(float x);
+float npy_log2f(float x);
+
+float npy_atan2f(float x, float y);
+float npy_hypotf(float x, float y);
+float npy_powf(float x, float y);
+float npy_fmodf(float x, float y);
+
+float npy_modff(float x, float* y);
+
+/*
+ * float C99 math functions
+ */
+
+npy_longdouble npy_sinl(npy_longdouble x);
+npy_longdouble npy_cosl(npy_longdouble x);
+npy_longdouble npy_tanl(npy_longdouble x);
+npy_longdouble npy_sinhl(npy_longdouble x);
+npy_longdouble npy_coshl(npy_longdouble x);
+npy_longdouble npy_tanhl(npy_longdouble x);
+npy_longdouble npy_fabsl(npy_longdouble x);
+npy_longdouble npy_floorl(npy_longdouble x);
+npy_longdouble npy_ceill(npy_longdouble x);
+npy_longdouble npy_rintl(npy_longdouble x);
+npy_longdouble npy_truncl(npy_longdouble x);
+npy_longdouble npy_sqrtl(npy_longdouble x);
+npy_longdouble npy_log10l(npy_longdouble x);
+npy_longdouble npy_logl(npy_longdouble x);
+npy_longdouble npy_expl(npy_longdouble x);
+npy_longdouble npy_expm1l(npy_longdouble x);
+npy_longdouble npy_asinl(npy_longdouble x);
+npy_longdouble npy_acosl(npy_longdouble x);
+npy_longdouble npy_atanl(npy_longdouble x);
+npy_longdouble npy_asinhl(npy_longdouble x);
+npy_longdouble npy_acoshl(npy_longdouble x);
+npy_longdouble npy_atanhl(npy_longdouble x);
+npy_longdouble npy_log1pl(npy_longdouble x);
+npy_longdouble npy_exp2l(npy_longdouble x);
+npy_longdouble npy_log2l(npy_longdouble x);
+
+npy_longdouble npy_atan2l(npy_longdouble x, npy_longdouble y);
+npy_longdouble npy_hypotl(npy_longdouble x, npy_longdouble y);
+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);
+
+#endif



More information about the Numpy-svn mailing list