1
2
3/*********************************************************************************************#
4#                       The PSLQ Integer Relation Algorithm                                   #
5#                                                                                             #
6# Aut.: Helaman R.P. Ferguson and David Bailey "A Polynomial Time, Numerically Stable         #
7#       Integer Relation Algorithm" (RNR Technical Report RNR-92-032)    helaman@super.org    #
8# Ref.: David Bailey and Simon Plouffe "Recognizing Numerical Constants" dbailey@nas.nasa.gov #
9# Cod.: Raymond Manzoni  raymman@club-internet.fr                                             #
10#*********************************************************************************************#
11# Creation:97/11                    #
12# New termination criteria:97/12/15 #
13# this code is free...              #
14
15Ported to Yacas 2000 Ayal Pinkus.
16
17Given a list of constants x find coefficients sol[i] such that
18      sum(sol[i]*x[i], i=1..n) = 0    (where n=Length(x))
19
20    x is the list of real expressions
21          N(x[i]) must evaluate to floating point numbers!
22    precision is the number of digits needed for completion;
23          must be greater or equal to log10(max(sol[i]))*n
24    returns the list of solutions with initial precision
25          and the confidence (the lower the better)
26
27    Example:
28
29    In> Pslq({2*Pi-4*Exp(1),Pi,Exp(1)},20)
30    Out> {1,-2,4};
31
32*/
33
34Pslq(x, precision) :=
35[
36  Local (ndigits, gam, A, B, H, n, i, j, k, s, y, tmp, t, m, maxi, gami,
37	 t0, t1, t2, t3, t4, mini, Confidence, norme,result);
38  n:=Length(x);
39  ndigits:=Builtin'Precision'Get();
40  Builtin'Precision'Set(precision+10); // 10 is chosen arbitrarily, but should always be enough. Perhaps we can optimize by lowering this number
41  Confidence:=10^(-MathFloor(N(Eval(precision/3))));
42//Echo("Confidence is ",Confidence);
43
44  gam:=N(Sqrt(4/3));
45  For (i:=1, i<=n,i++) x[i]:=N(Eval(x[i]));
46
47//Echo("1...");
48
49  A:=Identity(n); /*A and B are of Integer type*/
50  B:=Identity(n); /*but this doesn't speed up*/
51  s:=ZeroVector(n);
52  y:=ZeroVector(n);
53
54//Echo("2...");
55
56  For(k:=1,k<=n,k++)
57  [
58    tmp:=0;
59    For (j:=k,j<=n,j++) tmp:=tmp + N(x[j]^2);
60//tmp:=MathDivide(tmp,1.0);
61//Echo("tmp is ",tmp);
62//MathDebugInfo(tmp);
63/*If(Not IsPositiveNumber(tmp),
64  Echo("******** not a positive number: ",tmp)
65);
66If(Not IsNumber(tmp),
67  Echo("******** not a number: ",tmp)
68);
69If(LessThan(tmp,0),
70[
71  Echo("******** not positive: ",tmp);
72]
73);*/
74
75    s[k]:=MathSqrt(tmp);
76
77
78/*If(Not IsNumber(tmp),
79[
80Echo("************** tmp = ",tmp);
81]);
82If(Not IsNumber(s[k]),
83[
84Echo("************** s[k] = ",s[k]);
85]);*/
86
87  ];
88
89//Echo("3...");
90
91  tmp:=N(Eval(s[1]));
92/*If(Not IsNumber(tmp),
93[
94Echo("************** tmp = ",tmp);
95]);*/
96
97  For (k:= 1,k<= n,k++)
98  [
99    y[k]:=N(Eval(x[k]/tmp));
100    s[k]:=N(Eval(s[k]/tmp));
101
102//Echo("1..."," ",y[k]," ",s[k]);
103/*If(Not IsNumber(y[k]),
104[
105Echo("************** y[k] = ",y[k]);
106]);
107If(Not IsNumber(s[k]),
108[
109Echo("************** s[k] = ",s[k]);
110]);*/
111
112  ];
113  H:=ZeroMatrix(n, n-1);
114
115//Echo("4...",n);
116  For (i:=1,i<= n,i++)
117  [
118
119    if (i <= n-1)  [ H[i][i]:=N(s[i + 1]/s[i]); ];
120
121//Echo("4.1...");
122    For (j:= 1,j<=i-1,j++)
123    [
124//Echo("4.2...");
125      H[i][j]:= N(-(y[i]*y[j])/(s[j]*s[j + 1]));
126//Echo("4.3...");
127
128/*If(Not IsNumber(H[i][j]),
129[
130Echo("************** H[i][j] = ",H[i][j]);
131]
132);*/
133
134    ];
135  ];
136
137//Echo("5...");
138
139  For (i:=2,i<=n,i++)
140  [
141    For (j:=i-1,j>= 1,j--)
142    [
143//Echo("5.1...");
144      t:=Round(H[i][j]/H[j][j]);
145//Echo("5.2...");
146      y[j]:=y[j] + t*y[i];
147//Echo("2..."," ",y[j]);
148      For (k:=1,k<=j,k++) [ H[i][k]:=H[i][k]-t*H[j][k]; ];
149      For (k:=1,k<=n,k++)
150      [
151        A[i][k]:=A[i][k]-t*A[j][k];
152        B[k][j]:=B[k][j] + t*B[k][i];
153      ];
154    ];
155  ];
156  Local(found);
157  found:=False;
158
159//Echo("Enter loop");
160
161  While (Not(found))
162  [
163    m:=1;
164//Echo("maxi 1...",maxi);
165    maxi:=N(gam*Abs(H[1][1]));
166//Echo("maxi 2...",maxi);
167    gami:=gam;
168//Echo("3...");
169    For (i:= 2,i<= n-1,i++)
170    [
171      gami:=gami*gam;
172      tmp:=N(gami*Abs(H[i][i]));
173      if (maxi < tmp)
174      [
175        maxi:=tmp;
176//Echo("maxi 3...",maxi);
177        m:=i;
178      ];
179    ];
180//Echo("4...",maxi);
181    tmp:=y[m + 1];
182    y[m + 1]:=y[m];
183    y[m]:=tmp;
184//Echo("3..."," ",y[m]);
185//Echo("5...");
186    For (i:= 1,i<=n,i++)
187    [
188      tmp:=A[m + 1][ i];
189      A[m + 1][ i]:=A[m][ i];
190      A[m][ i]:=tmp;
191      tmp:=B[i][ m + 1];
192      B[i][ m + 1]:=B[i][ m];
193      B[i][ m]:=tmp;
194    ];
195    For (i:=1,i<=n-1,i++)
196    [
197      tmp:=H[m + 1][ i];
198
199      H[m + 1][ i]:=H[m][ i];
200      H[m][ i]:=tmp;
201    ];
202//Echo("7...");
203    if (m < n-1)
204    [
205      t0:=N(Eval(Sqrt(H[m][ m]^2 + H[m][ m + 1]^2)));
206
207      t1:=H[m][ m]/t0;
208      t2:=H[m][ m + 1]/t0;
209
210//      If(IsZero(t0),t0:=N(Confidence));
211//Echo("");
212//Echo("H[m][ m] = ",N(H[m][ m]));
213//Echo("H[m][ m+1] = ",N(H[m][ m+1]));
214
215//If(IsZero(t0),[t1:=Infinity;t2:=Infinity;]);
216//Echo("t0=",N(t0));
217//Echo("t1=",N(t1));
218//Echo("t2=",N(t2));
219
220      For (i:=m,i<=n,i++)
221      [
222        t3:=H[i][ m];
223        t4:=H[i][ m + 1];
224//Echo("    t1 = ",t1);
225//Echo("    t2 = ",t2);
226//Echo("    t3 = ",t3);
227//Echo("    t4 = ",t4);
228        H[i][ m]:=t1*t3 + t2*t4;
229//Echo("7.1... ",H[i][ m]);
230        H[i][ m + 1]:= -t2*t3 + t1*t4;
231//Echo("7.2... ",H[i][ m+1]);
232      ];
233    ];
234//Echo("8...");
235    For (i:= 1,i<= n,i++)
236    [
237      For (j := Min(i-1, m + 1),j>= 1,j--)
238      [
239        t:=Round(H[i][ j]/H[j][ j]);
240//Echo("MATRIX",H[i][ j]," ",H[j][ j]);
241//Echo("5... before"," ",y[j]," ",t," ",y[i]);
242        y[j]:=y[j] + t*y[i];
243//Echo("5... after"," ",y[j]);
244        For (k:=1,k<=j,k++) H[i][ k]:=H[i][ k]-t*H[j][ k];
245        For (k:= 1,k<=n,k++)
246        [
247          A[i][ k]:=A[i][ k]-t*A[j][ k];
248          B[k][ j]:=B[k][ j] + t*B[k][ i];
249        ];
250      ];
251    ];
252//Echo("9...",N(H[1],10));
253
254    /* Builtin'Precision'Set(10);*/ /*low precision*/
255//    maxi := N(H[1] .  H[1],10);
256    maxi := N(H[1] .  H[1]);
257//Echo("H[1] = ",H[1]);
258//Echo("N(H[1]) = ",N(H[1]));
259//Echo("N(H[1] . H[1]) = ",N(H[1] . H[1]));
260//Echo("maxi 4...",maxi);
261
262//Echo("9... maxi = ",maxi);
263
264    For (j:=2,j<=n,j++)
265    [
266//Echo("9.1...");
267      tmp:=N(H[j] . H[j],10);
268//Echo("9.2...");
269      if (maxi < tmp) [ maxi:=tmp; ];
270//Echo("maxi 5...",maxi);
271//Echo("9.3...");
272    ];
273//Echo("10...");
274    norme:=N(Eval(1/Sqrt(maxi)));
275    m:=1;
276    mini:=N(Eval(Abs(y[1])));
277//Echo("y[1] = ",y[1]," mini = ",mini);
278    maxi:=mini;
279
280//Echo("maxi 6...",maxi);
281//Echo("11...");
282    For (j:=2,j<=n,j++)
283    [
284      tmp:=N(Eval(Abs(y[j])));
285      if (tmp < mini)
286      [
287        mini:=tmp;
288        m:=j;
289      ];
290      if (tmp > maxi) [ maxi:=tmp; ];
291//Echo("maxi 7...",maxi);
292    ];
293    /* following line may be commented */
294//Echo({"Norm bound:",norme," Min=",mini," Conf=",mini/maxi," required ",Confidence});
295    if ((mini/maxi) < Confidence) /*prefered to : if mini < 10^(- precision) then*/
296    [
297    /* following line may be commented */
298/*      Echo({"Found with Confidence ",mini/maxi}); */
299      Builtin'Precision'Set(ndigits);
300      result:=Transpose(B)[m];
301      found:=True;
302    ]
303    else
304    [
305      maxi:=Abs(A[1][ 1]);
306      For (i:=1,i<=n,i++)
307      [
308//Echo("i = ",i," n = ",n);
309        For (j:=1,j<=n,j++)
310        [
311//Echo("j = ",j," n = ",n);
312          tmp:=Abs(A[i][ j]);
313          if (maxi < tmp) [ maxi:=tmp;];
314        ];
315      ];
316//Echo("maxi = ",maxi);
317      if (maxi > 10^(precision))
318      [
319        Builtin'Precision'Set(ndigits);
320        result:=Fail;
321        found:=True;
322      ];
323      Builtin'Precision'Set(precision+2);
324//Echo("CLOSE");
325    ];
326  ];
327  result;
328];
329
330/* end of file */
331
332
333
334
335