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
