math: Sin(Pi) and Cos(Pi/2) should return 0 but don't

[ijxy]

As the title says, the Sine and Cosine functions do not return zeros as expected.

Looking at the source, these values (which are surely special cases) are not checked in tests.

Could we add special cases to Sin and Cos (and Sincos) to return 0 identically at the required arguments, i.e. Sin(k*Pi) = Cos((2*k+1)*Pi/2) = 0, and add explicit tests for the first cases at least?

For example, a simple solution could be along the lines of:

func Sin(x float64) float64 {
	if r := x / math.Pi; r == float64(int(r)) {
		return 0
	}
        ...
}

func Cos(x float64) float64 {
	if r := x / (math.Pi / 2); r == float64(int(r)) {
		return 0
	}
	...
}

func Sincos(x float64) float64 {
	if r := x / math.Pi; r == float64(int(r)) {
                 // logic to determine whether Cos is + or -
		return 0, ±1
	}
	if r := x / (math.Pi / 2); r == float64(int(r)) {
                 // logic to determine whether Sin is + or -
		return ±1, 0
	}
        ...
}

What did you do?

https://play.golang.org/p/fqzS-dwzIFw

[randall77]

The number π isn't exactly representable in floating point. math.Pi is a close representation, but it isn't exact. So requiring that sin(N * math.Pi) be exactly zero isn't mathematically correct. And the larger N gets, the worse it gets, because the rounding error is being multiplied by N.

Your code won't even work, because the equality won't be satisfied for many multiples, because of the cited rounding error. For example, this code:

package main

import "math"
import "fmt"

func main() {
	for i := 0; i < 100; i++ {
		fmt.Printf("%d %g\n", i, (float64(i)*math.Pi)/math.Pi-float64(i))
	}
}

prints

0 0
1 0
2 0
3 0
4 0
5 0
6 0
7 0
8 0
9 0
10 0
11 -1.7763568394002505e-15
12 0
13 1.7763568394002505e-15
14 0
15 -1.7763568394002505e-15
16 0
17 0
18 0
19 0
20 0
21 0
22 -3.552713678800501e-15
23 0
24 0
25 0
26 3.552713678800501e-15
27 0
28 0
29 0
30 -3.552713678800501e-15
31 0
32 0
33 0
34 0
35 0
36 0
37 0
38 0
39 0
40 0
41 0
42 0
43 0
44 -7.105427357601002e-15
45 0
46 0
47 0
48 0
49 0
50 0
51 0
52 7.105427357601002e-15
53 0
54 0
55 0
56 0
57 0
58 0
59 0
60 -7.105427357601002e-15
61 0
62 0
63 0
64 0
65 0
66 0
67 0
68 0
69 0
70 0
71 0
72 0
73 0
74 0
75 0
76 0
77 0
78 0
79 0
80 0
81 0
82 0
83 1.4210854715202004e-14
84 0
85 0
86 0
87 0
88 -1.4210854715202004e-14
89 0
90 0
91 0
92 0
93 -1.4210854715202004e-14
94 0
95 0
96 0
97 0
98 0
99 1.4210854715202004e-14

Floats are weird. A lot of the things you expect from real numbers do not hold for floats.

[ijxy]

@randall77 I do understand that floating point arithmetic is not precise. However, would it not suffice to be within ε of an integer? e.g. Sin(x) = 0 for k - ε ≤ x/π ≤ k + ε where ε=5E-15 or something like it.

[ijxy]

Also, somewhat of an arbitrary and very subjective point, is that I think even if we cannot solve this for all k, we should at least solve it for the primary domain of the functions [-π, π].

[randall77]

I see two problems with solving this, even as you describe.
Consider the number float64(π). You got that number by doing var x float64 = π. But someone else may have gotten that number by some other mechanism, unrelated to π, say float64(z). When they compute the sin of that number, they want the most accurate result. 0 is not the correctly rounded sin(float64(z)) that they expect.
Second, doing this can lead to weird non-monotonicities in the output. If someone is computing sin(x) at values near π, they would expect that sin(x+ε) <= sin(x) for ε>0. I think truncating values to 0 destroys that property. (Caveat: I'm not sure our sin implementation is monotonic. I suspect it is, but I haven't checked.)

[ijxy]

All very fair points.

That said, unless it makes a noticeable difference in performance, I don't see why we would not add the simple special case checks. At best, it is getting the right answer and at worst it is giving the same (wrong) answer that we have at the moment. On my machine it adds about 5 ns:

BenchmarkMathSin-4   	180006896	        13.2 ns/op
BenchmarkSin-4       	132464058	        18.3 ns/op

package main_test

import (
    "testing"
    "math"
)

var GlobalF float64
var sin = math.Sin
var x = 1.23456789

func MathSin(x float64) float64 { return sin(x) }

func Sin(x float64) float64 {
    if r := x/math.Pi; r == math.Round(r) {
        return 0
    }
    return sin(x)
}

func BenchmarkMathSin(b *testing.B) {
    for i:= 0; i < b.N; i++ {
        GlobalF = MathSin(x)
    }
}

func BenchmarkSin(b *testing.B) {
    for i:= 0; i < b.N; i++ {
        GlobalF = Sin(x)
    }
}

[tandr]

On my machine it adds about 5 ns:

BenchmarkMathSin-4   	180006896	        13.2 ns/op
BenchmarkSin-4       	132464058	        18.3 ns/op

40% overhead