/* -*- c -*- */

#include "Python.h"
#include "numpy/arrayobject.h"
#define _UMATHMODULE
#include "numpy/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

#define _SIGN1(x) ((x) > 0 ? 1 : ((x) < 0 ? -1 : 0))
#define _SIGN2(x) ((x) == 0 ? 0 : 1)
#define _SIGNC(x) (((x).real > 0) ? 1 : ((x).real < 0 ? -1 : ((x).imag > 0 ? 1 : ((x).imag < 0) ? -1 : 0)))
/**begin repeat
#TYPE=BYTE,SHORT,INT,LONG,LONGLONG,FLOAT,DOUBLE,LONGDOUBLE,UBYTE,USHORT,UINT,ULONG,ULONGLONG#
#typ=byte,short,int,long,longlong,float,double,longdouble,ubyte,ushort,uint,ulong,ulonglong#
#func=_SIGN1*8,_SIGN2*5#
 */
static void
@TYPE@_sign(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];
        @typ@ t1;
	for(i=0; i<n; i++, i1+=is1, op+=os) {
                t1 = *((@typ@ *)i1);
		*((@typ@ *)op) = (@typ@) @func@(t1);
	}
}
/**end repeat**/

/**begin repeat
#TYPE=CFLOAT,CDOUBLE,CLONGDOUBLE#
#typ=cfloat,cdouble,clongdouble#
 */
static void
@TYPE@_sign(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];
        @typ@ t1;
	for(i=0; i<n; i++, i1+=is1, op+=os) {
                t1 = *((@typ@ *)i1);
		(*((@typ@ *)op)).real = _SIGNC(t1);
                (*((@typ@ *)op)).imag = 0;
	}
}
/**end repeat**/

#undef _SIGN1
#undef _SIGN2
#undef _SIGNC


static void
OBJECT_sign(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];
        PyObject *t1, *zero, *res;
        zero = PyInt_FromLong(0);
	for(i=0; i<n; i++, i1+=is1, op+=os) {
                t1 = *((PyObject **)i1);
                res = PyInt_FromLong((long) PyObject_Compare(t1, zero));
		*((PyObject **)op) = res;
	}
        Py_DECREF(zero);
}


/**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");
}