share/pslq/pslq_maxima.mac

   1 /* Pure Maxima implementation of PSLQ algorithm,
   2  * from: https://www.davidhbailey.com/dhbpapers/pslq-cse.pdf
   3  * Section 2, "The PSLQ Algorithm".
   4  *
   5  * <quote>
   6
   7    Let x be the n-long input real vector, and let nint denote the nearest integer function.
   8    Select γ ≥ √(4/3).
   9    Then perform the following operations:
  10
  11    Initialize:
  12
  13    1. Set the n × n matrices A and B to the identity.
  14
  15    2. Compute the n-long vector s as s[k] := sqrt(sum(x[j]^2, j, k, n)),
  16    and set y to the x vector, normalized by s[1].
  17
  18    3. Compute the initial n × (n − 1) matrix H as
  19    H[i,j] = 0 if i < j,
  20    H[j,j] := s[j + 1]/s[j], and
  21    H[i,j] := −y[i] y[j] /(s[j] s[j+1]) if i > j.
  22
  23    4. Reduce H:
  24    For i := 2 to n:
  25      for j := i − 1 to 1 step −1:
  26        set t := nint(H[i,j] / H[j,j]); and y[j] := y[j] + t y[i];
  27        for k := 1 to j:
  28          set H[i,k] := H[i,k] − t H[j,k];
  29        endfor;
  30        for k := 1 to n:
  31          set A[i,k] := A[i,k] − t A[j,k] and B[k,j] := B[k,j] + t B[k,i];
  32        endfor;
  33      endfor;
  34    endfor.
  35
  36    Iterate until an entry of y is within a reasonable tolerance of zero,
  37    or precision has been exhausted:
  38
  39    1. Select m such that γ^i |H[i,i]| is maximal when i = m.
  40
  41    2. Exchange the entries of y indexed m and m + 1,
  42       the corresponding rows of A and H,
  43       and the corresponding columns of B.
  44
  45    3. Remove the corner on H diagonal:
  46       If m ≤ n − 2
  47       then set t0 := sqrt(H[m,m]^2 + H[m,m + 1]^2),
  48                t1 := H[m,m]/t0
  49            and t2 := H[m,m + 1]/t0;
  50            for i := m to n:
  51              set t3 := H[i,m],
  52                  t4 := H[i,m + 1],
  53                  H[i,m] := t1 t3 + t2 t4
  54              and H[i,m + 1] := −t2 t3 + t1 t4;
  55            endfor;
  56       endif.
  57
  58    4. Reduce H:
  59       For i := m + 1 to n:
  60         for j := min(i − 1, m + 1) to 1 step −1:
  61           set t := nint(H[i,j] /H[j,j])
  62           and y[j] := y[j] + t y[i];
  63           for k := 1 to j:
  64             set H[i,k] := H[i,k] − t H[j, k];
  65           endfor;
  66           for k := 1 to n:
  67             set A[i,k] := A[i,k] − t A[j,k] and B[k,j] := B[k,j] + t B[k,i];
  68           endfor;
  69         endfor;
  70       endfor.
  71
  72    5. Norm bound: Compute M := 1/ max[j] |H[j,j]|.
  73       Then there can exist no relation vector whose Euclidean norm is less than M.
  74
  75    Upon completion, the desired relation is found in the column of B corresponding to the
  76    zero entry of y.
  77
  78  * </quote>
  79  *
  80  * NOTE 1: This implementation only implements the steps outlined in the paper;
  81  * there is no outer loop to control looping until a criterion is met.
  82  * The initialization step is pslq_init (no arguments, x is a global variable)
  83  * and pslq_iterate (no arguments, intermediate results are global variables).
  84  *
  85  * NOTE 2: I have translated the stated algorithm carefully, but it should
  86  * probably be reviewed for accuracy.
  87  */
  88
  89 /* Initialize: */
  90
  91 pslq_init () :=
  92    (kill (A, B, s, n, y, H),
  93     n: length (x),
  94
  95     /* 1. Set the n × n matrices A and B to the identity. */
  96     A: ident (n),
  97     B: ident (n),
  98
  99     /* 2. Compute the n-long vector s as s[k] := sqrt(sum(x[j]^2, j, k, n)), */
 100     s: makelist (sqrt (sum (x[j]^2, j, k, n)), k, 1, n),
 101
 102     /* and set y to the x vector, normalized by s[1]. */
 103     y: x / s[1],
 104
 105     /* 3. Compute the initial n × (n − 1) matrix H as
 106      *    H[i,j] = 0 if i < j,
 107      *    H[j,j] := s[j + 1]/s[j], and
 108      *    H[i,j] := −y[i] y[j] /(s[j] s[j+1]) if i > j.
 109      */
 110
 111     H: zeromatrix (n, n - 1),
 112     for j thru n - 1 do H[j, j]: s[j + 1]/s[j],
 113     for i: 2 thru n
 114         do for j thru i - 1
 115             do H[i, j]: - y[i] * y[j] / (s[j] * s[j + 1]),
 116
 117     /*  4. Reduce H:
 118      *     For i := 2 to n:
 119      *       for j := i − 1 to 1 step −1:
 120      *         set t := nint(H[i,j] / H[j,j]); and y[j] := y[j] + t y[i];
 121      *         for k := 1 to j:
 122      *           set H[i,k] := H[i,k] − t H[j,k];
 123      *         endfor;
 124      *         for k := 1 to n:
 125      *           set A[i,k] := A[i,k] − t A[j,k] and B[k,j] := B[k,j] + t B[k,i];
 126      *         endfor;
 127      *       endfor;
 128      *     endfor.
 129      */
 130
 131     (print ("pslq_init: just before step 4 (reduce H)."),
 132      display (A, B, s, y, H)),
 133
 134     for i: 2 thru n
 135         do for j: i - 1 thru 1 step -1
 136             do block ([t: nint (H[i, j] / H[j, j])],
 137                       y[j]: y[j] + t * y[i],
 138                       for k thru j
 139                           do H[i, k]: H[i, k] - t * H[j, k],
 140                       for k thru n
 141                           do (A[i, k]: A[i, k] - t * A[j, k],
 142                               B[k, j]: B[k, j] + t * B[k, i])));
 143
 144 nint (x) := round (x);
 145
 146 /* Iterate until an entry of y is within a reasonable tolerance of zero,
 147  * or precision has been exhausted:
 148  */
 149
 150 gamma: sqrt (4/3);
 151
 152 pslq_iterate () :=
 153    (/* 1. Select m such that γ^i |H[i,i]| is maximal when i = m. */
 154
 155     max_foo: gamma * abs (H[1, 1]),
 156     max_foo_i: 1,
 157     for i:2 thru n - 1
 158         do block ([foo: gamma^i * abs (H[i, i])],
 159                   if foo > max_foo
 160                       then (max_foo: foo, max_foo_i: i)),
 161     m: max_foo_i,
 162
 163     /* 2. Exchange the entries of y indexed m and m + 1,
 164      *    the corresponding rows of A and H,
 165      *    and the corresponding columns of B.
 166      */
 167
 168     exchange_elements (y, m, m + 1),
 169     exchange_rows (A, m, m + 1),
 170     exchange_rows (H, m, m + 1),
 171     exchange_columns (B, m, m + 1),
 172
 173     /*  3. Remove the corner on H diagonal:
 174      *     If m ≤ n − 2
 175      *     then set t0 := sqrt(H[m,m]^2 + H[m,m + 1]^2),
 176      *              t1 := H[m,m]/t0
 177      *          and t2 := H[m,m + 1]/t0;
 178      *          for i := m to n:
 179      *            set t3 := H[i,m],
 180      *                t4 := H[i,m + 1],
 181      *                H[i,m] := t1 t3 + t2 t4
 182      *            and H[i,m + 1] := −t2 t3 + t1 t4;
 183      *          endfor;
 184      *     endif.
 185      */
 186
 187     if m <= n - 2
 188         then block ([t0, t1, t2, t3, t4],
 189                     t0: sqrt (H[m, m]^2 + H[m, m + 1]^2),
 190                     t1: H[m, m] / t0,
 191                     t2: H[m, m + 1] / t0,
 192                     for i: m thru n
 193                         do (t3: H[i, m],
 194                             t4: H[i, m + 1],
 195                             H[i, m]: t1 * t3 + t2 * t4,
 196                             H[i, m + 1]: -t2 * t3 + t1 * t4)),
 197
 198     /* 4. Reduce H:
 199      *    For i := m + 1 to n:
 200      *      for j := min(i − 1, m + 1) to 1 step −1:
 201      *        set t := nint(H[i,j] /H[j,j])
 202      *        and y[j] := y[j] + t y[i];
 203      *        for k := 1 to j:
 204      *          set H[i,k] := H[i,k] − t H[j, k];
 205      *        endfor;
 206      *        for k := 1 to n:
 207      *          set A[i,k] := A[i,k] − t A[j,k] and B[k,j] := B[k,j] + t B[k,i];
 208      *        endfor;
 209      *      endfor;
 210      *    endfor.
 211      */
 212
 213     for i: m + 1 thru n
 214         do for j: min(i - 1, m + 1) thru 1 step -1
 215                do block ([t: nint (H[i, j] / H[j, j])],
 216                          y[j]: y[j] + t * y[i],
 217                          for k thru j
 218                              do H[i, k]: H[i, k] - t * H[j, k],
 219                          for k thru n
 220                              do (A[i, k]: A[i, k] - t * A[j, k],
 221                                  B[k, j]: B[k, j] + t * B[k, i])),
 222
 223     /* 5. Norm bound: Compute M := 1/ max[j] |H[j,j]|.
 224      *    Then there can exist no relation vector whose Euclidean norm is less than M.
 225      */
 226
 227     M: 1 / lmax (makelist (abs (H[j, j]), j, 1, n - 1)));
 228
 229 /* Upon completion, the desired relation is found in the column of B corresponding to the
 230  * zero entry of y.
 231  */
 232
 233 exchange_elements (l, i, j) := ([l[j], l[i]]: [l[i], l[j]], l);
 234 exchange_rows (M, i, j) := ([M[j], M[i]]: [M[i], M[j]], M);
 235 exchange_columns (M, i, j) := transpose (exchange_rows (transpose (M), i, j));
 236