I'm struggling to understand what is the actual problem described as "loss of significance" when a trigonometric range reduction argument is performed. I was trying to perform the analysis by myself but I'm not able to spot the problem. I know what the loss of significance is in floating point computation in general, but I don't understand in what terms is described in the algorithm I'm mentioning.

Taken from the paper "Radian Reduction for Trigonometric Function" (Hanek Payne Algorithm)

  1. Even if M is large enough to guarantee (p+K) valid bits, h' may have so many leading zero bits that there may be fewer than the (p+K) valid bits desired for the normalized equivalent argument.
  2. Similarly h' may have leading bits consisting entirely of ones, and if (1-h) is required for the function evaluation, (1-h') may contain fewer than the desired (p+K) valid bits.

Some notation... more or less from the paper

For any real number $x$ the range reduction consist in computing something like $$ \frac{2^i}{2\pi}x \; \text{mod} \; 2^i = I + h $$ where $i \in {0,1,2,3}$, $I$ is an integer that would allow to select the quadrant and finally $h \in [0,1)$ (real number in that interval). Actually the reduction consist in using $I$ as selector bits, and $h$ to compute an appropriate approximation of $sin$ or $cos$ in the sector indexed by $I$ using the argument $\frac{2\pi}{2^i} h$ as argument. For machine purpose it is assumed that

$$x = 2^k f ,\; \frac{1}{2} \leq f < 1$$

I.e. a floating point number (different from the IEEE standard, but the analysis would be similar), the number $f$ is represented using $p$ digits.

The authors of such paper later define $L$ as the largest integer such that the product of 2 numbers represented using $L$ digits are representable in a register. They also point out that the number $\frac{2}{\pi}$ has to be stored with an appropriate number of digits (virtually infinite). Again we can then write

$$ \begin{array}{l} x = 2^k \sum_{j=1}^{N-1} F_j 2^{-jL} \\ \frac{1}{2\pi} = D_0 + D_1 + D_2 \\ D_0 = \sum_{j=1}^{M} g_j 2^{-jL} \\ D_1 = g_{M+1}2^{-(M+1)L} \\ D_2 = \sum_{j=M+2}^{+\infty} g_j 2^{-jL} \end{array} $$

where $F_j$ is the $j$-th chunk of $L$ bits representing the mantissa $f$ and $g_j$ is the $j$-th chunk representing $\frac{1}{2\pi}$. What happens is that basically

$$ \frac{2^i}{2\pi}x \; \text{mod} \; 2^i \approx 2^i x D_0 \text{mod} \; 2^i $$

is being computed. In general we are interested in computing $h$ with $p+K$ digits, $p$ is the input precision while $K$ is a number of guard digits. What I don't really understand is why the points 7 and 8 I've highlighted are a problem for the reduction. I was trying to make an accurate analysis of the problem by myself but I'm not able to spot the problem. Can anyone explain to me?

  • $\begingroup$ The problem Payne-Hanek is solving is a particular case of subtractive cancellation, caused by the finite precision of floating-point formats. For an input $x$ we want to compute $x - i \frac{\pi}{2}$. However, in a naive implementation, the $\frac{\pi}{2}$ is represented as a floating-point number PIO2, and if $x$, and thus $i$, are sufficiently large, all leading bits cancel in the subtraction. Try running this experiment with single-precision arithmetic and $x \gt 2^{24}$, and compare to the mathematical result, it should become clear. $\endgroup$ – njuffa Dec 1 '16 at 15:27
  • $\begingroup$ @njuffa Wait wait... but in your comment I'm understanding that you're assuming the output would be represented in floating point format. But one could just want the reduced argument in fixed point, would that be a problem? (i mean the presence of many leading 0 or many leading 1?) $\endgroup$ – user8469759 Dec 1 '16 at 15:30
  • $\begingroup$ Payne-Hanek applies to floating-point computation, and the input is therefore a floating-point number, and the output either a floating-point number, or a pair of floating-point numbers normalized such that the "tail" is less than a half ulp of the "head". This number (or pair of numbers), representing the reduced argument, is then used in the core approximation of the trigonometric function. Internal to the implementation, Payne-Hanek reduction may use fixed-point (or integer) computation. $\endgroup$ – njuffa Dec 1 '16 at 15:41
  • $\begingroup$ The mechanics of Payne-Hanek (especially as described in the original paper) can be confusing. Maybe looking at an actual implementation can help; note that this is not a verbatim implementation of the algorithm in the paper, but uses the fundamental ideas of it. This commonly cited write-up may provide a clearer explanation than the original Payne-/Hanek paper. $\endgroup$ – njuffa Dec 1 '16 at 15:50
  • $\begingroup$ It's not the implementation itself that struggles me, it's the analysis of certain aspects that puzzles me. I don't understand the analysis of that specific point. I'm not sure what do you with the "pair" argument. I do understand that some "core implementation" of the trigonometric function can be implemented using the reduced argument in floating point. But in general this is not actually necessary, you could keep the output in fixed point (so you don't need to normalize it then) and implement the core approximation with such value, you wouldn't need in such a case the renormalization step. $\endgroup$ – user8469759 Dec 1 '16 at 16:02

Your Answer

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

Browse other questions tagged or ask your own question.