polylog(-v, z) wrong computation

[Jérôme Raulin <jerome@raulinfoissac.fr>]

Dear all,

polylog(-v, z) can be computed in different ways :
  sum_{k = 1}^{infinity} z^k k^v
Or
  z * An(z) / (1 - z)^(v + 1)
where An(z) is the nth Eulerian polynomial.

Comparing these formulae with Pari internal implementation (in default
precision on latest development version) :
> phi = (sqrt(5) - 1) / 2
%1 = 0.61803398874989484820458683436563811772
> polylog(-70, phi)
%2 = -8.286584156344620107 E117
  is obviously wrong as it should be > 0
> suminf(k = 1, phi^k * (k * 1.0)^70)
%3 = 4.2907022622606505026505656858632022776 E122
  is correct (the 1.0 is needed to force floating point computation and
avoid generating huge multiprecision integers)
> A(n) = sum(k = 1, n, k! * stirling(n, k, 2) * (x - 1)^(n - k));
> subst(A(70), x, phi) * phi / (1 - phi)^71
%5 = 4.2907022622606505026505656858632022811 E122
  is also correct.

For smaller v values, the result of polylog is correct.
Pari implementation uses the Eulerian formula (not with Stirling number but
with a polynomial recurrence).
A the end (1 - z)^(v + 1) is expanded symbolically to produce a rational
function. Then z is substituted by x.
The precision issue comes from the polynomial (1 - z)^(v + 1) which when
expanded has very large coefficients with alternate sign.
When evaluated with limited precision input massive cancellation happens.

> subst((1 - x)^71, x, phi)
%6 = -1.0912089477928565605 E-25
  is obviously wrong as it should be > 0
> (1 - phi)^71
%7 = 2.1074393480002467167065256550693471284 E-30
  is correct.

A first fix could be to do the polynomial evaluations before the expansion
of (1 - z)^(v + 1) and the polynomial division.
It should solve the precision issue.
Still when v gets larger we would start to face running time / complexity
issue as the Eulerian polynomial coefficients grow very large (the largest
coefficient of A70(z) is 100 digits long).

For larger v a better approach could be to fall back on the infinite sum
which still produce correct result very fast.

> suminf(k = 1, phi^k * (k * 1.0)^(10^6))
%8 = 8.7885110079460873612767697901006283040 E5883372
Which produces similar result to polylog approximation for negative v :
gamma(1 - v) * (-log(z))^(v - 1)
> gamma(1-(-10^6))*(-log(phi))^(-10^6-1)
%9 = 8.7885110079460873612767697901006026097 E5883372

I spotted the issue while solving the latest problem on a "famous
mathematical / computer programming problems" website which requires the
computation of log(polylog(-1234567, phi)).
Using the infinite summation formula, Pari produced the correct result in 5
seconds.

Cheers

Jerome