4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ If the cost of square root is not high then you can evaluate $\cos$ as $(1-\sin^2)^{1/2}$, that would improve it a bit further. $\endgroup$ – Maxim Umansky Sep 13 at 23:44
  • $\begingroup$ @MaximUmansky sqrt is about the same cost as sin, cos, which have fairly low cost here due to the limited range of the inputs. I will update the question with that information. $\endgroup$ – njuffa Sep 14 at 1:06
  • $\begingroup$ The taylor series expansion of $sin(y)/sin(x)$ at $y = x$ is surprisingly simple (requires $cot(x)$ and some integral exponentiation of $y-x$. Fortunately for you, the first term is 1! Not sure if this helps or not. wolframalpha.com/input/… $\endgroup$ – Charlie S Sep 14 at 2:55
  • $\begingroup$ @CharlieS cot is not available right now. I can synthesize it from div and sincos, or (given more time) create a dedicated implementation, since my current approach would need that anyhow. I would appreciate a more detailed answer along the lines of your comment. $\endgroup$ – njuffa Sep 14 at 5:16
  • $\begingroup$ Herbie is a nice tool if you often find yourself asking this kind of question. $\endgroup$ – Daniel Shapero Sep 14 at 22:27
2
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ That's a good idea. And if x is close to zero then we can use this expansion for sin(x)/sin(y) instead. $\endgroup$ – Maxim Umansky Sep 14 at 15:29
  • $\begingroup$ A nice way to write this answer could be: $\cot(x)*S_1(\delta) + S_2(\delta)$ where $S_1$ and $S_2$ are the two series. Perhaps the series can be summed up to something known? Not necessarily better for computing it but that would give an elegant form of the answer. $\endgroup$ – Maxim Umansky Sep 14 at 15:40
  • 1
    $\begingroup$ Actually, for $x \ll$1 and $\delta \ll$1 you just expand $\sin(x)$ and $\sin(x+\delta)$, and the leading term for $([\sin(x+\delta)/\sin(x)]-1)$ is $\delta/x$. $\endgroup$ – Maxim Umansky Sep 15 at 4:28
  • $\begingroup$ @MaximUmansky Simple is the best! It is evident there are many ways to approach this problem. $\endgroup$ – Charlie S Sep 15 at 14:50
  • $\begingroup$ Your first method is pretty good. I am just saying that for small $x$ you basically keep the first term there because it is dominant. It is not a different way to do it, it is just a small improvement to your first method. Your second method is not as useful as the first method, the convergence is slow there. $\endgroup$ – Maxim Umansky Sep 15 at 17:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.