Add transformations between harmonic numbers and polygamma functions
Test harmonic rewrites and expansions
Tests for rewriting of polygamma functions to harmonic numbers
Some more docs for polygamma functions
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)
Fix test failure about EulerGamma