cpx.sa
Generated by gen_html_sa_files from ICSI. Contact gomes@icsi.berkeley.edu for details
---------------------------> Sather 1.1 source file <--------------------------
-- Copyright (C) International Computer Science Institute, 1994. COPYRIGHT --
-- NOTICE: This code is provided "AS IS" WITHOUT ANY WARRANTY and is subject --
-- to the terms of the SATHER LIBRARY GENERAL PUBLIC LICENSE contained in --
-- the file "Doc/License" of the Sather distribution. The license is also --
-- available from ICSI, 1947 Center St., Suite 600, Berkeley CA 94704, USA. --
--------> Please email comments to sather-bugs@icsi.berkeley.edu. <----------
-- cpx.sa: Complex numbers.
-- Notes about deleted code:
-- The following require conversion from FLT to T and FLTD to T.
-- don't feel like putting those in $NUMBER{T}
-- create(c:CPX):SAME is r:SAME; return r.re(#T(c.re)).im(#T(c.im)) end;
-- create(c:CPXD):SAME is r:SAME; return r.re(c.re).im(#T(c.im)) end;
-- The following two functions do not make sense in this class, though
-- they could be defined in the concrete instantiations.
-- If we had them here, we'd also have to have create(c:CPX{FLT}) which
-- seems to defeat the whole purpose of instantiating the various versions
-- of cpx
-- cpx:CPX is return #CPX(self) end;
-- cpxd:CPXD is return #CPXD(self) end;
immutable class CPX{T < $REAL_NUMBER{T}} < $CPX_NUMBER{T,CPX{T}}
immutable class CPX{T < $REAL_NUMBER{T}} < $CPX_NUMBER{T,CPX{T}} is
-- Complex numbers.
--
-- Some of the algorithms are taken from:
-- Press, Flannery, Teukolsky, and Vettering, "Numerical Recipes in
-- C", second edition, Cambridge University Press, 1993.
--
-- Some of the choices of branch cut were chosen to be consistent
-- with:
-- Guy L. Steele, "Common Lisp, The Language", second edition,
-- 1990, Digital Press.
include COMPARABLE;
attr re,im:T; -- Real and imaginary parts.
create_real(x:T):SAME is return #(x,#T(0.0)) end;
create(re,im:T):SAME is
-- A complex number with real part `re' and imaginary part `im'.
r:SAME; return r.re(re).im(im)
end;
create(i: INT): SAME is
return create_real(T::create(i));
end;
create(f: FLT): SAME is
return create_real(T::create(f));
end;
create(f: FLTD): SAME is
return create_real(T::create(f));
end;
zero: SAME is return #SAME(T::zero,T::zero); end;
-- A zero valued complex number
one: SAME is return #SAME(T::one,T::zero); end;
-- A unit value (real part=1, imaginary part=0)
maxval: SAME is return #SAME(T::maxval,T::maxval); end;
-- Maximum values of real and imaginary part
is_eq(c: SAME): BOOL is return re=c.re and im=c.im end;
-- Return true if the two numbers are equal
is_nil:BOOL is return re.is_nil or im.is_nil end;
-- Return true if self is nil
nil:SAME is return #(re.nil,im.nil); end;
-- Return a nil valued complex number with real and imaginary
-- parts both nil
str:STR is
-- A string representation of self of the form "1.02+3.23i".
buf:FSTR;
if im>=#T(0.0) then
buf:=buf + re.str + ("+") + im.str + "i";
else buf:=buf + re.str + ("-") + (-im).str + "i" end;
return buf.str end;
plus(c:SAME):SAME is
-- The sum of self and `c'.
return #(re+c.re,im+c.im) end;
minus(c:SAME):SAME is
-- The difference of self and `c'.
return #(re-c.re,im-c.im) end;
negate:SAME is
-- The additive inverse of self.
return #(-re,-im) end;
conjugate:SAME is
-- The complex conjugate of self.
return #(re,-im) end;
times(c:SAME):SAME is
-- The product of self and `c'.
return #(re*c.re-im*c.im, re*c.im+im*c.re)
end;
div(c:SAME):SAME is
-- The ratio of self and `c'.
-- From Numerical Recipes in C, 2nd ed. p. 949.
if c.re.abs >=c.im.abs then
r::=c.im/c.re; den::=c.re+r*c.im;
return #((re+r*im)/den,(im-r*re)/den)
else r::=c.re/c.im; den::=c.im+r*c.re;
return #((re*r+im)/den,(im*r-re)/den)
end
end;
reciprocal:SAME is
-- The multiplicative inverse of self.
if re.abs>=im.abs then
r::=im/re; den::=re+r*im;
return #(#T(1.0)/den,-r/den)
else r::=re/im; den::=im+r*re;
return #(r/den,#T(-1.0)/den)
end
end;
exp:SAME is
-- The complex exponential `e^self'.
rp::=re.exp; return #(rp*im.cos,rp*im.sin)
end;
times(f:T):SAME is
-- Self times the floating point c.
return #(re*f,im*f)
end;
div(f:T):SAME is
-- Self div the floating point f.
return #(re/f,im/f)
end;
pow(c:SAME):SAME is
-- self^c = exp(c*log(self))
return (self.log*c).exp
end;
sqrt:SAME is
-- The square root of self.
-- From Numerical Recipes in C, 2nd ed. p. 949.
-- Steele, p. 302 chooses the branch cut by `e^((log z)/2)'
if re=#T(0.0) and im=#T(0.0) then return #(#T(0.0),#T(0.0)) end;
x::=re.abs; y::=im.abs;
w:T;
if x>=y then r::=y/x; w:=x.sqrt*(#T(0.5)*(#T(1.0)+(#T(1.0)+r*r).sqrt)).sqrt;
else r::=x/y; w:=y.sqrt*(#T(0.5)*(r+(#T(1.0)+r*r).sqrt)).sqrt end;
if re>=#T(0.0) then return #(w,im/(#T(2.0)*w));
elsif im>=#T(0.0) then return #(im/(#T(2.0)*w),w)
else return #(-im/(#T(2.0)*w),-w) end
end;
cube_root:SAME is
-- The cube root of self.
-- preliminary, but working.
return self.pow(create_real(T::one/#T(3.0)))
end;
square:SAME is
-- Self squared
return self*self
end;
cube:SAME is
-- Self cubed
return self*self*self
end;
abs: SAME is
-- For conformance with $NFE
return create_real(absolute);
end;
absolute:T is
-- The absolute magnitude of self.
-- Cannot overload abs
-- From Numerical Recipes in C, 2nd ed. p. 949.
x::=re.abs; y::=im.abs; temp:T;
if x=#T(0.0) then return y
elsif y=#T(0.0) then return x
elsif x>y then temp:=y/x; return x*(#T(1.0)+temp*temp).sqrt;
else temp:=x/y; return y*(#T(1.0)+temp*temp).sqrt end
end;
abs_squared:T is
-- The square of the absolute magnitude of self.
return re*re+im*im
end;
magnitude:T is
-- The absolute magnitude of self.
return absolute
end;
magnitude_squared:T is
-- The square of the absolute magnitude of self.
return abs_squared
end;
phase:T is
-- The angle part of the polar represenation of self.
-- `-pi < res <= pi'. Also get "-pi" from a negative real
-- part and a "-0.0" imaginary part.
-- They say 0+0i should be +0, 0-0i should be -0, -0+0i
-- should be +pi, and -0-0i should be -pi.
-- Same convention as Steele, p. 303.
return im.atan2(re)
end;
log:SAME is
-- The complex logarithm. The chosen branch is
-- `log |z| + i phase(z)'.
-- Same convention as Steele, p. 302.
return #((re*re+im*im).sqrt.log,phase)
end;
sign:SAME is
-- If not zero, a number with same phase as self but unit
-- magnitude. If it is, then returns self.
-- Steele, p. 304
--***
raise "CPX::sign:SAME not defined."
end;
is_within(r:T,c:SAME):BOOL is
-- self is within a circle around c with radius r
return (self-c).abs_squared <= r
end;
create_from_polar(mag,phase:T):SAME is
-- A complex number with magnitude `mag' and phase `phase'.
r:SAME; return r.re(mag * phase.cos).im(mag * phase.sin)
end;
cis(f:T):SAME is
-- Ignores self, e^i*f=cos f + i sin f .
-- Steele p. 304.
r:SAME; return r.re(f.cos).im(f.sin)
end;
-- preliminary, but working
sin:SAME is return (self.exp-(-self).exp)/#T(2.0) end;
cos:SAME is return (self.exp+(-self).exp)/#T(2.0) end;
tan:SAME is return self.sin/self.cos end;
asin:SAME is
-- -i log(iz + sqrt(1-z^2))
-- Steele p. 305.
--***
raise "CPX::asin:SAME missing."
end;
acos:SAME is
-- -i log(z + sqrt(1-z^2))
-- better to use: (pi/2)-asin(z)
-- Steele p. 305.
--***
raise "CPX::acos:SAME missing."
end;
atan:SAME is
-- (log(1+i*y)-log(1-i*y))/(2*i)
-- Steele p. 307.
--***
raise "CPX::atan:SAME missing."
end;
sinh:SAME is
-- (e^z-e^(-z))/2
-- Steele p. 308
--***
raise "CPX::sinh:SAME missing."
end;
cosh:SAME is
-- (e^z+e^(-z))/2
-- Steele p. 308
--***
raise "CPX::cosh:SAME missing."
end;
tanh:SAME is
-- (e^z-e^(-z))/(e^z+e^(-z))
-- Steele p. 308
--***
raise "CPX::tanh:SAME missing."
end;
asinh:SAME is
-- log(z+sqrt(1+z^2))
-- Steele p. 308
--***
raise "CPX::asinh:SAME missing."
end;
acosh:SAME is
-- log(z+(z+1)sqrt((z-1)/(z+1)))
-- Steele p. 308
--***
raise "CPX::acosh:SAME missing."
end;
atanh:SAME is
-- log((1+z)sqrt(1/(1-z^2)))
-- Steele p. 308
--***
raise "CPX::atanh:SAME missing."
end;
end; -- class CPX{T}
immutable class CPX < $CPX_NUMBER{FLT,CPX},$FMT
immutable class CPX < $CPX_NUMBER{FLT,CPX},$FMT is
include CPX{FLT};
fmt( fmt: STR ): STR
is
comma: INT;
fmt1: STR;
comma := fmt.search(';');
if comma >= 0 then
fmt1 := fmt.head(comma);
fmt := fmt.substring(comma+1)
else
fmt1 := fmt
end;
if ~ fmt1.is_prefix("polar") then
return re.fmt(fmt1) + '+' + im.fmt(fmt) + 'i'
else
fmt1 := fmt1.substring(5);
return absolute.fmt(fmt1) + "*e^" + phase.fmt(fmt) + 'i'
end;
-- return BASE_FORMAT::fmt_cpx( self, f )
end;
end; -- class CPX