[Python-checkins] CVS: python/dist/src/Objects complexobject.c,2.34,2.35

2001年3月18日 00:22:00 -0800

Update of /cvsroot/python/python/dist/src/Objects
In directory usw-pr-cvs1:/tmp/cvs-serv4744/python/dist/src/Objects
Modified Files:
	complexobject.c 
Log Message:
SF bug [ #409448 ] Complex division is braindead
http://sourceforge.net/tracker/?func=detail&aid=409448&group_id=5470&atid=105470
Now less braindead. Also added test_complex.py, which doesn't test much, but
fails without this patch.
Index: complexobject.c
===================================================================
RCS file: /cvsroot/python/python/dist/src/Objects/complexobject.c,v
retrieving revision 2.34
retrieving revision 2.35
diff -C2 -r2.34 -r2.35
*** complexobject.c	2001年03月11日 08:37:29	2.34
--- complexobject.c	2001年03月18日 08:21:57	2.35
***************
*** 30,34 ****
 static Py_complex c_1 = {1., 0.};
 
! Py_complex c_sum(Py_complex a, Py_complex b)
 {
 	Py_complex r;
--- 30,35 ----
 static Py_complex c_1 = {1., 0.};
 
! Py_complex
! c_sum(Py_complex a, Py_complex b)
 {
 	Py_complex r;
***************
*** 38,42 ****
 }
 
! Py_complex c_diff(Py_complex a, Py_complex b)
 {
 	Py_complex r;
--- 39,44 ----
 }
 
! Py_complex
! c_diff(Py_complex a, Py_complex b)
 {
 	Py_complex r;
***************
*** 46,50 ****
 }
 
! Py_complex c_neg(Py_complex a)
 {
 	Py_complex r;
--- 48,53 ----
 }
 
! Py_complex
! c_neg(Py_complex a)
 {
 	Py_complex r;
***************
*** 54,58 ****
 }
 
! Py_complex c_prod(Py_complex a, Py_complex b)
 {
 	Py_complex r;
--- 57,62 ----
 }
 
! Py_complex
! c_prod(Py_complex a, Py_complex b)
 {
 	Py_complex r;
***************
*** 62,67 ****
 }
 
! Py_complex c_quot(Py_complex a, Py_complex b)
 {
 	Py_complex r;
 	double d = b.real*b.real + b.imag*b.imag;
--- 66,79 ----
 }
 
! Py_complex
! c_quot(Py_complex a, Py_complex b)
 {
+ 	/******************************************************************
+ 	This was the original algorithm. It's grossly prone to spurious
+ 	overflow and underflow errors. It also merrily divides by 0 despite
+ 	checking for that(!). The code still serves a doc purpose here, as
+ 	the algorithm following is a simple by-cases transformation of this
+ 	one:
+ 
 	Py_complex r;
 	double d = b.real*b.real + b.imag*b.imag;
***************
*** 71,77 ****
 	r.imag = (a.imag*b.real - a.real*b.imag)/d;
 	return r;
 }
 
! Py_complex c_pow(Py_complex a, Py_complex b)
 {
 	Py_complex r;
--- 83,125 ----
 	r.imag = (a.imag*b.real - a.real*b.imag)/d;
 	return r;
+ 	******************************************************************/
+ 
+ 	/* This algorithm is better, and is pretty obvious: first divide the
+ 	 * numerators and denominator by whichever of {b.real, b.imag} has
+ 	 * larger magnitude. The earliest reference I found was to CACM
+ 	 * Algorithm 116 (Complex Division, Robert L. Smith, Stanford
+ 	 * University). As usual, though, we're still ignoring all IEEE
+ 	 * endcases.
+ 	 */
+ 	 Py_complex r;	/* the result */
+ 	 const double abs_breal = b.real < 0 ? -b.real : b.real;
+ 	 const double abs_bimag = b.imag < 0 ? -b.imag : b.imag;
+ 
+ 	 if (abs_breal >= abs_bimag) {
+ 		/* divide tops and bottom by b.real */
+ 	 	if (abs_breal == 0.0) {
+ 	 		errno = EDOM;
+ 	 		r.real = r.imag = 0.0;
+ 	 	}
+ 	 	else {
+ 	 		const double ratio = b.imag / b.real;
+ 	 		const double denom = b.real + b.imag * ratio;
+ 	 		r.real = (a.real + a.imag * ratio) / denom;
+ 	 		r.imag = (a.imag - a.real * ratio) / denom;
+ 	 	}
+ 	}
+ 	else {
+ 		/* divide tops and bottom by b.imag */
+ 		const double ratio = b.real / b.imag;
+ 		const double denom = b.real * ratio + b.imag;
+ 		assert(b.imag != 0.0);
+ 		r.real = (a.real * ratio + a.imag) / denom;
+ 		r.imag = (a.imag * ratio - a.real) / denom;
+ 	}
+ 	return r;
 }
 
! Py_complex
! c_pow(Py_complex a, Py_complex b)
 {
 	Py_complex r;
***************
*** 102,106 ****
 }
 
! static Py_complex c_powu(Py_complex x, long n)
 {
 	Py_complex r, p;
--- 150,155 ----
 }
 
! static Py_complex
! c_powu(Py_complex x, long n)
 {
 	Py_complex r, p;
***************
*** 117,121 ****
 }
 
! static Py_complex c_powi(Py_complex x, long n)
 {
 	Py_complex cn;
--- 166,171 ----
 }
 
! static Py_complex
! c_powi(Py_complex x, long n)
 {
 	Py_complex cn;