The ModulaTor
Erlangen's First Independent Modula-2 Journal! Nr. 1/Feb-1991
_____________________________________________________________
Modula-2 and the COMPLEX Type
An Evaluation of the Compiler's Complex Extension
by Guenter Dotzel, ModulaWare
The ISO Modula-2 Standardization Group decided in summer 1990 to extend the
language by the COMPLEX data type and associated operations. I was not happy with
the inclusion of complex in Modula-2, since complex arithmetic could be done easily with
functions returning complex values. This worked for years on the VAX for single and
double precision complex numbers. (c1+c2) it is not much shorter or easier to write or
read than cadd(c1,c2). But I see that procedure calls and parameter passing may cost
execution time. To see what is gained, I've therefore done some benchmarks to assure,
that complex is worth the effort. It took me a full week to extend the VAX/VMS Modula-2
compiler MVR to fully support complex including evaluation of constant complex
operations and VMS-debugger support for the new types.
100000 operations of each complex addition, multiplication and divison were performed.
The result of the benchmark is given in CPU seconds as measured on a VAXstationII
(figure 0). Column 1 is the time needed to perform the computation of c3 := c1+c2, c3 :=
c1*c2, and c3 := c1/c2 using the new complex type and it's associated operators.
Column 2 is c3 := CADD(c1,c2); c3 := CMUL(c1,c2); c3 := CDIV(c1,c2) using the new
complex type and function procedures programmed in Modula-2 as shown below. I did
use the old fashioned type transfer instead of ISO-Modula-2's SYSTEM.CAST. Old
habits die hard! Also INC is applied to real-type operands. This is a non-standard
Modula-2 language extension of our VAX/VMS implementation. Instead of c.r := c.r+d.r I
used the equivalent INC(c.r, d.r) to cosmetically get the best code possible (2-operand
instead of 3-operand ADDF-instruction is generated, see Appendix I), but this does by no
way influence the benchmark result. Column 3 is the speed-up gain, i.e. the value of
column 2 divided by that of column 1. The variables c1, c2, and c3 are of type COMPLEX
and were initialized to non-zero values. The data type COMPLEX is the new pervasive
type. The VAX machine code generated by MVR V3.01 for these procedures and their
calling sequence is listed in the Appendix I below):
TYPE reim=RECORD r,i: REAL END;
PROCEDURE CADD(a,b: COMPLEX): COMPLEX;
VAR c,d: reim;
BEGIN c:= reim(a); d:= reim(b);
INC(c.r, d.r); INC(c.i, d.i);
RETURN COMPLEX(c);
END CADD;
PROCEDURE CMUL (x,y: COMPLEX): COMPLEX;
VAR res,c1,c2: reim;
BEGIN c1 := reim(x); c2 := reim(y);
res.r := c1.r*c2.r - c1.i*c2.i;
res.i := c1.r*c2.i + c1.i*c2.r;
RETURN COMPLEX(res);
END CMUL;
PROCEDURE CDIV (x,y: COMPLEX): COMPLEX;
VAR res,c1,c2: reim;
BEGIN c1 := reim(x); c2 := reim(y);
res.i := c2.r*c2.r + c2.i*c2.i;
res.r := (c1.r*c2.r + c1.i*c2.i) / res.i;
res.i := (c1.i*c2.r - c1.r*c2.i) / res.i;
RETURN COMPLEX(res);
END CDIV;
COMPLEX arithmetic
CPU-times/s Speed-up gain
__________________________________________________
Column 1 2 3
__________________________________________________
Addition 2.07 4.34 2.1
Multiplication 4.35 6.49 1.5
Division 7.03 8.74 1.2
Figure 0: Build-in complex operators versus function procedures
The values of figure 0 show, that it is a clear advantage to have complex arithmetic
support in the language. It has to be considered, that most gain is with complex addition,
subtraction, negation and comparison operations. These operations are very simple and
where execution speed is a concern they can be programmed explicitely. For example all
calls of CADD in innermost loops can be replaced by the complex addition operation.
This saves the function procedure calling overhead. Under the provision that the
compiler generates good code (the quality can be inspected in Appendix I) the gain of a
build-in complex addition operation is zero.
Evaluation of Modula-2's Complex Extensions
The Standard has provisions for the data type COMPLEX, LONGCOMPLEX, two
functions to extract the real and imaginary part of a [long-]complex number, namely RE
and IM and a function CMPLX to construct a complex number from two [long-]real
numbers. Furthermore the operations +, -, *, /, negation and comparison (for equality
only) apply to complex numbers. Two library modules [Long-]ComplexMath provide the
functions Sqrt, Modulus, Argument and PolarToComplex.
After having used complex numbers in Modula-2 for a week, mainly to test the new
language features of the compiler, I'm concerned about the ISO-Modula-2 proposal on
the following details:
1. ComplexMath: The functions of module ComplexMath (see below) should be
included in RealMath. There is no need for the additional module. The module is too
simple and is easily implemented. I suggest a further extension of RealMath by some
real and complex functions. In contrast to the current position of the Standard, these
extensions can't be left for commercial support whatever that may be. A language
standard is needed and requested by the industry which by the way is commercial. The
Standard lists more than 10 pages text on threatening of the For-statement's control
variable, but complex is to be treated on 3 pages, including library. So if complex is being
included, then it has to be done in orthogonal fashion to be useful.
The following functions should be included in RealMath (see also paper from Charles R.
Bilbe, in Journal of Pascal, Ada and Modula-2, Nov/Dec-1985 and Scott Woodard's
complex number library on HP9000): arctan2, (* hyperbolic functions: *) arctanh, cosh,
sinh, tanh, log, log2, mod (* real modulus *), (* complex functions *) cabs (* absolute
value, renamed from Modulus *), csqrt (* square root, renamed from Sqrt *), conj (*
conjugate *), cinv (* inversion *), scmul (* multiplied by scalar real *), scdiv (* divided by
scalar real *), ccos, cexp, cln, csin, cpower.
Furthermore, the function names Argument, and PolarToComplex should be renamed to
obey the naming conventions of module RealMath (lower case letters only, short
names).
2. LongComplexMath: The module name is too long to be used in qualified mode. It's
boring to write LongComplexMath.PolarToComplex. The functions of module
LongComplexMath should be renamed and included in LongMath. See ComplexMath
above.
3. Complex number literals and constant expressions: It should be clarified what type
is constructed by CMPLX(const_real_expr, const_real_expr) and CMPLX(RR-type,
RR-type). For example:
CMPLX (1.0, 2.0) is of CC-type (complex number literal, compatible with both, complex-
and longcomplex-type and both number parts evaluated with at least longreal precision).
CMPLX(1.0, FLOAT(2)) is of complex-type and not of CC-type, since
FLOAT(const_expr) is always of real-type and never of RR-type. Note, that complex-
and longcomplex-type are incompatible types.
CMPLX(1.0, 2.0) / CMPLX(2.0, 1.0) is a constant complex expression of CC-type. Any
overflow that occurs when evaluating constant complex expressions at compile time
must be detected. Possible exceptions when evaluating complex constant expression
are: overflow, division by zero, and also reserved operand fault, which can for example
occur when declaring a complex constant number with a value of CMPLX( SYS-
TEM.CAST(REAL, 80000000H), 2.0), since the real-part is -0.0 which is the reserved
floating operand on the VAX.
No other notation for the complex numeric literal than CMPLX (..., ...) is needed. It seems
to me, that the paper of the ISO Working Group on Modula-2 and the Complex Type
(dated 16-Aug-1990) is not specific enough on that detail. Also specification for the
detection of mandatory and non-mandatory exceptions is missing.
4. Type of RE and IM: It should be clarified whether RE(const_complex_expr) and
IM(const_complex_expr) are of RR-type (real number literal) or constants of real- or
longreal-type, e.g.: in the declaration CONST c = CMPLX(1.0, FLOAT(2)); r=RE(c); is r of
RR- or real-type?
5. Conversion from complex to longcomplex and vice versa: A conversion from
values of data type COMPLEX to LONGCOMPLEX and vice versa can currently be
accomplished only by means of e.g.:
lc := CMPLX(LFLOAT(RE(c)), LFLOAT(IM(c));
c := CMPLX(FLOAT(RE(lc)), FLOAT(IM(lc));
where the variable lc is of type LONGCOMPLEX and c of type COMPLEX. It is
suggested, that VAL is applicable to complex data types too, e.g.:
lc := VAL(LONGCOMPLEX, c);
c := VAL(COMPLEX, lc);
To be orthogonal two new standard functions are needed which correspond to the VAL
operations, namely LCMPLX and RCMPLX respectively, e.g.:
lc := LCMPLX(c);
c:= RCMPLX(lc).
To extend the compiler by these two functions takes the compiler writer only a couple of
hours.
Appendix I: Code generated for CADD, CMUL CDIV:
16-NOV-1990 12:49:31 Test_Complex
16-NOV-1990 12:36:28 DUA0:[000000.LIDAS]TEST_COMPLEX.MOD;44
VAX/VMS Modula-2 MVR V3.0 (C) by ModulaWare GmbH; Machine Code
; Symbols for procedure CADD
FFFFFFF0 b = -16
FFFFFFF8 a = -8
FFFFFFE8 c = -24
FFFFFFE0 d = -32
00363 CADD:
0000 00363 .WORD ^M<> ;61
5E 20 C2 00365 SUBL2 #32, SP
F8 AD 04 BC 7D 00368 MOVQ @4(AP), a(FP)
F0 AD 08 BC 7D 0036D MOVQ @8(AP), b(FP)
E8 AD F8 AD 7D 00372 MOVQ a(FP), c(FP)
E0 AD F0 AD 7D 00377 MOVQ b(FP), d(FP)
E8 AD E0 AD 40 0037C ADDF2 d(FP), c(FP) ;62
EC AD E4 AD 40 00381 ADDF2 d+4(FP), c+4(FP) ;63
50 E8 AD 7D 00386 MOVQ c(FP), R0 ;64
04 0038A RET
08019F5C 8F DD 0038B PUSHL #134324060 ;65
00000000' EF 01 FB 00391 CALLS #1, LIB$SIGNAL
01 00398 NOP
; Symbols for procedure CMUL
FFFFFFF8 x = -8
FFFFFFF0 y = -16
FFFFFFE8 res = -24
FFFFFFE0 c1 = -32
FFFFFFD8 c2 = -40
00399 CMUL:
0600 00399 .WORD ^M<R9,R10> ;69
5E 28 C2 0039B SUBL2 #40, SP
F8 AD 04 BC 7D 0039E MOVQ @4(AP), x(FP)
F0 AD 08 BC 7D 003A3 MOVQ @8(AP), y(FP)
E0 AD F8 AD 7D 003A8 MOVQ x(FP), c1(FP)
D8 AD F0 AD 7D 003AD MOVQ y(FP), c2(FP)
5A E0 AD D8 AD 45 003B2 MULF3 c2(FP), c1(FP), R10 ;72
59 E4 AD DC AD 45 003B8 MULF3 c2+4(FP), c1+4(FP), R9
E8 AD 5A 59 43 003BE SUBF3 R9, R10, res(FP)
5A E0 AD DC AD 45 003C3 MULF3 c2+4(FP), c1(FP), R10
59 E4 AD D8 AD 45 003C9 MULF3 c2(FP), c1+4(FP), R9
EC AD 5A 59 41 003CF ADDF3 R9, R10, res+4(FP)
50 E8 AD 7D 003D4 MOVQ res(FP), R0
04 003D8 RET
08019F5C 8F DD 003D9 PUSHL #134324060 ;75
00000000' EF 01 FB 003DF CALLS #1, LIB$SIGNAL
01 003E6 NOP
; Symbols for procedure CDIV
FFFFFFF8 x = -8
FFFFFFF0 y = -16
FFFFFFE8 res = -24
FFFFFFE0 c1 = -32
FFFFFFD8 c2 = -40
003E7 CDIV:
0600 003E7 .WORD ^M<R9,R10> ;79
5E 28 C2 003E9 SUBL2 #40, SP
F8 AD 04 BC 7D 003EC MOVQ @4(AP), x(FP)
F0 AD 08 BC 7D 003F1 MOVQ @8(AP), y(FP)
E0 AD F8 AD 7D 003F6 MOVQ x(FP), c1(FP)
D8 AD F0 AD 7D 003FB MOVQ y(FP), c2(FP)
5A D8 AD D8 AD 45 00400 MULF3 c2(FP), c2(FP), R10 ;82
59 DC AD DC AD 45 00406 MULF3 c2+4(FP), c2+4(FP), R9
EC AD 5A 59 41 0040C ADDF3 R9, R10, res+4(FP)
5A E0 AD D8 AD 45 00411 MULF3 c2(FP), c1(FP), R10
59 E4 AD DC AD 45 00417 MULF3 c2+4(FP), c1+4(FP), R9
5A 59 40 0041D ADDF2 R9, R10
E8 AD 5A EC AD 47 00420 DIVF3 res+4(FP), R10, res(FP)
5A E4 AD D8 AD 45 00426 MULF3 c2(FP), c1+4(FP), R10
59 E0 AD DC AD 45 0042C MULF3 c2+4(FP), c1(FP), R9
5A 59 42 00432 SUBF2 R9, R10
EC AD 5A EC AD 47 00435 DIVF3 res+4(FP), R10, res+4(FP)
50 E8 AD 7D 0043B MOVQ res(FP), R0
04 0043F RET
08019F5C 8F DD 00440 PUSHL #134324060 ;86
00000000' EF 01 FB 00446 CALLS #1, LIB$SIGNAL
01 0044D NOP
Calling sequence for either procedure (c3 := CMUL(c1,c2) shown)
00000000' EF DF 01568 93$: PUSHAL c2
00000000' EF DF 0156E PUSHAL c1
00000000' EF 02 FB 01574 CALLS #2, CMUL
00000000' EF 50 7D 0157B MOVQ R0, c3
ComplexMath
__________________________________________________________________________________________________
DEFINITION MODULE ComplexMath;
(* Mathematical functions for the type COMPLEX *)
CONST
i = CMPLX (0.0, 1.0);
one = CMPLX (1.0, 0.0);
zero = CMPLX (0.0, 0.0);
PROCEDURE abs (z: COMPLEX): REAL;
(* Returns the length of z *)
PROCEDURE arg (z: COMPLEX): REAL;
(* Returns the angle that z subtends to the positive real axis *)
PROCEDURE conj (z: COMPLEX): COMPLEX;
(* Returns the complex conjugate of z *)
PROCEDURE power (base: COMPLEX; exponent: REAL): COMPLEX;
(* Returns the value of the number base raised to the power exponent *)
PROCEDURE sqrt (z: COMPLEX): COMPLEX;
(* Returns the principal square root of z *)
PROCEDURE exp (z: COMPLEX): COMPLEX;
(* Returns the complex exponential of z *)
PROCEDURE ln (z: COMPLEX): COMPLEX;
(* Returns the principal value of the natural logarithm of z *)
PROCEDURE sin (z: COMPLEX): COMPLEX;
(* Returns the sine of z *)
PROCEDURE cos (z: COMPLEX): COMPLEX;
(* Returns the cosine of z *)
PROCEDURE tan (z: COMPLEX): COMPLEX;
(* Returns the tangent of z *)
PROCEDURE arcsin (z: COMPLEX): COMPLEX;
(* Returns the arcsine of z *)
PROCEDURE arccos (z: COMPLEX): COMPLEX;
(* Returns the arccosine of z *)
PROCEDURE arctan (z: COMPLEX): COMPLEX;
(* Returns the arctangent of z *)
PROCEDURE polarToComplex (abs, arg: REAL): COMPLEX;
(* Returns the complex number with the specified polar coordinates *)
PROCEDURE scalarMult (scalar: REAL; z: COMPLEX): COMPLEX;
(* Returns the scalar product of scalar with z *)
PROCEDURE IsCMathException (): BOOLEAN;
(* Returns TRUE if the current coroutine is in the exceptional execution state
because of the raising of an exception in a routine from this module; otherwise
returns FALSE.
*)
(* not yet impl. on MVR | MAX *)
END ComplexMath.
LongComplexMath
__________________________________________________________________________________________________
DEFINITION MODULE LongComplexMath;
(* Mathematical functions for the type LONGCOMPLEX *)
CONST
i = CMPLX (0.0, 1.0);
one = CMPLX (1.0, 0.0);
zero = CMPLX (0.0, 0.0);
PROCEDURE abs (z: LONGCOMPLEX): LONGREAL;
(* Returns the length of z *)
PROCEDURE arg (z: LONGCOMPLEX): LONGREAL;
(* Returns the angle that z subtends to the positive real axis *)
PROCEDURE conj (z: LONGCOMPLEX): LONGCOMPLEX;
(* Returns the complex conjugate of z *)
PROCEDURE power (base: LONGCOMPLEX; exponent: LONGREAL): LONGCOMPLEX;
(* Returns the value of the number base raised to the power exponent *)
PROCEDURE sqrt (z: LONGCOMPLEX): LONGCOMPLEX;
(* Returns the principal square root of z *)
PROCEDURE exp (z: LONGCOMPLEX): LONGCOMPLEX;
(* Returns the complex exponential of z *)
PROCEDURE ln (z: LONGCOMPLEX): LONGCOMPLEX;
(* Returns the principal value of the natural logarithm of z *)
PROCEDURE sin (z: LONGCOMPLEX): LONGCOMPLEX;
(* Returns the sine of z *)
PROCEDURE cos (z: LONGCOMPLEX): LONGCOMPLEX;
(* Returns the cosine of z *)
PROCEDURE tan (z: LONGCOMPLEX): LONGCOMPLEX;
(* Returns the tangent of z *)
PROCEDURE arcsin (z: LONGCOMPLEX): LONGCOMPLEX;
(* Returns the arcsine of z *)
PROCEDURE arccos (z: LONGCOMPLEX): LONGCOMPLEX;
(* Returns the arccosine of z *)
PROCEDURE arctan (z: LONGCOMPLEX): LONGCOMPLEX;
(* Returns the arctangent of z *)
PROCEDURE polarToComplex (abs, arg: LONGREAL): LONGCOMPLEX;
(* Returns the complex number with the specified polar coordinates *)
PROCEDURE scalarMult (scalar: LONGREAL; z: LONGCOMPLEX): LONGCOMPLEX;
(* Returns the scalar product of scalar with z *)
PROCEDURE IsCMathException (): BOOLEAN;
(* Returns TRUE if the current coroutine is in the exceptional execution state
because of the raising of an exception in a routine from this module; otherwise
returns FALSE.
*)
(* not yet impl. on MVR | MAX *)
END LongComplexMath.
__________________________________________________________________________________________________
IMPRESSUM: The ModulaTor is an unrefereed journal. Technical papers are to be
taken as working papers and personal rather than organizational statements.
Items are printed at the discretion of the Editor based upon his judgement on
the interest and relevancy to the readership. Letters, announcements, and
other items of professional interest are selected on the same basis. Office of
publication: The Editor of The ModulaTor is Guenter Dotzel; he can be reached
by tel/fax: [removed due to abuse] or by
mailto:[email deleted due to spam]
ModulaWare home page
The ModulaTor download
![[Any browser]](starsbtn.gif){kind=small}