PSLQ
11-29-2017, 08:12 PM (This post was last modified: 12-05-2017 06:58 PM by Han.)
Post: #1
 Han Senior Member Posts: 1,775 Joined: Dec 2013
PSLQ
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}

Graph 3D | QPI | SolveSys
11-29-2017, 08:42 PM
Post: #2
 Han Senior Member Posts: 1,775 Joined: Dec 2013
RE: PSLQ
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);

Graph 3D | QPI | SolveSys
11-29-2017, 08:52 PM (This post was last modified: 11-29-2017 09:32 PM by Han.)
Post: #3
 Han Senior Member Posts: 1,775 Joined: Dec 2013
RE: PSLQ
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;

Graph 3D | QPI | SolveSys
11-29-2017, 08:58 PM
Post: #4
 salvomic Senior Member Posts: 1,364 Joined: Jan 2015
RE: PSLQ
Thank you!
very interesting, a great algorithm.

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
11-29-2017, 09:28 PM (This post was last modified: 11-29-2017 09:30 PM by salvomic.)
Post: #5
 salvomic Senior Member Posts: 1,364 Joined: Jan 2015
RE: PSLQ
(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 ...

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
11-29-2017, 09:33 PM (This post was last modified: 11-29-2017 09:33 PM by Han.)
Post: #6
 Han Senior Member Posts: 1,775 Joined: Dec 2013
RE: PSLQ
(11-29-2017 09:28 PM)salvomic Wrote:  Han,
QROOT() by me present a Syntax error maybe in (-b/2a)) at line 17 ...

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

Graph 3D | QPI | SolveSys
11-29-2017, 09:35 PM
Post: #7
 salvomic Senior Member Posts: 1,364 Joined: Jan 2015
RE: PSLQ
(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.

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
11-30-2017, 02:16 AM
Post: #8
 AlexFekken Member Posts: 151 Joined: May 2016
RE: PSLQ
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.
11-30-2017, 03:14 AM (This post was last modified: 11-30-2017 03:14 AM by Han.)
Post: #9
 Han Senior Member Posts: 1,775 Joined: Dec 2013
RE: PSLQ
(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

Graph 3D | QPI | SolveSys
11-30-2017, 03:56 AM
Post: #10
 AlexFekken Member Posts: 151 Joined: May 2016
RE: PSLQ
(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 :-)
11-30-2017, 11:05 AM
Post: #11
 ggauny@live.fr Senior Member Posts: 463 Joined: Nov 2014
RE: PSLQ
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.

Attached File(s) Thumbnail(s)

Gérard.
11-30-2017, 12:19 PM
Post: #12
 Han Senior Member Posts: 1,775 Joined: Dec 2013
RE: PSLQ
(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. [])

Graph 3D | QPI | SolveSys
11-30-2017, 12:29 PM
Post: #13
 ggauny@live.fr Senior Member Posts: 463 Joined: Nov 2014
RE: PSLQ
Thanks.

Gérard.
11-30-2017, 11:10 PM
Post: #14
 BruceH Member Posts: 175 Joined: Dec 2013
RE: PSLQ
(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.
12-02-2017, 11:55 PM
Post: #15
 AlexFekken Member Posts: 151 Joined: May 2016
RE: PSLQ
(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.
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)