The ModulaTor logo 7KB

The ModulaTor

Oberon-2 and Modula-2 Technical Publication

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]

Webdesign by www.otolo.com/webworx, 14-Jul-1998