I've been working recently on the AMPL interface for the GNU Scientific Library. Adding a function to AMPL normally requires providing first and second partial derivatives for it. Computing the derivatives is straightforward and there are plenty of tools, such as Wolfram Alpha, that automate this. However, transferring this to code is a tedious and error-prone process.

For example, consider the function
`gsl_sf_bessel_I0_scaled(double x)`

which is defined as
\(e^{-|x|} I_0(x)\), where \(I_0(x)\) is the
regular modified cylindrical Bessel function of zeroth order. It has a
simple derivative:

\( e^{-|x|} I_1(x)- \frac{x}{|x|} e^{-|x|} I_0(x) \)

The derivative can be expressed even in a more compact form using
\(I_0(x)\) and \(I_1(x)\) scaled by \(e^{-|x|}\),
`gsl_sf_bessel_I0_scaled`

and
`gsl_sf_bessel_I1_scaled`

.

Transferring this into C results in the following small program:

Let's compile and run the program:

```
$ gcc -pedantic -Wall -Wextra test.c -lgsl -lgslcblas -otest
$ ./test
-0.332111
```

It is
easy to check that the printed value is incorrect, `-0.148394`

should be printed instead. Can you spot the error?

`abs`

function operates on integers,
so the floating-point argument `1.5`

is silently converted into
the integer `1`

by discarding the fractional part. It can be
trivially fixed by replacing `abs`

with `fabs`

.
Unfortunately GCC doesn't give any warning even with `-Wall`

and
`-Wextra`

options, although there is an option called
`-Wconversion`

that makes GCC warn you of this kind of problems.
This small example demonstrates how important it is to test the code. Unfortunately it is rarely discussed in the context of OR. As Tim Hopper put it in his recent blog post:

Code testing was never mentioned in any of my classes in college, and I never hear operations researchers talk about it. Again, I have no doubt that the code behind much published work is full of mistakes. Operations researchers need good testing practices?

Returning to the question of derivatives, how do you actually test that the implementation is correct? One possibility is to use numerical differentiation. The simplest method is directly based on the definition of the derivative:

\( f'(x) \approx \frac{f(x + h) - f(x)}{h} \) for some small \(h\)

Although trivial to implement, this method suffers from several issues such as the difficulty of selecting \(h\) especially for \(x = 0\). So I used the Ridders' method described in the book Numerical Recipes in C: The Art of Scientific Computing. This methods starts from large \(h\) and successfully decreases it estimating error at each iteration.

Numerical differentiation proved to be very useful for testing. Although it took some time to implement the test infrastructure, it paid off by catching a number of bugs and corner cases. And now that the test infrastructure is in place, adding tests for new functions requires minimal amount of code.