home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / sa104os2.zip / SATHR104.ZIP / SATHER / LIBRARY / CPX.SA < prev    next >
Text File  |  1995-02-05  |  8KB  |  239 lines
  1. -- Copyright (C) International Computer Science Institute, 1994.  COPYRIGHT  --
  2. -- NOTICE: This code is provided "AS IS" WITHOUT ANY WARRANTY and is subject --
  3. -- to the terms of the SATHER LIBRARY GENERAL PUBLIC LICENSE contained in    --
  4. -- the file "Doc/License" of the Sather distribution.  The license is also   --
  5. -- available from ICSI, 1947 Center St., Suite 600, Berkeley CA 94704, USA.  --
  6. --------> Please email comments to "sather-bugs@icsi.berkeley.edu". <----------
  7.  
  8. -- cpx.sa: Complex numbers.
  9. -------------------------------------------------------------------
  10. value class CPX{T} is
  11.    -- Complex numbers.
  12.    -- 
  13.    -- Some of the algorithms are taken from:
  14.    -- Press, Flannery, Teukolsky, and Vettering, "Numerical Recipes in
  15.    -- C", second edition, Cambridge University Press, 1993.
  16.    -- 
  17.    -- Some of the choices of branch cut were chosen to be consistent
  18.    -- with:
  19.    -- Guy L. Steele, "Common Lisp, The Language", second edition,
  20.    -- 1990, Digital Press.
  21.  
  22.    attr re,im:T;        -- Real and imaginary parts.
  23.  
  24.    str:STR is
  25.       -- A string representation of self of the form "1.02+3.23i".
  26.       buf:FSTR;
  27.       if im>=#T(0.0) then 
  28.      buf:=buf + re.str + ("+") + im.str + "i";
  29.       else buf:=buf + re.str + ("-") + (-im).str + "i" end;
  30.       return buf.str end;
  31.    
  32.    
  33.    create(x:FLT):SAME is return #(#T(x),#T(0.0)) end;
  34.    
  35.    create(re,im:T):SAME is
  36.       -- A complex number with real part `re' and imaginary part `im'. 
  37.       r:SAME; return r.re(re).im(im) end;
  38.  
  39.    plus(c:SAME):SAME is
  40.       -- The sum of self and `c'.
  41.       return #(re+c.re,im+c.im) end;
  42.    
  43.    minus(c:SAME):SAME is
  44.       -- The difference of self and `c'.
  45.       return #(re-c.re,im-c.im) end;      
  46.  
  47.    negate:SAME is
  48.       -- The additive inverse of self.
  49.       return #(-re,-im) end;            
  50.  
  51.    conjugate:SAME is
  52.       -- The complex conjugate of self.
  53.       return #(re,-im) end;                  
  54.  
  55.    times(c:SAME):SAME is
  56.       -- The product of self and `c'.
  57.       return #(re*c.re-im*c.im, re*c.im+im*c.re) end;
  58.  
  59.    div(c:SAME):SAME is
  60.       -- The ratio of self and `c'.
  61.       -- From Numerical Recipes in C, 2nd ed. p. 949.
  62.       if c.re.abs>=c.im.abs then
  63.      r::=c.im/c.re; den::=c.re+r*c.im;
  64.      return #((re+r*im)/den,(im-r*re)/den)
  65. --    else r::=c.re/c.im; den::=c.im+r*c.re;                                    -- NLP
  66.       end; r::=c.re/c.im; den::=c.im+r*c.re;                                    -- NLP
  67. --       return #((re*r+im)/den,(im*r-re)/den) end end;                         -- NLP
  68.          return #((re*r+im)/den,(im*r-re)/den); end;                            -- NLP
  69.       
  70.    reciprocal:SAME is
  71.       -- The multiplicative inverse of self.
  72.       if re.abs>=im.abs then
  73.      r::=im/re; den::=re+r*im;
  74.      return #(#T(1.0)/den,-r/den)
  75. --    else r::=re/im; den::=im+r*re;                                            -- NLP
  76.       end; r::=re/im; den::=im+r*re;                                            -- NLP
  77. --       return #(r/den,#T(-1.0)/den) end end;                                  -- NLP
  78.          return #(r/den,#T(-1.0)/den); end;                                     -- NLP
  79.  
  80.    exp:SAME is
  81.       -- The complex exponential `e^self'.
  82.       rp::=re.exp; return #(rp*im.cos,rp*im.sin) end;       
  83.  
  84.    times(f:T):SAME is
  85.       -- Self times the floating point c.
  86.       return #(re*f,im*f) end;
  87.  
  88.    div(f:T):SAME is
  89.       -- Self div the floating point f.
  90.       return #(re/f,im/f) end;
  91.  
  92.    sqrt:SAME is
  93.       -- The square root of self.
  94.       -- From Numerical Recipes in C, 2nd ed. p. 949.
  95.       -- Steele, p. 302 chooses the branch cut by `e^((log z)/2)'
  96.       if re=#T(0.0) and im=#T(0.0) then return #(#T(0.0),#T(0.0)) end;
  97.       x::=re.abs; y::=im.abs;
  98.       w:T;
  99.       if x>=y then r::=y/x; w:=x.sqrt*(#T(0.5)*(#T(1.0)+(#T(1.0)+r*r).sqrt)).sqrt;
  100.       else r::=x/y; w:=y.sqrt*(#T(0.5)*(r+(#T(1.0)+r*r).sqrt)).sqrt end;
  101.       if re>=#T(0.0) then return #(w,im/(#T(2.0)*w));
  102.       elsif im>=#T(0.0) then return #(im/(#T(2.0)*w),w)
  103. --    else return #(-im/(#T(2.0)*w),-w) end end;                                -- NLP
  104.       end; return #(-im/(#T(2.0)*w),-w); end;                                   -- NLP
  105.      
  106.    abs:T is
  107.       -- The absolute magnitude of self.
  108.       -- From Numerical Recipes in C, 2nd ed. p. 949.
  109.       x::=re.abs; y::=im.abs; temp:T;
  110.       if x=#T(0.0) then return y
  111.       elsif y=#T(0.0) then return x
  112.       elsif x>y then temp:=y/x; return x*(#T(1.0)+temp*temp).sqrt;
  113. --    else temp:=x/y; return y*(#T(1.0)+temp*temp).sqrt end end;                -- NLP
  114.       end; temp:=x/y; return y*(#T(1.0)+temp*temp).sqrt; end;                   -- NLP
  115.  
  116.    abs_squared:T is
  117.       -- The square of the absolute magnitude of self.
  118.       return re*re+im*im end;
  119.    
  120.    magnitude:T is
  121.       -- The absolute magnitude of self.
  122.       return abs end;
  123.    
  124.    magnitude_squared:T is
  125.       -- The square of the absolute magnitude of self.
  126.       return abs_squared end;
  127.    
  128.    phase:T is
  129.       -- The angle part of the polar represenation of self.
  130.       -- `-pi < res <= pi'. Also get "-pi" from a negative real
  131.       -- part and a "-0.0" imaginary part.
  132.       -- They say 0+0i should be +0, 0-0i should be -0, -0+0i
  133.       -- should be +pi, and -0-0i should be -pi.
  134.       -- Same convention as Steele, p. 303.
  135.       return im.atan2(re) end;
  136.  
  137.    log:SAME is
  138.       -- The complex logarithm. The chosen branch is
  139.       -- `log |z| + i phase(z)'.
  140.       -- Same convention as Steele, p. 302.
  141.       return #((re*re+im*im).sqrt.log,phase) end;
  142.  
  143.    sign:SAME is
  144.       -- If not zero, a number with same sign as self but unit
  145.       -- magnitude. If it is, then returns self.
  146.       -- Steele, p. 304
  147.       --***
  148.       raise "CPX::sign:SAME not defined."
  149.    end;
  150.  
  151.    create_from_polar(mag,phase:T):SAME is
  152.       -- A complex number with magnitude `mag' and phase `phase'.
  153.       r:SAME; return r.re(mag * phase.cos).im(mag * phase.sin) end;
  154.  
  155.    cis(f:T):SAME is
  156.       -- Ignores self, e^i*f=cos f + i sin f . 
  157.       -- Steele p. 304.
  158.       r:SAME; return r.re(f.cos).im(f.sin) end;
  159.  
  160.    asin:SAME is
  161.       -- -i log(iz + sqrt(1-z^2))
  162.       -- Steele p. 305.
  163.       --***
  164.       raise "CPX::asin:SAME missing."             
  165.    end;
  166.  
  167.    acos:SAME is
  168.       -- -i log(z + sqrt(1-z^2))
  169.       -- better to use: (pi/2)-asin(z)
  170.       -- Steele p. 305.
  171.       --***
  172.       raise "CPX::acos:SAME missing."                   
  173.    end;
  174.  
  175.    atan:SAME is
  176.       -- (log(1+i*y)-log(1-i*y))/(2*i)
  177.       -- Steele p. 307.
  178.       --***
  179.       raise "CPX::atan:SAME missing."                         
  180.    end;
  181.  
  182.    sinh:SAME is
  183.       -- (e^z-e^(-z))/2
  184.       -- Steele p. 308
  185.       --***
  186.       raise "CPX::sinh:SAME missing."                               
  187.    end;
  188.  
  189.    cosh:SAME is
  190.       -- (e^z+e^(-z))/2
  191.       -- Steele p. 308
  192.       --***
  193.       raise "CPX::cosh:SAME missing."                               
  194.    end;      
  195.  
  196.    tanh:SAME is
  197.       -- (e^z-e^(-z))/(e^z+e^(-z))
  198.       -- Steele p. 308
  199.       --***
  200.       raise "CPX::tanh:SAME missing."     
  201.    end;
  202.  
  203.    asinh:SAME is
  204.       -- log(z+sqrt(1+z^2))
  205.       -- Steele p. 308
  206.       --***
  207.       raise "CPX::asinh:SAME missing."           
  208.    end;
  209.  
  210.    acosh:SAME is
  211.       -- log(z+(z+1)sqrt((z-1)/(z+1)))
  212.       -- Steele p. 308
  213.       --***
  214.       raise "CPX::acosh:SAME missing."                 
  215.    end;
  216.  
  217.    atanh:SAME is
  218.       -- log((1+z)sqrt(1/(1-z^2)))
  219.       -- Steele p. 308
  220.       --***  
  221.       raise "CPX::atanh:SAME missing."                       
  222.    end;      
  223.    
  224. end; -- class CPX{T}
  225.  
  226. -------------------------------------------------------------------
  227.    
  228. value class CPX is
  229.    include CPX{FLT}
  230. end; -- class CPX
  231.  
  232. -------------------------------------------------------------------
  233.  
  234. value class CPXD is
  235.    include CPX{FLTD}
  236. end; -- class CPXD
  237.  
  238. -------------------------------------------------------------------
  239.