2

I recently wrote How calculators do trigonometry where I wrote a simple program for computing $\sin(x)$. For completeness, I include a slightly modified version of the program here:

function y=mysin(x)
z=mod(x,2*pi);
if 3*pi/2<=z
  y=-sin1(2*pi-z);
elseif pi<=z
  y=-sin1(z-pi);
elseif pi/2<=z
  y=sin1(pi-z);
else
  y=sin1(z);
end
end

function y=sin1(x)
if x>pi/4
  y=cos1(pi/2-x);
else
  n=0:7;
  m=2*n+1;
  y=sum((-1).^n.*x.^(m)./factorial(m));
end
end

function y=cos1(x)
if x>pi/4
  y=sin1(pi/2-x);
else
  n=0:8;
  m=2*n;
  y=sum((-1).^n.*x.^(m)./factorial(m));
end
end

Provided all of the calculations were done in exact arithmetic, it would compute $\sin(x)$ for any real $x$ to within an accuracy of $10^{-16}$. I noticed that it has nevertheless has some error with large arguments; in particular, in Matlab I get 

abs(sin(50000)-mysin(50000))=3.4861e-14.

I have concluded that this is caused by the rounding in the first step, where I compute mod(x,2*pi), because in Matlab I get 

abs(mysin(50000)-sin(mod(50000,2*pi)))=0

How would I fix this? Would I need to perform this first step in higher precision somehow?

Edit: upon discussion with LutzL, I see that I can pose this problem in a more well-defined way as follows. Given any real number $x$, how can I calculate $x$ mod $2 \pi$ to within double precision, given a method for computing $\pi$ itself to arbitrary precision?

Community
  • 1
Ian
  • 101,645
  • Take a try to change into 'z=2pimod(x/(2*pi),1);' – Kay K. Dec 05 '15 at 23:41
  • @KayK. I think I see what you were going for (use a divisor that is floating-point representable in the call to mod), but that actually makes it slightly worse. – Ian Dec 05 '15 at 23:49
  • I see. Then what about z=2pi(x/(2pi)-int(x/(2pi)))? The idea is not to use the mod() function. – Kay K. Dec 05 '15 at 23:57
  • Would you try this too? z=mod(1000x,10002*pi)/1000; – Kay K. Dec 06 '15 at 00:00
  • @KayK. Your first two suggestions roughly double the error relative to my "naive" implementation. Your last one roughly triples it. (Let me stress that I see your line of thinking, and I do not think these are dumb ideas. I am just saying that as it turns out they do not work.) – Ian Dec 06 '15 at 00:03
  • Okay. Thanks :-). Then let's do an opposite one. z=mod(x/1000,2pi/1000)1000; – Kay K. Dec 06 '15 at 00:07
  • @KayK. With 1000 the error is roughly halved. But making it much larger increases the error again. – Ian Dec 06 '15 at 00:14
  • I see. I'd like to give a last try. What about mod( mod(x/100,2pi) 100 , 2*pi) ? – Kay K. Dec 06 '15 at 00:22
  • @KayK. The results are about the same as my original implementation for that version. – Ian Dec 06 '15 at 00:29
  • One suggestion is to do 'abs(mod(50000, 2pi) + int(50000/(2pi))2pi - 50000)' to actually see how much error mod() has, unless you've been already doing so.. – Kay K. Dec 06 '15 at 00:39
  • @KayK. It will depend on the exact order of operations. Surprisingly to me, mod(50000,2*pi)+2*pi*floor(50000/(2*pi))-50000 actually returns 0. But (2*pi*floor(50000/(2*pi))-50000)+mod(50000,2*pi) returns about $10^{-12}$. – Ian Dec 06 '15 at 00:43
  • I suggest reading Goldberg's paper What Every Computer Scientist Should Know About Floating-Point Arithmetic or some of the many StackOverflow questions about floating point computations. – shoover Dec 07 '15 at 17:09

1 Answers1

1

The floating point number $50000$ represents actually an interval of numbers contained in $[50000·(1-\mu),50000·(1+\mu)]$, with $\mu=2^{-53}$. The exact interval would require an analysis of the binary representation of $50000$.

Since the computer has no knowledge which of the numbers inside this interval is the given one, the most exact reduction by $2\pi$ still has an uncertainty of $50000·\mu\approx 5·10^{-12}$ which translates into an error interval of the same magnitude in the sine computation.

Thus the observed error between the two methods of computation of $3·10^{-14}$ is well within these margins.

Lutz Lehmann
  • 126,666
  • I understand the origin of the problem, as I said. My question was about possible solutions. – Ian Dec 07 '15 at 12:24
  • There is no problem, thus there are no solutions. What you are after is to duplicate the matlab algorithm or at least provably duplicate its behavior. For that you would need to specify that behavior. – Lutz Lehmann Dec 07 '15 at 16:28
  • I disagree with your perspective on the situation. I have reduced the problem to error in the result of mod(x,2*pi), since the results of my calculation are consistent with Matlab's built-in method if I precompose Matlab's method with mod(x,2*pi). I said this in the question. So assuming Matlab's method is itself floating point accurate, I could define my problem this way: I would like to devise a floating point accurate way to compute $x$ mod $2 \pi$ (not the floating point $\pi$ but the real $\pi$), given a method for computing $\pi$ itself to arbitrary precision. – Ian Dec 07 '15 at 16:54