// OEIS A259281
//
// "Decimal expansion of Sum'_{(x,y,z)=-infinity..infinity} 1/(x^2+y^2+z^2)^2,
//  where the 'prime' indicates that the term x=y=z=0 is to be left out."

DP:=90; // number of digits of precision for floating-point calculations
SetDefaultRealField(RealField(DP));
DPdisplay:=50; // number of digits of precision to display in output
M:=50; // number of convergence passes
nOutputValues:=18;

// The above settings will generate output correct to 50 significant digits;
//  if more digits are desired, increase DP, DPdisplay, and M.

nMax:=M+nOutputValues;
S:=0.0; // initialize the sum
u:=[];
for n in [1..nMax] do
   // Update S by adding the terms for all the points on the surface of the
   //  cube defined by |x| <= n, |y| <= n, |z| <= n:
   n2:=n*n;
   n4:=n2*n2;
   S+:=6.0/n4; // center of each of 6 faces of cube
   S+:=3.0/n4; // 4 corners in xy-, xz-, and yz-planes
   S+:=8.0/(9*n4); // 8 corners of cube
   for x in [1..n-1] do
      x2:=x*x;
      s2:=x2+n2;
      S+:=24.0/(s2*s2); // one coordinate is 0, one is n, other is in [1..n-1]
      s2:=2*x2+n2;
      S+:=24.0/(s2*s2); // two coordinates are equal and in [1..n-1]
      s2:=x2+2*n2;
      S+:=24.0/(s2*s2); // two coordinates both equal n; other is in [1..n-1]
   end for;
   for x in [1..n-2] do
      x2plusn2:=x*x+n2;
      for y in [x+1..n-1] do
         s2:=x2plusn2+y*y;
         S+:=48.0/(s2*s2); // one coordinate is equal to n;
                           //  the other two are distinct and in [1..n-1]
      end for;
   end for;
   u[#u+1]:=S*n^M;
end for;

// The code below has the effect of taking the M-th differences of the values
//  (S(n)*n^M) and dividing by M!; the resulting real-valued sequence will
//  rapidly approach the limiting value of S (and will do so more rapidly at
//  larger values of M).

for j in [1..M] do
   for i in [1..#u-j] do
      u[i]:=(u[i+1]-u[i])/j;
      if j eq M then
         ChangePrecision(u[i],DPdisplay);
      end if;
   end for;
end for;