------------------------->  GNU Sather - sourcefile  <-------------------------
-- Copyright (C) 1994 by International Computer Science Institute            --
-- This file is part of the GNU Sather library. It is free software; you may --
-- redistribute  and/or modify it under the terms of the GNU Library General --
-- Public  License (LGPL)  as published  by the  Free  Software  Foundation; --
-- either version 2 of the license, or (at your option) any later version.   --
-- This  library  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 Doc/LGPL for more details.       --
-- The license text is also available from:  Free Software Foundation, Inc., --
-- 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA                     --
-------------->  Please email comments to <bug-sather@gnu.org>  <--------------

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

-------------------------------------------------------------------
