Moved scipy directory to numpy

1 files changed, 1847 insertions, 0 deletions

diff --git a/numpy/base/src/umathmodule.c.src b/numpy/base/src/umathmodule.c.src
new file mode 100644
index 000000000..5096f3361
--- /dev/null
+++ b/numpy/base/src/umathmodule.c.src
@@ -0,0 +1,1847 @@
+/* -*- c -*- */
+
+#include "Python.h"
+#include "scipy/arrayobject.h"
+#define _UMATHMODULE
+#include "scipy/ufuncobject.h"
+#include "abstract.h"
+#include <math.h>
+
+
+/* A whole slew of basic math functions are provided originally by Konrad Hinsen. */
+
+#if !defined(__STDC__) && !defined(_MSC_VER)
+extern double fmod (double, double);
+extern double frexp (double, int *);
+extern double ldexp (double, int);
+extern double modf (double, double *);
+#endif
+#ifndef M_PI
+#define M_PI 3.14159265358979323846264338328
+#endif
+
+#ifndef HAVE_INVERSE_HYPERBOLIC
+static double acosh(double x)
+{
+	return log(x + sqrt((x-1.0)*(x+1.0)));
+}
+
+static double asinh(double xx)
+{
+	double x;
+	int sign;
+	if (xx < 0.0) {
+		sign = -1;
+		x = -xx;
+	}
+	else {
+		sign = 1;
+		x = xx;
+	}
+	return sign*log(x + sqrt(x*x+1.0));
+}
+
+static double atanh(double x)
+{
+	return 0.5*log((1.0+x)/(1.0-x));
+}
+#endif
+
+#ifdef HAVE_HYPOT
+#if !defined(NeXT) && !defined(_MSC_VER)
+extern double hypot(double, double);
+#endif
+#else
+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
+
+
+
+/* Define isnan, isinf, isfinite, signbit if needed */
+/* Use fpclassify if possible */
+/* isnan, isinf --
+       these will use macros and then fpclassify if available before
+       defaulting to a dumb convert-to-double version...
+
+   isfinite -- define a macro if not already available
+   signbit -- if macro available use it, otherwise define a function
+               and a dumb convert-to-double version for other types.
+*/
+
+#if defined(fpclassify)
+
+#if !defined(isnan)
+#define isnan(x) (fpclassify(x) == FP_NAN)
+#endif
+#if !defined(isinf)
+#define isinf(x) (fpclassify(x) == FP_INFINITE)
+#endif
+
+#else  /* check to see if already have a function like this */
+
+#if !defined(HAVE_ISNAN)
+
+#if !defined(isnan)
+#include "_isnan.c"
+#endif
+#endif /* HAVE_ISNAN */
+
+#if !defined(HAVE_ISINF)
+#if !defined(isinf)
+#define isinf(x) (!isnan((x)) && isnan((x)-(x)))
+#endif
+#endif /* HAVE_ISINF */
+
+#endif /* defined(fpclassify) */
+
+
+/* Define signbit if needed */
+#if !defined(signbit)
+#include "_signbit.c"
+#endif
+
+
+/* Now defined the extended type macros */
+
+#if !defined(isnan)
+
+#if !defined(HAVE_LONGDOUBLE_FUNCS) || !defined(HAVE_ISNAN)
+#define isnanl(x) isnan((double)(x))
+#endif
+
+#if !defined(HAVE_FLOAT_FUNCS) || !defined(HAVE_ISNAN)
+#define isnanf(x) isnan((double)(x))
+#endif
+
+#else /* !defined(isnan) */
+
+#define isnanl(x) isnan((x))
+#define isnanf(x) isnan((x))
+
+#endif /* !defined(isnan) */
+
+
+#if !defined(isinf)
+
+#if !defined(HAVE_LONGDOUBLE_FUNCS) || !defined(HAVE_ISINF)
+#define isinfl(x) (!isnanl((x)) && isnanl((x)-(x)))
+#endif
+
+#if !defined(HAVE_FLOAT_FUNCS) || !defined(HAVE_ISINF)
+#define isinff(x) (!isnanf((x)) && isnanf((x)-(x)))
+#endif
+
+#else /* !defined(isinf) */
+
+#define isinfl(x) isinf((x))
+#define isinff(x) isinf((x))
+
+#endif /* !defined(isinf) */
+
+
+#if !defined(signbit)
+#define signbitl(x) ((longdouble) signbit((double)(x)))
+#define signbitf(x) ((float) signbit((double) (x)))
+#else
+#define signbitl(x) signbit((x))
+#define signbitf(x) signbit((x))
+#endif
+
+#if !defined(isfinite)
+#define isfinite(x) (!(isinf((x)) || isnan((x))))
+#endif
+#define isfinitef(x) (!(isinff((x)) || isnanf((x))))
+#define isfinitel(x) (!(isinfl((x)) || isnanl((x))))
+
+
+/* First, the C functions that do the real work */
+
+/* if C99 extensions not availble
+
+then define dummy functions that use the double versions for
+
+sin, cos, tan
+sinh, cosh, tanh,
+fabs, floor, ceil, fmod, sqrt, log10, log, exp, fabs
+asin, acos, atan,
+asinh, acosh, atanh
+
+hypot, atan2, pow
+
+*/
+
+/**begin repeat
+
+#kind=(sin,cos,tan,sinh,cosh,tanh,fabs,floor,ceil,sqrt,log10,log,exp,asin,acos,atan)*2#
+#typ=longdouble*16, float*16#
+#c=l*16,f*16#
+#TYPE=LONGDOUBLE*16, FLOAT*16#
+*/
+#ifndef HAVE_@TYPE@_FUNCS
+@typ@ @kind@@c@(@typ@ x) {
+	return (@typ@) @kind@((double)x);
+}
+#endif
+/**end repeat**/
+
+/**begin repeat
+
+#kind=(atan2,hypot,pow,fmod)*2#
+#typ=longdouble*4, float*4#
+#c=l*4,f*4#
+#TYPE=LONGDOUBLE*4,FLOAT*4#
+*/
+#ifndef HAVE_@TYPE@_FUNCS
+@typ@ @kind@@c@(@typ@ x, @typ@ y) {
+	return (@typ@) @kind@((double)x, (double) y);
+}
+#endif
+/**end repeat**/
+
+/**begin repeat
+#kind=modf*2#
+#typ=longdouble, float#
+#c=l,f#
+#TYPE=LONGDOUBLE, FLOAT#
+*/
+#ifndef HAVE_@TYPE@_FUNCS
+@typ@ modf@c@(@typ@ x, @typ@ *iptr) {
+        double nx, niptr, y;
+        nx = (double) x;
+        y = modf(nx, &niptr);
+        *iptr = (@typ@) niptr;
+        return (@typ@) y;
+}
+#endif
+/**end repeat**/
+
+
+#if !defined(HAVE_INVERSE_HYPERBOLIC_FLOAT)
+#ifdef HAVE_FLOAT_FUNCS
+static float acoshf(float x)
+{
+	return logf(x + sqrtf((x-1.0)*(x+1.0)));
+}
+
+static float asinhf(float xx)
+{
+	float x;
+	int sign;
+	if (xx < 0.0) {
+		sign = -1;
+		x = -xx;
+	}
+	else {
+		sign = 1;
+		x = xx;
+	}
+	return sign*logf(x + sqrtf(x*x+1.0));
+}
+
+static float atanhf(float x)
+{
+	return 0.5*logf((1.0+x)/(1.0-x));
+}
+#else
+static float acoshf(float x)
+{
+  return (float)acosh((double)(x));
+}
+
+static float asinhf(float x)
+{
+  return (float)asinh((double)(x));
+}
+
+static float atanhf(float x)
+{
+  return (float)atanh((double)(x));
+}
+#endif
+#endif
+
+
+#if !defined(HAVE_INVERSE_HYPERBOLIC_LONGDOUBLE)
+#ifdef HAVE_LONGDOUBLE_FUNCS
+static longdouble acoshl(longdouble x)
+{
+	return logl(x + sqrtl((x-1.0)*(x+1.0)));
+}
+
+static longdouble asinhl(longdouble xx)
+{
+	longdouble x;
+	int sign;
+	if (xx < 0.0) {
+		sign = -1;
+		x = -xx;
+	}
+	else {
+		sign = 1;
+		x = xx;
+	}
+	return sign*logl(x + sqrtl(x*x+1.0));
+}
+
+static longdouble atanhl(longdouble x)
+{
+	return 0.5*logl((1.0+x)/(1.0-x));
+}
+#else
+static longdouble acoshl(longdouble x)
+{
+  return (longdouble)acosh((double)(x));
+}
+
+static longdouble asinhl(longdouble x)
+{
+  return (longdouble)asinh((double)(x));
+}
+
+static longdouble atanhl(longdouble x)
+{
+  return (longdouble)atanh((double)(x));
+}
+#endif
+#endif
+
+
+
+
+/* Don't pass structures between functions (only pointers) because how
+   structures are passed is compiler dependent and could cause
+   segfaults if ufuncobject.c is compiled with a different compiler
+   than an extension that makes use of the UFUNC API
+*/
+
+/**begin repeat
+
+#typ=float, double, longdouble#
+#c=f,,l#
+*/
+
+/* constants */
+static c@typ@ nc_1@c@ = {1., 0.};
+static c@typ@ nc_half@c@ = {0.5, 0.};
+static c@typ@ nc_i@c@ = {0., 1.};
+static c@typ@ nc_i2@c@ = {0., 0.5};
+/*
+static c@typ@ nc_mi@c@ = {0., -1.};
+static c@typ@ nc_pi2@c@ = {M_PI/2., 0.};
+*/
+
+static void
+nc_sum@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
+{
+	r->real = a->real + b->real;
+	r->imag = a->imag + b->imag;
+	return;
+}
+
+static void
+nc_diff@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
+{
+	r->real = a->real - b->real;
+	r->imag = a->imag - b->imag;
+	return;
+}
+
+static void
+nc_neg@c@(c@typ@ *a, c@typ@ *r)
+{
+	r->real = -a->real;
+	r->imag = -a->imag;
+	return;
+}
+
+static void
+nc_prod@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
+{
+	@typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
+	r->real = ar*br - ai*bi;
+	r->imag = ar*bi + ai*br;
+	return;
+}
+
+static void
+nc_quot@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
+{
+
+	@typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
+	@typ@ d = br*br + bi*bi;
+	r->real = (ar*br + ai*bi)/d;
+	r->imag = (ai*br - ar*bi)/d;
+	return;
+}
+
+static void
+nc_floor_quot@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
+{
+	@typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
+	@typ@ d = br*br + bi*bi;
+	r->real = floor@c@((ar*br + ai*bi)/d);
+	r->imag = 0;
+	return;
+}
+
+static void
+nc_sqrt@c@(c@typ@ *x, c@typ@ *r)
+{
+	@typ@ s,d;
+	if (x->real == 0. && x->imag == 0.)
+		*r = *x;
+	else {
+		s = sqrt@c@(0.5*(fabs@c@(x->real) + hypot@c@(x->real,x->imag)));
+		d = 0.5*x->imag/s;
+		if (x->real > 0.) {
+			r->real = s;
+			r->imag = d;
+		}
+		else if (x->imag >= 0.) {
+			r->real = d;
+			r->imag = s;
+		}
+		else {
+			r->real = -d;
+			r->imag = -s;
+		}
+	}
+	return;
+}
+
+static void
+nc_log@c@(c@typ@ *x, c@typ@ *r)
+{
+	@typ@ l = hypot@c@(x->real,x->imag);
+	r->imag = atan2@c@(x->imag, x->real);
+	r->real = log@c@(l);
+	return;
+}
+
+static void
+nc_exp@c@(c@typ@ *x, c@typ@ *r)
+{
+	@typ@ a = exp@c@(x->real);
+	r->real = a*cos@c@(x->imag);
+	r->imag = a*sin@c@(x->imag);
+	return;
+}
+
+static void
+nc_pow@c@(c@typ@ *a, c@typ@ *b, c@typ@ *r)
+{
+	@typ@ ar=a->real, br=b->real, ai=a->imag, bi=b->imag;
+
+	if (br == 0. && bi == 0.) {
+		r->real = 1.;
+		r->imag = 0.;
+	}
+	else if (ar == 0. && ai == 0.) {
+		r->real = 0.;
+		r->imag = 0.;
+	}
+	else {
+		nc_log@c@(a, r);
+		nc_prod@c@(r, b, r);
+		nc_exp@c@(r, r);
+	}
+	return;
+}
+
+
+static void
+nc_prodi@c@(c@typ@ *x, c@typ@ *r)
+{
+	r->real = -x->imag;
+	r->imag = x->real;
+	return;
+}
+
+
+static void
+nc_acos@c@(c@typ@ *x, c@typ@ *r)
+{
+	nc_prod@c@(x,x,r);
+	nc_diff@c@(&nc_1@c@, r, r);
+	nc_sqrt@c@(r, r);
+	nc_prodi@c@(r, r);
+	nc_sum@c@(x, r, r);
+	nc_log@c@(r, r);
+	nc_prodi@c@(r, r);
+	nc_neg@c@(r, r);
+	return;
+	/* return nc_neg(nc_prodi(nc_log(nc_sum(x,nc_prod(nc_i,
+	   nc_sqrt(nc_diff(nc_1,nc_prod(x,x))))))));
+	*/
+}
+
+static void
+nc_acosh@c@(c@typ@ *x, c@typ@ *r)
+{
+	nc_prod@c@(x, x, r);
+	nc_diff@c@(&nc_1@c@, r, r);
+	nc_sqrt@c@(r, r);
+	nc_prodi@c@(r, r);
+	nc_sum@c@(x, r, r);
+	nc_log@c@(r, r);
+	return;
+	/*
+	  return nc_log(nc_sum(x,nc_prod(nc_i,
+	  nc_sqrt(nc_diff(nc_1,nc_prod(x,x))))));
+	*/
+}
+
+static void
+nc_asin@c@(c@typ@ *x, c@typ@ *r)
+{
+	c@typ@ a, *pa=&a;
+	nc_prod@c@(x, x, r);
+	nc_diff@c@(&nc_1@c@, r, r);
+	nc_sqrt@c@(r, r);
+	nc_prodi@c@(x, pa);
+	nc_sum@c@(pa, r, r);
+	nc_log@c@(r, r);
+	nc_prodi@c@(r, r);
+	nc_neg@c@(r, r);
+	return;
+	/*
+	      return nc_neg(nc_prodi(nc_log(nc_sum(nc_prod(nc_i,x),
+	      nc_sqrt(nc_diff(nc_1,nc_prod(x,x)))))));
+	*/
+}
+
+
+static void
+nc_asinh@c@(c@typ@ *x, c@typ@ *r)
+{
+	nc_prod@c@(x, x, r);
+	nc_sum@c@(&nc_1@c@, r, r);
+	nc_sqrt@c@(r, r);
+	nc_diff@c@(r, x, r);
+	nc_log@c@(r, r);
+	nc_neg@c@(r, r);
+	return;
+	/*
+	  return nc_neg(nc_log(nc_diff(nc_sqrt(nc_sum(nc_1,nc_prod(x,x))),x)));
+	*/
+}
+
+static void
+nc_atan@c@(c@typ@ *x, c@typ@ *r)
+{
+	c@typ@ a, *pa=&a;
+	nc_diff@c@(&nc_i@c@, x, pa);
+	nc_sum@c@(&nc_i@c@, x, r);
+	nc_quot@c@(r, pa, r);
+	nc_log@c@(r,r);
+	nc_prod@c@(&nc_i2@c@, r, r);
+	return;
+	/*
+	  return nc_prod(nc_i2,nc_log(nc_quot(nc_sum(nc_i,x),nc_diff(nc_i,x))));
+	*/
+}
+
+static void
+nc_atanh@c@(c@typ@ *x, c@typ@ *r)
+{
+	c@typ@ a, *pa=&a;
+	nc_diff@c@(&nc_1@c@, x, r);
+	nc_sum@c@(&nc_1@c@, x, pa);
+	nc_quot@c@(pa, r, r);
+	nc_log@c@(r, r);
+	nc_prod@c@(&nc_half@c@, r, r);
+	return;
+	/*
+	  return nc_prod(nc_half,nc_log(nc_quot(nc_sum(nc_1,x),nc_diff(nc_1,x))));
+	*/
+}
+
+static void
+nc_cos@c@(c@typ@ *x, c@typ@ *r)
+{
+	@typ@ xr=x->real, xi=x->imag;
+	r->real = cos@c@(xr)*cosh@c@(xi);
+	r->imag = -sin@c@(xr)*sinh@c@(xi);
+	return;
+}
+
+static void
+nc_cosh@c@(c@typ@ *x, c@typ@ *r)
+{
+	@typ@ xr=x->real, xi=x->imag;
+	r->real = cos(xi)*cosh(xr);
+	r->imag = sin(xi)*sinh(xr);
+	return;
+}
+
+
+#define M_LOG10_E 0.434294481903251827651128918916605082294397
+
+static void
+nc_log10@c@(c@typ@ *x, c@typ@ *r)
+{
+	nc_log@c@(x, r);
+	r->real *= M_LOG10_E;
+	r->imag *= M_LOG10_E;
+	return;
+}
+
+static void
+nc_sin@c@(c@typ@ *x, c@typ@ *r)
+{
+	@typ@ xr=x->real, xi=x->imag;
+	r->real = sin@c@(xr)*cosh@c@(xi);
+	r->imag = cos@c@(xr)*sinh@c@(xi);
+	return;
+}
+
+static void
+nc_sinh@c@(c@typ@ *x, c@typ@ *r)
+{
+	@typ@ xr=x->real, xi=x->imag;
+	r->real = cos@c@(xi)*sinh@c@(xr);
+	r->imag = sin@c@(xi)*cosh@c@(xr);
+	return;
+}
+
+static void
+nc_tan@c@(c@typ@ *x, c@typ@ *r)
+{
+	@typ@ sr,cr,shi,chi;
+	@typ@ rs,is,rc,ic;
+	@typ@ d;
+	@typ@ xr=x->real, xi=x->imag;
+	sr = sin@c@(xr);
+	cr = cos@c@(xr);
+	shi = sinh(xi);
+	chi = cosh(xi);
+	rs = sr*chi;
+	is = cr*shi;
+	rc = cr*chi;
+	ic = -sr*shi;
+	d = rc*rc + ic*ic;
+	r->real = (rs*rc+is*ic)/d;
+	r->imag = (is*rc-rs*ic)/d;
+	return;
+}
+
+static void
+nc_tanh@c@(c@typ@ *x, c@typ@ *r)
+{
+	@typ@ si,ci,shr,chr;
+	@typ@ rs,is,rc,ic;
+	@typ@ d;
+	@typ@ xr=x->real, xi=x->imag;
+	si = sin@c@(xi);
+	ci = cos@c@(xi);
+	shr = sinh@c@(xr);
+	chr = cosh@c@(xr);
+	rs = ci*shr;
+	is = si*chr;
+	rc = ci*chr;
+	ic = si*shr;
+	d = rc*rc + ic*ic;
+	r->real = (rs*rc+is*ic)/d;
+	r->imag = (is*rc-rs*ic)/d;
+	return;
+}
+
+/**end repeat**/
+
+
+/**begin repeat
+
+#TYPE=(BOOL, BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE)*2#
+#OP=||, +*13, ^, -*13#
+#kind=add*14, subtract*14#
+#typ=(Bool, byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble)*2#
+*/
+
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		*((@typ@ *)op)=*((@typ@ *)i1) @OP@ *((@typ@ *)i2);
+	}
+}
+
+/**end repeat**/
+
+/**begin repeat
+
+#TYPE=(CFLOAT, CDOUBLE, CLONGDOUBLE)*2#
+#OP=+*3,-*3#
+#kind=add*3,subtract*3#
+#typ=(float, double, longdouble)*2#
+
+*/
+
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		((@typ@ *)op)[0]=((@typ@ *)i1)[0] @OP@ ((@typ@ *)i2)[0];
+		((@typ@ *)op)[1]=((@typ@ *)i1)[1] @OP@ ((@typ@ *)i2)[1];
+	}
+}
+
+/**end repeat**/
+
+
+/** Routines borrowed from numarray **/
+
+/* The following routine is used in the event of a detected integer *
+** divide by zero so that a floating divide by zero is generated.   *
+** This is done since Numeric uses the floating point exception    *
+** sticky bits to detect errors. The last bit is an attempt to      *
+** prevent optimization of the divide by zero away, the output value *
+** should always be 0                                               *
+*/
+
+/* These should really be altered to just set the corresponding bit
+   in the floating point status flag.  Need to figure out how to do that
+   on all the platforms...
+*/
+
+static int numeric_zero = 0.0;
+
+#if !defined(generate_divbyzero_error)
+static void generate_divbyzero_error(void) {
+	double dummy;
+	dummy = 1./numeric_zero;
+	return;
+}
+#endif
+
+#if !defined(generate_overflow_error)
+static double numeric_two = 2.0;
+static void generate_overflow_error(void) {
+	double dummy;
+	dummy = pow(numeric_two,1000);
+	return;
+}
+#endif
+
+
+static int ulonglong_overflow(ulonglong a, ulonglong b)
+{
+        ulonglong ah, al, bh, bl, w, x, y, z;
+
+#if SIZEOF_LONGLONG == 64
+	ah = (a >> 32);
+	al = (a & 0xFFFFFFFFL);
+	bh = (b >> 32);
+	bl = (b & 0xFFFFFFFFL);
+#elif SIZEOF_LONGLONG == 128
+	ah = (a >> 64);
+	al = (a & 0xFFFFFFFFFFFFFFFFL);
+	bh = (b >> 64);
+	bl = (b & 0xFFFFFFFFFFFFFFFFL);
+#else
+	ah = al = bh = bl = 0;
+#endif
+
+        /* 128-bit product:  z*2**64 + (x+y)*2**32 + w  */
+	w = al*bl;
+	x = bh*al;
+	y = ah*bl;
+	z = ah*bh;
+
+	/* *c = ((x + y)<<32) + w; */
+#if SIZEOF_LONGLONG == 64
+	return z || (x>>32) || (y>>32) ||
+		(((x & 0xFFFFFFFFL) + (y & 0xFFFFFFFFL) + (w >> 32)) >> 32);
+#elif SIZEOF_LONGLONG == 128
+	return z || (x>>64) || (y>>64) ||
+		(((x & 0xFFFFFFFFFFFFFFFFL) + (y & 0xFFFFFFFFFFFFFFFFL) + (w >> 64)) >> 64);
+#else
+	return 0;
+#endif
+
+}
+
+static int slonglong_overflow(longlong a0, longlong b0)
+{
+	ulonglong a, b;
+        ulonglong ah, al, bh, bl, w, x, y, z;
+
+        /* Convert to non-negative quantities */
+	if (a0 < 0) { a = -a0; } else { a = a0; }
+	if (b0 < 0) { b = -b0; } else { b = b0; }
+
+
+#if SIZEOF_LONGLONG == 64
+	ah = (a >> 32);
+	al = (a & 0xFFFFFFFFL);
+	bh = (b >> 32);
+	bl = (b & 0xFFFFFFFFL);
+#elif SIZEOF_LONGLONG == 128
+	ah = (a >> 64);
+	al = (a & 0xFFFFFFFFFFFFFFFFL);
+	bh = (b >> 64);
+	bl = (b & 0xFFFFFFFFFFFFFFFFL);
+#else
+	ah = al = bh = bl = 0;
+#endif
+
+	w = al*bl;
+	x = bh*al;
+	y = ah*bl;
+	z = ah*bh;
+
+        /*
+	   ulonglong c = ((x + y)<<32) + w;
+	   if ((a0 < 0) ^ (b0 < 0))
+	   *c = -c;
+	   else
+	   *c = c
+	   */
+
+#if SIZEOF_LONGLONG == 64
+	return z || (x>>31) || (y>>31) ||
+		(((x & 0xFFFFFFFFL) + (y & 0xFFFFFFFFL) + (w >> 32)) >> 31);
+#elif SIZEOF_LONGLONG == 128
+	return z || (x>>63) || (y>>63) ||
+		(((x & 0xFFFFFFFFFFFFFFFFL) + (y & 0xFFFFFFFFFFFFFFFFL) + (w >> 64)) >> 63);
+#else
+	return 0;
+#endif
+}
+
+/** end direct numarray code **/
+
+static void
+BOOL_multiply(char **args, intp *dimensions, intp *steps, void *func) {
+	register intp i;
+	intp is1=steps[0], is2=steps[1], os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for (i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		*((Bool *)op) = *((Bool *)i1) && *((Bool *)i2);
+	}
+}
+
+/**begin repeat
+
+#TYP= UBYTE,USHORT,UINT, ULONG#
+#typ= ubyte, ushort, uint, ulong#
+#bigtyp= int, int, double, double#
+*/
+
+static void
+@TYP@_multiply(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0], is2=steps[1], os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	@bigtyp@ temp;
+	for (i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		temp = (@bigtyp@)(*((@typ@ *)i1)) * (@bigtyp@)(*((@typ@ *)i2));
+		if (temp > MAX_@TYP@)
+			generate_overflow_error();
+		*((@typ@ *)op) = temp;
+	}
+}
+
+/**end repeat**/
+
+static void
+ULONGLONG_multiply(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0], is2=steps[1], os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	ulonglong temp;
+	for (i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		temp = *((ulonglong *)i1) * *((ulonglong *)i2);
+		if (ulonglong_overflow(*((ulonglong *)i1), *((ulonglong *)i2)))
+			generate_overflow_error();
+		*((ulonglong *)op) = temp;
+	}
+}
+
+/**begin repeat
+
+#TYP= BYTE,SHORT,INT, LONG#
+#typ= byte, short, int, long#
+#bigtyp= int, int, double, double#
+*/
+
+static void
+@TYP@_multiply(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0], is2=steps[1], os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	@bigtyp@ temp;
+	for (i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		temp = (@bigtyp@)*((@typ@ *)i1) * (@bigtyp@)*((@typ@ *)i2);
+		if (temp > MAX_@TYP@)
+			generate_overflow_error();
+		else if (temp < MIN_@TYP@)
+			generate_overflow_error();
+		*((@typ@ *)op) = temp;
+	}
+}
+
+/**end repeat**/
+
+static void
+LONGLONG_multiply(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0], is2=steps[1], os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	longlong temp;
+	for (i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		temp = *((longlong *)i1) * *((longlong *)i2);
+		if (slonglong_overflow(*((longlong *)i1), *((longlong *)i2)))
+			generate_overflow_error();
+		*((longlong *)op) = temp;
+	}
+}
+
+
+/**begin repeat
+
+#TYP=FLOAT,DOUBLE,LONGDOUBLE#
+#typ=float,double,longdouble#
+*/
+static void
+@TYP@_multiply(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	/* fprintf(stderr, "Multiplying %d elements of type @typ@\n", n);
+	fprintf(stderr, "args= %p, %p, %p\n", i1, i2, op);
+	fprintf(stderr, "steps=%d, %d, %d\n", is1, is2, os); */
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		*((@typ@ *)op)=*((@typ@ *)i1) * *((@typ@ *)i2);
+	}
+}
+/**end repeat**/
+
+
+/**begin repeat
+
+#TYP=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG#
+#typ=char, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong#
+#otyp=float*4, double*6#
+*/
+static void
+@TYP@_divide(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i, is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		if (*((@typ@ *)i2)==0) {
+			generate_divbyzero_error();
+			*((@typ@ *)op)=0;
+		}
+		else {
+			*((@typ@ *)op)= *((@typ@ *)i1) / *((@typ@ *)i2);
+		}
+	}
+}
+static void
+@TYP@_true_divide(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i, is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		if (*((@typ@ *)i2)==0) {
+			generate_divbyzero_error();
+			*((@otyp@ *)op)=0;
+		}
+		else {
+			*((@otyp@ *)op)=				\
+				*((@typ@ *)i1) / (double)*((@typ@ *)i2);
+		}
+	}
+}
+#define @TYP@_floor_divide @TYP@_divide
+/**end repeat**/
+
+/**begin repeat
+
+#TYP=(FLOAT,DOUBLE,LONGDOUBLE)*2#
+#typ=(float,double,longdouble)*2#
+#kind=divide*3, true_divide*3#
+*/
+static void
+@TYP@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i, is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		*((@typ@ *)op)=*((@typ@ *)i1) / *((@typ@ *)i2);
+	}
+}
+/**end repeat**/
+
+/**begin repeat
+
+#TYP=FLOAT,DOUBLE,LONGDOUBLE#
+#typ=float,double,longdouble#
+#c=f,,l#
+*/
+static void
+@TYP@_floor_divide(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i, is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		*((@typ@ *)op)=floor@c@(*((@typ@ *)i1) / *((@typ@ *)i2));
+	}
+}
+/**end repeat**/
+
+
+/**begin repeat
+
+#TYP=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG#
+#typ=char, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong#
+#btyp=float*4, double*6#
+*/
+static void
+@TYP@_power(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i, is1=steps[0],is2=steps[1];
+	register intp os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	@btyp@ x, y, v;
+        @typ@ z;
+
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		x = *((@typ@ *)i1);
+		y = *((@typ@ *)i2);
+		z = (@typ@) y;
+		if ((x < 0.0) && (y != z))  v = 1.0/numeric_zero;
+		else v = pow(x,y);
+		*((@typ@ *)op) = (@typ@) v;
+	}
+}
+/**end repeat**/
+
+/**begin repeat
+
+#TYP=UBYTE, BYTE, SHORT, USHORT, INT, UINT, LONG, ULONG, LONGLONG, ULONGLONG, FLOAT, DOUBLE, LONGDOUBLE#
+#typ=ubyte, char, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble#
+*/
+static void
+@TYP@_conjugate(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i, is1=steps[0], os=steps[1], n=dimensions[0];
+	char *i1=args[0], *op=args[1];
+	for(i=0; i<n; i++, i1+=is1, op+=os) {
+		*((@typ@ *)op)=*((@typ@ *)i1);
+	}
+}
+/**end repeat**/
+
+/**begin repeat
+
+#TYP=CFLOAT, CDOUBLE, CLONGDOUBLE#
+#typ=float, double, longdouble#
+*/
+static void
+@TYP@_conjugate(char **args, intp *dimensions, intp *steps, void *func) {
+	register intp i, is1=steps[0], os=steps[1], n=dimensions[0];
+	char *i1=args[0], *op=args[1];
+
+	for(i=0; i<n; i++, i1+=is1, op+=os) {
+		((@typ@ *)op)[0]=((@typ@ *)i1)[0];
+		((@typ@ *)op)[1]=-(((@typ@ *)i1)[1]);
+	}
+}
+/**end repeat**/
+
+
+/**begin repeat
+
+#TYPE=BOOL,UBYTE,USHORT,UINT,ULONG,ULONGLONG#
+#typ=Bool, ubyte, ushort, uint, ulong, ulonglong#
+*/
+static void
+@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i, n;
+	intp is1=steps[0], os=steps[1];
+	char *i1=args[0], *op=args[1];
+
+    n=dimensions[0];
+
+	for(i=0; i<n; i++, i1+=is1, op+=os) {
+		*((@typ@ *)op) = *((@typ@*)i1);
+	}
+}
+/**end repeat**/
+
+/**begin repeat
+
+#TYPE=BYTE,SHORT,INT,LONG,LONGLONG,FLOAT,DOUBLE,LONGDOUBLE#
+#typ=byte, short, int, long, longlong, float, double, longdouble#
+*/
+static void
+@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i, n;
+	intp is1=steps[0], os=steps[1];
+	char *i1=args[0], *op=args[1];
+
+	n=dimensions[0];
+	for(i=0; i<n; i++, i1+=is1, op+=os) {
+		*((@typ@ *)op) = *((@typ@ *)i1) < 0 ? -*((@typ@ *)i1) : *((@typ@ *)i1);
+	}
+}
+/**end repeat**/
+
+/**begin repeat
+   #TYPE=CFLOAT,CDOUBLE,CLONGDOUBLE#
+   #typ= float, double, longdouble#
+   #c= f,,l#
+*/
+static void
+@TYPE@_absolute(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i, n;
+	register intp is1=steps[0], os=steps[1];
+	char *i1=args[0], *op=args[1];
+	n=dimensions[0];
+	for(i=0; i<n; i++, i1+=is1, op+=os) {
+		*((@typ@ *)op) = (@typ@)sqrt@c@(((@typ@ *)i1)[0]*((@typ@ *)i1)[0] + ((@typ@ *)i1)[1]*((@typ@ *)i1)[1]);
+	}
+}
+/**end repeat**/
+
+/**begin repeat
+
+#kind=greater, greater_equal, less, less_equal, equal, not_equal, logical_and, logical_or, bitwise_and, bitwise_or, bitwise_xor#
+#OP=>, >=, <, <=, ==, !=, &&, ||, &, |, ^#
+**/
+static void
+BOOL_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	Bool in1, in2;
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		in1 = (*((Bool *)i1) != 0);
+		in2 = (*((Bool *)i2) != 0);
+		*((Bool *)op)= in1 @OP@ in2;
+	}
+}
+/**end repeat**/
+
+/**begin repeat
+
+#TYPE=(BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE)*4#
+#OP= >*13, >=*13, <*13, <=*13#
+#typ=(byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble)*4#
+#kind= greater*13, greater_equal*13, less*13, less_equal*13#
+*/
+
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		*((Bool *)op)=*((@typ@ *)i1) @OP@ *((@typ@ *)i2);
+	}
+}
+/**end repeat**/
+
+
+/**begin repeat
+#TYPE=(CFLOAT,CDOUBLE,CLONGDOUBLE)*4#
+#OP= >*3, >=*3, <*3, <=*3#
+#typ=(cfloat, cdouble, clongdouble)*4#
+#kind= greater*3, greater_equal*3, less*3, less_equal*3#
+*/
+
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		if (((@typ@ *)i1)->real == ((@typ@ *)i2)->real)
+			*((Bool *)op)=((@typ@ *)i1)->imag @OP@ \
+				((@typ@ *)i2)->imag;
+		else
+			*((Bool *)op)=((@typ@ *)i1)->real @OP@ \
+				((@typ@ *)i2)->real;
+	}
+}
+/**end repeat**/
+
+
+/**begin repeat
+#TYPE=(BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE)*4#
+#typ=(byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble)*4#
+#OP= ==*13, !=*13, &&*13, ||*13#
+#kind=equal*13, not_equal*13, logical_and*13, logical_or*13#
+*/
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		*((Bool *)op) = *((@typ@ *)i1) @OP@ *((@typ@ *)i2);
+	}
+}
+/**end repeat**/
+
+
+/**begin repeat
+
+#TYPE=(CFLOAT, CDOUBLE, CLONGDOUBLE)*4#
+#typ=(float, double, longdouble)*4#
+#OP= ==*3, !=*3, &&*3, ||*3#
+#OP2= &&*3, ||*3, &&*3, ||*3#
+#kind=equal*3, not_equal*3, logical_and*3, logical_or*3#
+*/
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		*((Bool *)op) = (*((@typ@ *)i1) @OP@ *((@typ@ *)i2)) @OP2@ (*((@typ@ *)i1+1) @OP@ *((@typ@ *)i2+1));
+	}
+}
+/**end repeat**/
+
+
+/** OBJECT comparison for OBJECT arrays **/
+
+/**begin repeat
+
+#kind=greater, greater_equal, less, less_equal, equal, not_equal#
+#op=GT, GE, LT, LE, EQ, NE#
+*/
+static void
+OBJECT_@kind@(char **args, intp *dimensions, intp *steps, void *func) {
+	register intp i, is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		*((Bool *)op)=PyObject_RichCompareBool(*((PyObject **)i1),
+						       *((PyObject **)i2),
+						       Py_@op@);
+	}
+}
+/**end repeat**/
+
+/**begin repeat
+
+#TYPE=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE#
+#typ=byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble#
+*/
+static void
+@TYPE@_negative(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],os=steps[1], n=dimensions[0];
+	char *i1=args[0], *op=args[1];
+	for(i=0; i<n; i++, i1+=is1, op+=os) {
+		*((@typ@ *)op) = - *((@typ@ *)i1);
+	}
+}
+/**end repeat**/
+
+#define BOOL_negative BOOL_logical_not
+
+
+/**begin repeat
+#TYPE=BOOL,BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE#
+#typ=Bool, byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble#
+*/
+static void
+@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],os=steps[1], n=dimensions[0];
+	char *i1=args[0], *op=args[1];
+	for(i=0; i<n; i++, i1+=is1, op+=os) {
+		*((Bool *)op) = ! *((@typ@ *)i1);
+	}
+}
+/**end repeat**/
+
+/**begin repeat
+#TYPE=CFLOAT,CDOUBLE,CLONGDOUBLE#
+#typ=cfloat, cdouble, clongdouble#
+*/
+static void
+@TYPE@_logical_not(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],os=steps[1], n=dimensions[0];
+	char *i1=args[0], *op=args[1];
+	for(i=0; i<n; i++, i1+=is1, op+=os) {
+		*((Bool *)op) = ! (((@typ@ *)i1)->real ||	\
+				   ((@typ@ *)i1)->imag);
+	}
+}
+/**end repeat**/
+
+
+
+
+/**begin repeat
+#TYPE=BYTE,SHORT,INT,LONG,LONGLONG#
+#typ=byte, short, int, long, longlong#
+#ftyp=float*2,double*2,longdouble*1#
+#c=f*2,,,l*1#
+*/
+static void
+@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	register @typ@ ix,iy, tmp;
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		ix = *((@typ@ *)i1);
+		iy = *((@typ@ *)i2);
+		if (iy == 0 || ix == 0) {
+			if (iy == 0) generate_divbyzero_error();
+			*((@typ@ *)op) = 0;
+		}
+		else if ((ix > 0) == (iy > 0)) {
+			*((@typ@ *)op) = ix % iy;
+		}
+		else {  /* handle mixed case the way Python does */
+			tmp = ix % iy;
+			if (tmp) tmp += iy;
+			*((@typ@ *)op)= tmp;
+		}
+	}
+}
+/**end repeat**/
+
+/**begin repeat
+#TYPE=UBYTE,USHORT,UINT,ULONG,ULONGLONG#
+#typ=ubyte, ushort, uint, ulong, ulonglong#
+*/
+static void
+@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	register @typ@ ix,iy;
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		ix = *((@typ@ *)i1);
+		iy = *((@typ@ *)i2);
+		if (iy == 0) {
+			generate_divbyzero_error();
+			*((@typ@ *)op) = 0;
+		}
+		*((@typ@ *)op) = ix % iy;
+	}
+}
+/**end repeat**/
+
+/**begin repeat
+#TYPE=FLOAT,DOUBLE,LONGDOUBLE#
+#typ=float,double,longdouble#
+#c=f,,l#
+*/
+static void
+@TYPE@_remainder(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	@typ@ x, y, res;
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		x = *((@typ@ *)i1);
+		y = *((@typ@ *)i2);
+		res = x - floor@c@(x/y)*y;
+		*((@typ@ *)op)= res;
+	}
+}
+/**end repeat**/
+
+
+/**begin repeat
+
+#TYPE=(BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG)*6#
+#typ=(byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong)*6#
+#OP= %*10, &*10, |*10, ^*10, <<*10, >>*10#
+#kind=fmod*10, bitwise_and*10, bitwise_or*10, bitwise_xor*10, left_shift*10, right_shift*10#
+
+*/
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		*((@typ@ *)op)=*((@typ@ *)i1) @OP@ *((@typ@ *)i2);
+	}
+}
+/**end repeat**/
+
+
+/**begin repeat
+   #TYPE=BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG#
+   #typ=byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong#
+*/
+static void
+@TYPE@_invert(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0], os=steps[1], n=dimensions[0];
+	char *i1=args[0], *op=args[1];
+	for(i=0; i<n; i++, i1+=is1, op+=os) {
+		*((@typ@ *)op) = ~ *((@typ@*)i1);
+	}
+}
+/**end repeat**/
+
+static void
+BOOL_invert(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0], os=steps[1], n=dimensions[0];
+	char *i1=args[0], *op=args[1];
+	for(i=0; i<n; i++, i1+=is1, op+=os) {
+		*((Bool *)op) = (*((Bool *)i1) ? FALSE : TRUE);
+	}
+}
+
+
+/**begin repeat
+#TYPE=BOOL,BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE#
+#typ=Bool, byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble#
+
+*/
+static void
+@TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		*((Bool *)op)=(*((@typ@ *)i1) || *((@typ@ *)i2)) && !(*((@typ@ *)i1) && *((@typ@ *)i2));
+	}
+}
+/**end repeat**/
+
+
+/**begin repeat
+#TYPE=CFLOAT,CDOUBLE,CLONGDOUBLE#
+#typ=cfloat, cdouble, clongdouble#
+*/
+static void
+@TYPE@_logical_xor(char **args, intp *dimensions, intp *steps, void *func)
+{
+	Bool p1, p2;
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		p1 = ((@typ@ *)i1)->real || ((@typ@ *)i1)->imag;
+		p2 = ((@typ@ *)i2)->real || ((@typ@ *)i2)->imag;
+		*((Bool *)op)= (p1 || p2) && !(p1 && p2);
+	}
+}
+/**end repeat**/
+
+
+
+/**begin repeat
+
+#TYPE=(BOOL,BYTE,UBYTE,SHORT,USHORT,INT,UINT,LONG,ULONG,LONGLONG,ULONGLONG,FLOAT,DOUBLE,LONGDOUBLE)*2#
+#OP= >*14, <*14#
+#typ=(Bool, byte, ubyte, short, ushort, int, uint, long, ulong, longlong, ulonglong, float, double, longdouble)*2#
+#kind= maximum*14, minimum*14#
+*/
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+		*((@typ@ *)op)=*((@typ@ *)i1) @OP@ *((@typ@ *)i2) ? *((@typ@ *)i1) : *((@typ@ *)i2);
+	}
+}
+/**end repeat**/
+
+/**begin repeat
+
+#TYPE=(CFLOAT,CDOUBLE,CLONGDOUBLE)*2#
+#OP= >*3, <*3#
+#typ=(cfloat, cdouble, clongdouble)*2#
+#kind= maximum*3, minimum*3#
+*/
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+	register intp i;
+	intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
+	char *i1=args[0], *i2=args[1], *op=args[2];
+        @typ@ *i1c, *i2c;
+	for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+                i1c = (@typ@ *)i1;
+                i2c = (@typ@ *)i2;
+                if ((i1c->real @OP@ i2c->real) ||                       \
+                    ((i1c->real==i2c->real) && (i1c->imag @OP@ i2c->imag)))
+                        memcpy(op, i1, sizeof(@typ@));
+                else
+                        memcpy(op, i2, sizeof(@typ@));
+	}
+}
+/**end repeat**/
+
+
+
+/*** isinf, isinf, isfinite, signbit ***/
+/**begin repeat
+#kind=isnan*3, isinf*3, isfinite*3, signbit*3#
+#TYPE=(FLOAT, DOUBLE, LONGDOUBLE)*4#
+#typ=(float, double, longdouble)*4#
+#c=(f,,l)*4#
+*/
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+        register intp i;
+        intp is=steps[0], os=steps[1], n=dimensions[0];
+        char *ip=args[0], *op=args[1];
+        for(i=0; i<n; i++, ip+=is, op+=os) {
+                *((Bool *)op) = (Bool) (@kind@@c@(*((@typ@ *)ip)) != 0);
+        }
+}
+/**end repeat**/
+
+
+/**begin repeat
+#kind=isnan*3, isinf*3, isfinite*3#
+#TYPE=(CFLOAT, CDOUBLE, CLONGDOUBLE)*3#
+#typ=(float, double, longdouble)*3#
+#c=(f,,l)*3#
+#OP=||*6,&&*3#
+*/
+static void
+@TYPE@_@kind@(char **args, intp *dimensions, intp *steps, void *func)
+{
+        register intp i;
+        intp is=steps[0], os=steps[1], n=dimensions[0];
+        char *ip=args[0], *op=args[1];
+        for(i=0; i<n; i++, ip+=is, op+=os) {
+                *((Bool *)op) = @kind@@c@(((@typ@ *)ip)[0]) @OP@        \
+                        @kind@@c@(((@typ@ *)ip)[1]);
+        }
+}
+/**end repeat**/
+
+
+
+
+/****** modf ****/
+
+/**begin repeat
+#TYPE=FLOAT, DOUBLE, LONGDOUBLE#
+#typ=float, double, longdouble#
+#c=f,,l#
+*/
+static void
+@TYPE@_modf(char **args, intp *dimensions, intp *steps, void *func)
+{
+        register intp i;
+        intp is1=steps[0],os1=steps[1],os2=steps[2],n=dimensions[0];
+        char *i1=args[0], *op1=args[1], *op2=args[2];
+        @typ@ x1, y1, y2;
+        for (i=0; i<n; i++, i1+=is1, op1+=os1, op2+=os2) {
+                x1 = *((@typ@ *)i1);
+                y1 = modf@c@(x1, &y2);
+                *((@typ@ *)op1) = y1;
+                *((@typ@ *)op2) = y2;
+        }
+}
+/**end repeat**/
+
+#define HAVE_DOUBLE_FUNCS
+/**begin repeat
+#TYPE=FLOAT, DOUBLE, LONGDOUBLE#
+#typ=float, double, longdouble#
+#c=f,,l#
+*/
+#ifdef HAVE_@TYPE@_FUNCS
+static void
+@TYPE@_frexp(char **args, intp *dimensions, intp *steps, void *func)
+{
+        register intp i;
+        intp is1=steps[0],os1=steps[1],os2=steps[2],n=dimensions[0];
+        char *i1=args[0], *op1=args[1], *op2=args[2];
+        @typ@ x1, y1;
+        int y2;
+        for (i=0; i<n; i++, i1+=is1, op1+=os1, op2+=os2) {
+                x1 = *((@typ@ *)i1);
+                y1 = frexp@c@(x1, &y2);
+                *((@typ@ *)op1) = y1;
+                *((int *) op2) = y2;
+        }
+}
+
+static void
+@TYPE@_ldexp(char **args, intp *dimensions, intp *steps, void *func)
+{
+        register intp i;
+        intp is1=steps[0],is2=steps[1],os=steps[2],n=dimensions[0];
+        char *i1=args[0], *i2=args[1], *op=args[2];
+        @typ@ x1, y1;
+        int x2;
+        for (i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
+                x1 = *((@typ@ *)i1);
+                x2 = *((int *)i2);
+                y1 = ldexp@c@(x1, x2);
+                *((@typ@ *)op) = y1;
+        }
+}
+#endif
+/**end repeat**/
+#undef HAVE_DOUBLE_FUNCS
+
+
+static PyUFuncGenericFunction frexp_functions[] = {
+#ifdef HAVE_FLOAT_FUNCS
+        FLOAT_frexp,
+#endif
+        DOUBLE_frexp
+#ifdef HAVE_LONGDOUBLE_FUNCS
+        ,LONGDOUBLE_frexp
+#endif
+};
+
+static void * blank3_data[] = { (void *)NULL, (void *)NULL, (void *)NULL};
+static char frexp_signatures[] = {
+#ifdef HAVE_FLOAT_FUNCS
+        PyArray_FLOAT, PyArray_FLOAT, PyArray_INT,
+#endif
+        PyArray_DOUBLE, PyArray_DOUBLE, PyArray_INT
+#ifdef HAVE_LONGDOUBLE_FUNCS
+        ,PyArray_LONGDOUBLE, PyArray_LONGDOUBLE, PyArray_INT
+#endif
+};
+
+
+static PyUFuncGenericFunction ldexp_functions[] = {
+#ifdef HAVE_FLOAT_FUNCS
+        FLOAT_ldexp,
+#endif
+        DOUBLE_ldexp
+#ifdef HAVE_LONGDOUBLE_FUNCS
+        ,LONGDOUBLE_ldexp
+#endif
+};
+
+static char ldexp_signatures[] = {
+#ifdef HAVE_FLOAT_FUNCS
+        PyArray_FLOAT, PyArray_INT, PyArray_FLOAT,
+#endif
+        PyArray_DOUBLE, PyArray_INT, PyArray_DOUBLE
+#ifdef HAVE_LONGDOUBLE_FUNCS
+        ,PyArray_LONGDOUBLE, PyArray_INT, PyArray_LONGDOUBLE
+#endif
+};
+
+
+
+#include "__umath_generated.c"
+
+
+#include "ufuncobject.c"
+
+#include "__ufunc_api.c"
+
+static double
+pinf_init(void)
+{
+    double mul = 1e10;
+    double tmp = 0.0;
+    double pinf;
+
+    pinf = mul;
+    for (;;) {
+        pinf *= mul;
+        if (pinf == tmp) break;
+        tmp = pinf;
+    }
+    return pinf;
+}
+
+static double
+pzero_init(void)
+{
+    double div = 1e10;
+    double tmp = 0.0;
+    double pinf;
+
+    pinf = div;
+    for (;;) {
+        pinf /= div;
+        if (pinf == tmp) break;
+        tmp = pinf;
+    }
+    return pinf;
+}
+
+/* Less automated additions to the ufuncs */
+
+static void
+InitOtherOperators(PyObject *dictionary) {
+        PyObject *f;
+        int num=1;
+
+#ifdef HAVE_LONGDOUBLE_FUNCS
+        num += 1;
+#endif
+#ifdef HAVE_FLOAT_FUNCS
+        num += 1;
+#endif
+        f = PyUFunc_FromFuncAndData(frexp_functions, blank3_data,
+				    frexp_signatures, num,
+                                    1, 2, PyUFunc_None, "frexp",
+                                    "Split the number, x, into a normalized"\
+				    " fraction (y1) and exponent (y2)",0);
+        PyDict_SetItemString(dictionary, "frexp", f);
+        Py_DECREF(f);
+
+        f = PyUFunc_FromFuncAndData(ldexp_functions, blank3_data, ldexp_signatures, num,
+                                    2, 1, PyUFunc_None, "ldexp",
+                                    "Compute y = x1 * 2**x2.",0);
+        PyDict_SetItemString(dictionary, "ldexp", f);
+        Py_DECREF(f);
+        return;
+}
+
+static struct PyMethodDef methods[] = {
+        {"frompyfunc", (PyCFunction) ufunc_frompyfunc,
+         METH_VARARGS | METH_KEYWORDS, doc_frompyfunc},
+	{"update_use_defaults", (PyCFunction) ufunc_update_use_defaults,
+	 METH_VARARGS , NULL},
+	{NULL,		NULL, 0}		/* sentinel */
+};
+
+DL_EXPORT(void) initumath(void) {
+	PyObject *m, *d, *s, *s2, *c_api;
+	double pinf, pzero, mynan;
+
+	/* Create the module and add the functions */
+	m = Py_InitModule("umath", methods);
+
+	/* Import the array */
+	if (import_array() < 0) return;
+
+	/* Initialize the types */
+	if (PyType_Ready(&PyUFunc_Type) < 0)
+                return;
+
+	/* Add some symbolic constants to the module */
+	d = PyModule_GetDict(m);
+
+	c_api = PyCObject_FromVoidPtr((void *)PyUFunc_API, NULL);
+	if (PyErr_Occurred()) goto err;
+	PyDict_SetItemString(d, "_UFUNC_API", c_api);
+	Py_DECREF(c_api);
+	if (PyErr_Occurred()) goto err;
+
+	s = PyString_FromString("0.4.0");
+	PyDict_SetItemString(d, "__version__", s);
+	Py_DECREF(s);
+
+	/* Load the ufunc operators into the array module's namespace */
+	InitOperators(d);
+
+        InitOtherOperators(d);
+
+	PyDict_SetItemString(d, "pi", s = PyFloat_FromDouble(M_PI));
+	Py_DECREF(s);
+	PyDict_SetItemString(d, "e", s = PyFloat_FromDouble(exp(1.0)));
+	Py_DECREF(s);
+
+#define ADDCONST(str) PyModule_AddIntConstant(m, #str, UFUNC_##str)
+#define ADDSCONST(str) PyModule_AddStringConstant(m, "UFUNC_" #str, UFUNC_##str)
+
+        ADDCONST(ERR_IGNORE);
+        ADDCONST(ERR_WARN);
+        ADDCONST(ERR_CALL);
+        ADDCONST(ERR_RAISE);
+        ADDCONST(ERR_DEFAULT);
+
+        ADDCONST(SHIFT_DIVIDEBYZERO);
+        ADDCONST(SHIFT_OVERFLOW);
+        ADDCONST(SHIFT_UNDERFLOW);
+        ADDCONST(SHIFT_INVALID);
+
+        ADDCONST(FPE_DIVIDEBYZERO);
+        ADDCONST(FPE_OVERFLOW);
+        ADDCONST(FPE_UNDERFLOW);
+        ADDCONST(FPE_INVALID);
+
+        ADDSCONST(PYVALS_NAME);
+
+#undef ADDCONST
+#undef ADDSCONST
+        PyModule_AddIntConstant(m, "UFUNC_BUFSIZE_DEFAULT", (long)PyArray_BUFSIZE);
+
+        pinf = pinf_init();
+        pzero = pzero_init();
+        mynan = pinf / pinf;
+
+        PyModule_AddObject(m, "PINF", PyFloat_FromDouble(pinf));
+        PyModule_AddObject(m, "NINF", PyFloat_FromDouble(-pinf));
+        PyModule_AddObject(m, "PZERO", PyFloat_FromDouble(pzero));
+        PyModule_AddObject(m, "NZERO", PyFloat_FromDouble(-pzero));
+        PyModule_AddObject(m, "NAN", PyFloat_FromDouble(mynan));
+
+	s = PyDict_GetItemString(d, "conjugate");
+	s2 = PyDict_GetItemString(d, "remainder");
+	/* Setup the array object's numerical structures with appropriate
+	   ufuncs in d*/
+	PyArray_SetNumericOps(d);
+
+	PyDict_SetItemString(d, "conj", s);
+	PyDict_SetItemString(d, "mod", s2);
+
+ err:
+	/* Check for errors */
+	if (PyErr_Occurred())
+		Py_FatalError("can't initialize module umath");
+}