## hilbert-matrix.a68

```
1  COMMENT
2
3  @section Synopsis
4
5  Compute the determinant of a Hilbert matrix using fractions.
6
7  A comprehensive implementation of fractions is in the MC Algol 68 test set (appl08).
8
9  Order   Result
10   1      1 / 1
11   2      1 / 12
12   3      1 / 2,160
13   4      1 / 6,048,000
14   5      1 / 266,716,800,000
15   6      1 / 186,313,420,339,200,000
16   7      1 / 2,067,909,047,925,770,649,600,000
17   8      1 / 365,356,847,125,734,485,878,112,256,000,000
18   9      1 / 1,028,781,784,378,569,697,887,052,962,909,388,800,000,000
19  10      1 / 46,206,893,947,914,691,316,295,628,839,036,278,726,983,680,000,000,000
20
21  COMMENT
22
23  BEGIN
24
25  # Data structure. #
26
27        MODE FRAC = STRUCT (LINT n, d), LINT = LONG LONG INT;
28
29        OP NOM = (FRAC u) LINT: n OF u,
30           DEN = (FRAC u) LINT: d OF u;
31
32        PR precision 101 PR # Now LINT holds a googol. #
33
34  # Basic operations. #
35
36        OP RECIPROCAL = (LINT i) FRAC: IF i >= 0 THEN (1, i) ELSE (-1, -i) FI;
37
38        OP - = (FRAC u) FRAC: (-NOM u, DEN u);
39
40        OP + = (FRAC u, FRAC v) FRAC:
41           BEGIN LINT k = DEN u GCD DEN v;
42                 LINT du = DEN u OVER k, dv = DEN v OVER k;
43                 LINT n = NOM u * dv + NOM v * du;
44                 LINT l = n GCD k, d = dv * du;
45                 (n OVER l, k OVER l * d)
46           END;
47
48        OP +:= = (REF FRAC u, FRAC v) REF FRAC: u := u + v;
49
50        OP - = (FRAC u, FRAC v) FRAC:
51           BEGIN LINT k = DEN u GCD DEN v;
52                 LINT du = DEN u OVER k, dv = DEN v OVER k;
53                 LINT n = NOM u * dv - NOM v * du;
54                 LINT l = n GCD k, d = dv * du;
55                 (n OVER l, k OVER l * d)
56           END;
57
58        OP -:= = (REF FRAC u, FRAC v) REF FRAC: u := u - v;
59
60        OP * = (FRAC u, v) FRAC:
61           BEGIN LINT i = NOM u GCD DEN v, j = NOM v GCD DEN u;
62                 ((NOM u OVER i) * (NOM v OVER j), (DEN u OVER j) * (DEN v OVER i))
63           END;
64
65        OP *:= = (REF FRAC u, FRAC v) REF FRAC: u := u * v;
66
67        OP / = (FRAC u, FRAC v) FRAC:
68           IF LINT i = NOM u GCD NOM v, j = DEN v GCD DEN u;
69              NOM v >= 0
70           THEN ((NOM u OVER i) * (DEN v OVER j), (DEN u OVER j) * (NOM v OVER i))
71           ELSE (- (NOM u OVER i) * (DEN v OVER j), - (DEN u OVER j) * (NOM v OVER i))
72           FI;
73
74        OP /:= = (REF FRAC u, FRAC v) REF FRAC: u := u / v;
75
76  # Comparing rationals with integrals. #
77
78        OP = = (FRAC u, LINT i) BOOL: NOM u = i ANDF DEN u = 1;
79
80        OP /= = (FRAC u, LINT i) BOOL: NOT (u = i);
81
82  # Matrix algebra. #
83
84        OP INNER = ([] FRAC u, v) FRAC:
85           # Innerproduct of two arrays of rationals #
86           BEGIN FRAC s := (0, 1);
87                 FOR i TO UPB u
88                 DO s +:= u[i] * v[i]
89                 OD;
90                 s
91           END;
92
93        PRIO INNER = 8;
94
95        PROC lu decomposition = (REF [, ] FRAC a, REF [] INT p) VOID:
96             # LU decomposition by Crout's method, of a matrix of rationals. #
97             BEGIN INT n = 1 UPB a;
98                   FOR k TO n
99                   DO FRAC piv := (0, 1), INT k1 := k - 1;
100                      REF INT pk = p[k];
101                      REF [] FRAC aik = a[, k], aki = a[k,];
102                      FOR i FROM k TO n
103                      DO aik[i] -:= a[i, 1 : k1] INNER aik[1 : k1];
104                         IF piv = LINT (0) AND aik[i] /= LINT (0)
105                         THEN piv := aik[i];
106                              pk := i
107                         FI
108                      OD;
109                      IF piv = LINT (0)
110                      THEN print((newline, newline, "Singular matrix"));
111                           stop
112                      FI;
113                      IF pk /= k
114                      THEN FOR i TO n
115                           DO FRAC r = aki[i];
116                              aki[i] := a[pk, i];
117                              a[pk, i] := -r
118                           OD
119                      FI;
120                      FOR i FROM k + 1 TO n
121                      DO aki[i] -:= aki[1 : k1] INNER a[1 : k1, i] /:= piv
122                      OD
123                   OD
124             END;
125
126        PROC determinant = ([,] FRAC a) FRAC:
127             # Determinant of a decomposed matrix is its trace. #
128             BEGIN FRAC d := (1, 1);
129                   FOR i TO 1 UPB a
130                   DO d *:= a[i, i]
131                   OD;
132                   d
133             END;
134
135  # Recursive definition of greatest common divisor. #
136
137        OP GCD = (LINT a, b) LINT:
138           IF b = 0
139           THEN ABS a
140           ELSE b GCD (a MOD b)
141           FI;
142
143        PRIO GCD = 8;
144
145  # Table of required results. #
146
147        [] LINT table = BEGIN
148           LONG LONG 1,
149           LONG LONG 12,
150           LONG LONG 2160,
151           LONG LONG 6048000,
152           LONG LONG 266716800000,
153           LONG LONG 186313420339200000,
154           LONG LONG 2067909047925770649600000,
155           LONG LONG 365356847125734485878112256000000,
156           LONG LONG 1028781784378569697887052962909388800000000,
157           LONG LONG 46206893947914691316295628839036278726983680000000000
158        END;
159
160  # Compute determinant of Hilbert matrix of increasing rank. #
161
162        printf((\$"Determinant of the Hilbert matrix - LONG LONG INT"\$));
163        FOR n TO UPB table
164        DO [1 : n, 1 : n] FRAC a;
165           FOR i TO n
166           DO a[i,i] := RECIPROCAL LINT (i * 2 - 1);
167              FOR j FROM i + 1 TO n
168              DO a[i, j] := a[j, i] := RECIPROCAL LINT (i + j - 1)
169              OD
170           OD;
171           lu decomposition(a, LOC [1 : n] INT);
172
173           OP PRETTY = (LINT n) VOID:
174              IF n < 0
175              THEN print ("-");
176                   PRETTY -n
177              ELIF n < 1000
178              THEN print (whole (n, 0))
179              ELSE PRETTY (n % 1000);
180                   printf ((\$","3d\$, n %* 1000))
181              FI;
182
183           FRAC det = determinant (a);
184           newline (standout);
185           PRETTY NOM det;
186           print (" / ");
187           PRETTY DEN det;
188           newline (standout);
189           PRETTY LONG LONG 1;
190           print (" / ");
191           PRETTY table[n];
192           newline (standout)
193       OD;
194       new line (standout)
195  END

```

 © 2002-2024 J.M. van der Veer (jmvdveer@xs4all.nl)