[PYTHON MATRIX-SIG] Bug in LinearAlgebra

[Jim Hugunin]

> only to uncover
> a problem with multiply.reduce. Check this out:
...

This is a really terrible bug.  It only shows up when doing multiply (or
divide) as a reduction on an array of complex numbers, but still this is
ugly.  The problem lies in the low-level C functions in fast_umath and
umath modules.  You should replace the following two functions with the
forms given below to fix it (in fast_umathmodule.c and umathmodule.c). 
This one's bad enough I thought I'd post a patch right away.

A gold star to anyone who can explain what stupid mistake I made when
writing these functions the first time.

static void CFLOAT_multiply(char **args, int *dimensions, int *steps, void
*func) {
  int i, is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
  char *i1=args[0], *i2=args[1], *op=args[2];
  float tmp;
  for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
    tmp=((float *)i1)[0] * ((float *)i2)[0] - ((float *)i1)[1] * ((float
*)i2)[1]; ((float *)op)[1]=((float *)i1)[1] * ((float *)i2)[0] + ((float
*)i1)[0] * ((float *)i2)[1];
	((float *)op)[0]=tmp;
  }
}

static void CDOUBLE_multiply(char **args, int *dimensions, int *steps, void
*func) {
  int i, is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
  char *i1=args[0], *i2=args[1], *op=args[2];
  double tmp;
  for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
    tmp=((double *)i1)[0] * ((double *)i2)[0] - ((double *)i1)[1] *
((double *)i2)[1]; ((double *)op)[1]=((double *)i1)[1] * ((double *)i2)[0]
+ ((double *)i1)[0] * ((double *)i2)[1];
	((double *)op)[0]=tmp;
  }
}


> Here's the fixed (but not really tested) version of determinant:
> 
> def determinate(a):
>     _assertRank2(a)
>     _assertSquareness(a)
>     ev = eigenvalues(a)
>     return Numeric.product(ev)

If somebody more LA familiar than myself can confirm and test this I'll
make the change in the standard distribution.

-Jim

=================
MATRIX-SIG  - SIG on Matrix Math for Python

send messages to: matrix-sig@python.org
administrivia to: matrix-sig-request@python.org
=================