Re: PSLQ
- To: mathgroup at smc.vnet.net
- Subject: [mg80819] Re: [mg80785] PSLQ
- From: danl at wolfram.com
- Date: Mon, 3 Sep 2007 06:11:03 -0400 (EDT)
- References: <200709010435.AAA26690@smc.twtelecom.net>
> I am trying to use Peter Bertok's PSLQ algorithm at > > http://library.wolfram.com/infocenter/MathSource/4263/ > > to check for linear relations between logarithms, but it is not > working properly. > > If I do > > PSLQ[N[{Log[3/2], -Log[2], -Log[4/3]}, 200], 200] > > , the program correctly outputs > {1, 1, -1}. > > However, if I try > > PSLQ[N[{Log[3/2], -Log[5], -Log[4/3]}, 200], 200] > > , instead of outputting a set of three large integers which would give > a linear relation equalling zero correct to 10^(-200) > (which I could interpret to mean there was no true linear > relationship, because of the size of the integers), > it outputs a string of error messages. > > Can anyone tell me how to tweak this program to eliminate these > errors? > > Failing that, is there another Mathematica implementation of PSLQ > which would do the same thing? > > I have searched the archive and there was a link to a Mathematica > program by David Bailey, but the link no longer works. > > Thanks for any help with this, > > Jimmy Mc Laughlin > > PS Please reply also to my e-mail address > Clashton at gmail.com I suspect that Bertok code is doing roughly as it should. In general PSLQ just "runs out of precision", so to speak, and has to give up, when there is no reasonable candidate relation. If you in fact want to see huge multipliers, you can get that effect using lattice reduction methods. But I suspect you just want a graceful failure mode. Around 9-10 years ago I wrote a PSLQ implementation that I think is a bit, well, not so good. I give the code below, with the caveat that it is "as is". It might or might not find your desired integer relations. Don't expect it to find you amorous relations, or wealthy distant relations. I simply offer it as an alternative to what you already have, on the off chance that it might be of use. Feel free to improve on it in any way you like (and if it starts finding wealthy distant relations, I'd like to hear about that). I no longer recall how this stuff works. My only comment on the actual code is that the value 'const' can be made into a parameter for this algorithm; it is equivalent to lambda^2 from PSLQ literature. badNum[n_? NumberQ] := (n == 0 && Accuracy[n] < 3) || (n != 0 && Precision[n] < 4) vectorNumberQ[vec_?VectorQ] /; Length[vec] != 0 := Apply[Or, Map[NumberQ, vec]] vectorNumberQ[_] := False innerProduct[v1_, v2_, False] := v1.v2 innerProduct[v1_, v2_, True] := v1.Conjugate[v2] normSquare[vec_, cmplx_: False] := innerProduct[vec, vec, cmplx] norm[vec_, cmplx_: False] := Sqrt[normSquare[vec, cmplx]] reduceH[ohh_, len_, start_, mat_, inv_, iter_] := Module[{j, k, r, hh = ohh, newmat = mat, newinv = inv, bad = False}, For[j = start, j <= len, j++, For[k = If[iter == 0, j - 1, Min[j - 1, start]], k > 0, k--, r = If[hh[[k, k]] != 0, Round[hh[[j, k]]/hh[[k, k]]], 0]; If[r != 0, For[m = 1, m <= k, m++, hh[[j, m]] -= r*hh[[k, m]]; If[hh[[j, m]] == 0, hh[[j, m]] = 0]; bad = bad || badNum[hh[[j, m]]];]; newmat[[j]] -= r*newmat[[k]]; newinv[[k]] += r*newinv[[j]];];];]; {newmat, newinv, hh, bad}] Corner[ohh_, j_, len_] := Module[{hh = ohh, a = ohh[[j + 1, j]], b = ohh[[j, j]], c = ohh[[j, j + 1]], bstar, cstar, d, k, m1, m2, m3, m4, tmp}, bstar = Conjugate[b]; cstar = Conjugate[c]; d = Sqrt[b*bstar + c*cstar]; m1 = bstar/d; m2 = -c/d; m3 = cstar/d; m4 = b/d; hh[[{j, j + 1}, {j, j + 1}]] = {{d, 0}, {a*m1, a*m2}}; For[k = j + 2, k <= len, k++, tmp = hh[[k, j]]; hh[[k, j]] = hh[[k, j]]*m1 + hh[[k, j + 1]]*m3; hh[[k, j + 1]] = tmp*m2 + hh[[k, j + 1]]*m4;]; hh] PSLQ[ovec_?VectorQ] /; vectorNumberQ[ovec] := Module[{cmplx = ! FreeQ[ovec, Complex], len = Length[ovec], mat, hh, savevec, inv, s, prec, j, k, eps, tau, best, bestj = 1, value, usevec, vecs, vec, bad, const, qq, rr, iter = 0, orighh, newmat, oldmat, oldinv}, prec = Internal`EffectivePrecision[ovec]; savevec = SetPrecision[ovec, prec]; savevec = savevec/norm[savevec, cmplx]; eps = 10^(-1.5*prec/len); mat = IdentityMatrix[len]; inv = IdentityMatrix[len]; {qq, rr} = QRDecomposition[Transpose[Prepend[mat, savevec]]]; hh = SetPrecision[Drop[Transpose[Drop[rr, 1]], 1], prec]; orighh = hh; const = If[cmplx, 4, 2]; While[True, {mat, inv, hh, bad} = reduceH[hh, len, bestj + 1, mat, inv, iter]; iter++; For[s = 1, s < len, s++, bad = bad || badNum[hh[[s, s]]]]; If[bad, bad = False; newmat = mat.orighh; If[Internal`EffectivePrecision[newmat] <= $MachinePrecisi on, newmat = SetPrecision[newmat, 2*$MachinePrecisio n];]; {qq, rr} = QRDecomposition[Transpose[newmat]]; hh = Transpose[rr]; bad = Apply[Or, Map[badNum, Flatten[hh]]]; If[bad, Return[$Failed]]; If[Dimensions[hh] =!= {len, len - 1}, Return[$Failed]]; ]; For[s = 1, s < len, s++, If[hh[[s, s]] == 0, Break[]]]; If[s < len, Return[Last[inv]]]; For[j = 1, j < len, j++, If[Abs[hh[[j, j]]] > eps, Break[]]]; If[j == len, Return[$F ailed]]; best = 0; tau = 1; For[j = 1, j < s, j++, value = tau*Abs[hh[[j, j]]]; If[value > best, best = value; bestj = j]; tau *= const;]; mat[[{bestj, bestj + 1}]] = mat[[{bestj + 1, bestj}]]; inv[[{bestj, bestj + 1}]] = inv[[{bestj + 1, bestj}]]; hh[[{bestj, bestj + 1}]] = hh[[{bestj + 1, bestj}]]; If[bestj < len - 1, hh = Corner[hh, bestj, len]]; ]; ] Your examples: In[159]:= Timing[PSLQ[N[{Log[3/2], -Log[2], -Log[4/3]}, 200]]] Out[159]= {0.01, {1, 1, -1}} In[161]:= Timing[PSLQ[N[{Log[3/2], -Log[5], -Log[4/3]}, 200]]] Out[161]= {0.431, $Failed} Daniel Lichtblau Wolfram Research