PSLQ.frink

Download or view PSLQ.frink in plain text format


/** 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. */

}


Download or view PSLQ.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 19944 days, 10 hours, 30 minutes ago.