Accurate and efficient computation of the logarithm of the ratio of two sines

Question

For exploratory work related to special function implementations, I need to compute $\log \frac{\sin y}{\sin x} $, where $0 \le x \le y \le 2x < \frac{\pi}{2}$. Cases with $x \approx y$ in particular are critical to overall accuracy.

Given that the ratio of the sines is often close to unity, I want to use the log1p function to compute the logarithm as accurately as possible, which means I need to find a way to compute $\frac{\sin y}{\sin x}-1$ accurately.

Given the pre-conditions, based on the Sterbenz lemma, $\delta = y - x$ can be computed exactly with binary floating-point arithmetic. With the help of the angle-sum and half-angle formulas, I then get

$$ \sin y = \sin(x+\delta) = \sin(x) + \left(\sin(\delta) \cos(x) - 2 \sin^{2}\left(\frac{\delta}{2}\right) \sin x \right) $$

from which follows immediately

$$\frac{\sin y}{\sin x} - 1 = \sin(\delta) \frac{\cos x}{\sin x} - 2 \sin^{2} \left(\frac{\delta}{2}\right) $$

Given the pre-conditions, there is no risk of cancellation in the subtraction, since the minuend is at least twice as large as the subtrahend, and usually much larger than that. This computation is performance sensitive, and since a function sincos is available that computes $\sin$ and $\cos$ in one go, I have also considered rewriting the above as follows to reduce the cost of computing all transcendentals to just two sincos calls (presumably trading-off with a small increase in round-off error)

$$2\sin\left(\frac{\delta}{2}\right) \cos\left(\frac{\delta}{2}\right)\frac{\cos x}{\sin x} - 2\sin^{2}\left(\frac{\delta}{2}\right)$$

This could be further transformed into the following but I have not checked yet whether this is actually advantageous

$$2\sin\left(\frac{\delta}{2}\right) \left(\cos\left(\frac{\delta}{2}\right)\frac{\cos x}{\sin x} - \sin\left(\frac{\delta}{2}\right)\right)$$

Is there an alternate arrangement of this computation that also maintains full accuracy and further minimizes computational cost? The availability of fused-multiply add (FMA) can be assumed. Abstract operational costs are as follows: add, sub, mul, fma = 1; div, sqrt, sin, cos = 10; log, log1p, sincos = 15; tan = 20.

Answer 1

Consider the following taylor series expansion of $sin(y)/sin(x)-1$ at $y=x$, with $δ=y-x$: $$sin(y)/sin(x)-1=δcot(x)-\frac12δ^2-\frac16δ^3cot(x)+\frac1{24}δ^4+\frac1{120}δ^5cot(x)...$$ Thanks Wolfram! https://www.wolframalpha.com/input/?i=series+sin%28y%29%2Fsin%28x%29+y+%3D+x

This only requires you to compute a single trigonometric operation cot(x) and terms $δ^n/n!$. If you can compute cot(x) accurately and δ is small, its easy to see that this converges pretty quickly. You could even re-use values $δ^n/n!$ for subsequent iterations.

If x is close to zero such that cot(x) is garbage, then you may have to try something else, maybe L'Hopital's Rule?

EDIT:

A less "clever" approach is to instead consider the following taylor series at $x=0$ where $y/x=a$: $$sin(ax)/sin(x)=a+\frac12a(1-a^2)x^2+...nasty\ terms$$ For your domain, $a$ is between 1 and 2. As long as you can calculate $a$ accurately, then it will never a problem. Evaluating the function at $x=0$ is also very well behaved. Perhaps you can switch between the two forms as needed.