kshitij12345 raised an issue two years ago concerning the fact that pytorch's implementation of the 1st order polygamma function incorrectly produces finite values for negative integers. The issue was closed, but the issue still seems to be present in the latest nightlies and the latest release of pytorch.

n = 1

>>> t = torch.tensor([-1.])
>>> torch.polygamma(1, t)
tensor([1.2914e+15])
>>> scipy.special.polygamma(1, t.numpy())
array([inf], dtype=float32)

>>> t = torch.tensor([-501.])
>>> torch.polygamma(1, t)
tensor([2.0831e+09])
>>> scipy.special.polygamma(1, t.numpy())
array([inf], dtype=float32)

>>> t = torch.tensor([-float('inf')])
>>> torch.polygamma(1, t)
tensor([nan])
>>> scipy.special.polygamma(1, t.numpy())
array([inf], dtype=float32)
>>> 

n > 1

>>> t = torch.tensor([float('inf')])
>>> torch.polygamma(2, t)
tensor([nan])
>>> scipy.special.polygamma(2, t.numpy())
array([-0.], dtype=float32)
>>> torch.polygamma(3, t)
tensor([nan])
>>> scipy.special.polygamma(3, t.numpy())
array([0.], dtype=float32)
>>> torch.polygamma(4, t)
tensor([nan])
>>> scipy.special.polygamma(4, t.numpy())
array([-0.], dtype=float32)

cc @mruberry @rgommers @kshitij12345

_Originally posted in #55357