Inverse Harmonic and Inverse Digamma Functions

[Albert Chan]

Something unexpected when I try to test invPsi() accuracy.
Note: Psi(1) = -euler_gamma = -0.57721566490153286 ...

XCas> 1 - invPsi(-0.57721566490153286)       → 5.55111512313e-16
XCas> 1 - invPsi(-euler_gamma)                     → 1.713074127e-13 ???

Above were done on XCas 1.5.0-63 (win32), with 12-digits display float format.
Raising to 17-digits, expressions now agreed, error = 5.5511151231258e-16

Rounding constant with changing float display format is likely a XCas bug.

This bug seems to affect only euler_gamma.
Example, setting float display format to 6 digits:

XCas> 3.1415926535897931 - float(pi)                          → 0.0
XCas> 0.57721566490153286 - float(euler_gamma)       → 0.470199e-6 ???