r"""
An implementation of the PSLQ algorithm developed by H. R. P. Ferguson, D. H. Bailey and S. Arno.
AUTHORS:
-- Agnes Jany (2009-05-10)
EXAMPLES:
sage: t=vector([11,27,31])
sage: p=pslq(t)
Relation found:
(-8, 9, -5)
sage: p=pslq(t,iteration=true)
Iteration no. 0
Iteration no. 1
Iteration no. 2
Iteration no. 3
Relation found:
(-8, 9, -5)
sage: p=pslq(t,digits=1000,iteration=true)
Iteration no. 0
Iteration no. 1
Iteration no. 2
Iteration no. 3
Iteration no. 4
Relation found:
(-6, -1, 3)
sage: p=pslq(t,tau=1.01,iteration=true)
Iteration no. 0
Iteration no. 1
Iteration no. 2
Iteration no. 3
Iteration no. 4
Iteration no. 5
Relation found:
(-8, 9, -5)
sage: p=pslq(t,tau=1.99,iteration=true)
Iteration no. 0
Iteration no. 1
Relation found:
(0, -31, 27)
sage: p=pslq(t,tau=1,iteration=true)
Choose tau such that 1<tau<2 !
sage: p=pslq(t,digits=-5,iteration=true)
Set digits to a positive integer!
sage: t=vector(RDF,[11,27,31])
sage: p=pslq(t)
Relation found:
(-6, -1, 3)
sage: t=vector(RealField(213),[11,27,31])
sage: p=pslq(t)
Relation found:
(-1, 5, -4)
sage: alp=vector(RealField(212),17)
sage: for i in range(0,17,1):
sage: alp[i]=((3^(1/4)-2^(1/4))^i)
sage: p=pslq(alp,digits=62,tau=1.000006145)
Relation found:
(1, 0, 0, 0, -3860, 0, 0, 0, -666, 0, 0, 0, -20, 0, 0, 0, 1)
sage:p=pslq(alp,tau=1.000006145)
Relation found:
(192, -1536, 162, 85, 287, 274, 814, -503, 1148, 786, -206, -464, -279,
83, 144, -316, -725)
sage: p=pslq(alp,digits=100,tau=1.000006145)
/usr/local/sage/local/bin/sage-sage: line 347: 3995 Killed
python "$@"
sage: t=vector([1,(pi),(sqrt(2)),(sqrt(3)),(ln(1+sqrt(2))),(ln(2+sqrt(3))),
sage: 0.926390055174046729218163586547779014444960190107335046732521921271418504594036683829313473075349968212])
sage: p=pslq(t)
Relation found:
(4, -7, 17, -6, 21, 42, -75)
sage: t=vector([gamma(1/2),gamma(7/2),gamma(13/2)])
sage: p=pslq(t)
Relation found:
(15, -8, 0)
sage: t=vector([gamma(1/3),1,sqrt(5)])
sage: p=pslq(t,30)
Relation found:
(727989880380837, -562243113417127, -620731142335100)
sage: p=pslq(t,100)
/usr/local/sage/local/bin/sage-sage: line 347: 4092 Segmentation fault
python "$@"
"""
#*****************************************************************************
# Copyright (C) 2009 Agnes Jany <agnes.jany@yahoo.de>
#
# Distributed under the terms of the GNU General Public License (GPL)
# http://www.gnu.org/licenses/
#*****************************************************************************
class pslq(object):
def __init__(self, vec ,digits=50,tau=(4/3)^0.5,iteration=false,printing=true):
r"""
Initialization step of the PSLQ algorithm.
INPUT:
vec -- input vector
digits -- integer
tau -- real number
iteration -- boolean value
NOTES:
Define vectors x and s. Define matrices H,A and B. Perform first reduction. See [FBA] p.8.
REFERENCES:
[FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno :
"Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
if tau<=1 or tau>=2:
print "Choose tau such that 1<tau<2 !"
else:
if digits <=0:
print "Set digits to a positive integer!"
else:
self.printing=printing
sys.setrecursionlimit(100000)
self.it=iteration
self.d=digits
self.gamma=sqrt((4*tau^2)/(4-tau^2))
self.count=-1
n=len(vec)
RF=RealField(212)
x=vector( RF,(1/norm(vec))*(vec))
s=self._set_s_(x,n)
H=self._define_H_(s,x,n)
A1=MatrixSpace(ZZ,n)
A=A1.identity_matrix()
B1=MatrixSpace(ZZ,n)
B=B1.identity_matrix()
H,A,B,x=self._Reduce_(H,A,B,x,n,false)
self._iterate_(H,A,B,x,n)
def _iterate_(self,H,A,B,x,n):
r"""
This function executes the iteration steps of PSLQ until it is halted by _terminate_().
INPUT:
H -- matrix
A,B -- matrices with integer values and A=B^-1
x -- vector
n -- length of the input vector
NOTES:
This function uses the algorithm of [FBA] to find integer relations.
REFERENCES:
[FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno :
"Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
if self._terminate_(H,A,B,x,n)==true:
if self.printing:
print self.reason[0]
print self.reason[1]
return self.reason
else:
self.count=self.count+1
if self.it:
print "Iteration no.", self.count
r=self._compare_gamma_(self.gamma,H,n)
if r<n-2:
self._predef_(H,n,r)
self.R=self._define_R_(r,n)
x=x*self.R
H=self.R*H
A=self.R*A
B=B*self.R
if r<n-2:
Q=self._define_Q_(H,n,r)
H=H*Q
H,A,B,x = self._Reduce_(H,A,B,x,n,r)
self._iterate_(H,A,B,x,n)
def _set_s_(self,x,n):
r"""
Define a vector of partial sums.
INPUT:
x -- vector
n -- length of the input vector
OUTPUT:
s -- vector of partial sums.
NOTES:
See [FBA] p.4.
REFERENCES:
[FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno :
"Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
RF=RealField(212)
j=0
s=vector(RF,copy(x))
while j<n:
k=j
sum=0
while k<n:
sum=sum+x[k]^2
k=k+1
s[j]=sqrt(sum)
j=j+1
return s
def _nint_(self,x):
r"""
Nearest integer function.
INPUT:
x -- real number.
OUTPUT:
xnew -- nearest integer to x. For exact half-integer values: integer with greater absolute value.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
if x/0.5 not in ZZ or x in ZZ:
xnew=round(x)
if x/0.5 in ZZ:
if x>0:
xnew=ceil(x)
else:
xnew=floor(x)
return xnew
def _define_H_(self,s,x,n):
r"""
Define a $n\times n-1$ lower trapezodial matrix H.
INPUT:
s -- vector of partial sums
x -- vector
n -- length of the input vector
OUTPUT:
H -- matrix
NOTES:
See [FBA] p.5.
REFERENCES:
[FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno :
"Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
RF=RealField(212)
H=matrix(RF,n,n-1)
for i in range(0,n,1):
for j in range(0,n-1,1):
if i<j:
H[i,j]=0
if i==j:
H[i,j]=s[i+1]/s[i]
if j<i:
H[i,j]=(-x[i]*x[j])/(s[j]*s[j+1])
return H
def _compare_gamma_(self,gamma,H,n):
r"""
Choose an integer j such that $\gamma^j |h_{jj}|\geq \gamma^i |h_{ii}|$ for all $1\leq i \leq n-1$.
INPUT:
gamma -- constant real number
H -- matrix
n -- length of the input vector
OUTPUT:
j -- integer
NOTES:
See [FBA] p.9.
REFERENCES:
[FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno :
"Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
j=0
while j<=n-2:
i=0
while i<=n-2:
g1=(gamma^j)*abs(H[j,j])
g2=(gamma^i)*abs(H[i,i])
if g1<g2:
break
else:
i=i+1
if i==n-1:
return j
break
else:
j=j+1
def _define_R_(self,j,n):
r"""
Define a permutation matrix R.
INPUT:
j -- integer
n -- length of the input vector
OUTPUT:
R -- matrix
NOTES:
See [FBA] p.9.
REFERENCES:
[FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno :
"Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
R1=MatrixSpace(ZZ,n)
R=R1.identity_matrix()
R[j,j]=0
R[j,j+1]=1
R[j+1,j+1]=0
R[j+1,j]=1
return R
def _predef_(self,H,n,r):
r"""
Remember certain entries of H before defining Q.
INPUT:
H -- matrix
r -- integer
n -- length of the input vector
NOTES:
See [FBA] p.8-9.
REFERENCES:
[FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno :
"Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
self.Qa=H[r,r]
self.Qb=H[r+1,r]
self.Qc=H[r+1,r+1]
self.Qd=sqrt(self.Qc^2+self.Qb^2)
def _define_Q_(self,H,n,r):
r"""
Define a $n-1\times n-1$ matrix Q.
INPUT:
H -- matrix
r -- integer
n -- length of the input vector
OUTPUT:
Q -- matrix
NOTES:
See [FBA] p.9.
REFERENCES:
[FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno :
"Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
RF=RealField(500)
b=self.Qb
c=self.Qc
d=self.Qd
Q1=MatrixSpace(RF,n-1)
Q=Q1.identity_matrix()
if not (b/d).is_infinity():
Q[r,r]=b/d
if not (-c/d).is_infinity():
Q[r,r+1]=-c/d
if not (c/d).is_infinity():
Q[r+1,r]=c/d
if not (b/d).is_infinity():
Q[r+1,r+1]=b/d
return Q
def _Reduce_(self,H,A,B,x,n,r):
r"""
Perform Hermite reduction.
INPUT:
H,A,B -- matrices
x -- vector
r -- integer
n -- length of the input vector
OUTPUT:
H,A,B -- matrices
x -- vector
NOTES:
See [FBA] p.8-9 and p.17-18.
REFERENCES:
[FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno :
"Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
if self.count==-1:
for i in range(1,n,1):
for j in range(i-1,-1,-1):
t=Integer(self._nint_((H[i,j]/H[j,j])))
x[j]=(copy(x[j])+t*x[i])
for k in range(0,j+1,1):
H[i,k]=(copy(H[i,k])-t*H[j,k])
for k in range(0,n,1):
A[i,k]=(copy(A[i,k])-t*A[j,k])
B[k,j]=(copy(B[k,j])+t*B[k,i])
return H,A,B,x
else:
for i in range(r+1,n,1):
z=min(i-1,r+1)
for j in range(z,-1,-1):
t=Integer(self._nint_((H[i,j]/H[j,j])))
x[j]=(copy(x[j])+t*x[i])
for k in range(0,j+1,1):
H[i,k]=(copy(H[i,k])-t*H[j,k])
for k in range(0,n,1):
A[i,k]=(copy(A[i,k])-t*A[j,k])
B[k,j]=(copy(B[k,j])+t*B[k,i])
return H,A,B,x
def _zero_H_(self,H,n):
r"""
Tests for zeros in the diagonal elements of H.
INPUT:
H -- matrices
n -- length of the input vector
OUTPUT:
boolean value
NOTES:
See [FBA] p.9.
REFERENCES:
[FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno :
"Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
i=0
while i<=n-2:
if (abs(H[i,i]))<10^(-1000):
break
else:
i=i+1
if i<n-1:
return true
else:
return false
def _ratio_x_(self,x,n):
r"""
Check if the ratio of the smallest and largest x vector entries is very small.
INPUT:
x -- vector
n -- length of the input vector
OUTPUT:
boolean value
NOTES:
See [B] p.3.
REFERENCES:
[B] D. H. Bailey:
"Finding New Mathematical Identities via Numerical Computations" ACM SIGNUM, vol. 33, no. 1 (Jan. 1998), pp. 17-22.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
xmax=max(abs(x[i]) for i in range(0,n,1))
xmin=min((abs(x[i]) for i in range(0,n,1)))
xratio=xmin/xmax
if xratio<10^(-20):
return true
else:
return false
def _zero_x_(self,x,n,d):
r"""
Test for zero entries in x.
INPUT:
x -- vector
n -- length of the input vector
d -- integer
OUTPUT:
integer
NOTES:
Set precision to $10^{-d}$. If x has a zero-entry return its index else return -1. See [FBA] p.9.
REFERENCES:
[FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno :
"Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
i=0
while i<=n-2:
if abs(x[i])<10^(-d):
break
else:
i=i+1
if i==n-1:
return -1
else:
return i
def _check_A_(self,A):
r"""
Test for too large entries in A.
INPUT:
A -- matrix
OUTPUT:
boolean value
NOTES:
See [FBA] p.18 .
REFERENCES:
[FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno :
"Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
n=A.nrows()
m=A.ncols()
Amax=0
for i in range (0,n,1):
z=max(abs(A[i,j]) for j in range(0,m,1))
if z>Amax:
Amax=z
if Amax>10^(10000):
return true
else:
return false
def _terminate_(self,H,A,B,x,n):
r"""
Check termination conditions.
INPUT:
H,A,B -- matrices
x -- vector
n -- length of the input vector
OUTPUT:
boolean value
NOTES:
See [FBA] p.9 and [B] p.3 .
REFERENCES:
[FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno :
"Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999.
[B] D. H. Bailey:
"Finding New Mathematical Identities via Numerical Computations" ACM SIGNUM, vol. 33, no. 1 (Jan. 1998), pp. 17-22.
AUTHORS:
- Agnes Jany (2009-05-10)
"""
if self._zero_H_(H,n)==true: #wieso zuerst das und dann erst zero x..
self.tex="\text{No relation found:}\quad H_{ii}=0."
self.reason="No relation found:","H[i,i]=0"
return true
else:
zerox=self._zero_x_(x,n,self.d)
if zerox!=-1:
self.tex=latex(B.column(zerox))
self.reason="Relation found:",B.column(zerox)
return true
else:
if self._ratio_x_(x,n):
zerox2=self._zero_x_(x,n,floor(self.d/4))
if zerox2!=-1:
self.tex=latex(B.column(zerox2))
self.reason="Relation found:",B.column(zerox2)
return true
else:
if self._check_A_(A):
self.tex="\text{No relation found:}\quad A_{ij}\quad \text{too large.}"
self.reason="No relation found:","A[i,j] too large"
return true
else:
return false
def _latex_(self):
r"""
Return \Latex representation of pslq output.
EXAMPLES:
sage: t=vector([11,27,31])
sage: a=pslq(t,digits=20)
Relation found:
(-1, 5, -4)
sage: latex(a)
\left(-1,5,-4\right)
"""
return self.tex