Search code examples
numpyprecisiontrigonometry

Precision issue when using numpy.sin

I am trying to generate a square wave starting from a sine wave using numpy.
I ran into a problem where:

  • using the numpy.pi constant gives inaccurate results
  • using the 3.14 constant gives me the correct result

To create the square wave the sign of the sine wave is taken and then the -1s are transformed into 0s using the maximum function.

N = 15
x = np.arange(1, N)
s = np.maximum(np.sign(np.sin(3.14 * (x - 1))), np.zeros_like(x))
print(s)

>> [0. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 1.]

While using the constant for pi in numpy gives:

N = 15
x = np.arange(1, N)
s = np.maximum(np.sign(np.sin(np.pi * (x - 1))), np.zeros_like(x))
print(s)

>> [0. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 0.]

where the sequence is not alternating 1s and 0s. With more elements you can observe patches of more than two 0s next to each other.

[0. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 1. 0. 0. 0. 1. 0. 0. 0.]

I am mostly trying to understand how using more precise values can lead to imprecise results.

Solution

  • Your computation is numerically unstable. Using a higher precision will not solve the problem but just hide it a little bit more. The main issue is that sin(PI * x) should be theoretically always 0, but no floating-point function is perfect, so there is an error (typically 1 ULP). This error change the sign of the result. This is what you observe. This problem is independent of the precision. Here is a geometrical illustration of what you compute:

    enter image description here

    If you just want alternated 0-1 values, you can add a shift to the sin function (eg. np.maximum(np.sign(np.sin(np.pi * (x - 1) + np.pi/2)), np.zeros_like(x))). If the shift is pi/2, then you can use a cos function.