/* complex.q: complex numbers
   $Id: complex.q,v 1.2 2004/01/28 16:54:43 agraef Exp $ */

/* This file is part of the Q programming system.

   The Q programming system is free software; you can redistribute it and/or
   modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2, or (at your option)
   any later version.

   The Q programming system is distributed in the hope that it will be
   useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */

/* Complex numbers are represented as pairs of Num's (i.e., "reals"). All the
   usual basic operations (abs, arg, re, im) are supported, as well as complex
   arithmetic and complex sqrt, exp/ln, trig and hyperbolic functions. The
   conj function computes the complex conjugate. */

include stdlib, math;

/* Basic operations: */

public iscomplex Z, arg Z, re Z, im Z, conj Z;

iscomplex (_:Num,_:Num) = true;
iscomplex _:Num		= true;
iscomplex _		= false otherwise;

abs (X:Num,Y:Num)	= sqrt (X*X+Y*Y);

arg (X:Num,Y:Num)	= atan2 Y X;
arg X:Num		= atan2 0 X;

re (X:Num,_:Num)	= X;
re X:Num		= X;

im (_:Num,Y:Num)	= Y;
im _:Int		= 0;
im _:Float		= 0.;

conj (X:Num,Y:Num)	= (X,-Y);
conj X:Num		= X;

/* The following is a more or less literal transcription from Bronstein/
   Semendjajew. There might be more efficient or precise ways to do these
   things, so if you have any suggestions, please let me know. */

/* Complex sqrt: */

sqrt (X:Num,Y:Num)	= sqrt (abs (X,Y)) * (cos (arg (X,Y)/2),
				sin (arg (X,Y)/2));
sqrt X:Num		= (0.,sqrt (-X)) if X<0;

/* Exp/logarithms: */

exp (X:Num,Y:Num)	= exp X*(cos Y,sin Y);
ln (X:Num,Y:Num)	= (ln (abs (X,Y)),arg (X,Y))
				if (X<>0) or else (Y<>0);
ln X:Num		= (ln (abs X),arg X)
				if X<0;

lg (X:Num,Y:Num)	= ln (X,Y) / ln 2;
log (X:Num,Y:Num)	= ln (X,Y) / ln 10;

/* Trig functions: */

sin (X:Num,Y:Num)	= (sin X*cosh Y,cos X*sinh Y);
cos (X:Num,Y:Num)	= (cos X*cosh Y,-sin X*sinh Y);
tan (X:Num,Y:Num)	= (sin (2*X)/(cos (2*X)+cosh (2*Y)),
				sinh (2*Y)/(cos (2*X)+cosh (2*Y)));

asin (X:Num,Y:Num)	= (0,-1)*ln ((-Y,X)+sqrt (1-(X,Y)*(X,Y)));
acos (X:Num,Y:Num)	= (0,-1)*ln ((X,Y)+sqrt ((X,Y)*(X,Y)-1));
atan (X:Num,Y:Num)	= (0,-1)*ln ((1-Y,X)/(1+Y,-X))/2
				if (X<>0) or else (1+Y<>0);

/* Hyperbolic functions: */

sinh (X:Num,Y:Num)	= (sinh X*cos Y,cosh X*sin Y);
cosh (X:Num,Y:Num)	= (cosh X*cos Y,sinh X*sin Y);
tanh (X:Num,Y:Num)	= (sinh (2*X)/(cosh (2*X)+cos (2*Y)),
				sin (2*Y)/(cosh (2*X)+cos (2*Y)));

asinh (X:Num,Y:Num)	= ln ((X,Y)+sqrt ((X,Y)*(X,Y)+1));
acosh (X:Num,Y:Num)	= ln ((X,Y)+sqrt ((X,Y)*(X,Y)-1));
atanh (X:Num,Y:Num)	= ln ((1+X,Y)/(1-X,-Y))/2
				if (Y<>0) or else (1-X<>0);

/* Arithmetic: */

-(X:Num,Y:Num)		= (-X,-Y);

(X1:Num,Y1:Num)+(X2:Num,Y2:Num)
			= (X1+X2,Y1+Y2);
X1:Num+(X2:Num,Y2:Num)	= (X1+X2,Y2);
(X1:Num,Y1:Num)+X2:Num
			= (X1+X2,Y1);

(X1:Num,Y1:Num)-(X2:Num,Y2:Num)
			= (X1-X2,Y1-Y2);
X1:Num-(X2:Num,Y2:Num)	= (X1-X2,-Y2);
(X1:Num,Y1:Num)-X2:Num	= (X1-X2,Y1);

(X1:Num,Y1:Num)*(X2:Num,Y2:Num)
			= (X1*X2-Y1*Y2,X1*Y2+Y1*X2);
X1:Num*(X2:Num,Y2:Num)	= (X1*X2,X1*Y2);
(X1:Num,Y1:Num)*X2:Num	= (X1*X2,Y1*X2);

(X1:Num,Y1:Num)/(X2:Num,Y2:Num)
			= ((X1*X2+Y1*Y2)/(X2*X2+Y2*Y2),
				(Y1*X2-X1*Y2)/(X2*X2+Y2*Y2))
				if (X2<>0) or else (Y2<>0);
X1:Num/(X2:Num,Y2:Num)	= ((X1*X2)/(X2*X2+Y2*Y2),(-X1*Y2)/(X2*X2+Y2*Y2))
				if (X2<>0) or else (Y2<>0);
(X1:Num,Y1:Num)/X2:Num	= (X1/X2,Y1/X2)
				if X2<>0;

/* General exponentiation: */

(X1:Num,Y1:Num)^(X2:Num,Y2:Num)
			= exp (ln (X1,Y1)*(X2,Y2))
				if (X1<>0) or else (Y1<>0);
			= (0.,0.)
				if (X2<>0) or else (Y2<>0);
X1:Num^(X2:Num,Y2:Num)	= exp (ln X1*(X2,Y2))
				if X1<>0;
			= (0.,0.)
				if (X2<>0) or else (Y2<>0);
(X1:Num,Y1:Num)^X2:Num	= exp (ln (X1,Y1)*X2)
				if (X1<>0) or else (Y1<>0);
			= (0.,0.)
				if X2<>0;
X1:Num^X2:Num		= exp (ln X1*X2)
				if X1<0;

/* Equality: */

/* Note: since stdlib.q provides tuple equality/inequality, we only have to
   define comparisons between complex and real numbers here. */

(X1:Num=(X2:Num,Y2:Num))
			= (X1=X2) and then (Y2=0);
((X1:Num,Y1:Num)=X2:Num)
			= (X1=X2) and then (Y1=0);

X1:Num<>(X2:Num,Y2:Num)	= (X1<>X2) or else (Y2<>0);
(X1:Num,Y1:Num)<>X2:Num	= (X1<>X2) or else (Y1<>0);