Re: Diophantine matters


Richard Fateman schrieb:

D Herring wrote:

Looking at the Mathematica notebook available at
http://library.wolfram.com/infocenter/MathSource/4263/
it appears that this algorithm is not terribly complicated.

It does require some basic linear algebra operations.

- Daniel

look at
http://www.cecm.sfu.ca/organics/papers/bailey/paper/html/node3.html
which may help


The Bailey-Plouffe paper is exactly what I needed; I find the PSLQ
algorithm beautifully, and probably unbeatably, simple. I have
implemented the Bailey-Plouffe version in Derive with the following
modifications:

1. For easier handling I have transposed the matrix B of potential
relation vectors, since matrices in Derive are stored as vectors of
rows.

2. The termination criterion in the paper is incorrect or incomplete:
not |y_min| < f*eps should be tested for, where y_min is the y element
of smallest absolute value and f is some constant fudge factor (the
Mathematica implementation does it this way too), but |y_min| <
|b_min|*f*eps, where |b_min| is the Euclidean norm of the relation
vector belonging to y_min. This corresponds to testing b.x < |b|*|x|*eps
with the input data vector x. (In fact one could form y anew after each
iteration without accumulating rounding errors as y = B.x/|x|.)

3. I have dispensed with the accumulation of the inverse A = B^(-1) of
the matrix B of potential relation vectors, as I found no use for this
matrix (the Mathematica implementation ignores it too): testing for
exhaustion of numeric precision as suggested in the paper again fails
for normalization reasons: a relation exhausting the given precision is
generally found well before max_(i,j)|a_(i,j)| exceeds 1/eps. So far I
haven't seen an example where such an addditional exit is needed to
terminate the iteration; can anybody come up with one?

The Mathematica implementation differs from the Bailey-Plouffe version
of the algorithm in that full (rather than block) matrix reduction of H
is employed during the iteration, and in that the termination criterion
involves the diagonal elements of H along with y (a constant fudge
factor as large as f = 10^5 is applied). Perhaps this corresponds to an
older version of the algorithm?

As not to waste resources, I am systematically truncating the iterated H
and y to the specified working precision plus some safety margin. Along
with the integer solution I am outputting the bound M = 1/max_j |H_j| on
the norm of any possible relation vector, as this is certainly useful. I
have assumed here that |H_j| in the paper denotes the Euclidian norm of
the matrix row H_j, else the bound will be incorrect.

This is my current Derive 6.10 code:

plsq_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(
y SUB j := 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^(2-PrecisionDigits),d:=PrecisionDigits+2,
n,a,b,y,h,i,j,m,s,t) := PROG(
IF(3 > 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),
plsq_aux(n,2,n-1,y,h,b,d),
LOOP(
m := FIRST(FIRST(ITERATE([IF(j_ SUB 2 < ABS(h SUB i_ SUB i_),
[i_,ABS(h SUB i_ SUB i_)*SQRT(3)/2],
REPLACE(j_ SUB 2*SQRT(3)/2,j_,2)),i_+1],
[j_,i_],[[1,0],1],n-1))),
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(m < n-1, PROG(
s := SIGN(h SUB m SUB j),
i := m,
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)))),
plsq_aux(n,m+1,m+1,y,h,b,d),
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(ABS(y SUB m) < eps*APPROX(ABS(b SUB m),d),
RETURN([b SUB m, APPROX(1/MAX(VECTOR(ABS(t_),t_,h)),d)]))))

[PrecisionDigits:=30, NotationDigits:=30, Notation:=Scientific]

pslq([1, 2.78838211157540371710574518532,
2.42559455216300418672182398066, 2.42928300559800284992107700363])

[[-380,-611,169,689], 305.226562795605888861851499927]

To really obtain the precision specified by the parameter
PrecisionDigits, it seems necessary to increase the parameter
NotationDigits as well, as shown. I haven't investigated this effect.
You may activate the Display instruction in order to watch the iteration
progress. Can the PSLQ code be made to look more elegant in Maxima?

Martin.
.

Quantcast