Fix rewrite of polygamma in terms of zeta functions
Make harmonic numbers tractable for Gruntz algorithm.

This adds nothing really new but consistency.

In [11]: gruntz(harmonic(n), n, oo)
Out[11]: oo

In [12]: gruntz(harmonic(n, 2), n, oo)
Out[12]: pi**2/6

In [13]: gruntz(harmonic(n, 3), n, oo)
Out[13]: -polygamma(2, 1)/2

In [14]: gruntz(harmonic(n, 4), n, oo)
Out[14]: pi**4/90

In [15]: gruntz(harmonic(n, 5), n, oo)
Out[15]: -polygamma(4, 1)/24

In [16]: gruntz(harmonic(n, 6), n, oo)
Out[16]: pi**6/945

We can already do this by direct computation:

In [40]: harmonic(oo)
Out[40]: zoo

In [41]: harmonic(oo, 2)
Out[41]: pi**2/6

In [42]: harmonic(oo, 3)
Out[42]: zeta(3)

In [43]: harmonic(oo, 4)
Out[43]: pi**4/90

In [44]: harmonic(oo, 5)
Out[44]: zeta(5)

In [45]: harmonic(oo, 6)
Out[45]: pi**6/945

Note how the polygamma correspond to the zeta functions:

In [14]: gruntz(harmonic(n, 7), n, oo)
Out[14]: -polygamma(6, 1)/720

In [15]: gruntz(harmonic(n, 7), n, oo).rewrite(zeta)
Out[15]: zeta(7)

In [16]: harmonic(oo, 7)
Out[16]: zeta(7)