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