# Strange pattern in rounding errors?

This will look at first like a posting about trigonometry, then maybe about statistics, then finally about peculiarities of either

• a certain random process; or
• the pseudorandom number generator that I'm using; or
• other (specify!).

My question is: Which is it? And what's actually going on? I suspect it's the second alternative, but I'm not at all confident about that.

In the course of doing a bit of amateur cartography, I derived this little trigonometric relation:

If $(\cos\alpha,\sin\alpha)$ is in the right half of the unit circle (in other words, $\cos\alpha>0$), and $$\tan\gamma=\dfrac{\sin\alpha\sin\beta}{\cos\alpha+\cos\beta},$$ and $\cos\gamma$ is also positive, then $$\tan\dfrac\gamma2=\tan\dfrac\alpha2\cdot\tan\dfrac\beta2.$$ Numerical evidence bore out what I had derived, so now I should live happily ever after.

(And I was moderately intrigued by the resemblance to the simpler and more familiar tangent half-angle formula $\dfrac{\sin\alpha+\sin\beta}{\cos\alpha+\cos\beta}=\tan\dfrac{\alpha+\beta}{2}$.)

But then I asked what happens in the left half of the circle, where the cosine is negative. The answer turns out to be $$-\cot\frac\gamma2 = \tan\frac\alpha2\cdot\tan\frac\beta2.$$ [But see the "later note" below.]

But instead of deriving this by massaging trigonometric identities I got lazy and did some "experimental mathematics". Using R, I entered these commands:

a <- pi/180*(runif(1000)*(177 - 93) +93)
b <- pi/180*(runif(1000)*(177 - 93) +93)
c <- atan( sin(a)*sin(b)/(cos(a)+cos(b) ))
u <- -1/(tan(a/2)*tan(b/2))
coefficients(lm(tan(c/2) ~ u))


(Intercept) u
0 1

anova(lm(tan(c/2) ~ u))


Analysis of Variance Table

Response: tan(c/2)
Df Sum Sq Mean Sq F value Pr(>F)
u 1 29.222 29.222 2.1747e+34 < 2.2e-16 *
Residuals 998 0.000 0.000
---
Signif. codes: 0 ‘
’ 0.001 ‘*’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Warning message:
In anova.lm(lm(tan(c/2) ~ u)) :
ANOVA F-tests on an essentially perfect fit are unreliable

plot(u,residuals(lm(tan(c/2) ~ u)))


So this puts $a$ and $b$ between $\pi/2$ and $\pi$, or more precisely, between $93\cdot\frac{\pi}{180}$ and $177\cdot\frac{\pi}{180}$, and chooses 1000 such pairs $(a,b)$ independently, and they're uniformly distributed in that region. Then it sets $c=\arctan(\sin a\sin b/(\cos a+\cos b))$, and $u=-1/(\tan a\tan b)$. Then we look at coefficients from a simple linear regression of $\tan(c/2)$ on $u$, and the software reports $0$ for the intercept and $1$ for the slope. An analysis of variance gives $0$ as the sum of squares of residuals, so it seems we have a perfect fit.

Finally, I plotted $u$ on the horizontal axis and the residuals on the vertical axis, and I got the following! If real numbers rather than approximations could be used, they would of course all be $0$, so this is about rounding errors, but still I wouldn't expect to see a pattern like this. I tried it a dozen or so times with pretty much the same result, and I tried it with the angles in the first quadrant and the identity that holds in that quadrant, with the same result again.

LATER NOTE: Well, haste makes waste, I guess. I should have let $\gamma$ be the "other value of" the arctangent once I moved into the second quadrant, i.e. $\gamma=\arctan(\cdots\cdots\cdots)+\pi$ as soon as the argument to the arctangent function was more than $\pi/2$. That way we still have the identity $\tan\frac\gamma2=\tan\frac\alpha2\cdot\tan\frac\beta2$. However, this doesn't upset the main point of this question. As I said, this already works in the first quadrant; I simply hadn't yet noticed it because at that point I was still doing things intelligently rather than numerically.

-
I see this got 6 views while I was still correcting some dumb typos. Where should $\gamma$ appear without getting divided by $2$ and where should $\gamma/2$ appear. I hope I've got all of them right by now, and I apologize for the distress that the incorrect versions must have caused. –  Michael Hardy Dec 12 '12 at 22:51
I don't have a good explanation for this, but it's presumably an artifact of using double precision. A double precision number is a 52-bit approximation to something between 1 and 2 times a power of 2 (and a sign). This means if $u$ is between $0.125$ and $0.25$, you'll get errors in multiples of $2^{-55}$, which is about $2.8 \cdot 10^{-17}$. And that roughly matches the spacing here. As $u$ increases by factors of $2$, the spacing should also double, which fits the graph (it changes when $u$ passes roughly $0.5$, $0.25$, and maybe $0.125$). –  Henry Cohn Dec 13 '12 at 0:23
However, what I don't see how to explain is why the lines aren't perfectly flat (since double precision error should come in exact multiples in principle). Am I missing something? –  Henry Cohn Dec 13 '12 at 0:24
There should also be a good explanation of why the lines are most flat near $-0.2$ (but I don't see what it is offhand). –  Henry Cohn Dec 13 '12 at 0:25
For the R-illiterate amongst us, could you confirm that this is what the plot is of? You're sampling random vectors of alphas and betas in the range 93deg to 177deg (converted to radians). Next you compute the vector of gammas = arctan ((sin alpha.sin beta)/(cos alpha+cos beta)) Finally you're plotting tan(gamma/2)+cot(alpha/2)cot(beta/2) vs -cot(alpha/2)cot(beta/2)? –  Anthony Quas Dec 13 '12 at 0:27

For what it's worth, if you do the same computation on Mathematica, you see something quite different (and very much like what @Henry Cohn predicted):

-
That's great! Michael Hardy's picture shows more or less the same thing, except the horizontal lines are much fuzzier in the R version. I'm still puzzled by why that should be. I believe the Mathematica picture shows the true round-off error in double precision. One possibility is that the R code is computing something slightly different; for example, the slope of the regression presumably is not 1 exactly, but just extremely close to 1. However, it seems like this should introduce a gradual drift in the horizontal lines, not apparently random noise. (Maybe I'm missing something?) –  Henry Cohn Dec 13 '12 at 13:25
Or R could be computing differently in another way, for example using standardized or studentized residuals, but I don't think it is for these commands. I don't see how this could be a plotting issue, so it's got to have something to do with the underlying computation, but I'm not sure what. –  Henry Cohn Dec 13 '12 at 13:38
I'd be interested in seeing the results of a simpler computation (which in fact is what I did in Mathematica): plotting $\tan(\gamma/2)+\cot(\alpha/2)\cot(\beta/2)$ vs $-\cot(\alpha/2)\cot(\beta/2)$. I have some suspicions (without knowing anything about R) about what the lm and residuals functions are actually doing. –  Anthony Quas Dec 13 '12 at 13:59

Apparently it's something to do with the "residuals" function in R.

If you do this

h <- fitted.values(lm(tan(c/2)~u))
plot(u,h-tan(c/2),ylim=c(-2e-16,2e-16),cex=0.1)

instead (with the previous code being the same) you get the following:

-
Case closed (at least as far as I'm concerned) –  Anthony Quas Dec 13 '12 at 17:57
It would still be nice to know what R is doing to get the messier pattern Michael Hardy observed, but it's really more of an R question. –  Henry Cohn Dec 13 '12 at 21:30