Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2007
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

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






  • Prev by Date: Re: StringCases kills Kernel
  • Next by Date: Re: is there a better way to do constraint logic
  • Previous by thread: Re: PSLQ
  • Next by thread: Re: PSLQ