Radovan Omorjan pointed me to this page of the „Kitchin Group“. It is all about numerical stuff done with Python. One trick is very interesting. For numerical differentiation, we can use the simple differential quotient as an estimate. Using the Taylor formula, we can get the error as follows.

\(\dfrac{f(x+h)-f(x)}{h}  = f'(x) + \dfrac{h}{2} f'(\xi)\)

So it is in the order of |h|. The problem, however, is the evaluation of f(x+h)-f(x). If f(x) is not exactly zero, the two numbers have many equal digits and the subtraction cancels these digits leaving much less correct digits. Here is an example computed in EMT.

>format(20,16); ...
>h=1e-8; exp(1+h), exp(1), exp(1+h)-exp(1)
  2.7182818556418633 
  2.7182818284590451 
  0.0000000271828182 
>(exp(1+h)-exp(1))/h, E,
  2.7182818218562943 
  2.7182818284590451

In the result, only 8 digits are correct. If we take a smaller h, then the result is much worse.

>h=1e-14; exp(1+h), exp(1), exp(1+h)-exp(1)
  2.7182818284590722 
  2.7182818284590451 
  0.0000000000000271 
>(exp(1+h)-exp(1))/h, E,
  2.7089441800853820 
  2.7182818284590451

The  optimal h must keep a subtle balance between the accuracy of the formula and the accuracy of the subtraction.

The difference formula can be improved by taking two other points, or even taken more points. The Taylor expansion shows

\(\dfrac{f(x+h)-f(x-h)}{2h} = f'(x) + \dfrac{h^2}{6} f“(\xi)\)

The optimal h for this is the third root of the machine epsilon.

>h=1e-5; (exp(1+h)-exp(1-h))/(2*h), E,
  2.7182818285176320 
  2.7182818284590451

The function diff() of EMT uses four points, which makes things a bit better, but not much.

>diff("exp(x)",1), E,
  2.7182818284604378 
  2.7182818284590451

For analytic functions, there is as simple trick to avoid the cancellation. We take a complex value for h! This is just as good an approximation as the simple differential with real h. Indeed a purely imaginary step yields

\(\text{re}\, \dfrac{f(x+ih)-f(x)}{ih} = \text{im}\, \dfrac{f(x+ih)-f(x)}{h} = \dfrac{\text{im}\,f(x+ih)}{h}\)

since the imaginary part of f(x) is zero. Now, there is no subtraction to be done, and we can use very small h.

>h=1e-16; im(exp(1+I*h))/h, E
  2.7182818284590451 
  2.7182818284590451

The trick is only good for analytic functions, and we assume that we can compute complex values easily. If the function involves complicated real computations we cannot use the trick at all.