
/* complex.q: complex numbers
   $Id: complex.q,v 1.17 2007/09/27 10:09:04 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. */

/* This module supports all the usual basic operations on complex numbers
   (abs, arg, re, im, conj), as well as complex and mixed (complex/real)
   arithmetic and complex sqrt, exp/ln, trig and hyperbolic functions. Please
   note that, as of Q 7.7, Complex is now implemented as an abstract data
   type, very much like Rational (cf. rational.q). The virtual constructor
   (+:) takes two Real values X and Y and returns the complex number X+:Y
   which has X as its real and Y as its imaginary part; this representation is
   also used as Complex's canonical view. An alternative way to write the same
   number using an arithmetic expression is X+i*Y, where the constant i
   defined by this module denotes the imaginary unit, 0+:1. */

include math, stdlib;

public type Complex : Num = virtual (+:) X Y @(+) | private const complex X Y;

@-0x80000000
view (complex X:Real Y:Real)	= '(X +: Y);
@0

/* This alias is provided for Haskell compatibility. */

public virtual (complex::+:) X Y as :+;

/* The imaginary unit: */

public var const i = complex 0 1;

/* Basic operations: */

public arg Z, re Z, im Z, re_im Z, conj Z, cis X, polar R X;

abs (complex X:Real Y:Real)	= sqrt (X*X+Y*Y);

arg (complex X:Real Y:Real)	= atan2 Y X;
arg X:Real			= atan2 0 X;

re (complex X:Real _:Real)	= X;
re X:Real			= X;

im (complex _:Real Y:Real)	= Y;
im _:Real			= 0;

re_im (complex X:Real Y:Real)	= (X,Y);
re_im X:Real			= (X,0);

conj (complex X:Real Y:Real)	= complex X (-Y);
conj X:Real			= X;

cis X:Real			= complex (cos X) (sin X);

polar R:Real X:Real		= R*cis 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 Z:Complex			= sqrt (abs Z) * complex (cos A) (sin A)
				    where A = arg Z/2;
sqrt X:Real			= complex 0.0 (sqrt (-X)) if X<0;

/* Exp/logarithms: */

exp (complex X:Real Y:Real)	= exp X*complex (cos Y) (sin Y);
ln Z:Complex			= complex (ln (abs Z)) (arg Z);
ln X:Real			= complex (ln (abs X)) (arg X) if X<0;

lg Z:Complex			= ln Z / ln 2;
log Z:Complex			= ln Z / ln 10;

/* Trig functions: */

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

asin Z:Complex			= -i*ln (i*Z+sqrt (1-Z*Z));
acos Z:Complex			= -i*ln (Z+sqrt (Z*Z-1));
atan Z:Complex			= complex 0.0 inf if Z=i;
				= complex 0.0 (-inf) if Z=-i;
				= -i*0.5*ln ((1+i*Z)/(1-i*Z)) where Z = 1.0*Z;

/* Hyperbolic functions: */

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

asinh Z:Complex			= ln (Z+sqrt (Z*Z+1));
acosh Z:Complex			= ln (Z+sqrt (Z*Z-1));
atanh Z:Complex			= complex inf 0.0 if Z=1;
				= complex (-inf) 0.0 if Z=-1;
				= ln ((1+Z)/(1-Z))/2 where Z = 1.0*Z;
atanh X:Real			= atanh (complex X 0);

/* Virtual constructor: */

complex X1:Real Y1:Real +: complex X2:Real Y2:Real
				= complex (X1-Y2) (Y1+X2);
X1:Real +: complex X2:Real Y2:Real
				= complex (X1-Y2) X2;
complex X1:Real Y1:Real +: X2:Real
				= complex X1 (Y1+X2);
X:Real +: Y:Real		= complex X Y;

/* Arithmetic: */

-complex X:Real Y:Real		= complex (-X) (-Y);

complex X1:Real Y1:Real+complex X2:Real Y2:Real
				= complex (X1+X2) (Y1+Y2);
X1:Real+complex X2:Real Y2:Real	= complex (X1+X2) Y2;
complex X1:Real Y1:Real+X2:Real	= complex (X1+X2) Y1;

complex X1:Real Y1:Real-complex X2:Real Y2:Real
				= complex (X1-X2) (Y1-Y2);
X1:Real-complex X2:Real Y2:Real	= complex (X1-X2) (-Y2);
complex X1:Real Y1:Real-X2:Real	= complex (X1-X2) Y1;

complex X1:Real Y1:Real*complex X2:Real Y2:Real
				= complex (X1*X2-Y1*Y2) (X1*Y2+Y1*X2);
X1:Real*complex X2:Real Y2:Real	= complex (X1*X2) (X1*Y2);
complex X1:Real Y1:Real*X2:Real	= complex (X1*X2) (Y1*X2);

complex X1:Real Y1:Real/complex X2:Real Y2:Real
				= complex ((X1*X2+Y1*Y2)/(X2*X2+Y2*Y2))
				  ((Y1*X2-X1*Y2)/(X2*X2+Y2*Y2));
X1:Real/Z2:Complex		= complex X1 0/Z2;
Z1:Complex/X2:Real		= Z1/complex X2 0;

/* General exponentiation: */

Z1:Complex^Z2:Complex		= exp (ln Z1*Z2) if Z1<>0;
				= complex 0. 0. if Z2<>0;
X1:Real^Z2:Complex		= exp (ln X1*Z2) if X1<>0;
				= complex 0. 0. if Z2<>0;
Z1:Complex^X2:Real		= exp (ln Z1*X2) if Z1<>0;
				= complex 0. 0. if X2<>0;
X1:Real^X2:Real			= exp (ln X1*X2) if X1<0;

/* Equality: */

(complex X1:Real Y1:Real=complex X2:Real Y2:Real)
				= (X1=X2) and then (Y1=Y2);
(X1:Real=complex X2:Real Y2:Real)
				= (X1=X2) and then (Y2=0);
(complex X1:Real Y1:Real=X2:Real)
				= (X1=X2) and then (Y1=0);

complex X1:Real Y1:Real<>complex X2:Real Y2:Real
				= (X1<>X2) or else (Y1<>Y2);
X1:Real<>complex X2:Real Y2:Real
				= (X1<>X2) or else (Y2<>0);
complex X1:Real Y1:Real<>X2:Real
				= (X1<>X2) or else (Y1<>0);
