PSLQ in Yacas and MuPAD

-------------
PSLQ in Yacas
-------------

Reference : pslq
proteus-yacas-windows-21feb2001.zip
<http://www.xs4all.nl/~apinkus/backups/>

See <http://www.xs4all.nl/~apinkus/infoindex.html>

-- See attached forward pslq.mu
(pslq algorithm source written out in MuPAD code) --

/**************************************************************************#
# The PSLQ Integer Relation Algorithm #
# #
# Aut.: Helaman R.P. Ferguson and David Bailey "A Polynomial #
# Time, Numerically Stable Integer Relation Algorithm" #
# (RNR Technical Report RNR-92-032) helaman@... #
# Ref.: David Bailey and Simon Plouffe "Recognizing Numerical Constants" #
# dbailey@... #
# Cod.: Raymond Manzoni raymman@... #
#**************************************************************************#
# Creation:97/11 #
# New termination criteria:97/12/15 #
# this code is free... #

Ported to Yacas 2000 Ayal Pinkus.

Given a list of constants x find coefficients sol[i] such that
sum(sol[i]*x[i], i=1..n) = 0 (where n=Length(x))

x is the list of real expressions
N(x[i]) must evaluate to floating point numbers!
precision is the number of digits needed for completion;
must be greater or equal to log10(max(sol[i]))*n
returns the list of solutions with initial precision
and the confidence (the lower the better)

Example:

In> Pslq({2*Pi-4*Exp(1),Pi,Exp(1)},20)
Out> {1,-2,4};

*/

Pslq(x, precision) :=
[
Local (ndigits, gam, A, B, H, n, i, j, k, s, y, tmp, t, m, maxi, gami,
t0, t1, t2, t3, t4, mini, Confidence, norme,result);
n:=Length(x);
ndigits:=GetPrecision();
Precision(precision+2);
Confidence:=10^(-MathFloor(N(Eval(precision/3))));
gam:=N(Sqrt(4/3));
For (i:=1, i<=n,i++) x[i]:=N(x[i]);

A:=Identity(n); /*A and B are of Integer type*/
B:=Identity(n); /*but this doesn't speed up*/
s:=ZeroVector(n);
y:=ZeroVector(n);

For(k:=1,k<=n,k++)
[
tmp:=0;
For (j:=k,j<=n,j++) tmp:=tmp + x[j]^2;
s[k]:=N(Eval(Sqrt(tmp)));
];
tmp:=s[1];
For (k:= 1,k<= n,k++)
[
y[k]:=N(Eval(x[k]/tmp));
s[k]:=N(Eval(s[k]/tmp));
];
H:=ZeroMatrix(n, n-1);

For (i:=1,i<= n,i++)
[

if (i <= n-1) [ H[i][i]:=s[i + 1]/s[i]; ];
For (j:= 1,j<=i-1,j++)
[
H[i][j]:= -(y[i]*y[j])/(s[j]*s[j + 1]);
];
];

For (i:=2,i<=n,i++)
[
For (j:=i-1,j>= 1,j--)
[
t:=Round(H[i][j]/H[j][j]);
y[j]:=y[j] + t*y[i];
For (k:=1,k<=j,k++) [ H[i][k]:=H[i][k]-t*H[j][k]; ];
For (k:=1,k<=n,k++)
[
A[i][k]:=A[i][k]-t*A[j][k];
B[k][j]:=B[k][j] + t*B[k][i];
];
];
];
Local(found);
found:=False;

While (Not(found))
[
m:=1;
maxi:=N(gam*Abs(H[1][1]));
gami:=gam;
For (i:= 2,i<= n-1,i++)
[
gami:=gami*gam;
tmp:=N(gami*Abs(H[i][i]));
if (maxi < tmp)
[
maxi:=tmp;
m:=i;
];
];
tmp:=y[m + 1];
y[m + 1]:=y[m];
y[m]:=tmp;
For (i:= 1,i<=n,i++)
[
tmp:=A[m + 1][ i];
A[m + 1][ i]:=A[m][ i];
A[m][ i]:=tmp;
tmp:=B[i][ m + 1];
B[i][ m + 1]:=B[i][ m];
B[i][ m]:=tmp;
];
For (i:=1,i<=n-1,i++)
[
tmp:=H[m + 1][ i];
H[m + 1][ i]:=H[m][ i];
H[m][ i]:=tmp;
];
if (m < n-1)
[
t0:=N(Eval(Sqrt(H[m][ m]^2 + H[m][ m + 1]^2)));
t1:=H[m][ m]/t0;
t2:=H[m][ m + 1]/t0;
For (i:=m,i<=n,i++)
[
t3:=H[i][ m];
t4:=H[i][ m + 1];
H[i][ m]:=t1*t3 + t2*t4;
H[i][ m + 1]:= -t2*t3 + t1*t4;
];
];
For (i:= 1,i<= n,i++)
[
For (j := Min(i-1, m + 1),j>= 1,j--)
[
t:=Round(H[i][ j]/H[j][ j]);
y[j]:=y[j] + t*y[i];
For (k:=1,k<=j,k++) H[i][ k]:=H[i][ k]-t*H[j][ k];
For (k:= 1,k<=n,k++)
[
A[i][ k]:=A[i][ k]-t*A[j][ k];
B[k][ j]:=B[k][ j] + t*B[k][ i];
];
];
];
/* Precision(10);*/ /*low precision*/
maxi := (H[1] . H[1]);
For (j:=2,j<=n,j++)
[
tmp:=(H[j] . H[j]);
if (maxi < tmp) [ maxi:=tmp; ];
];
norme:=N(Eval(1/Sqrt(maxi)));
m:=1;
mini:=N(Eval(Abs(y[1])));
maxi:=mini;
For (j:=2,j<=n,j++)
[
tmp:=N(Eval(Abs(y[j])));
if (tmp < mini)
[
mini:=tmp;
m:=j;
];
if (tmp > maxi) [ maxi:=tmp; ];
];
/* following line may be commented */
/* Echo({"Norm bound:",norme," Min=",mini," Conf=",mini/maxi}); */
if ((mini/maxi) < Confidence) /*prefered to : if mini < 10^(-
precision) then*/
[
/* following line may be commented */
/* Echo({"Found with Confidence ",mini/maxi}); */
Precision(ndigits);
result:=Transpose(B)[m];
found:=True;
]
else
[
maxi:=Abs(A[1][ 1]);
For (i:=1,i<=n,i++)
[
For (j:=1,j<=n,j++)
[
tmp:=Abs(A[i][ j]);
if (maxi < tmp) [ maxi:=tmp;];
];
];
if (maxi > 10^(precision))
[
Precision(ndigits);
result:=Fail;
found:=True;
];
Precision(precision+2);
];
];
result;
];

/* end of file */

-------------
PSLQ in MuPAD
-------------

Ayal Pinkus wrote (private e-mail)

<< I found this piece of code in the MuPAD scripts. Please
find attached the relevant file, pslq.mu, which is the
pslq algorithm written out in MuPAD code. It states
explicitly that it is free >>

file pslq.mu

//***************************************************************************
// The PSLQ Integer Relation Algorithm

//

// Aut.: Helaman R.P. Ferguson and David Bailey "A Polynomial Time,
// Numerically Stable
// Integer Relation Algorithm" (RNR Technical Report RNR-92-032)
// helaman@...
// Ref.: David Bailey and Simon Plouffe "Recognizing Numerical Constants"
// dbailey@...
// Cod.: Raymond Manzoni raymman@...

//****************************************************************************
// Creation:97/11
// New termination criteria:97/12/15
// this code is free...

/*++ Given a list of constants x find coefficients sol[i] such that
sum(sol[i]*x[i], i=1..n) = 0 (where n=nops(x))

x is the list of real expressions
float(x[i]) must evaluate to numeric type!
precision is the number of digits needed for completion;
must be greater or equal to log10(max(sol[i]))*n
returns the list of solutions with initial precision
and the confidence (the lower the better)

Example:

>> pslq([2*PI+E,PI,E],20);

Found with Confidence , 1.4780013e-28
+- -+
| 1, -2, -1 |
+- -+
++*/

misc::pslq:=
proc(x:DOM_LIST, precision:Type::PosInt)
local ndigits, gam, A, B, H, n, i, j, k, s, y, tmp, t, m, maxi, gami,
t0, t1, t2, t3, t4, mini, Confidence, norme, MExpr;
save DIGITS;
begin
MExpr:=Dom::Matrix();
n:=nops(x);
ndigits:=DIGITS;
DIGITS:=precision+2;
Confidence:=10^(-floor(precision/3));
gam:=float(sqrt(4/3));
for i from 1 to n do
x[i]:=float(x[i]);
if not testtype(x[i], Type::Real)=TRUE then
error("List must consist of real numbers")
end_if
end_for;
A:=MExpr(n, n, 1, Diagonal); //A and B are of Integer type
B:=MExpr(n, n, 1, Diagonal); //but this doesn't speed up
s:=MExpr(n, 1);
y:=MExpr(n, 1);
for k from 1 to n do
tmp:=0;
for j from k to n do
tmp:=tmp + x[j]^2
end_for;
s[k]:=sqrt(tmp)
end_for;
tmp:=s[1];
for k from 1 to n do
y[k]:=x[k]/tmp;
s[k]:=s[k]/tmp
end_for;
H:=MExpr(n, n-1);
for i from 1 to n do
if i <= n-1 then
H[i, i]:=s[i + 1]/s[i]
end_if;
for j from 1 to i-1 do
H[i, j]:=-(y[i]*y[j])/(s[j]*s[j + 1])
end_for
end_for;
for i from 2 to n do
for j from i-1 downto 1 do
t:=round(H[i, j]/H[j, j]);
y[j]:=y[j] + t*y[i];
for k from 1 to j do
H[i, k]:=H[i, k]-t*H[j, k]
end_for;
for k from 1 to n do
A[i, k]:=A[i, k]-t*A[j, k];
B[k, j]:=B[k, j] + t*B[k, i]
end_for
end_for
end_for;
while TRUE do
m:=1;
maxi:=float(gam*abs(H[1, 1]));
gami:=gam;
for i from 2 to n-1 do
gami:=gami*gam;
tmp:=float(gami*abs(H[i, i]));
if maxi < tmp then
maxi:=tmp;
m:=i
end_if
end_for;
tmp:=y[m + 1];
y[m + 1]:=y[m];
y[m]:=tmp;
for i from 1 to n do
tmp:=A[m + 1, i];
A[m + 1, i]:=A[m, i];
A[m, i]:=tmp;
tmp:=B[i, m + 1];
B[i, m + 1]:=B[i, m];
B[i, m]:=tmp
end_for;
for i from 1 to n-1 do
tmp:=H[m + 1, i];
H[m + 1, i]:=H[m, i];
H[m, i]:=tmp
end_for;
if m < n-1 then
t0:=sqrt(H[m, m]^2 + H[m, m + 1]^2);
t1:=H[m, m]/t0;
t2:=H[m, m + 1]/t0;
for i from m to n do
t3:=H[i, m];
t4:=H[i, m + 1];
H[i, m]:=t1*t3 + t2*t4;
H[i, m + 1]:=-t2*t3 + t1*t4
end_for
end_if;
for i from 1 to n do
for j from min(i-1, m + 1) downto 1 do
t:=round(H[i, j]/H[j, j]);
y[j]:=y[j] + t*y[i];
for k from 1 to j do
H[i, k]:=H[i, k]-t*H[j, k]
end_for;
for k from 1 to n do
A[i, k]:=A[i, k]-t*A[j, k];
B[k, j]:=B[k, j] + t*B[k, i]
end_for
end_for
end_for;
DIGITS:=8; //low precision
maxi := linalg::scalarProduct(linalg::row(H, 1), linalg::row(H, 1));
for j from 2 to n do
tmp:=linalg::scalarProduct(linalg::row(H, j), linalg::row(H, j));
if maxi < tmp then
maxi:=tmp
end_if
end_for;
norme:=1/sqrt(maxi);
m:=1;
mini:=float(abs(y[1]));
maxi:=mini;
for j from 2 to n do
tmp:=float(abs(y[j]));
if tmp < mini then
mini:=tmp;
m:=j
end_if;
if tmp > maxi then
maxi:=tmp;
end_if
end_for;
// following line may be commented
userinfo(5,"Norm bound:".expr2text(norme)." Min=".expr2text(mini).
" Conf=".expr2text(mini/maxi));
if (mini/maxi) < Confidence then //prefered to : if mini < 10^(-
precision) then
userinfo(5,"Found with Confidence ".expr2text(mini/maxi));
DIGITS:=ndigits;
linalg::transpose(linalg::col(B, m));
break
end_if;
maxi:=abs(A[1, 1]);
for i from 1 to n do
for j from 1 to n do
tmp:=abs(A[i, j]);
if maxi < tmp then
maxi:=tmp
end_if
end_for
end_for;
if maxi > 10^(precision) then
DIGITS:=ndigits;
return(FAIL)
end_if;
DIGITS:=precision+2;
end_while;
end_proc:

// end of file