PSLQ

[Han]

PSLQ (Partial Sums LQ Decomposition) can be used to determine if there exists a set of integer coefficients \( \mathbf{a} = (a_1, a_2, \dotsm, a_n) \) for a given vector \( \mathbf{x} = (x_1, x_2, \dotsm, x_n) \) such that their dot product is 0 (i.e. \( \sum a_i x_i = 0 \) ).


Please note that this is beta; I have not implemented a full set of exit conditions!

Code:

NINT(x)
BEGIN
  local n:=IP(x);
  if x > n+.5 then return(n+1); end;
  if x < n-.5 then return(n-1); end;
  return(n);
END;


EXPORT PSLQ(lst)
BEGIN

local MA,MB,MH;
local j,k,n,m,r;
local s:={}, y:={};
local t;
local g:=2/sqrt(3);
local M1,M2;
local t0,t1,t2,t3,t4;
local run:=0;


n:=size(lst);

MA:=idenmat(n);
MB:=MA;
MH:=MA;
MH:=delcol(MH,n);

for k from 1 to n do
  s(k):=sqrt(sum(lst(j)^2,j,k,n));
end;

t:=1/s(1);

//for k from 1 to n do
//  y(k):=t*lst(k);
//  s(k):=t*s(k);
//end;

y:=t*lst;
s:=t*s;

for j from 1 to n-1 do
  MH(j,j):=s(j+1)/s(j);
  for k from j+1 to n do
    MH(k,j):=-y(k)*y(j)/(s(j)*s(j+1));
  end;
end;


for k from 2 to n do
  for j from k-1 downto 1 do
    t:=NINT(MH(k,j)/MH(j,j));
    y(j):=y(j)+t*y(k);

    for m from 1 to j do
      MH(k,m):=MH(k,m)-t*MH(j,m);
    end;

    for m from 1 to n do
      MA(k,m):=MA(k,m)-t*MA(j,m);
      MB(m,j):=MB(m,j)+t*MB(m,k);
    end;

  end;
end;


// iterate
repeat

m:=1; M1:=g*abs(MH(1,1));

for k from 2 to n-1 do
  M2:=g^k*abs(MH(k,k));
  if M2>M1 then
    m:=k; M1:=M2;
  end;
end;


M2:=y(m);
y(m):=y(m+1);
y(m+1):=M2;

MH:=swaprow(MH,m,m+1);
MA:=swaprow(MA,m,m+1);
MB:=swapcol(MB,m,m+1);

if m<=n-2 then 
  t0:=sqrt(MH(m,m)^2+MH(m,m+1)^2);
  t1:=MH(m,m)/t0;
  t2:=MH(m,m+1)/t0;
  for j from m to n do
    t3:=MH(j,m);
    t4:=MH(j,m+1);
    MH(j,m):=t1*t3+t2*t4;
    MH(j,m+1):=-t2*t3+t1*t4;
  end;
end;

for k from m+1 to n do
  for j from min(k-1,m+1) downto 1 do
    t:=NINT(MH(k,j)/MH(j,j));
    y(j):=y(j)+t*y(k);
    for r from 1 to j do
      MH(k,r):=MH(k,r)-t*MH(j,r);
    end;
    for r from 1 to n do
      MA(k,r):=MA(k,r)-t*MA(j,r);
      MB(r,j):=MB(r,j)+t*MB(r,k);
    end;

  end;
end;

m:=1; M1:=abs(y(1));
for j from 2 to n do
  M2:=abs(y(j));
  if M2<M1 then M1:=M2; m:=j; end;
end;

run:=(M1<1E-10);
until run end;

return(col(MB,m));

END;

This can be used to solve a number of interesting "algebraic form" problems. For example, the number 3.96811878507 is the decimal approximation of \( x=\sqrt{3} + \sqrt{5} \). The value \( x \) is a root of the polynomial \( x^4 - 16x^2 + 4 \). Moreover, numbers of the form \( \sqrt{a} + \sqrt{b} \) are roots of \( x^4 - 2(a+b)x^2 + (a-b)^2 \). We can use PSLQ on the vector \( ( x^4, x^3, x^2, x, 1 ) \) to obtain the result: \( ( 1, 0, -16, 0, 4 ) \). Additionally solving \( -2(a+b) = -16 \) and \( (a-b)^2=4 \) simultaneously results in the values a=3 and b=5.

Code:


X:=SQRT(3)+SQRT(5);
PSLQ({X^4,X^3,X^2,X,1); // returns {1,0,-16,0,4}

Another example may be to rewrite \( x \approx 7.09439510239 \) in the preferred form \( 2\pi/3 + 5 \). Notice that \( x \) satisfies: \( 0 = 2\pi -3x + 15 \). So if we apply PLSQ to the vector \( (\pi, 7.09439510239, 1 ) \) we get the vector \( (2,-3,15) \)

Code:


X:=2*PI/3+5;
PSLQ({PI,X,1}); // returns {2,-3,15}

Also, if we have a number of the form \( x= r_1 + \sqrt{r_2} \) where the \( r_i \)'s are rational numbers, then simply apply PSLQ to \( ( x^2, x, 1) \) since such numbers are roots of some quadratic polynomial. The result from PSLQ -- \( ( a,b,c ) \) -- are the coefficient of the quadratic polynomial, and can then be used to rewrite \( x\) into a nice algebraic form using \( x = \frac{-b \pm \sqrt{b^2-4ac}}{2a} \)

Code:


X:=2/3+4*SQRT(3);
PSLQ({X^2,X,1}); // returns {9,-12,-428}