# Astoundingly large difference when evaulating trigonometric identity with NumPy

According to Wolfram Alpha and the Sage computer algebra system, the following identity holds: $$\cos\left(\arctan\left(\frac{l_1-l_2}{d}\right)\right) = \frac{1}{\sqrt{1 + \frac{(l_1-l_2)^2}{d^2}}}$$

However, when I tried to verify it with an arbitrary example in NumPy, I noticed a rather large difference in the actual values computed by both sides of the idenitity. I have used the following code:

    l1 = 10; l2 = 8; d = 17
from numpy import arctan2, cos, sin, sqrt

alpha = arctan2((l1-l2),d)
left = cos(alpha)

right = sqrt(1 + ((l1-l2)**2)/(d**2))


Evaluating the results for left and right yielded the following:

    left = 0.99315060432287616
right = 1.0


It is tempting to write this off as simply being a numerical error, but since I have very little experience in how large numerical errors can get, I am not so sure. Is this possible or am I missing something (obvious)?

• your right is incorrectly entered. it should be right = 1/sqrt() When I enter the formulas into my Ti-89 I get a match out to 12 digits at 0.99315... – Godric Seer Apr 26 '13 at 15:39
• Yeah, I think he would have caught that if the square root expression had been evaluated correctly. – Michael Grant Apr 26 '13 at 15:40
• yeah, when I entered his expression I got 1.007, which definitely would have been visible. – Godric Seer Apr 26 '13 at 15:42
• Oh, really. That is embarrasing. Thank you for noticing. – Daniel Eberts Apr 26 '13 at 15:42
• Made this mistake before myself. – Michael Grant Apr 26 '13 at 17:15

I have a suspicion that your Python expression for the right-hand side is doing integer division, not floating-point division. As a result, ((l1-l2)**2)/(d**2) is being evaluated as zero, and the term inside the square root is one.

In fact, you forgot the reciprocal in your right-hand expression as well, but that's not the first problem...

In MATLAB:

>> l1 = 10; l2 = 8; d = 17;
>> alpha = atan2((l1-l2),d)
alpha =
0.1171
>> left = cos(alpha)
left =
0.9932
>> right = 1/sqrt(1 + ((l1-l2)^2)/(d^2))
right =
0.9932

• I've gotten into the habit of using from __future__ import division at the top of all my python scripts for this purpose. – Geoffrey Irving Apr 26 '13 at 16:28
• @GeoffreyIrving, turn that into an answer and I'll bet you get some serious rep for that. Good advice. – Michael Grant Apr 27 '13 at 15:52

As Michael C. Grant pointed out, the problem is that the division operation is integer division, not floating point division. In versions of Python before 3, dividing two integers using / rounds down to an integer. This behavior can be changed by adding

from __future__ import division


at the top of your script. Note that this declaration changes the behavior of / only for the current file. If you want the original behavior occasionally, use // (double slash). This is discussed in more detail here:

PEP 238 -- Changing the Division Operator

Note that while most languages including C have similar semantics for division, the semantics are more confusing in Python because integers are not usually promoted to floats when passed to functions. Thus, a function written as if it was operating on floats may in fact be operating on integers.