# Avoiding numerical cancellation question for $\sin x -\sin y$ for $x \approx y$

When trying to avoid cancellation, one tries to reformulate the equation in order to avoid subtraction between almost equal terms.

In $\sin (x) - \sin (y), x \approx y$ the suggested solution is to reformulate it to $$2\cos\left(\frac{x+y}{2}\right)\sin\left(\frac{x-y}{2}\right)$$

But I don't understand how it is any better? The subtraction between $x, y$ remains. Is it because $\lvert\,\sin(x)-\sin(y)\,\rvert\leq \lvert\,x-y\,\rvert$, so cancellation is less likely to happen between $x,y$ than the sines?

It really depends on how exactly $x$ and $y$ are given. Frequently what we really have is not $x$ and $y$ but $x$ and $y-x$. We can then gain precision if we can analytically write the subtraction of nearly equal numbers exclusively in terms of $y-x$, like you did here, because we are already given $y-x$ accurately (more accurately than we would get if we computed it directly).

There are actually some standard library functions out there specialized to this exact purpose, for example "log1p", which is used to compute $\log(1+x)$ for small $x$.

• So assuming $x=y+h$ it would become $2\cos (x+\frac{h}{2})\sin(\frac{h}{2})$. Can it be assumed that $y-x$ is given though? As a beginner in numerics I am never sure when the exact values are available. Oct 24, 2017 at 12:46
• @B.Swan Right. Then at least the $\sin(h/2)$ can be accurately computed if $h$ is already given accurately.
– Ian
Oct 24, 2017 at 12:47
• @B.Swan It really depends on how the problem is given, as I said. But if $y-x$ is not already given accurately, then this trick gives a modest benefit at best (because you will already commit severe error inside the sine). Let me put it another way: this approach does not really help you define a function delta_sin(x,y). It helps you define a function delta_sin(x,h). If your caller needs delta_sin(x,y), then you will need to do something else to help them. But if they need delta_sin(x,h), then you can help them. Does that make sense?
– Ian
Oct 24, 2017 at 12:48

While a function $$f : A \to B$$ is a triple, consisting of a domain $$A$$, a codomain $$B$$ and a rule which assigns to each element $$x \in A$$ exactly one element $$f(x) \in B$$, too many focus exclusively on rule and forget to carefully specify the domain and the codomain.

In this case the function in question is $$f : \mathcal F \times \mathcal F \rightarrow \mathbb R$$, where $$f(x,y) = \sin(x) - \sin(y),$$ and $$\mathcal F$$ is the set of machine numbers, say, double precision floating point numbers. I will explain below why $$\mathcal F$$ this is the correct domain.

The problem of computing a difference $$d = a - b$$ between two real numbers $$a$$ and $$b$$ is ill-conditioned when $$a \approx b$$. Indeed if $$\hat{a} = a(1+\delta_a)$$ and $$\hat{b} = b(1 + \delta_b)$$ are the best available approximations of $$a$$ and $$b$$, then we cannot hope to compute a better approximation of $$d$$ than $$\hat{d} = \hat{a} - \hat{b}.$$ The relative error $$r = \frac{d - \hat{d}}{d},$$ satisfies the bound $$|r| \leq \frac{|a| + |b|}{|a-b|} \max\{|\delta_a|,|\delta_b|\}.$$ When $$a \approx b$$, we cannot guarantee that the difference $$d$$ is computed with a small relative error. In practice, the relative error is large. We say that the subtraction magnifies the error committed when replacing $$a$$ with $$\hat{a}$$ and $$b$$ with $$\hat{b}$$.

In your situation $$a = \sin(x)$$ and $$b = \sin(y)$$. Errors are committed when computing the sine function. No matter how skilled we are, the best we can hope for is to obtain the floating point representation of $$a$$, i.e. $$\text{fl}(a) = \sin(x)(1 + \delta)$$, where $$|\delta| \leq u$$ and $$u$$ is the unit roundoff. Why? The computer may well have extra wide registers for internal use, but eventually, the result has to be rounded to, say, double precision, so that the result can be stored in memory. It follows, that if we compute $$f$$ using the definition and $$x \approx y$$, then the computed result will have relative error which is many times the unit roundoff.

In order to avoid the offending subtraction, we turn to the function $$g : \mathcal F \times \mathcal F \to \mathbb R$$ given by $$g(x,y) = 2 \cos \left( \frac{x+y}{2} \right) \sin \left(\frac{x-y}{2} \right)$$ In absence of rounding errors $$f(x,y) = g(x,y)$$, but in floating point arithmetic they behave quite differently. The subtraction of two floating point numbers $$x$$ and $$y$$ is perfectly safe. In fact, if $$y/2 \leq x \leq 2y$$, then subtraction is done with one guard digit, then $$x-y$$ is computed exactly.

We are not entirely in the clear, as $$x + y$$ need not be a floating point number, but is computed with a relative error bounded by the unit roundoff. In the event that $$(x+y)/2 \approx (\frac{1}{2} + k) \pi$$ where $$k \in \mathbb Z$$, the calculation of $$g$$ suffers from the fact that cosine is ill-conditioned near a root.

Using a conditional to pick the correct expressions allows us to cover a larger subset of the domain.

In general, why $$\mathcal F$$ rather than $$\mathbb R$$? Consider the simpler problem of computing $$f : \mathbb R \rightarrow \mathbb R$$. In general, you do not know the exact value of $$x$$, and the best you can hope for is $$\hat{x}$$, the floating point representation of $$x$$. The impact of this error is controlled by the condition number of $$f$$. There is nothing you can do about large condition numbers, except switch to better hardware or simulate a smaller unit roundoff $$u'$$. This leaves you with the task of computing $$f(\hat{x})$$, where $$\hat{x} \in \mathcal F$$ is a machine number. That is why $$\mathcal F$$ is the natural domain during this the second stage of designing an algorithm for computing approximations of $$f : \mathbb R \to \mathbb R$$.

• The book I am using ("Algorithmic mathematics" by Vygen/Hougardy, I think it's only available in German) actually does define a computation problem as a triple with two sets $A, B$ and a relation $R$, which is a subset of $A \times B$, and in numerical computational problems it actually sets $A,B$ as subsets of $\mathbb{R}$, not necessarily as the machine numbers. It does that in the definition of the condition, too. I will take my time to think about your answer though and discuss it with my professor, thank you. Oct 24, 2017 at 20:38
• @B.Swan Conditioning should certainly be considered for subsets of $\mathbb R$. Once this first stage is complete, one must consider how the target function is evaluated for floating point numbers. This is frequently difficult as seen above, and a rewrite or an approximation must be sought. Then comes the third and final stage which consist evalutating the chosen approximation. We want $T(x)$, but can only get $\hat{A}(\hat{x})$. The error is can be expressed as $T(x) - \hat{A}(\hat{x}) = T(x) - T(\hat{x}) + T(\hat{x}) - A(\hat{x}) + A(\hat{x}) - \hat{A}(\hat{x})$ Oct 24, 2017 at 21:55

One possibility: To determine the sine you use the Maclaurin series, and the faster this converges the fewer ill-comditioned operations you need to perform for the factor $\sin(\frac{x-y}{2})$. The small argument in that factor gets you there.

• Unless $x,y$ are large and their sines are nearly equal, I don't think this has anything to do with it.
– Ian
Oct 24, 2017 at 12:33
• If $x$ and $y$ are large the Maclaurin series themselves can become ill-conditioned. You can use periodicity and symmetry to reduce them if you have a sufficiently accurate rendering of $\pi$. Oct 24, 2017 at 12:37
• That periodicity trick is actually very badly conditioned, cf. math.stackexchange.com/questions/1561713/… What I'm really saying here is that if $x,y$ are, say, confined to $[-3\pi/4,3\pi/4]$, the series themselves are not badly conditioned but the difference is badly conditioned. In these cases accelerating the convergence of the series does you rather little good. (Indeed, you actually pass the problem off to the cosine function without making it any better.)
– Ian
Oct 24, 2017 at 12:40
• Periodicity and symmetry. You can reduce arguments to an absolute value less than or equal to $\pi/4$, actually. Oct 24, 2017 at 12:42
• Confining to $[-\pi/4,\pi/4]$ doesn't improve matters much more than confining to $[-3\pi/4,3\pi/4]$. Indeed the problem with keeping the series well-conditioned is to avoid being close to a zero of sine other than $x=0$. But again, provided you stick to a range where the series themselves are well-conditioned (or you use another method entirely to compute the trig functions themselves), the problem of computing the difference persists.
– Ian
Oct 24, 2017 at 12:43