GENERIC MODULEMatrixSupport (R);

Arithmetic for Modula-3, see doc for details

<* INLINE *> PROCEDUREGiven the matrix of all weights of clows of length l compute the weight matrix for all clows of length (l+1). Take the result of 'newClow' as diagonal and the result of 'extendClow' as lower triangle of the weight matrix.AssertEqualWidth (n, m: CARDINAL; ) = BEGIN <* ASSERT n = m, "Width or height of operands don't match." *> END AssertEqualWidth; PROCEDURENewZero (m, n: CARDINAL; ): T = VAR z := NEW(T, m, n); BEGIN FOR i := FIRST(z^) TO LAST(z^) DO FOR j := FIRST(z[i]) TO LAST(z[i]) DO z[i, j] := R.Zero; END; END; RETURN z; END NewZero; PROCEDURENewOne (n: CARDINAL; ): T = VAR z := NEW(T, n, n); BEGIN FOR i := FIRST(z^) TO LAST(z^) DO z[i, i] := R.One; FOR j := FIRST(z[i]) TO i - 1 DO z[i, j] := R.Zero; z[j, i] := R.Zero; END; END; RETURN z; END NewOne; PROCEDURETranspose (x: T; ): T = VAR z := NEW(T, NUMBER(x[0]), NUMBER(x^)); BEGIN FOR i := FIRST(x[0]) TO LAST(x[0]) DO FOR j := FIRST(x^) TO LAST(x^) DO z[i, j] := x[j, i]; END; END; RETURN z; END Transpose;

Note that only the lower triangle of 'cl0' is used. It is somewhat waste of space, but saving space would lead to more complicated (time inefficient) array handling

PROCEDURECompute the determinant with about n^4 multiplications without division according to the clow decomposition algorithm of Mahajan and Vinay, and Berkowitz presented by Günter Rote:LongerClow (x: T; cl0: T; ): T = VAR cl1 := NEW(T, NUMBER(x^), NUMBER(x[0])); sum := R.Zero; BEGIN (* Compute the weights of all clow sequences where the last clow is closed and a new one is started. *) cl1[0, 0] := R.Zero; FOR i := FIRST(x^) TO LAST(x^) - 1 DO FOR j := i TO LAST(x^) DO sum := R.Sub(sum, R.Mul(x[i, j], cl0[j, i])); END; cl1[i + 1, i + 1] := sum; END; (* Compute the weights of all clow sequences where the last (open) clow is extended by a new arc. This is essentially a multiplication of the matrix 'x' with the lower triangular matrix 'cl0' where the result is restricted to a lower triangular matrix without the diagonal. *) FOR i := FIRST(x^) + 1 TO LAST(x^) DO FOR j := 0 TO i - 1 DO sum := R.Zero; FOR k := j TO LAST(x^) DO sum := R.Add(sum, R.Mul(x[i, k], cl0[k, j])); END; cl1[i, j] := sum; END; END; RETURN cl1; END LongerClow;

```
Division-Free Algorithms for the
Determinant and the Pfaffian: Algebraic and Combinatorial
Approaches.
```

PROCEDUREDeterminant (x: T; ): R.T = VAR y := NewOne(NUMBER(x^)); sum := R.Zero; BEGIN AssertEqualWidth(NUMBER(x^), NUMBER(x[0])); FOR i := 1 TO LAST(x^) DO y := LongerClow(x, y); END; (* This is equal to the first part of LongerClow. *) FOR i := FIRST(x^) TO LAST(x^) DO FOR j := i TO LAST(x^) DO sum := R.Sub(sum, R.Mul(x[i, j], y[j, i])); END; END; IF NUMBER(x^) MOD 2 = 0 THEN RETURN sum; ELSE RETURN R.Neg(sum); END; END Determinant; BEGIN END MatrixSupport.