Best way to implement numpy.sin(x) / x where x might contain 0

Question

What I am doing now is:

import numpy as np

eps = np.finfo(float).eps

def sindiv(x):
    x = np.abs(x)
    return np.maximum(eps, np.sin(x)) / np.maximum(eps, x)

But there is quite a lot of additional array operation. Is there a better way?

Answer 1

As others have said, numpy.sinc() is the easiest.

I want to include a copy of its current implementation in NumPy 1.21.2 (link) to show there's no special tricks:

y = pi * where(x == 0, 1.0e-20, x)
return sin(y)/y

It's basically just sin(x)/x. Note that in creating y: multiplication by pi, where(), and x == 0 will create at least 2 intermediate arrays plus the final array for y. And then sin(y)/y creates two more arrays. In total at least 5 arrays are created by numpy.sinc(); and by my count your sindiv() also creates at least 5 arrays, so it's not actually that wasteful.

---

Here is another implementation:

TINY = np.finfo(float).tiny  # ≈ 2e-308 (smallest 'normal' float)
def mysinc(x):
    y = np.abs(np.pi*x) + TINY
    return np.sin(y)/y

I'm pretty sure this returns identical values to numpy.sinc(). The reason being sin(x) == x for relatively 'large' values of x:

x = np.ldexp(1, -26, dtype=np.double)  # x = 2**-26 ≈ 1.5e-8
print(np.sin(x) == x)  # True

x = np.ldexp(1, -32, dtype=np.longdouble)  # x = 2**-32 ≈ 2.3e-10
print(np.sin(x) == x)  # True

So for small enough x (ignore pi factors), mysinc(x) = (x+TINY)/(x+TINY) = x/x = np.sinc(x). The exact threshold this happens does not matter too much so long as TINY < np.spacing(x) when it occurs so that x + TINY = x in this regime.

(The cutoff is around the square-root of the machine epsilon as can be understood from the Taylor series sin(x) = x - x**3/6 + ... = x(1-x**2/6) + .... So TINY is always small enough to not matter.)

Timings

import numpy as np

eps = np.finfo(float).eps
tiny = np.finfo(float).tiny

def npsinc(x):
    y = np.pi * np.where(x == 0, 1.0e-20, x)
    return np.sin(y)/y

def sindiv(x):
    x = np.pi * np.abs(x)
    return np.maximum(eps, np.sin(x)) / np.maximum(eps, x)

def mysinc(x):
    y = np.abs(np.pi*x) + tiny
    return np.sin(y)/y

def mysinc2(x):
    y = np.abs(np.pi*x)
    y += tiny  # in-place addition
    return np.sin(y)/y

# Test data
x = np.random.rand(100)
x[np.random.randint(100, size=10)] = 0

%timeit npsinc(x)
# 10.9 µs ± 18.3 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit sindiv(x)
# 9.4 µs ± 12.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit mysinc(x)
# 7.38 µs ± 15.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit mysinc2(x)
# 8.64 µs ± 20.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Curiously using mysinc2() with in-place addition seems to be slower, and using in-place numpy.abs() and in-place numpy.sin() is even slower. Not entirely sure why, but see this related question.

Regardless, if you really need performance, you can try using Cython to generate C code and do things properly instead of playing tricks with NumPy:

%%cython

from libc.math cimport M_PI, sin
cimport cython
cimport numpy as np
import numpy as np

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cdef _cysinc(double[:] x, double[:] out):
    cdef size_t i
    for i in range(x.shape[0]):
        if x[i] == 0:
            out[i] = 1
        else:
            out[i] = sin(M_PI*x[i])/(M_PI*x[i])

def cysinc(np.ndarray x):
    out = np.empty_like(x)
    _cysinc(x.ravel(), out.ravel())
    return out

%timeit cysinc(x)
# 4.38 µs ± 11.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

As always, don't prematurely optimize, just use numpy.sinc() to begin with.

---

Side note

There's a question Is boost::math::sinc_pi unnecessarily complicated? that asks about the benefits of using a Taylor expansion about x=0. In summary, almost none, but maybe they are doing it for other reasons.

To emphasise, there is nothing unstable about floating point division, or dividing a small number by a small number since you're just dividing the significands and subtracting the exponents.

If you calculate sinc(x) as sin(x)/x, instead of a direct Taylor series or other method that sums to convergence beyond the machine epsilon np.spacing(sinc(x)), you will be off by at most np.spacing(sinc(x)) coming from the round-off error in division /, just as you'd get with multiplication *. (Assuming no subnormal business, which even here does not matter in the treatment of sin(x)/x.)

Answer 2

You could use numpy.sinc, which computes sin(pi x)/(pi x):

In [20]: x = 2.4

In [21]: np.sin(x)/x
Out[21]: 0.28144299189631289

In [22]: x_over_pi = x / np.pi

In [23]: np.sinc(x_over_pi)
Out[23]: 0.28144299189631289

In [24]: np.sinc(0)
Out[24]: 1.0

Answer 3

What about allowing div by zero and replace NaNs later?

import numpy as np

def sindiv(x):
    a = np.sin(x)/x
    a = np.nan_to_num(a)
    return a

If you don't want warnings, supress them via seterr

Of course, using a could be eliminated:

def sindiv(x):
    return np.nan_to_num(np.sin(x)/x)

Answer 4

In numpy array notation (so you get back a np array):

def sindiv(x):
    return np.where(np.abs(x) < 0.01, 1.0 - x*x/6.0, np.sin(x)/x)

Here I've made "epsilon" fairly large for testing and used the first two terms of the taylor series for the approximation. In practice, I'd change 0.01 to some small multiple of your eps (machine epsilon).

xx = np.arange(-0.1, 0.1, 0.001)
yy = sinxdiv(xx)
type(yy)

outputs numpy.ndarray and the values are continuous (and differentiable, if that's important) near the origin.

If you don't want the double evaluation (i.e. both branches are evaluated in the above), then I think you have to go with a loop as I don't believe there is any sort of "lazy where" option.

def sindiv(x):
    sox = np.zeros(x.size)
    for i in xrange(x.size):
        xv = x[i]
        if np.abs(xv) < 0.001: # For testing, use a small multiple of machine epsilon
            sox[i] = 1.0 - xv * xv / 6.0
        else:
            sox[i] = np.sin(xv) / xv
    return sox

To make this really pythonic though it would be best to check the type of x and just do the non-array version if it is not an array.