Wrong values calculated by self-defined function

[German]

For a fit I need a fit function utilizing the digamma function. Unfortunately, gnuplot does not
know the digamma function, so I tried to implement my own definition of a digamma function using a (known) approximation:

digamma(x) = (x<7) ? digamma(x+1.0)-1.0/(x*1.0) : (x=x-0.5, xx=1.0/x, xx2=xx**2, xx4=xx2**2.0, log(x)+(1./24.)*xx2-(7.0/960.0)*xx4+(31.0/8064.0)*xx4*xx2-(127.0/30720.0)*xx4**2.0)

In python my self-defined function and the built-in function correlate nicely:

from math import log
from scipy.special import digamma

def digam(x):
    if x<7:
        return digam(x+1)-1.0/x
    x=x-0.5
    xx=1.0/x
    xx2=xx*xx
    xx4=xx2*xx2
    return log(x)+(1./24.)*xx2-(7.0/960.0)*xx4+(31.0/8064.0)*xx4*xx2-(127.0/30720.0)*xx4*xx4

But the values I retrieve from my self-defined gnuplot function are not very presice:

for example:

self-defined digamma(1) gnuplot: -0.503242515779268151
self-defined digam(1) python: -0.5772156649542977
built-in digamma pyton: -0.5772156649015329

My question: Is is possible to tell gnuplot to be more accurate in calculations?

Many thanks in advance
German

[Hans-Bernhard Broeker]

It's possible, and even obvious. The solution is that you've been trying too hard. Almost all of those temporary variables are just not needed. Just write down the formula as you most likely found it in your source (formatted for legibitliy here):

digamma(x) = (x<7) \
  ? digamma(x+1.0)-1.0/x \
  : (  x1=x-0.5 \
          , log(x1) \
                + (1./24.)/x1**2 \
                - (7.0/960.0)/x1**4 \
                + (31.0/8064.0)/x1**6 \
                - (127.0/30720.0)/x1**8 \
        )

The primary key to maintain sufficient precision is to avoid writing integer exponents as floating-point numbers. :-)