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, eliasen@mindspring.com