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}

[Han]

Code for identifying \( x = \sqrt{a} + \sqrt{b} \):

Code:


EXPORT SSQR(x)
BEGIN
  local a:=PSLQ({x^4,x^3,x^2,x,1});
  local r:=sqrt(abs(a(5)));

  if (a(2) OR a(4)) then return(x); end;
  if abs(a(1))<>1 then return(x); end;
  if a(1)<0 then a:=-a; end;
  if a(5)<0 then return(x); end;
  if IP(r)<>r then return(x); end;
  if a(3) MOD 2 then return(x); end;
  a:=-a(3)/2;
  if a<r then
    local t:=r;
    r:=a; a:=t;
  end;

  return(expr("'sqrt(" + (a-r)/2 + ")+sqrt(" + (a+r)/2 + ")'"));

END;

Requires PSLQ above. Example:

Code:

X:=SQRT(3)+SQRT(5);
SSQR(X);

[Han]

Code for identifying \( x = p_1/q_1 + \sqrt{p_2/q_2} \):

Code:

EXPORT QROOT(x)
BEGIN
  local a:=PSLQ({x^2,x,1});
  if a(1)<0 then a:=-a; end;

  local b:=a(2);
  local c:=a(3);
  a:=a(1);

  if IP(a)<>a then return(x); end;
  if IP(b)<>b then return(x); end;
  if IP(c)<>c then return(x); end;

  local t:=(-b+sqrt(b^2-4*a*c))/(2*a);

  if x==t then
    return(expr("'" + string(QPI(-b/(2*a))) + "+" + string(QPI(sqrt(b^2-4*a*c)/(2*a))) + "'"));
  else
    return(expr("'" + string(QPI(-b/(2*a))) + "-" + string(QPI(sqrt(b^2-4*a*c)/(2*a))) + "'"));
  end;

END;

[salvomic]

Thank you!
very interesting, a great algorithm.

[salvomic]

(11-29-2017 08:52 PM)Han Wrote:  ...

Code:

...
  if x==t then
    return(expr("'" + string(QPI(-b/2a)) + "+" + string(QPI(sqrt(b^2-4*a*c)/(2*a))) + "'"));
  else
    return(expr("'" + string(QPI(-b/2a)) + "-" + string(QPI(sqrt(b^2-4*a*c)/(2*a))) + "'"));
  end;

END;

Han,
QROOT() by me present a Syntax error maybe in (-b/2a)) at line 17 ...
please, check it.

[Han]

(11-29-2017 09:28 PM)salvomic Wrote:  Han,
QROOT() by me present a Syntax error maybe in (-b/2a)) at line 17 ...
please, check it.

I must have copied the wrong version of the code (one on calc vs one typed out onto notepad). Should be correct now.

[salvomic]

(11-29-2017 09:33 PM)Han Wrote:  I must have copied the wrong version of the code (one on calc vs one typed out onto notepad). Should be correct now.

yes, it is! now it's ok.

[AlexFekken]

I like the irony, or if you like: paradox, of these functions. So sorry, if you are not interested, for the caveat rant for those expecting silver bullets wherever they go.

The motivation for PSLQ (and QPI) seems to be that "arbitrary numbers" are likely to be of a special type: rational, semi-rational (in QPI sense), or algebraic (or in the same field extension of Q as some other numbers). This is also illustrated by the various naming suggestions for QPI, especially those involving the word "exact" or otherwise suggesting that you are (always, usually?) getting a "better" representation of a number. And in situations where that is actually the case these functions are of course very useful.

But, as a sanity check, from a mathematical perspective these special types of numbers have measure zero, i.e. the probability that a randomly picked real number (within a given range, say between 0 and 1) belongs to one of these special types is 0, rather than close to 1. So mathematically your expectation that these functions will be useful when working with arbirtrary numbers should be 0 too.
And at the other end of the scale, i.e. when you take into account that machine numbers are not "arbitrary (real) numbers" to begin with: all machine floats are already rational, so why bother about finding "better" (?) rational, semi-rational or algebraic representations? This is borderline pseudo-science, depending on the context.

And on top of this it is all relative as well. "3.141592653" coming out of a long calculation on the Prime may be "better" represented as pi, whereas coming out of a similar calculation using Free32 it is extremely unlikely that pi is a better representation (unless there is an issue with the calculation itself). Who is going to check whether the accuracy justifies or invalidates the "better" representation (not your average student, I would expect)?

So it seems that we are dealing with some shadowy area between Q and R, representing it all with just a subset of Q in a way that actually seems to make some sense. But only as long as you know what you are doing.

[Han]

(11-30-2017 02:16 AM)AlexFekken Wrote:  [b]But, as a sanity check, from a mathematical perspective these special types of numbers have measure zero, i.e. the probability that a randomly picked real number (within a given range, say between 0 and 1) belongs to one of these special types is 0, rather than close to 1. So mathematically your expectation that these functions will be useful when working with arbirtrary numbers should be 0 too.

I suppose it depends on whether you are looking at pure measure, or "relative measure" -- in the sense that most students are given problems in which the answers are generally not from the set of reals, but instead from the set of algebraic numbers. (I.e. they are guaranteed with probability 1.)

But even for folks such as physicists are interested in these "nice" numbers because they show up as constants in various "important" physics equations. In fact, there is an entire field of study devoted to finding the "identity" of various constants all due to the advent of more efficient algorithms for integer coefficients such as PSLQ.

http://www2.lbl.gov/Science-Articles/Arc...rithm.html

[AlexFekken]

(11-30-2017 03:14 AM)Han Wrote:  I suppose it depends on whether you are looking at pure measure, or "relative measure" -- in the sense that most students are given problems in which the answers are generally not from the set of reals, but instead from the set of algebraic numbers. (I.e. they are guaranteed with probability 1.)

It has been a long time since I was a student and I agree that this was certainly the case in "my days". Really sad if this is still the case though.
It also begs the question of why you would allow the use of calculators for such problems (or why you would give them such problems when you allow the use of calculators)...

(11-30-2017 03:14 AM)Han Wrote:  But even for folks such as physicists are interested in these "nice" numbers because they show up as constants in various "important" physics equations. In fact, there is an entire field of study devoted to finding the "identity" of various constants all due to the advent of more efficient algorithms for integer coefficients such as PSLQ.

I am sure there is value in that in certain areas of physics.
But as the article points out: "Very high precision arithmetic is needed by PSLQ, or else nonsense results are obtained" which, if you add the obvious requirement that you need to relate that to the accuracy of your input, is sort of the main point that I was trying to make.

Better finish by saying that I do appreciate your efforts in implementating these algorithms. Just worried about some of the people who will use them with their eyes closed :-)

[ggauny@live.fr]

Hello Han,

Really amazing, as usual from you !!!

Only a question, Instead getting {1,0,-16,0,4}, I get here this [1 0 −16 0 4]

Am-I wrong ?

Thanks.

[Han]

(11-30-2017 11:05 AM)ggauny@live.fr Wrote:  Hello Han,

Really amazing, as usual from you !!!

Only a question, Instead getting {1,0,-16,0,4}, I get here this [1 0 −16 0 4]

Am-I wrong ?

Thanks.

No, I just typed with the wrong symbols ({} vs. [])

[ggauny@live.fr]

Thanks.

[BruceH]

(11-30-2017 02:16 AM)AlexFekken Wrote:  But, as a sanity check, from a mathematical perspective these special types of numbers have measure zero [...]

This is very true.

However, there are good, non-mathematical reasons for having such a function: I refer you to this excellent post from Valentin Albillo in which he determines the underlying cause of an HP-35S ROM bug by analysing the faulty values returned.

[AlexFekken]

(11-30-2017 11:10 PM)BruceH Wrote:  However, there are good, non-mathematical reasons for having such a function: I refer you to this excellent post from Valentin Albillo in which he determines the underlying cause of an HP-35S ROM bug by analysing the faulty values returned.

Nice reference, thanks for that.