PSLQ revisited
clicl...@freenet.de
Recently I posted an implementation in Derive 6.10 of Ferguson's
integer-relation algorithm PSLQ. Meanwhile I noticed and corrected some
minor deficiencies in the code: My input checking inadvertantly excluded
vectors of length 2, and the termination test has to be moved to the
beginning of the main loop to preclude failures for some rational input
vectors. A few other changes are meant to inhance readability or
efficiency. The PSLQ parameter gamma remains at SQRT(4/3).
pslq_aux(n,i0,j0,y,h,b,d,i,j,k,t) := PROG(
i := i0,
LOOP(
j := MIN(i-1,j0),
LOOP(
IF(0 /= t:=ROUND(h SUB i SUB j/h SUB j SUB j),PROG(
APPROX(y SUB j :+ t*y SUB i,d),
k := [1,...,j],
h SUB i SUB k := APPROX(h SUB i SUB k-t*h SUB j SUB k,d),
b SUB j :+ t*b SUB i)),
IF(1 > j:=j-1,exit)),
IF(n < i:=i+1,exit)))
pslq(x,eps:=10^(1-PrecisionDigits),d:=PrecisionDigits+2,
n,a,b,y,h,i,j,m,s,t) := PROG(
IF(SOME(NOT(RATIONAL?(t_)) OR t_ = 0,t_,VECTOR(APPROX(t_),t_,x)) OR
2 > n:=DIM(x),RETURN("inapplicable")),
s := VECTOR(SQRT(t_),t_,FIRST(ITERATE(
[REPLACE(t_ SUB i_^2+t_ SUB (i_+1),t_,i_),i_-1],[t_,i_],
[REPLACE(x SUB n^2,x,n),n-1],n-1))),
y := APPROX(x/FIRST(s),d),
h := VECTOR(VECTOR(IF(j_ < i_,
APPROX(-x SUB i_*x SUB j_/(s SUB j_*s SUB (j_+1)),d),
IF(j_ = i_,APPROX(s SUB (j_+1)/s SUB j_,d),0)),j_,1,n-1),i_,1,n),
b := IDENTITY_MATRIX(n),
pslq_aux(n,2,n-1,y,h,b,d),
LOOP(
m := FIRST(ITERATE([IF(ABS(y SUB j_) > ABS(y SUB i_),i_,j_),
i_+1],[j_,i_],[1,2],n-1)),
"DISPLAY([b SUB m,ABS(y SUB m)])",
IF(y SUB m^2 < eps^2*b SUB m^2,RETURN([b SUB m,
APPROX(ABS(y SUB m)/MAX(VECTOR(ABS(t_),t_,y)),d),
APPROX(SQRT(1/MAX(VECTOR(t_^2,t_,h))),d)])),
m := FIRST(ITERATE([IF(h SUB j_ SUB j_^2 < (4/3)^(i_-j_)*
h SUB i_ SUB i_^2,i_,j_),i_+1],[j_,i_],[1,2],n-2)),
j := [m,m+1],
i := REVERSE(j),
y SUB i := y SUB j,
h SUB i := h SUB j,
b SUB i := b SUB j,
IF(n-1 > i:=m,PROG(
s := SIGN(h SUB m SUB j),
LOOP(
t := h SUB i SUB j,
h SUB i SUB j := APPROX([s*t,CROSS(s,t)],d),
IF(n < i:=i+1,exit)))),
pslq_aux(n,m+1,m+1,y,h,b,d)))
The first two of the following examples would have failed previously:
[PrecisionDigits:=30,NotationDigits:=30,Notation:=Scientific]
pslq([1,2])
[[2,-1],0,2.23606797749978969640917366873]
pslq([1,1,1])
[[1,-1,0],0,1.22474487139158904909864203735]
pslq([pi,#e,LN(2),euler_gamma])
[[9.439541*10^6,-1.0360028*10^7,-3.287929*10^6,1.3605*10^6],
0.115344673313855049732308815137,
8.23224397166993621056247068053*10^6]
The first number after the relation vector is the "confidence level",
which should be below 10^-10 or so, the second is a bound on the
Euclidean norm of any possible relation vector, as before. Thus each of
the first pair of relations is highly significant, while a confidence
level of order unity reveals that there is no relation between pi, #e,
LN(2), and euler_gamma (its Euclidean norm would have to be > 8.2*10^6).
Here is the same thing for complex input vectors. This works for real
vectors too, but does not necessarily give the same answer as the real
version, because Ferguson's parameter gamma is set to SQRT(2) here.
cround(x) := ROUND(RE(x))+#i*ROUND(IM(x))
cpslq_aux(n,i0,j0,y,h,b,d,i,j,k,t) := PROG(
i := i0,
LOOP(
j := MIN(i-1,j0),
LOOP(
IF(0 /= t:=cround(h SUB i SUB j/h SUB j SUB j),PROG(
APPROX(y SUB j :+ t*y SUB i,d),
k := [1,...,j],
h SUB i SUB k := APPROX(h SUB i SUB k-t*h SUB j SUB k,d),
b SUB j :+ t*b SUB i)),
IF(1 > j:=j-1,exit)),
IF(n < i:=i+1,exit)))
cabs2(x) := RE(x)^2+IM(x)^2
cpslq(x,eps:=10^(1-PrecisionDigits),d:=PrecisionDigits+2,
n,a,b,y,h,i,j,m,s,t) := PROG(
IF(SOME(NOT(NUMBER?(t_)) OR t_ = 0,t_,VECTOR(APPROX(t_),t_,x)) OR
2 > n:=DIM(x),RETURN("inapplicable")),
s := VECTOR(SQRT(t_),t_,FIRST(ITERATE(
[REPLACE(cabs2(t_ SUB i_)+t_ SUB (i_+1),t_,i_),i_-1],[t_,i_],
[REPLACE(cabs2(x SUB n),x,n),n-1],n-1))),
y := APPROX(x/FIRST(s),d),
h := VECTOR(VECTOR(IF(j_ < i_,
APPROX(-CONJ(x SUB i_)*x SUB j_/(s SUB j_*s SUB (j_+1)),d),
IF(j_ = i_,APPROX(s SUB (j_+1)/s SUB j_,d),0)),j_,1,n-1),i_,1,n),
b := IDENTITY_MATRIX(n),
cpslq_aux(n,2,n-1,y,h,b,d),
LOOP(
m := FIRST(ITERATE([IF(cabs2(y SUB j_) > cabs2(y SUB i_),i_,j_),
i_+1],[j_,i_],[1,2],n-1)),
"DISPLAY([b SUB m,ABS(y SUB m)])",
IF(cabs2(y SUB m) < eps^2*cabs2(b SUB m),RETURN([b SUB m,
APPROX(SQRT(cabs2(y SUB m)/MAX(VECTOR(cabs2(t_),t_,y))),d),
APPROX(SQRT(1/MAX(VECTOR(cabs2(t_),t_,h))),d)])),
m := FIRST(ITERATE([IF(cabs2(h SUB j_ SUB j_) < 2^(i_-j_)*
cabs2(h SUB i_ SUB i_),i_,j_),i_+1],[j_,i_],[1,2],n-2)),
j := [m,m+1],
i := REVERSE(j),
y SUB i := y SUB j,
h SUB i := h SUB j,
b SUB i := b SUB j,
IF(n-1 > i:=m,PROG(
s := SIGN(h SUB m SUB j),
LOOP(
t := h SUB i SUB j,
h SUB i SUB j := APPROX([CONJ(s)*t,CROSS(s,t)],d),
IF(n < i:=i+1,exit)))),
cpslq_aux(n,m+1,m+1,y,h,b,d)))
Again a simple example:
[PrecisionDigits:=30,NotationDigits:=30,Notation:=Scientific]
cpslq([1,1+#i,1-#i])
[[0,1,-#i],0,1.11803398874989484820458683436]
The real and complex PSLQ algorithms can be put to use for the rational
factorization of univariate polynomials; in fact, the speed of this
procedure compares favorably with the built-in algorithm of Derive 6.10
for real polynomials. Here is a complex factorization example, taken
from the recent thread "factoring over the complex rationals":
[PrecisionDigits:=30,NotationDigits:=30,Notation:=Scientific]
SORT(NSOLUTIONS(29*x^6+2*x^5-10*x^4-60*x^3+10*x^2+24*x+37,x))
[-0.616855149136108389349551112943-
0.960173605641725364006512446171*#i,
-0.616855149136108389349551112943+
0.960173605641725364006512446171*#i,
-0.516109317059961663101235546354-
0.694169983821462170219007915817*#i,
-0.516109317059961663101235546354+
0.694169983821462170219007915817*#i,
1.09848170757538039727837286619-
0.320203274731460944143529952405*#i,
1.09848170757538039727837286619+
0.320203274731460944143529952405*#i]
croot_poly(r,x,m,d:=PrecisionDigits+2,p) := PROG(
p := cpslq(ITERATES(APPROX(t_*r,d),t_,1,m)),
LOOP(IF(FIRST(p SUB 1) = 0,p SUB 1 := REST(p SUB 1),exit)),
RETURN([p SUB 1*VECTOR(x^i_,i_,0,DIM(p SUB 1)-1),p SUB 2,p SUB 3]))
croot_poly(-0.616855149136108389349551112943-
0.960173605641725364006512446171*#i,x,5)
[-2*x^3-3*x^2+1+#i*(5*x^3-x^2-2*x-6),
1.59570517785194071425804385836*10^(-25),
2.28314364426950302068647354558]
This works by taking one zero of the given polynomial and constructing a
smaller polynomial with integer coefficients from it. Note that the root
must be approximated accurately up to the numeric precision employed in
PSLQ, which in turn must be large enough not to miss any factor. In
general, though not in the present case, the result may contain spurious
factors absent in the original polynomial; these have to be removed by a
subsequent GCD calculation.
How does the efficiency of PSLQ-based polynomial factorization compare
with the native factorization algorithms on the other systems?
Martin.
Daniel Lichtblau
Depends on which PSLQ factorization method you have in mind, and what
class of polynomials.
Generally speaking, the knapsack method by van Hoeij (with
improvements by Novocin and perhaps others) is, I believe, best in
terms of optimality. The method you show above was first discussed in
the original LLL paper (though using LLL reduction rather than PSLQ).
It is better than the recombination method for those polynomials that
are "bad", i.e. have many factors modulo any prime, but far fewer
factors over the rationals.
Daniel Lichtblau
Wolfram Research