The function jinc(x) that I wrote about yesterday is almost trivial to implement, but not quite. I’ll explain why it’s not quite as easy as it looks and how one might implement it in C and Python.

The function jinc(x) is defined as J₁(x) / x, so if you have code to compute J₁ then it ought to be a no-brainer. For example, why not use the following C code?

    #include <math.h>
    double jinc(double x) {
        return j1(x) / x;
    }

The problem is that if you pass in 0, the code will divide by 0 and return a NaN. The function jinc(x) is defined to be 1/2 at x = 0 because that’s the limit of J_1(x)(x) / x as x goes to 0. So we try again:

    #include <math.h>
    double jinc(double x) {
        return (x == 0.0) ? 0.5 : j1(x) / x;
    }

Does that work? Technically, it could still fail — we’ll come back to that at the end — but we’ll assume for now that it’s OK.

We could write the analogous Python code, and it would be adequate as long as we’re only calling the function with scalars and not NumPy arrays.

    from scipy.special import j1
    def jinc(x):
        if x == 0.0:
            return 0.5
        return j1(x) / x

Now suppose you want to plot this function. You create an array of points, say

    x = np.linspace(-1, 1, 25)

and plot jinc(x). You’ll get a warning: “ValueError: The truth value of an array with one element is ambiguous. Use a.any() or a.all().” Incidentally, if we called linspace with an even integer in the last argument, our array of points would avoid zero and the naive implementation of jinc would work.

When Python tries to apply jinc to an array, it doesn’t know how to interpret the test x == 0. The warning suggests “Do you mean if any component of x is 0? Or if all components of x are 0?” Neither option is what we want. We want to apply jinc as written to each element of x. We could do this by calling the vectorize function.

    jinc = np.vectorize(jinc)

This replaces our original jinc function with one that handles NumPy arrays correctly.

There is an extremely unlikely scenario in which the code above could fail. The value of J₁(x) is approximately x/2 for small values of x. If the floating point value x is so small that 0.5*x returns 0, our function will return 0, even though it should return 0.5. The C code above works for values of x as small as DBL_MIN and even values much smaller. (DBL_MIN is not the smallest value of a double, only the smallest normalized double.) But if you set

    x = DBL_MIN / pow(2.0, 52);

then jinc(x) will return 0. If you want to be absolutely safe, you could change the implementation to

    #include <math.h>
    double jinc(double x) {
        return (fabs(x) < 1e-8) ? 0.5 : j1(x) / x;
    }

Why test for whether the absolute value is less than 10^-8 rather than a much smaller number? For small x, the error in approximating jinc(x) with 1/2 is on the order of x²/16. So for x as large as 10^-8, the approximation error is below the resolution of a double. As a bonus, the function jinc(x) will be more efficient for |x| < 10^-8 since it avoids a call to j1.

Related posts

• Jinc function
• Sine approximation for small angles
• Functions in math.h that seem unnecessary