[pure-lang-svn] SF.net SVN: pure-lang: [346] pure/trunk/lib/math.pure

[<ag...@us...>]

Revision: 346
          http://pure-lang.svn.sourceforge.net/pure-lang/?rev=346&view=rev
Author:   agraef
Date:     2008-07-01 03:43:04 -0700 (Tue, 01 Jul 2008)

Log Message:
-----------
Add complex numbers.

Modified Paths:
--------------
    pure/trunk/lib/math.pure

Modified: pure/trunk/lib/math.pure
===================================================================
--- pure/trunk/lib/math.pure	2008-06-30 20:00:55 UTC (rev 345)
+++ pure/trunk/lib/math.pure	2008-07-01 10:43:04 UTC (rev 346)
@@ -1,5 +1,5 @@
 
-/* Pure basic math routines. */
+/* Pure math routines. Also defines the complex numbers. */
 
 /* Copyright (c) 2008 by Albert Graef <Dr....@t-...>.
 
@@ -93,3 +93,262 @@
 asinh x::int | asinh x::bigint = asinh (double x);
 acosh x::int | acosh x::bigint = acosh (double x);
 atanh x::int | atanh x::bigint = atanh (double x);
+
+/* Complex numbers. We provide both rectangular (x+:y) and polar (r<:a)
+   representations, where (x,y) are the Cartesian coordinates and (r,t) the
+   radius (absolute value) and angle (in radians) of a complex number,
+   respectively. The +: and <: constructors bind weaker than all other
+   arithmetic operators and are non-associative. */
+
+infix 5 +: <: ;
+
+/* Constructor equations to normalize the polar representation r<:t so that r
+   is always nonnegative and t falls in the range -pi<t<=pi. */
+
+r::int<:t	|
+r::bigint<:t	|
+r::double<:t	= -r <: t+pi if r<0;
+
+r<:t::int	|
+r<:t::bigint	|
+r<:t::double	= r <: atan2 (sin t) (cos t) if t<-pi || t>pi;
+		= r <: pi if t==-pi;
+
+/* The imaginary unit. */
+
+def i = 0+:1;
+
+/* The following operations all work with both the rectangular and the polar
+   representation, promoting real inputs to complex where appropriate. When
+   the result of an operation is again a complex number, it generally uses the
+   same representation as the input (except for explicit conversions). */
+
+/* We mostly follow Bronstein/Semendjajew here. Please mail me if you know
+   better formulas for any of these. */
+
+/* Convert any kind of number to a complex value. */
+
+complex z@(x+:y) | complex z@(r<:t) = z;
+complex x::int | complex x::bigint = x+:0;
+complex x::double = x+:0.0;
+
+/* Convert between polar and rectangular representations. */
+
+polar (x+:y)	= sqrt (x*x+y*y) <: atan2 y x;
+rect  (r<:t)	= r*cos t +: r*sin t;
+
+/* For convenience, make these work with real values, too. */
+
+polar x::int | polar x::bigint = x<:0;
+polar x::double = x<:0.0;
+
+rect x::int | rect x::bigint = x+:0;
+rect x::double = x+:0.0;
+
+/* Create complex values on the unit circle. Note: To quickly compute
+   exp (x+:y) in polar form, use exp x*cis y. */
+
+cis t		= rect (1<:t);
+
+/* Modulus (absolute value) and argument (angle). Note that you can also find
+   both of these in one go by converting to polar form. */
+
+abs (x+:y)	= sqrt (x*x+y*y);
+abs (r<:t)	= r;
+
+arg (x+:y)	= atan2 y x;
+arg (r<:t)	= t;
+
+arg x::int	|
+arg x::bigint	|
+arg x::double	= atan2 0 x;
+
+/* Real and imaginary part. */
+
+re (x+:y)	= x;
+re (r<:t)	= r*sin t;
+
+re x::int	|
+re x::bigint	|
+re x::double	= x;
+
+im (x+:y)	= y;
+im (r<:t)	= r*cos t;
+
+im x::int	= 0;
+im x::bigint	= 0L;
+im x::double	= 0.0;
+
+/* Complex conjugate. */
+
+conj (x+:y)	= x +: -y;
+conj (r<:t)	= r <: -t;
+
+conj x::int	|
+conj x::bigint	|
+conj x::double	= x;
+
+/* Complex sqrt. */
+
+sqrt (x+:y)	= sqrt (x*x+y*y) * (cos (t/2) +: sin (t/2));
+sqrt (r<:t)	= sqrt r <: t/2;
+
+// Complex square roots of negative reals.
+sqrt x::int	|
+sqrt x::bigint	|
+sqrt x::double	= 0.0 +: sqrt (-x) if x<0;
+
+/* Complex exponential and logarithms. */
+
+exp (x+:y)	= exp x * (cos y +: sin y);
+exp (r<:t)	= exp (r*cos t) <: r*sin t;
+
+ln z@(x+:y)	= ln (abs z) +: arg z;
+ln (r<:t)	= polar (ln r +: t);
+
+log z@(x+:y)	|
+log z@(r<:t)	= ln z / ln 10;
+
+// Complex logarithms of negative reals.
+ln x::double	= ln (abs x) +: arg x if x<0;
+log x::double	= ln x / ln 10 if x<0;
+
+/* Complex trig functions. */
+
+sin (x+:y)	= sin x*cosh y +: cos x*sinh y;
+cos (x+:y)	= cos x*cosh y +: -sin x*sinh y;
+tan (x+:y)	= (sin (2*x) +: sinh (2*y)) / (cos (2*x)+cosh (2*y));
+
+// These are best computed in rect and then converted back to polar.
+sin z@(r<:t)	= polar $ sin $ rect z;
+cos z@(r<:t)	= polar $ cos $ rect z;
+tan z@(r<:t)	= polar $ tan $ rect z;
+
+// Use complex logarithms for the inverses.
+asin z@(x+:y)	|
+asin z@(r<:t)	= -i*ln (i*z+sqrt (1-z*z));
+acos z@(x+:y)	|
+acos z@(r<:t)	= -i*ln (z+sqrt (z*z-1));
+atan z@(x+:y)	|
+atan z@(r<:t)	= 0.0 +: inf if z==i;
+		= 0.0 +: -inf if z==-i;
+		= -i*0.5*ln ((1+i*z)/(1-i*z));
+
+/* Complex hyperbolic functions. */
+
+sinh (x+:y)	= sinh x*cos y +: cosh x*sin y;
+cosh (x+:y)	= cosh x*cos y +: sinh x*sin y;
+tanh (x+:y)	= (sinh (2*x) +: sin (2*y)) / (cosh (2*x)+cos (2*y));
+
+sinh z@(r<:t)	= polar $ sinh $ rect z;
+cosh z@(r<:t)	= polar $ cosh $ rect z;
+tanh z@(r<:t)	= polar $ tanh $ rect z;
+
+asinh z@(x+:y)	|
+asinh z@(r<:t)	= ln (z+sqrt (z*z+1));
+acosh z@(x+:y)	|
+acosh z@(r<:t)	= ln (z+sqrt (z*z-1));
+atanh z@(x+:y)	|
+atanh z@(r<:t)	= inf +: 0.0 if z==1;
+		= -inf +: 0.0 if z==-1;
+		= ln ((1+z)/(1-z))/2;
+
+// These inverse hyperbolic trigs have complex results for some reals.
+acosh x::double	= acosh (x+:0);
+atanh x::double	= atanh (x+:0);
+
+/* Complex arithmetic. */
+
+-(x+:y)		= -x +: -y;
+-(r<:t)		= r <: t+pi;
+
+(x1+:y1) + (x2+:y2)	= x1+y2 +: y1+y2;
+z1@(r1<:t1)+z2@(r2<:t2)	= polar $ rect z1 + rect z2;
+
+(x1+:y1) - (x2+:y2)	= x1-y2 +: y1-y2;
+z1@(r1<:t1)-z2@(r2<:t2)	= polar $ rect z1 - rect z2;
+
+(x1+:y1) * (x2+:y2)	= x1*x2-y1*y2 +: x1*y2+y1*x2;
+(r1<:t1) * (r2<:t2)	= r1*r2 <: t1+t2;
+
+(x1+:y1) / (x2+:y2)	= (x1*x2+y1*y2 +: y1*x2-x1*y2) / (x2*x2+y2*y2);
+(r1<:t1) / (r2<:t2)	= r1/r2 <: t1-t2;
+
+z1@(x1+:y1)^z2@(x2+:y2)	= exp (ln z1*z2) if z1!=0;
+			= 0.0+:0.0 if z2!=0;
+			= 1.0+:0.0 otherwise;
+z1@(r1<:t1)^z2@(r2<:t2)	= exp (ln z1*z2) if z1!=0;
+			= 0.0<:0.0 if z2!=0;
+			= 1.0<:0.0 otherwise;
+
+// Complex powers of negative reals.
+x1::double^x2::double	= exp (ln x1*x2) if x1<0;
+
+/* Mixed rect/polar and polar/rect forms always return a rect result. */
+
+z1@(x1+:y1)+z2@(r2<:t2)	= z1 + rect z2;
+z1@(r1<:t1)+z2@(x2+:y2)	= rect z1 + z2;
+
+z1@(x1+:y1)-z2@(r2<:t2)	= z1 - rect z2;
+z1@(r1<:t1)-z2@(x2+:y2)	= rect z1 - z2;
+
+z1@(x1+:y1)*z2@(r2<:t2)	= z1 * rect z2;
+z1@(r1<:t1)*z2@(x2+:y2)	= rect z1 * z2;
+
+z1@(x1+:y1)/z2@(r2<:t2)	= z1 / rect z2;
+z1@(r1<:t1)/z2@(x2+:y2)	= rect z1 / z2;
+
+z1@(x1+:y1)^z2@(r2<:t2)	= z1 ^ rect z2;
+z1@(r1<:t1)^z2@(x2+:y2)	= rect z1 ^ z2;
+
+/* Mixed complex/real and real/complex forms yield a rect or polar result,
+   depending on what the complex input was. */
+
+(x1+:y1)+x2	= x1+x2 +: y1;
+x1+(x2+:y2)	= x1+x2 +: y2;
+z1@(r1<:t1)+x2	= rect z1 + x2;
+x1+z2@(r2<:t2)	= x1 + rect z2;
+
+(x1+:y1)-x2	= x1-x2 +: y1;
+x1-(x2+:y2)	= x1-x2 +: -y2;
+z1@(r1<:t1)-x2	= rect z1 - x2;
+x1-z2@(r2<:t2)	= x1 - rect z2;
+
+(x1+:y1)*x2	= x1*x2 +: y1*x2;
+x1*(x2+:y2)	= x1*x2 +: x1*y2;
+(r1<:t1)*x2	= r1*x2 <: t1;
+x1*(r2<:t2)	= x1*r2 <: t2;
+
+(x1+:y1)/x2	= x1/x2 +: y1/x2;
+x1/z2@(x2+:y2)	= (x1*x2 +: -x1*y2) / (x2*x2+y2*y2);
+(r1<:t1)/x2	= r1/x2 <: t1;
+x1/(r2<:t2)	= x1/r2 <: -t2;
+
+z1@(x1+:y1)^x2	= z1 ^ (x2+:0);
+x1^z2@(x2+:y2)	= (x1+:0) ^ z2;
+(r1<:t1)^x2	= r1^x2 <: t1*x2;
+x1^z2@(r2<:t2)	= (x1<:0) ^ z2;
+
+/* Equality. */
+
+(x1+:y1) == (x2+:y2)	= x1==y2 && y1==y2;
+(r1<:t1) == (r2<:t2)	= r1==r2 && t1==t2;
+
+(x1+:y1) != (x2+:y2)	= x1!=y2 || y1!=y2;
+(r1<:t1) != (r2<:t2)	= r1!=r2 || t1!=t2;
+
+z1@(_+:_)==z2@(_<:_)	= z1 == rect z2;
+z1@(_<:_)==z2@(_+:_)	= rect z1 == z2;
+
+z1@(_+:_)!=z2@(_<:_)	= z1 != rect z2;
+z1@(_<:_)!=z2@(_+:_)	= rect z1 != z2;
+
+(x1+:y1)==x2	= x1==x2 && y1==0;
+x1==(x2+:y2)	= x1==x2 && y2==0;
+z1@(r1<:t1)==x2	= z1 == (x2<:0);
+x1==z2@(r2<:t2)	= (x1<:0) == z2;
+
+(x1+:y1)!=x2	= x1!=x2 || y1!=0;
+x1!=(x2+:y2)	= x1!=x2 || y2!=0;
+z1@(r1<:t1)!=x2	= z1 != (x2<:0);
+x1!=z2@(r2<:t2)	= (x1<:0) != z2;


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.