I had try following code:
from mpmath import *
mp.dps = 30
rf16 = lambda x,y: fsub(loggamma(x+y),loggamma(x))
a = 0.00016798648139376704
b = 8.082147903385792
print rf16(a,b)
print diff(rf16,(a,b),(0,0))
The outputs were:
-3.898517388468276439663575182e-17
7.35104632426086107405258486621e-17
However, the results should be equal for the two functions, i.e. rf16(a,b) = diff(rf16,(a,b),(0,0))
I don't know why the error happens.