Path: utzoo!mnetor!uunet!husc6!cmcl2!brl-adm!adm!4526P%NAVPGS.BITNET@CUNYVM.CUNY.EDU
From: 4526P%NAVPGS.BITNET@CUNYVM.CUNY.EDU (LT Scott A. Norton, USN)
Newsgroups: comp.lang.pascal
Subject: Complex arithmetic
Message-ID: <11219@brl-adm.ARPA>
Date: 12 Jan 88 02:42:44 GMT
Sender: news@brl-adm.ARPA
Lines: 205

As an electrical engineer, I do a lot of complex math.  My current
complex arithmetic package for Pascal is a bit cumbersome, and I wonder
if anyone can see any way to clean it up.

My first gripe is that Pascal doesn't provide, or let me define, infix
arithmetic operators for complex numbers.  Fortran has complex math
built in, and Ada lets you define "+", for example, as an operator
for user-defined types.  I can't see any way around this, so my
complex math is based on functions, with the result of a LISP-like
next of parenthesis.  add( times( a,b ), times( c,d ) );

The other restriction is the fact that a function can not return a
record.  So I have to return pointers to records.  The records
themselves I have to allocate with new() inside the function.  Then,
to prevent the growth of garbage, I needed some way of disposing of these
allocated records.  My solution was that every complex math function
disposes of its arguments.  Finally, I needed conversion routines
to transform complex variables into allocated records, and an assignment
procedure to store results in variables.

Is there any cleaner way to implement complex numbers in Pascal?  My
only other idea was to mark() the heap, but that is not standard Pascal.

The declarations for my complex numbers follows.  If you can see a way
to improve them, please let me know.  Conversely, if you want too use
them, feel free.  Since they were developed by a federal employee( me ),
they are public domain.
-------------------------------------------------------------------
const pi = 3.14159265358979;
type
     complex = record
        re, im : real;
        end; { record }
     cxptr = @complex;

{ These operations return pointers to complex numbers.  To prevent
garbage accumulation, they dispose of their arguments, which are
also @complex.  }

function v( a : complex ) : cxptr;  { 'value' function.  Converts  }
                                    { complex to cxptr.            }
     var result : cxptr;
     begin
        new(result);
        result@ := a;
        v := result;
        end; { v }

function cmplx(a,b : real ) : cxptr;  { a + b*j }
     var result : cxptr;
     begin
        new(result);
        result@.re := a;  result@.im := b;
        cmplx := result;
        end; { cmplx }

procedure assign( var a : complex; b : cxptr ); { Assign a := b@ }
     begin
        a := b@;
        dispose( b );
        end; { assign }

function plus (a,b : cxptr) : cxptr;
 external; function plus;
     var result : cxptr;
     begin
        new(result);
        result@ .re := a@ .re + b@ .re;
        result@ .im := a@ .im + b@ .im;
        plus := result;
        dispose (a); dispose (b);
     end; { plus }

function minus(a,b : cxptr) : cxptr;
 external; function minus;
     var result : cxptr;
     begin
        new(result);
        result@ .re := a@ .re - b@ .re;
        result@ .im := a@ .im - b@ .im;
        minus := result;
        dispose (a); dispose (b);
     end; { minus }

function times (a,b : cxptr) : cxptr;
 external; function times;
     var result : cxptr;
     begin
        new(result);
        result@ .re := a@ .re * b@ .re - a@ .im * b@ .im;
        result@ .im := a@ .im * b@ .re + a@ .re * b@ .im;
        times := result;
        dispose (a); dispose (b);
     end; { times }

function magsq (a : cxptr) : real;
 external; function magsq;
         { Calculates the magnitude squared of the argument a }
     begin
        magsq := sqr(a@.re) + sqr(a@.im);
     end; { magsq }

function divide (a,b : cxptr) : cxptr;
 external; function divide;
     var result : cxptr;
         m : real;
     begin
        m := 1/ magsq( b );
        new(result);
        result@ .re := ( a@.re * b@.re + a@.im * b@.im ) * m;
        result@ .im := ( a@.im * b@.re - a@.re * b@.im ) * m;
        divide := result;
        dispose (a); dispose (b);
     end; { divide }

function scale (a : real; b : cxptr) : cxptr;
 external; function scale;
     var result : cxptr;
     begin
        new(result);
        result@ .re := a * b@ .re;
        result@ .im := a * b@ .im;
        scale := result;
        dispose (b);
     end; { scale }

function rplus (a : cxptr; b : real) : cxptr;   { complex + real }
 external; function rplus;
     var result : cxptr;
     begin
        new(result);
        result@ .re := a@ .re + b;
        result@ .im := a@ .im;
        rplus := result;
        dispose (a);
     end; { rplus }

function jexp (a : real) : cxptr;     { e ** (j a) }
 external; function jexp;
     var result : cxptr;
     begin
        new(result);
        result@ .re := cos(a);
        result@ .im := sin(a);
        jexp := result;
        end; { jexp }

function angle (x : cxptr) : real;  { -270 <= angle < 90 degrees }
 external; function angle;
     var a : real;
     begin
        if x@.re = 0 then if x@.im > 0 then a:= -270 else a:= -90
                     else a:= arctan( x@.im / x@.re ) * 180 / pi;
        if x@.re < 0 then a := a - 180;
        angle := a;
     end; { angle }

function csqrt(a : cxptr) : cxptr;   { Complex square root }
 external; function csqrt;
 { This function is based on an algorithm describes in
   Numerical Recipes by William H. Press, et al, 1986, Cambridge
   University Press. }
     const tiny = 1e-15;
     var result : cxptr;
         Q, R : real;
     begin
        new(result);
        with a@ do begin
           if ( abs(re)<tiny ) and ( abs(im)<tiny )
             then begin
              result@.re := 0; result@.im := 0;
              end
            else begin
              if abs(re)>abs(im)
               then R := abs(re) * sqrt( 1 + sqr( im / re ) )
               else R := abs(im) * sqrt( 1 + sqr( re / im ) );
              Q := abs(re) + R;
              if re > 0
                then begin
                 result@.re :=sqrt(0.5*Q);
                 result@.im :=abs( im ) * sqrt(0.5/Q);
                 end { then }
                else begin
                 result@.re :=abs( im ) * sqrt(0.5/Q);
                 result@.im :=sqrt(0.5*Q);
                 end; { else }
              if im < 0
                then result@.im := -result@.im;
              end; { else ( csqrt non-zero ) }
           end; { with }
        csqrt := result;
        dispose (a);
     end; { csqrt }

function power( a : cxptr; b : real ) : cxptr;
{ Raises a to the power b  }
   begin { function power }
      power := scale( exp( 0.5 * b * ln(magsq( a )) ),
                      jexp( b * angle( a ) ) );
      end; { function power }
----------------------------------------------------------------------
LT Scott A. Norton, USN     | From Internet, if you need a gateway, use
Naval Postgraduate School   |    4526p%navpgs.bitnet@jade.berkeley.edu
Monterey, CA 93943-5018     | or 4526p%navpgs.bitnet@ucscc.ucsc.edu
4526P@NavPGS.BITNET         | The WISCVM gateway closed 15 Dec 87. )