[Caml-list] Re: Complex Arithmetic

[wester@ilt.fhg.de]

Hi David,

as far as I understood Kahan says that the complex arithmetic in 
FORTRAN and many C++ libraries is wrong. I tested:

sqrt( (0.0 , -1.0)**2 )  vs. sqrt( (0.0 , 1.0)**2 ) 

on a DEC Alpha with FORTRAN and C++ (source and output below). 
They both give the correct result. A naive approach abviously gives
the same result (0.0, 1.0) for both expressions. Can you tell us how you 
implemented the correct behaviour in OCaml (or Nml)?

BTW what about sqrt( (-1.0 , -1.0)**2 ) and sqrt( (-1.0 , 1.0)**2 ), which
result in (1.0, 1.0) and (1.0, -1.0) on the DEC Alpha. What should the
results be in this case?

Rolf

// c++ on DEC Alpha
#include <complex>
#include <iostream>

main()
{
	complex<double> a,b,c,d,e,f;
	a = complex<double>(0.0,-1.0);
	b = complex<double>(0.0, 1.0);
	c = a*a;
	d = b*b;
	e = sqrt(c);
	f = sqrt(d);
	cout << a << " , " << b << endl; 
	cout << c << " , " << d << endl; 
	cout << e << " , " << f << endl; 
}
Output:
(0,-1) , (0,1)
(-1,-0) , (-1,0)
(6.12323e-17,-1) , (6.12323e-17,1)


//f77 on DEC Alpha
	complex*16 a,b,c,d,e,f
	a = (0.0,-1.0)
	b = (0.0, 1.0)
	c = a**2
	d = b**2
	e = sqrt(c)
	f = sqrt(d)
	write(6,*) a,b
	write(6,*) c,d
	write(6,*) e,f
	stop
	end
Output:	
(0.000000000000000E+000,-1.00000000000000)
(0.000000000000000E+000,1.00000000000000)
(-1.00000000000000,0.000000000000000E+000)
(-1.00000000000000,0.000000000000000E+000)
(0.000000000000000E+000,-1.00000000000000)
(0.000000000000000E+000,1.00000000000000)
-------------------------------------
Rolf Wester
wester@ilt.fhg.de
-------------------
To unsubscribe, mail caml-list-request@inria.fr.  Archives: http://caml.inria.fr