function y=gpolylog(n,a)
%GPOLYLOG Generalized polylogarithm. 
%   Y=GPOLYLOG(N,A) computes in double-precision arithmetic
%   the generalized polylogarithm function S_{N,2)(A),
%   0 < A <= 1. 
%Reference: Eq.(2.9) in K.S. Koelbig, Nielsen's generalized 
%polylogarithms, SIAM J. Math. Anal. 17 (1986), 1232-1258.

eps0=.5*10^(-15);
err=1; nu=0; h=1; y=0;
while err>eps0
  y0=y;
  nu=nu+1;
  t=h*a^(nu+1)/(nu+1)^(n+1);
  y=y0+t;
  err=abs(t/y);
  h=h+1/(nu+1);
end
%nu
