# PSLQ.frink

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