% from 20.3.1989 hk, version 0.03
% revised 5.11.1989 hk, version 0.30
% arithmetic part 26.1.1990 hk


%\setcounter{chapter}{4}
\chapter{Basic Arithmetic}

The MAS system consists of a kernel with
list processing, input / output facilities and
an interaction part. The MAS language of
the interaction part supports only 'primitive' 
arithmetic operations i.e. the
operators \verb/+/, \verb/-/, \verb/*/, \verb./., etc.
are only valid on $\beta$--integers (atoms).
All further arithmetic functions must be imported 
from Modula--2 libraries. 

Available basic arithmetic libraries are:
\index{basic arithmetic}\index{arithmetic}
\begin{quote}
\begin{obeylines} 
  ALDES/SAC--2 Digit Arithmetic System,
  ALDES/SAC--2 Integer Arithmetic System,
  ALDES/SAC--2 Rational Number System,
  ALDES/SAC--2 Modular Integer Arithmetic System, 
  ALDES/SAC--2 Integer Factorization System, 
  ALDES/SAC--2 Set of Integers System, 
  ALDES/SAC--2 Combinatorical System, 
  MAS          Arbitrary Precision Floating Point System, 
\end{obeylines}
\index{digit}\index{modular integer}
\index{integer}\index{rational number}
\index{integer factorization}\index{set of integer}
\index{floating point}
\end{quote}

Which functions are accessible from this libraries
can be determined by the MAS \verb/HELP/ function.\index{HELP} 
The listings of the definition modules of these libraries
is contained in the folder \verb/\HELP/ on the MAS
distribution disk. 

In the following we will discuss the 
libraries for {\em arbitrary precision integers}, 
{\em arbitrary precision rational numbers} and 
{\em arbitrary precision floating point numbers}.
Polynomial arithmetic will be discussed later. 

For further reading the following documents are recommended:
\begin{enumerate}
\item D. E. Knuth: {\em The Art of Computer Programming},
      Vol. II: {\em Seminumerical Algorithms},
      Chap. 4: {\em Arithmetic}. \\
      This book contains an in--depth treatment of most
      basic arithmetic algorithms.
\item G. E. Collins, M. Mignotte, F. Winkler: 
      {\em Arithmetic in Basic Algebraic Domains}.
      In B. Buchberger, G. E. Collins, R. Loos:
      {\em Computer Algebra}, Computing Supplement 4.
      This article contains an overview on 
      the computational best algebraic algorithms.  
\end{enumerate}
It is further recommended to look into the
original ALDES / SAC--2 algorithm documentation  
or into the translated MAS libraries in Modula--2.


\section{Integer Arithmetic}

The arbitrary precision integers are in the following 
called {\bf integers}. The small integers are explicitly
denoted as $\beta$--integers. 

We will first discuss the representation of integers by
lists. Although it is not required to understand 
the representation to use the functions on integers, 
it is required to understand the algorithms.
Further the representation has important influence 
on the computing time, i.e. the complexity of
the algorithm. 

The base for the number representation of integers 
is $\beta$.

{\bf Note:} 
Let $A \in \Z$, $A \geq 0$ then there
exist unique numbers $n$, $a_0, \ldots, a_{n-1}$ $\in \Z$,
with $0 \leq a_i < \beta$ $(i=0, \ldots, n-1)$ 
and $a_{n-1} \neq 0$ such that:
\begin{displaymath}
       A = \sum_{i=0}^{n-1} a_i \beta^i.
\end{displaymath}
For $A < 0$ the same representation
holds under the condition 
$-\beta < a_i \leq 0$ $(i=0, \ldots, n-1)$ 
and $a_{n-1} \neq 0$.

%{\bf Proof:} Induction on $\vert A \vert$ by division with remainder.

{\bf Definition:} List representation of integers.
\begin{quote}
       $A = 0$ is represented by the atom $0$.
       \\
       $A \neq 0$, $\vert A \vert < \beta$, is represented 
       by the atom $a_0$, and 
       \\
       $A \neq 0$, $\vert A \vert \geq \beta$, is represented 
       by the list 
       $(a_0, \ldots, a_{n-1})$.
\end{quote}
\index{list representation}\index{integer representation}
\index{representation}
The `least significant' $\beta$--digit is the first atom 
in the list representation. The `most significant' 
digit appears as last element in the list representation.
So e.g. $A$ is even iff $a_0$ (the first element in the list)
is even. 

Example: 
\begin{quote}
      $\beta = 0 \cdot \beta^0 + 1 \cdot \beta^1$, so $n = 2$ and 
      the representation is $(0,1)$. 
      \\
      $-\beta = 0 \cdot \beta^0 + (-1) \cdot \beta^1$, 
      the list representation is $(0,-1)$. 
\end{quote}


\subsection{Algorithms}

The programs of the most important integer algorithms 
and their complexity are summarized next.
Therefore let
$\A$ be the set of atoms,
$\L$ be the set of lists,
$\O = \A \cup \L$ be the set of objects,
$ \I = \{ x \in \O: \ x$ represents an element of ${\rm \bf Z} \}$ 
be the set of integers 
and $\L(\I)$ be the set of lists over integers.
\index{integer algorithms}\index{integer}
Further let $L(a)$ denote the number of $\beta$--digits 
of $a$ 
(if $a \geq \beta$ then this is equal to the list 
length of the representation of $a$).
We will also write $L(a)$ for $O(L(a))$, i.e. we will not count for 
constant factors.
The computing time functions $t, \tma, \tmi, \ta$ 
are defined as before in section \ref{secCOMPL}.
The time analyses are due to G. E. Collins.

\begin{deflist}{$c \gets IPRODK(a,b)$}
\item[$b \gets INEG(a)$] $a, b \in \I$. 
     $b = -a$ is the negative of $a$. 
     Every $\beta$--digit of $a$ must be processed
     to change the sign, so
     $t = L(a)$, $c = L(a)$.
     \index{INEG}
\item[$s \gets ISIGN(a)$] $a \in \I$, $s \in \{ -1, 0, 1 \}$. 
     $s$ is the sign of $a$. 
     As long as $\beta$--digits of $a$ are zero, the list must 
     be processed, so
     $\tma = L(a)$, $\tmi = 1$, $\ta = 1$, $c = 0$.
     \index{ISIGN}
\item[$b \gets IABS(a)$] $a, b \in \I$. 
     $b = \vert a \vert$ is the absolute value of $a$. 
     If ${\rm sign}(a) = -1$ then $a$ must change the sign, so
     $\tma = L(a)$, $\tmi = 1$, $\cma = L(a)$, $\cmi = 0$.
     \index{IABS}
\item[$s \gets ICOMP(a,b)$] $a, b \in \I$, $s \in \{ -1, 0, 1 \}$. 
     $s = {\rm sign}(a-b)$ is the sign of $a-b$. 
     At most the smaller number of digits of $a$ and $b$ must be
     compared, so
     $\tma = {\rm min}\{ L(a), L(b) \}$, 
     $\tmi = 1$, $c = 0$.
     \index{ICOMP}
\item[$c \gets ISUM(a,b)$] $a, b, c \in \I$.
     $c = a + b$ is the sum of $a$ and $b$. 
     At least the smaller number of digits of $a$ and $b$ must be
     added, so
     $\tmi = {\rm min}\{ L(a), L(b) \}$, 
     $\cmi = {\rm min}\{ L(a), L(b) \}$. 
     If always carrys occur, the maximal computing time is 
     proportional to
     the greater number of digits of $a$ and $b$, so
     $\tma = {\rm max}\{ L(a), L(b) \}$, 
     $\cma = {\rm max}\{ L(a), L(b) \}$. 
     However carrys occur `seldom', so
     $\ta = \tmi$, $\ca = \cmi$. 
     \index{ISUM}
\item[$c \gets IDIF(a,b)$] $a, b, c \in \I$.
     $c = a - b$ is the difference of $a$ and $b$. 
     $b$ is negated, then the sum is computed, so
     $\tmi = L(b) $, $\cmi = L(b) $. 
     The maximal computing time is proportional to 
     the maximal computing time of $ISUM$, so 
     $\tma = {\rm max}\{ L(a), L(b) \}$, 
     $\cma = {\rm max}\{ L(a), L(b) \}$. 
     \index{IDIF}
\item[$c \gets IPROD(a,b)$] $a, b, c \in \I$.
     $c = a \cdot b$ is the product of $a$ and $b$.
     The product is computed by the classical method, multiplying
     each $\beta$--digit of $a$ by each $\beta$--digit of $b$, so 
     $\tma = L(a) \cdot L(b)$ $= L(a)^2$ if $L(a) = L(b)$.   
     $\cma = L(a) + L(b)$ since the cells of the result are reused
     and not taken from the available cells during summation. 
     \index{IPROD}
\item[$c \gets IPRODK(a,b)$] $a, b, c \in \I$.
     $c = a \cdot b$ is the product of $a$ and $b$.
     The product is computed by Karatsubas method, 
     splitting each integer in two halves and then using only
     3 multiplications, 4 summations and some shifting 
     to obtain the product.
     More precisely
     if $a = a_1 \lambda + a_0$, $b = b_1 \lambda + b_0$ 
     where $\lambda$ such that $L(a_1) \approx L(a_0)$, then
     $$a \cdot b = 
            a_1 b_1 \lambda^2
          + ( a_1 b_1 
              + (a_1 - a_0) (b_0 - b_1) 
              +  a_0 b_0 ) \lambda
          + a_0 b_0.$$
     If $L(a) = L(b)$ then    
     $\tma = L(a)^{\log_2(3)} = L(a)^{1.585\ldots}$. 
     However $\cma = L(a) \cdot L(b)$, since the summation algorithms
     are used explicitly.
     \\
     {\em Trade--off:} $IPRODK$ is only superior to $IPROD$ if 
     the length of the integers exceed $55$ $\beta$--digits,
     i.e. the integers exceed about $500$ decimal digits 
     ($=$ $8$ lines full of decimal digits on the screen).   
     \index{IPRODK}
\item[$IQR(a,b, q,r)$] $a, b, q, r \in \I$.
     $b \neq 0$, $q = \lfloor a/b \rfloor$, $r = a - q \cdot b$.
     The classical division method (dividing the leading
     $\beta$--digits and subtracting $b$ times the 
     trial quotient digit)
     has been refined to the division of the two leading 
     digits and subtracting $b$ times the trial quotient digit.
     The advantage of this refinement is that the trial
     quotient digit is in nearly all cases the true quotient digit 
     of the quotient, so seldom adjustments of the remainder 
     are necessary. 
     The computing time is proportional to 
     the time for the computation of the product of $b$ and $q$, so  
     $\tma = L(b) \cdot L(q)$ and $\cma = L(b) \cdot L(q)$.
     \index{IQR}
\item[$q \gets IQ(a,b)$] $a, b, q \in \I$.
     $b \neq 0$, $q = \lfloor a/b \rfloor$. 
     Internally $IQR$ is called, see $IQR$ for the
     computing time.
     \index{IQ}
\item[$r \gets IREM(a,b)$] $a, b, r \in \I$.
     $b \neq 0$, $r = a - \lfloor a/b \rfloor \cdot b$.
     Internally $IQR$ is called, see $IQR$ for the
     computing time.
     \index{IREM}
\item[$c \gets IEXP(a,n)$] $a, c \in \I$, $0 \leq n \in \A$.
     $c = a^n$ is the $n$--th power of $a$ (exponentiation).
     $c$ is computed by a binary exponentiation method,
     i.e. $c$ is recursively computed as 
     $(a^{\lfloor n/2 \rfloor})^2 \cdot a^{n - \lfloor n/2 \rfloor 2}$.
     The computing time is proportional to 
     the time for the computation of the `biggest' product, so  
     $\tma = (n \cdot L(a))^2$, $\cma = 2 \cdot n \cdot L(a)$, 
     $2$ because of the classical multiplication algorithm.
     \index{IEXP}
\item[$c \gets IGCD(a,b)$] $a, b, c \in \I$.
     $c = {\rm gcd}(a,b)$ is the greatest common divisor of $a$ and $b$.
     The classical Euclidean algorithm uses the invariant 
     ${\rm gcd}(a,b) = {\rm gcd}(b,{\rm rem}(a,b))$ for 
     the computation of the gcd until $b = 0$. 
     Assume w.l.o.g. $L(a) \geq L(b)$, 
     then the maximal computing time is  
     $\tma = L(b) \cdot ( L(a) - L( {\rm gcd}(a,b) ) + 1 )$.
     Thus if ${\rm gcd}(a,b) = 1$ and $L(a) = L(b)$ then 
     $\tma = L(a)^2$ and $\cma = L(a)^2$.
     The number of division steps is bounded by 
     $\lceil \log_{\phi}(\sqrt{5} a) \rceil - 2$  
     $\approx 4.8 \log_{10}(a) - 0.32$. 
     The average number of division steps is 
     given by 
     $\frac{12 \ln(2)}{\pi^2} \ln(a)$ 
     $\approx 1.9 \log_{10}(a)$. 
     \index{IGCD}
\item[$c \gets ILCM(a,b)$] $a, b, c \in \I$.
     $c = {\rm lcm}(a,b)$ is the least common multiple of $a$ and $b$.
     $c$ is computed as $\frac{a}{{\rm gcd}(a,b)} \cdot b$.
     So the maximal computing time is  
     $\tma = 2 L(a)^2$.  
     \index{ILCM}
\item[$a \gets IRAND(n)$] $0 \leq n \in \A$, $a \in \I$.
     $a$ is a random integer with random sign and 
     $\vert a \vert < 2^n$. $t = L(a)$ and $c = L(a)$.
     \index{IRAND}
\item[$a \gets IREAD()$] $a \in \I$.
     An integer $a$ is read from the actual input stream.
     A conversion from decimal representation to 
     $\beta$ representation is done, so
     $\tma = L(a)^2$ and $\cma = L(a)^2$. \\
     The syntax for integers is:
\begin{verbatim}
int = [ "+" | "-" ] 
      [ "0" | ... | "9" ] { "0" | ... | "9" }
\end{verbatim}
     \index{syntax}\index{integer syntax}
     {\bf Note:} No blanks may appear within an integer !
     \index{IREAD}
\item[$IWRITE(a)$] $a \in \I$.
     The integer $a$ is written to the actual output stream.
     A conversion from $\beta$ representation to 
     decimal representation is done, so
     $\tma = L(a)^2$ and $\cma = L(a)^2$.
     The syntax is the same as for $IREAD$.
     \index{IWRITE}
\item[$ILWRIT(a)$] $a \in \L(\I)$.
     $a$ is a list of integers. Each element of $a$ is written
     to the output stream by $IWRITE$.
     So $\tma = t_{\rm IWRITE} \cdot {\rm length}(a)$
     and $\cma = c_{\rm IWRITE} \cdot {\rm length}(a)$.
     \index{ILWRIT}
\item[$l \gets IFACT(a)$]  $a \in \I$, $l \in \L(\I)$.
     $l$ is the list of the integer prime factors of $a$, 
     multiple factors occur multiple times in $l$.
     (e.g. $IFACT(8) = (2,2,2)$.)
     For large $L(a)$ the computing time has
     exponential growth: $\tma = c^{L(a)}$.
     \index{IFACT}
\item[$l \gets SMPRM()$] $l \in \L(\A)$.
     $l$ is the list of the prime numbers between $2$ and $1000$.
     \index{SMPRM}
\item[$l \gets DPGEN(m,k)$] $l \in \L(\A)$, 
     $m, k \in \A$, $m$ and $k$ positive, $k \leq 1000$, 
     $m + 2 \cdot k < \beta$.
     $l$ is the list of all prime numbers $p$, such that
     $m \leq p < m + 2 \cdot k$.    
     \index{DPGEN}
\end{deflist}

This concludes the summary of 
integer arithmetic functions.

Examples:
\begin{verbatim}
      a:=LIST(0,7,3).    --> (0 7 3)
      IWRITE(a).         --> 864691132213231616
      ISIGN(a).          --> 1
      INEG(a).           --> (0 -7 -3)
      ISUM(a,a).         --> (0 14 6)

      b:=IEXP(2,29).     --> (0 1) 
      IREM(a,b).         --> 0
      IQ(a,b).           --> (7 3) 
      IPROD(a,b).        --> (0 0 7 3)
      ICOMP(a,b).        --> 1

      IFACT(a).          --> (2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 
                              2 2 2 2 2 2 2 2 2 2 2 23 70026641)
\end{verbatim}
In the example
an integer $a$ is set to the list $(0,7,3)$.
That means $a$ represents the integer
$0 \cdot \beta^0 + 7 \cdot \beta^1 + 3 \cdot \beta^2$. 
The decimal representation of $a$ is generated in the 
second line.
Since ${\rm sign}(a) = {\rm sign}(7) = 1$ we have $a > 0$.
$-a$ has the representation  $(0,-7,-3)$.
$a+a$ has the representation $(0,14,6)$.

Next $b$ is set to $\beta = 2^{29}$, 
represented by the list $(0,1)$.
The remainder by division by $b$ is zero, $a$ rem $b = 0$,
($=$ the first element of the list representation of $a$).
The quotient $a/b$ is the rest list of the representation of $a$
i.e. $(7,3)$.
$a \cdot b$ shifts $a$ one $\beta$--digit to the left 
$(0,0,7,3)$.
Since $7 > 1$ and 
the list representation of $a$ is 
longer than that of $b$ we have 
${\rm sign}(a-b) = 1$, i.e. $a > b$.

Finally the prime factors of $a$ are computed,
since $a$ is divisible by $\beta$ we must have $29$ 
occurrences of $2$ in the prime factors.
And $7 + 3 \beta$ seems to have the prime factors
$23$ and $70026641$. 

For illustration we list 
the algorithms \verb/ICOMP/ and \verb/IPROD/ in Modula--2
in MAS.
The function of the algorithms should be clear from the 
step comments and the integer representation 
discussed before.

{\footnotesize
\begin{verbatim}
PROCEDURE ICOMP(A,B: LIST): LIST;
(*Integer comparison.  A and B are integers.  s=SIGN(A-B).*)
VAR   AL, AP, BL, BP, DL, SL, UL, VL: LIST;
BEGIN
(*1*) (*A and B single-precision.*)
      IF (A < BETA) AND (B < BETA) THEN SL:=MASSIGN(A-B);
         RETURN(SL); END;
(*2*) (*A single-precision.*)
      IF A < BETA THEN SL:=-ISIGNF(B); RETURN(SL); END;
(*3*) (*B single-precision.*)
      IF B < BETA THEN SL:=ISIGNF(A); RETURN(SL); END;
(*4*) (*compare corresponding digits.*) SL:=0; AP:=A; BP:=B;
      REPEAT ADV(AP,AL,AP); ADV(BP,BL,BP); UL:=MASSIGN(AL); 
             VL:=MASSIGN(BL);
             IF UL*VL = -1 THEN SL:=UL; RETURN(SL); END;
             DL:=AL-BL;
             IF DL <> 0 THEN SL:=MASSIGN(DL); END;
             UNTIL (AP = SIL) OR (BP = SIL);
(*5*) (*same length*)
      IF (AP = SIL) AND (BP = SIL) THEN RETURN(SL); END;
(*6*) (*use sign of longer input.*)
      IF AP = SIL THEN SL:=-ISIGNF(BP); ELSE
         SL:=ISIGNF(AP); END;
      RETURN(SL);
(*9*) END ICOMP;
\end{verbatim}
}
\index{ICOMP}
\verb/BETA/ denotes $\beta$,
\verb/MASSIGN/ denotes the sign function for $\beta$--integers.
\verb/ISIGNF/ denotes the sign function for integers in
the libraries (for compatibility with the ALDES / SAC--2 libraries).
\verb/SIL/ denotes the empty list.

{\footnotesize
\begin{verbatim}
PROCEDURE IPROD(A,B: LIST): LIST;
(*Integer product.  A and B are integers.  C=A*B.*)
VAR   AL, AP, APP, BL, BP, C, C2, CL, CLP, CP, CPP, EL, FL, 
      I, ML, NL, TL: LIST;
BEGIN
(*1*) (*A or B zero.*)
      IF (A = 0) OR (B = 0) THEN C:=0; RETURN(C); END;
(*2*) (*A and B single-precision.*)
      IF (A < BETA) AND (B < BETA) THEN DPR(A,B,CLP,CL);
         IF CLP = 0 THEN C:=CL; ELSE C:=LIST2(CL,CLP); END;
         RETURN(C); END;
(*3*) (*A or B single-precision.*)
      IF A < BETA THEN C:=IDPR(B,A); RETURN(C); END;
      IF B < BETA THEN C:=IDPR(A,B); RETURN(C); END;
(*4*) (*interchange if B is longer.*) ML:=LENGTH(A); NL:=LENGTH(B);
      IF ML >= NL THEN AP:=A; BP:=B; ELSE AP:=B; BP:=A; END;
(*5*) (*set product to zero.*) C2:=LIST2(0,0); C:=C2;
      FOR I:=1 TO ML+NL-2 DO C:=COMP(0,C); END;
(*6*) (*multiply digits and add products.*) CP:=C;
      REPEAT APP:=AP; ADV(BP,BL,BP);
             IF BL <> 0 THEN CPP:=CP; CLP:=0;
                REPEAT ADV(APP,AL,APP); DPR(AL,BL,EL,FL);
                       CL:=FIRST(CPP); CL:=CL+FL; CL:=CL+CLP;
                       CLP:=CL DIV BETA; TL:=CLP*BETA; CL:=CL-TL;
                       SFIRST(CPP,CL); CLP:=EL+CLP; CPP:=RED(CPP);
                       UNTIL APP = SIL;
                SFIRST(CPP,CLP); END;
             CP:=RED(CP);
             UNTIL BP = SIL;
(*7*) (*leading digit zero*)
      IF CLP = 0 THEN SRED(C2,SIL); END;
      RETURN(C);
(*9*) END IPROD;
\end{verbatim}
}
\index{IPROD}
\verb/DPR/ denotes the $\beta$--digit product,
\verb/IDPR/ denotes the integer product with a $\beta$--digit,
\verb/SFIRST/ (Set FIRST) 
and \verb/SRED/ (Set RED)
are list modifying
functions not available for interactive use in MAS.


\subsection{Exercises}

\begin{enumerate}
\item Compute the gcd of $156562431911123$ and 
      $442677773754356$ in three ways:
      \begin{enumerate}
      \item Use the the built--in algorithm for the gcd.
      \item Write an recursive algorithm in MAS and
            use it for the computation.
      \item Write an iterative algorithm in MAS and
            use it for the computation.
      \item $*$ Write an iterative algorithm for the 
            extended gcd.
      \end{enumerate}
      \index{greatest common divisor}
\item Write algorithms to compute the determinant of 
      a square matrix over the integers. 
      Assume that the matrices are represented as lists of lists
      of integers. Use the Laplace method for the expansion 
      of the determinant.
      Therefore write the following sub--algorithms: 
      \begin{enumerate}
      \item Deletion of an element of a vector (represented as list).
      \item Deletion of a column of a matrix (represented 
            as list of vectors).
      \item Determinant expansion using Laplace's method.
      \item An algorithm to generate matrices with
            respect to a function of the matrix entries.
      \item $*$ An algorithm with Gauss's method for the 
            computation of the determinant.
      \item $*$ Write the determinant algorithm for other
            coefficient domains, e.g. rational numbers.
      \end{enumerate}
      \index{determinant}
\end{enumerate}

The solutions to the exercises are discussed in the sequel.

{\bf Exercise 1.} 
The built--in gcd algorithm has the name \verb/IGCD/.

For the MAS algorithms we use the euclidean method:
\begin{displaymath}
      {\rm gcd}(a,b) = \left\{
                       \begin{array}{ll}
                              a & {\rm if} \ \  b = 0, \\
  {\rm gcd}(b, {\rm rem}(a,b) ) & {\rm else}.
                       \end{array}
                       \right.
\end{displaymath}

With this definition the
recursive gcd algorithm can be formulated as follows:
\begin{verbatim}
       PROCEDURE ggt(a,b);
       IF b = 0 THEN RETURN(a) 
                ELSE RETURN(ggt(b,IREM(a,b))) 
          END ggt.
\end{verbatim}
The integer remainder is called \verb/IREM/.

The iterative gcd algorithm can be formulated 
using the WHILE--statement and an additional variable 
\verb/d/:
\begin{verbatim}
       PROCEDURE GGT(a,b);
       VAR   d: ANY;
       BEGIN 
       WHILE b # 0 DO d:=b; 
             b:=IREM(a,b); a:=d END;
       RETURN(a) END GGT.
\end{verbatim}

The numbers can be converted to internal representation by \verb/IREAD/.
The sample output with computing times follows:
\begin{verbatim}
       a:=IREAD(). 156562431911123 
       b:=IREAD(). 442677773754356 

       c:=IGCD(a,b).
       {0 sec} ANS: 7

       c:=ggt(a,b).
       {4 sec} ANS: 7

       c:=GGT(a,b).
       {2 sec} ANS: 7
\end{verbatim}
The gcd is always $7$. The built--in algorithm is 
the fastest: it needs unmeasurable many seconds.
The recursive algorithm is the slowest: it needs 
about $4$ seconds and the iterative algorithm 
need $2$ seconds. 


{\bf Exercise 2.} 
Let $A = ( a_{ij} )$ be a $n \times n$ matrix over 
a commutative ring without zero divisors.
Then the determinant has the following expansion:
\begin{displaymath}  
       {\rm det}(A) = \sum_{i=1}^{n} (-1)^{i+j} \cdot a_{ij} \cdot
                                     {\rm det}(A_{ij})
       \ \ \ \mbox{\rm for some} \ \  j, (1 \leq j \leq n).
\end{displaymath}
Where $A_{ij}$ is the matrix $A$ without the $i$--th row
and the $j$--th column. 

The computing time can be estimated by the following 
considerations:
The determinant expansion produces $n !$ summands of products
of $n$ factors. 
Assume $L(A) = {\rm max}\{ L(a_{ij}), 1 \leq i, j \leq n \}$.
So $\tma = n ! \cdot L(A)^n$ and 
$\tmi = L(A)^n$, e.g. if the matrix is in triangular form.

In contrast it needs about $n^3$ field operations 
and $n^5$ ring operations 
to transform a matrix to upper triangular from by the 
Gauss respective the Bareiss algorithm.
I.e. the Laplace expansion is very slow.

The determinant expansion algorithm can be formulated as follows 
\begin{verbatim}
   PROCEDURE det(M);
   (*Determinant. M is a matrix (i.e a list 
   of row lists). The determinant of M is computed. *)
   VAR   i, d, dp, s, N, MP, V, VP, v: LIST;
   BEGIN 
   (*1*) d:=0;
         IF M = NIL THEN RETURN(d) END;
   (*2*) ADV(M,V,MP); 
         IF MP = NIL THEN RETURN(FIRST(V)) END;
   (*3*) s:=1; i:=0;
         WHILE V # NIL DO ADV(V,v,V); i:=i+1;
               IF v # 0 THEN 
                  N:=delcolumn(MP,i); dp:=det(N); 
                  dp:=IPROD(v,dp);
                  IF s < 0 THEN dp:=INEG(dp) END; 
                  d:=ISUM(d,dp);
                  END;
               s:=-s; END;
   (*4*) RETURN(d);
   (*9*) END det.
\end{verbatim}
\index{det}
Step 2 handles the recursion base of a $1 \times 1$ matrix.
In step 3 the determinant is expanded w.r.t. the first 
row (that is $j=1$ in our formula).
The test if some $a_{ij} = 0$ 
(denoted by \verb/v/) is important, since in 
this case the expansion of the sub--matrix is useless 
and much computing time can be saved.
Since the first row of the matrix has been removed in step 2,
it is only necessary to remove the $i$--th column of the
matrix. The matrices $A_{i1}$ are denoted by \verb/N/.
The factor $(-1)^{i+1}$ is handled by a sign flag 
denoted by \verb/s/.

The vector element deletion algorithm is
\begin{verbatim}
   PROCEDURE delelem(V,i);
   (*Delete element in list. V is a list. The i-th 
   element of V is deleted. 0 <= i <= length(V). *)
   VAR   U, VP, v, j: ANY;
   BEGIN 
   (*1*) IF i <= 0  THEN RETURN(V) END;
         IF V = NIL THEN RETURN(V) END;
   (*2*) VP:=V; j:=0; U:=NIL;
         REPEAT j:=j+1; 
                IF VP = NIL THEN RETURN(V) END; 
                ADV(VP,v,VP); U:=COMP(v,U); 
                UNTIL j = i;
   (*3*) U:=RED(U); U:=INV(U); U:=CONC(U,VP);
         RETURN(U);
   (*9*) END delelem.
\end{verbatim}
Step 1 does some precondition checks.
In step 2 the beginning of the list is copied to avoid 
the modification of existing data. 
In step 3 the (new) beginning is reversed and 
concatenated with the rest of the vector.

The column delete algorithm is straight forward: 
for each row the previous algorithm is called.
\begin{verbatim}
   PROCEDURE delcolumn(M,i);
   (*Delete column in matrix. M is a matrix (i.e a list 
   of row lists). In each row the i-th element is 
   deleted. *)
   VAR   N, MP, V: ANY;
   BEGIN 
   (*1*) IF i <= 0  THEN RETURN(M) END;
         IF M = NIL THEN RETURN(M) END;
   (*2*) MP:=M; N:=NIL;
         WHILE MP # NIL DO ADV(MP,V,MP); 
               V:=delelem(V,i); N:=COMP(V,N); 
               END;
   (*3*) N:=INV(N); RETURN(N);
   (*9*) END delcolumn.
\end{verbatim}

The sample matrices are:
\begin{enumerate}
\item the unit matrix $E = ( e_{ij} )$ with:
      \begin{math}
             e_{ij} = \left\{
                       \begin{array}{ll}
                              1 & {\rm if} \ \  i = j, \\
                              0 & {\rm else},
                       \end{array}
                       \right.
      \end{math}
\item an upper triangular matrix $U = ( u_{ij} )$ with:
      \begin{math}
             u_{ij} = \left\{
                       \begin{array}{ll}
                              j & {\rm if} \ \ i \leq j, \\
                              0 & {\rm else},
                       \end{array}
                       \right.
      \end{math}
\item a lower triangular matrix $L = ( l_{ij} )$ with:
      \begin{math}
             l_{ij} = \left\{
                       \begin{array}{ll}
                              j & {\rm if} \ \ i \geq j, \\
                              0 & {\rm else}.
                       \end{array}
                       \right.
      \end{math}
\end{enumerate}

The generating functions for the matrices are as follows:
\begin{verbatim}
   PROCEDURE eh(i,j);
   IF i = j THEN RETURN(1) ELSE RETURN(0) END eh.

   PROCEDURE ut(i,j);
   IF i <= j THEN RETURN(j) ELSE RETURN(0) END ut.

   PROCEDURE lt(i,j);
   IF i >= j THEN RETURN(j) ELSE RETURN(0) END lt.
\end{verbatim}
\verb/eh/ means `Einheitsmatrix' = unit matrix,
\verb/ut/ means upper triangular matrix and
\verb/lt/ means lower triangular matrix.

The matrix generation procedure takes the 
generating function \verb/f/ and the dimension \verb/k/
as input and delivers the appropriate matrix. 
\begin{verbatim}
   PROCEDURE mat(f,k);
   VAR   i, j, V, M: LIST;
   BEGIN
   (*1*) M:=NIL;
         IF k <= 0 THEN RETURN(M) END;
   (*2*) i:=0;
         WHILE i < k DO i:=i+1; j:=0; V:=NIL; 
               WHILE j < k DO j:=j+1;
                     V:=COMP(f(i,j),V) END;
               V:=INV(V); M:=COMP(V,M);
               END;
   (*4*) M:=INV(M); RETURN(M);
   (*9*) END mat.
\end{verbatim}

A sample output follows:
\begin{verbatim}
   A:=mat(ut,5).
   {4 sec} ANS: ((1 2 3 4 5) (0 2 3 4 5) (0 0 3 4 5) 
                (0 0 0 4 5) (0 0 0 0 5))
 
   d:=det(A). {246 sec} ANS: 120

   A:=mat(lt,10).
   {14 sec} ANS: ((1 0 0 0 0 0 0 0 0 0) (1 2 0 0 0 0 
   0 0 0 0) (1 2 3 0 0 0 0 0 0 0) (1 2 3 4 0 0 0 0 0 0) 
   (1 2 3 4 5 0 0 0 0 0) (1 2 3 4 5 6 0 0 0 0) (1 2 3 
   4 5 6 7 0 0 0) (1 2 3 4 5 6 7 8 0 0) (1 2 3 4 5 6 7 
   8 9 0) (1 2 3 4 5 6 7 8 9 10))
 
   d:=det(A). {64 sec} ANS: 3628800
\end{verbatim}
In the first example a $5 \times 5$ upper triangular
matrix is generated. The determinant is $120$
and the computation takes $246$ seconds. 

In the second example a $10 \times 10$ lower triangular
matrix is generated. The determinant is $3628800$ and
the computing time is $64$ seconds. 

The second computation is much faster, since the 
determinant is expanded with respect to the 
first row. And the lower triangular matrix has 
most times zero in these rows.

