arithmetic/src/basictypes/complex/ComplexBasic.mg

GENERIC MODULE ComplexBasic(R);
Arithmetic for Modula-3, see doc for details

Abstract: Complex numbers and basic operations 


FROM Arithmetic IMPORT Error;

<* UNUSED *>
CONST
  Module = "ComplexBasic.";

PROCEDURE Add (READONLY x, y: T; ): T =
  VAR z: T;
  BEGIN
    z.re := R.Add(x.re, y.re);
    z.im := R.Add(x.im, y.im);
    RETURN z;
  END Add;

PROCEDURE Sub (READONLY x, y: T; ): T =
  VAR z: T;
  BEGIN
    z.re := R.Sub(x.re, y.re);
    z.im := R.Sub(x.im, y.im);
    RETURN z;
  END Sub;

PROCEDURE Neg (READONLY x: T; ): T =
  VAR z: T;
  BEGIN
    z.re := R.Neg(x.re);
    z.im := R.Neg(x.im);
    RETURN z;
  END Neg;

PROCEDURE Conj (READONLY x: T; ): T =
  VAR z: T;
  BEGIN
    z.re := x.re;
    z.im := R.Neg(x.im);
    RETURN z;
  END Conj;

PROCEDURE IsZero (READONLY x: T; ): BOOLEAN =
  BEGIN
    RETURN R.IsZero(x.re) AND R.IsZero(x.im);
  END IsZero;

PROCEDURE Equal (READONLY x, y: T; ): BOOLEAN =
  BEGIN
    RETURN R.Equal(x.re, y.re) AND R.Equal(x.im, y.im);
  END Equal;

PROCEDURE Mul (READONLY x, y: T; ): T =
  VAR z: T;
  BEGIN
    z.re := R.Sub(R.Mul(x.re, y.re), R.Mul(x.im, y.im));
    z.im := R.Add(R.Mul(x.im, y.re), R.Mul(x.re, y.im));
    RETURN z;
  END Mul;
multiply x by the conjugate of y 
PROCEDURE MulConj (READONLY x, y: T; ): T =
  VAR z: T;
  BEGIN
    z.re := R.Add(R.Mul(x.re, y.re), R.Mul(x.im, y.im));
    z.im := R.Sub(R.Mul(x.im, y.re), R.Mul(x.re, y.im));
    RETURN z;
  END MulConj;

Use Karatsuba's trick tend to produce truncation errors

xr yr - xi yi xr yi + xi yr

(xr - xi) (yr + yi) - xr yi + xi yr 

<* UNUSED *>
PROCEDURE FastMul (READONLY x, y: T; ): T =
  VAR
    xyri    := R.Mul(R.Sub(x.re, x.im), R.Add(y.re, y.im));
    xryi    := R.Mul(x.re, y.im);
    xiyr    := R.Mul(x.im, y.re);
    z   : T;
  BEGIN
    z.re := R.Add(R.Sub(xyri, xryi), xiyr);
    z.im := R.Add(xryi, xiyr);
    RETURN z;
  END FastMul;

PROCEDURE DivScale (READONLY x, y: T; ): T RAISES {Error} =
  VAR denom := R.Add(R.Mul(y.re, y.re), R.Mul(y.im, y.im));
  BEGIN
    (* Err.divide_by_zero will be thrown by Div *)
    RETURN T{R.Div(x.re, denom), R.Div(x.im, denom)};
  END DivScale;

PROCEDURE Div (READONLY x, y: T; ): T RAISES {Error} =
  BEGIN
    RETURN DivScale(Mul(x, Conj(y)), y);
  END Div;

PROCEDURE Rec (READONLY x: T; ): T RAISES {Error} =
  BEGIN
    RETURN DivScale(Conj(x), x);
  END Rec;
I have not found a method which a) is generic (works without comparisons) b) matches the needs of the Euclidean algorithm for determining the greatest common divisor (that is, there must be a magnitude indicator which is reduced by a modulo operation). It work for some Gaussian numbers though, but not for all, e.g. GCD(-3-3i,-1-4i) (the order is important, the opposite would work!). The last example lead to a cycle: -3-3i,-1-4i,1-4i,3-3i,4-i,4+i,3+3i,... In this case the Euclidean algorithm may be rewritten to process the pair (x*u,y) for several units u in parallel. 
PROCEDURE Mod (READONLY x, y: T; ): T RAISES {Error} =
  VAR
    denom := R.Add(R.Mul(y.re, y.re), R.Mul(y.im, y.im));
    (* Err.divide_by_zero will be thrown by Mod *)
    xbig    := MulConj(x, y);
    r   : T;
  BEGIN
    r.re := R.Mod(xbig.re, denom);
    r.im := R.Mod(xbig.im, denom);
    r := Mul(r, y);              (* in fact, r is now AbsSqr(y) as big as
                                    before *)
    r.re := R.Div(r.re, denom);  (* is always divisible *)
    r.im := R.Div(r.im, denom);
    RETURN r;
  END Mod;

PROCEDURE DivMod (READONLY x, y: T; ): QuotRem RAISES {Error} =
  VAR
    denom := R.Add(R.Mul(y.re, y.re), R.Mul(y.im, y.im));
    (* Err.divide_by_zero will be thrown by Div *)
    xbig := MulConj(x, y);
    re   := R.DivMod(xbig.re, denom);
    im   := R.DivMod(xbig.im, denom);
    r := Mul(T{re.rem, im.rem}, y); (* in fact, r is now AbsSqr(y) as big
                                       as before *)
  BEGIN
    r.re := R.Div(r.re, denom);  (* is always divisible *)
    r.im := R.Div(r.im, denom);
    <* ASSERT Equal(r, Mod(x, y)) *>
    RETURN QuotRem{T{re.quot, im.quot}, r};
  END DivMod;

PROCEDURE Square (READONLY x: T; ): T =
  VAR
    xri: R.T;
    z  : T;
  BEGIN
    z.re := R.Sub(R.Mul(x.re, x.re), R.Mul(x.im, x.im));
    xri := R.Mul(x.re, x.im);
    z.im := R.Add(xri, xri);     (* double *)
    RETURN z;
  END Square;

PROCEDURE Scale (READONLY x: T; y: R.T; ): T =
  VAR z: T;
  BEGIN
    z.re := R.Mul(x.re, y);
    z.im := R.Mul(x.im, y);
    RETURN z;
  END Scale;

BEGIN
  Zero := T{re := R.Zero, im := R.Zero};
  One := T{re := R.One, im := R.Zero};
  I := T{re := R.Zero, im := R.One};
  MinusOne := T{re := R.MinusOne, im := R.Zero};
END ComplexBasic.