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