[Caml-list] Re: Complex Arithmetic

[Caml-list] Re: Complex Arithmetic

[wester]

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

Re: [Caml-list] Re: Complex Arithmetic

[David McClain]

> 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)?

Well, well... it appears that DEC Alpha has taken Kahan to heart! I have
never used an Alpha, and perhaps I should. I understood that their
implementation of FP traps is a bit strange, using trap gates interspersed
with the code so you don't actually get the trap at the location of the
error unless you pepper your code to death with trap gates... But I have
heard many great things about the speed of these boxes!

> 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?

These cases again illustrate Kahan's principal about the limitations of the
language of [rectangular] pairs... (-1,-1) when squared exceeds the
principal Riemann sheet and so all attempts at representing complex
arithmetic with rectilinear coordinates of the principal sheet are doomed to
failure. Same for (-1,1).

So in these cases, even NML will fail like all the others. Only polar
representation which can display arbitrary levels of the Riemann surface
could possibly succeed in such cases.

My NML implementation [correct for the principal sheet arithmetic only] is
available to anyone, and I simply (ha!) made sure that all the algebraic
relations between 0+ and 0- were properly observed. For example, IIRC,

    x + 0+ != x  (this fails in the case where x = 0-. It succeeds
everywhere else).
    x + 0- == x

    x * 0- = 0- unless x = 0- in which case 0- * 0- = 0+ (sheesh!!)

(of course I may be mistaken here in my mental recall. I have to reexamine
my code to be certain)

When implementing NML, I actually wrote a simple-minded algebraic theorem
prover in OCaml to assist me in this process. That's another terrific thing
about OCaml - you can perform symbolic calculations very easily, perhaps
even more easily than in Lisp (although I also love Lisp). The whole prover
took less than a hundred lines of OCaml or so... and it ran correctly the
first time out!

- DM


-------------------
To unsubscribe, mail caml-list-request@inria.fr.  Archives: http://caml.inria.fr