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