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 20217 days, 23 hours, 34 minutes ago.