4526P%NAVPGS.BITNET@CUNYVM.CUNY.EDU (LT Scott A. Norton, USN) (01/12/88)
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. )shankar@hpclscu.HP.COM (Shankar Unni) (01/14/88)
/ hpclscu:comp.lang.pascal / 4526P%NAVPGS.BITNET@CUNYVM.CUNY.EDU (LT Scott A. Norton, USN) / 6:42 pm Jan 11, 1988 / > 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. This is one among several big hits against Pascal, at least as it is defined by the ANSI folks. Lots of other features (seperate compilation, type coercion, pointer manipulation, extensibility, etc) are also missing. There is an evolving standard for Extended Pascal (also by the ANSI committee) which has most of these features, and a whole lot more (to the point where it looks like a union of Pascal, FORTRAN and COBOL (shudder!!)). But the most important "feature", i.e. extensibility of the language (ability to define and overload operators, add "predefines" (inline procedures)) is still absent. C++ (or objective C) is much better in this regard, giving the programmer the ability to define classes of objects, and overloaded functions and operators (a lot like ADA) to generate inline code for these types. As a starter, the compiler comes with a few such types implemented as add-on include files, to give the programmer a good example to base his/her modules on. One of these, of course, is a "complex" module. I realise, of course, that this is no consolation to you. As Pascal stands today (or even as it stands tomorrow, for that matter), you may have to define "LISP-like" functions to do your stuff. BTW, most Pascals allow users to return structures as function return values, so at least you should be able to get rid of most of the new/dispose stuff: simply define functions like: function add (var c1, c2 : complex) : complex; ... ----------------------------------------------------------------------------- Shankar Unni allegra) Unix Mail: shankar@hpda.hp.com UUCP: ucbvax)!hplabs!hpda!shankar sun) (AT&T: (408) 447-5797) decwrl) -----------------------------------------------------------------------------
bobdi@omepd (Bob Dietrich) (01/15/88)
I realize this is not an immediate help to you, but the draft proposed Extended Pascal Standard has complex numbers built-in as a "simple" type (e.g., infix operators are supported), along with functions operating on complex numbers. You can also have parameters and functions of type complex, as with any data type. Incidently, Extended Pascal allows functions to be of most any structured type (except, of course, files, since they are not an assignable type). Bob Dietrich Intel Corporation, Hillsboro, Oregon (503) 696-2092 (messages x2111) usenet: tektronix!ogcvax!omepd!ihf1!bobd or tektronix!psu-cs!omepd!ihf1!bobd or ihnp4!verdix!omepd!ihf1!bobd