/** This is an implementation of the PSLQ algorithm as defined in

    Bailey, D. H., & Broadhurst, D. J. (2000).
    Parallel integer relation detection: Techniques and applications.
    Mathematics of Computation, 70(236), 1719–1737.
    doi:10.1090/s0025-5718-00-01278-3 

    https://sci-hub.se/10.1090/S0025-5718-00-01278-3 

    Also see:
    Bailey, D. H. (2000). Integer relation detection.
     Computing in Science & Engineering, 2(1), 24–28. doi:10.1109/5992.814653 

    https://sci-hub.se/10.1109/5992.814653 
*/

use Matrix.frink

/** This performs the PSLQ algorithm on the real-valued input vector x.  x
    is assumed to have its values in elements 1 to n to match matrix notation.
*/
PSLQ[x, prec=getPrecision[]] :=
{
   setPrecision[prec]
   bestY = googol
   
   // Initialize
   // Step 1
   n = length[x]-1

   A = Matrix.makeIdentity[n]
   B = Matrix.makeIdentity[n]

   // Step 2
   s = new array[[n+1], undef]
   for k = 1 to n
   {
      sum = 0
      for j = k to n
         sum = sum + (x@j)^2
      s@k = sqrt[sum, prec+5]
   }

   t = 1/(s@1)

   y = new array[[n+1], undef]
   for k = 1 to n
   {
      y@k = t x@k
      s@k = t s@k
   }

   H = new Matrix[n, n]
   // Step 3:  Initial H
   for j = 1 to n-1
   {
      for i = 1 to j-1
      {
         H.set[i,j,0]
      }
      H.set[j,j, s@(j+1) / (s@j)]

      for i = j+1 to n
      {
         H.set[i,j, -y@i y@j / (s@j s@j+1)]
      }
   }

   // Step 4: Reduce H
   for i = 2 to n
   {
      for j = i-1 to 1 step -1
      {
         t = round[H.get[i,j] / H.get[j,j]]
         y@j = y@j + t y@i
         
         for k = 1 to j
         {
            H.set[i, k, H.get[i,k] - t H.get[j, k]]
         }

         for k = 1 to n
         {
            A.set[i,k, A.get[i,k] - t A.get[j, k]]
            B.set[k,j, B.get[k,j] + t B.get[k, i]]
         }
      }
   }

   println["A is \n" + A.formatMatrix[]]
   println["B is \n" + B.formatMatrix[]]
   println["H is \n" + H.formatMatrix[]]

   // The algorithm says "select gamma > sqrt[4/3]" but with no more guidance
   // than that.   That means something greater than 1.1547005383792515914.
   gamma = 12ee-1

   // Iteration.

   do
   {
      // Step 1.  Select m such that gamma^i |Hii| is maximal when i = m
      maxm = 0
      maxval = 0
      for i = 1 to n
      {
         val = gamma^i abs[H.get[i,i]]
         if val > maxval
         {
            maxm = i
            maxval = val
         }
      }

      m = maxm
      // println["Maximum m is $m"]

      // Step 2.  Exchange the entries of y indexed m and m+1, the corresponding
      // rows of A and H, and the corresponding columns of B.

      // Swap y@m and y@(m+1)
      temp = y@(m+1)
      y@(m+1) = y@m
      y@m = temp

      A.exchangeRows[m, m+1]
      H.exchangeRows[m, m+1]
      B.exchangeCols[m, m+1]

      // Step 3.
      // Remove corner on H diagonal
      if m <= n-2
      {
         t0 = sqrt[H.get[m,m]^2 + H.get[m, m+1]^2, prec+5]
         t1 = H.get[m, m] / t0
         t2 = H.get[m, m+1] / t0

         for i = m to n
         {
            t3 = H.get[i, m]
            t4 = H.get[i, m+1]
            H.set[i, m,    t1 t3 + t2 t4]
            H.set[i, m+1, -t2 t3 + t1 t4]
         }
      }

      // Step 4.
      // Reduce H
      for i = m + 1 to n
      {
         for j = min[i-1, m+1] to 1 step -1
         {
            t = round[H.get[i,j] / H.get[j,j]]
            y@j = y@j + t y@i

            for k = 1 to j
            {
               H.set[i,k, H.get[i,k] - t H.get[j, k]]
            }

            for k = 1 to n
            {
               A.set[i, k, A.get[i,k] - t A.get[j, k]]
               B.set[k, j, B.get[k,j] + t B.get[k, i]]
            }
         }
      }

      // Step 5.
      // Norm Bound:  Compute M := 1/max_j |Hjj|
      // The notation max_j is unclear but maybe the next page clarifies it
      
      maxval = 0
      for j = 1 to n-1  // Should this be n?
      {
         val = abs[H.get[j,j]]
         if val > maxval
            maxval = val
      }

      M = 1 / maxval

      println[]
      println["M is $M"]
      println["y is " + rest[y]]
      smallestY = min[abs[rest[y]]]
      col = 0
      for i = 1 to length[y]-1
         if abs[y@i] == smallestY
            col = i

       println["Smallest y is $smallestY"]
       println["Relation is found in column $col of B"]
       if smallestY < bestY
       {
          bestY = smallestY
          println["B is \n" + B.formatMatrix[]]
       }
   } while smallestY > 10^-(prec-5)

   println["A is \n" + A.formatMatrix[]]
   println["B is \n" + B.formatMatrix[]]
   println["H is \n" + H.formatMatrix[]]

   /** The integer relation is now found in the column of B
       corresponding to the smallest entry of y. */
}