(*
<< c:/job2007/PSLQ/pslq.m;
x0 =   7.97640656766789629481666285130155418456;
x[1] = N[Zeta[6], 200] ;
x[2] = N[ Zeta[3]^2, 200] ;
PSLQ[{x[1], x[2], x0}, 10];
PSLQ[{x[1], x[2], x0}, 20];
%/96;
(83*Pi^6)/10080 + Zeta[3]^2/24;
% - x0
*)
PSLQ[inx_List, prec_] := Block[{x, n = Length[inx], \[Gamma] = 2/Sqrt[3], A,
      B, H, D, Dinv, t, i, j, k, l, iter, \[Alpha], \[Beta], \[Lambda],
      \[Delta], r, R}, x = N[inx/Sqrt[inx . inx], prec];
      s = Sqrt[MapIndexed[Plus @@ Drop[x^2, First[#2] - 1] & , x]];
      A = B = IdentityMatrix[n];
      H = Table[Which[i > j, ((-x[[i]])*x[[j]])/(s[[j]]*s[[j + 1]]), i == j,
         s[[i + 1]]/s[[i]], i < j, 0], {i, 1, n}, {j, 1, n - 1}];
      t = HermiteReduce[H]; D = First[t]; Dinv = Inverse[D]; H = Last[t];
      x = x . Dinv; A = D . A; B = B . Dinv; For[iter = 0,
       iter < $IterationLimit, ++iter,
       r = MaxIndex[MapIndexed[\[Gamma]^First[#2]*Abs[#1] & , Tr[H, List]]];
        If[r < n - 1, \[Alpha] = H[[r,r]]; \[Beta] = H[[r + 1,r]];
          \[Lambda] = H[[r + 1,r + 1]]; \[Delta] =
           Sqrt[\[Beta]^2 + \[Lambda]^2]]; R = IdentityMatrix[n]; t = R[[r]];
        R[[r]] = R[[r + 1]]; R[[r + 1]] = t; x = x . R; H = R . H; A = R . A;
        B = B . R; If[r < n - 1, H = H . Table[Which[i == r && j == r,
             \[Beta]/\[Delta], i == r && j == r + 1, -\[Lambda]/\[Delta],
             i == r + 1 && j == r, \[Lambda]/\[Delta], i == r + 1 &&
              j == r + 1, \[Beta]/\[Delta], (i == j && j != r) ||
              (i == j && j != r + 1), 1, True, 0], {i, 1, n - 1},
            {j, 1, n - 1}]]; t = HermiteReduce[H]; D = First[t];
        Dinv = Inverse[D]; H = Last[t]; x = x . Dinv; A = D . A;
        B = B . Dinv; If[Min[Abs[Union[x, Tr[H, List]]]] <= 10^(-prec + 5),
         Break[]]]; Return[Transpose[B][[MaxIndex[-Abs[x]]]]]]

Recognize2[n_, po_, v_] := PSLQ[N[Table[n^i, {i, 0, po}], 100], 100] .
     Table[v^i, {i, 0, po}]

MaxIndex[x_List] := Block[{i = 1, j}, Do[If[x[[j]] > x[[i]], i = j],
       {j, 2, Length[x]}]; i]

HermiteReduce[H_] := Block[{n = Length[H], i, j, k, q, H2 = H,
      D = IdentityMatrix[Length[H]]},
     Do[q = Round[H[[i,j]]/H[[j,j]]]; Do[H2[[i,k]] -= q*H2[[j,k]],
         {k, 1, j}]; Do[D[[i,k]] -= q*D[[j,k]], {k, 1, n}]; , {i, 2, n},
       {j, i - 1, 1, -1}]; {D, H2}]
