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