Omnimaga
General Discussion => Other Discussions => Math and Science => Topic started by: Xeda112358 on November 20, 2014, 10:07:48 am

Hi again! This time I was fooling around with tangent (who knows why? I don't .__.) Anyways, I was having some fun and I figured out that Horner's method can be abused to approximate tangent by a similar logic I used to approximate arctangent: all you have to do is approximate infinitely many terms of the Maclaurin series with one simple function!
So before the fun, I'll show what Horner's method is (it's a super simple idea):
If you've never seen the identity that ##sin(x)=x\frac{x^{3}}{3!}+\frac{x^{5}}{5!}\frac{x^{7}}{7!}+...##, that's fine, just trust me (you'll get to that in calculus and it's crazy cool stuff). Horner's method is like what a true 31337 programmer would come up with if you want to efficiently approximate that polynomial. First, factor out the first term in the series (##x##):
##\sin(x)=x(1\frac{x^{2}}{3!}+\frac{x^{4}}{5!}\frac{x^{6}}{7!}+...##
Now in that inside term, factor out the ##\frac{x^{2}}{3!}##:
##\sin(x)=x(1\frac{x^{2}}{3!}(1\frac{x^{2}}{4\cdot 5}+\frac{x^{4}}{4\cdot 5\cdot 6\cdot 7\cdot }\frac{x^{6}}{4\cdot 5\cdot 6\cdot 7\cdot 8\cdot 9}+...##
Rinse, repeat:
##\sin(x)=x(1\frac{x^{2}}{3!}(1\frac{x^{2}}{4\cdot 5}(1\frac{x^{2}}{6\cdot 7}(1\frac{x^{2}}{8\cdot 9}...##
So the programmer sees, "aw, shweet, I only need to compute x^{2} once and reuse it a bunch of times with some constant multiplications!" So yeah, Horner's method is snazztacular, and I urge you to use it when possible.
Now let's look at tan(x). The Maclaurin series expansion for tan(x) is uuuugly. Each term uses Bernoulli numbers which are recursively computed and are intricately linked with some really cool functions like the Riemann Zeta function. For comparison (thanks to Mathworld for the table):
##\cos(x)=\sum\limits_{k=0}^{\infty}{\frac{(1)^{k}}{(2k)!}x^{2k}}## (simple, cute, compact)
##\sin(x)=\sum\limits_{k=0}^{\infty}{\frac{(1)^{k}}{(2k+1)!}x^{2k+1}}## (simple, cute, compact)
##\tan(x)=\sum\limits_{k=0}^{\infty}{\frac{(1)^{k}4^{k+1}\left(4^{k+1}1\right)B_{2k+2}}{(2k+2)!}x^{2k+1}}## (wtf are math)
That's what it looks like to divide sine by cosine. Not pretty (well, okay, it is actually pretty gorgeous, but whatevs). So the first few terms of tan(x) are:
##\tan(x)=x+\frac{x^{3}}{3}+\frac{2x^{5}}{15}+\frac{17x^{7}}{315}+\frac{62x^{9}}{2835}+\frac{1382x^{11}}{155925}+...##
Nasty gorgeous, but applying Horner's Method, we end up with:
##\tan(x)=x(1+\frac{x^{2}}{3}(1+\frac{2x^{2}}{5}(1+\frac{17x^{2}}{42}(1+\frac{62x^{2}}{153}(1+...##
So seeing this, I realized, "wait a minute, those coefficients sure don't look like they are converging to zero or one, but something in between!" And that made me think about how ##\frac{1}{1ax^{2}}=(1+ax^{2}(1+ax^{2}(1+ax^{2}(1+ax^{2}(1+...##. Now, I am a mathematician, but I am also a programmer so what I am about to say is blasphemous and I am sorry, but I decided to numerically evaluate and assume those constants converge. I ended up with approximately .405285. What I am saying is:
##\tan(x)\approx x(1+\frac{x^{2}}{3}(1+\frac{2x^{2}}{5}(1+\frac{17x^{2}}{42}(1+\frac{62x^{2}}{153}(1+.405285x^{2}(1+.405285x^{2}(1+.405285x^{2}(1+...##
So then I ended up with:
##\tan(x)\approx x(1+\frac{x^{2}}{3}(1+\frac{2x^{2}}{5}(1+\frac{17x^{2}}{42}(1+\frac{62x^{2}}{153}\frac{1}{1.405285x^{2}}))))##
##\tan(x)\approx x(1+\frac{x^{2}}{3}(1+\frac{2x^{2}}{5}(1+\frac{17x^{2}}{42}(1+\frac{62x^{2}}{15362x^{2}}))))##
##\tan(x)\approx x(1+\frac{x^{2}}{3}(1+\frac{2x^{2}}{5}(1+\frac{17x^{2}}{42}\frac{153}{15362x^{2}})))##It's still a bit of work to compute, but I had fun! No to do the mathy thing and actually prove the constants converge and maybe even find it's exact value!

Good job! Do you think this can have applications in z80 programming (is it faster, more memory efficient, etc.)?

I think that if you really want speed, you should use a LUT for the sine function, calculate the cosine as sine+90°, and calculate the tangent as sine/cosine. This does require a bit of memory and isn't really related to math though.
Related: I previously had no idea how the sine, cosine and tangent were actually calculated, and good job on simplifying it so much!
EDIT: while the z80 will probably not be able to calculate that at high speed for use in games, I guess it might be usefull for CASrelated programs.

This is really cool! I found using ries (http://mrob.com/pub/ries/ries.php?target=.405285&rst=) that the constant is arctan(pi)/pi (if it turns out to be a constant, which seems to be the case).
EDIT: that is probably not correct. I just saw that first and it had a tangent in it. Looking at the others, 4/(pi^2) looks like it could be promising. But those are all approximations anyway.

Wow this topic is really interesting!
It's pretty funny that I'm reading this now, because it's Mathday at school today :)

Good job! Do you think this can have applications in z80 programming (is it faster, more memory efficient, etc.)?
This would not be the best way to compute tangent (though it is decent). For 21 bits of accuracy, there are probably much better ways :P It's just neat and if you've never seen Horner's method, it can get you thinking about math programming in the future.
I think that if you really want speed, you should use a LUT for the sine function, calculate the cosine as sine+90°, and calculate the tangent as sine/cosine. This does require a bit of memory and isn't really related to math though.
Related: I previously had no idea how the sine, cosine and tangent were actually calculated, and good job on simplifying it so much!
EDIT: while the z80 will probably not be able to calculate that at high speed for use in games, I guess it might be usefull for CASrelated programs.
Yeah, LUT's are the fastest method. There are fairly efficient ways to use small LUTs to compute the fun functions (trig and exponentials, among others). Glad to show you something new!
This is really cool! I found using ries (http://mrob.com/pub/ries/ries.php?target=.405285&rst=) that the constant is arctan(pi)/pi (if it turns out to be a constant, which seems to be the case).
EDIT: that is probably not correct. I just saw that first and it had a tangent in it. Looking at the others, 4/(pi^2) looks like it could be promising. But those are all approximations anyway.
Yeah, 4/pi^{2} looks most promising, but it might not converge to anything nice. Then again, I have lots of identities that I've derived about the Bernoulli numbers, so who knows! For example, if you sum them up to get 'C', then do c/(c1), you get euler's number, and if you sum the absolute value of the terms, you get a function involving cotangent!Wow this topic is really interesting!
It's pretty funny that I'm reading this now, because it's Mathday at school today :)
Math day is best day.
I worked on this a little last night, and I will probably work on this some more. I might actually figure out the constant and I would be more inclined to believe that it has something to do with e.

The best part about being friends with mathematicians a gazillion times smarter than me? They look at a problem and are like, "why don't you just use stuff you learned in Calc 2?" .__. Anyways, I figured out the main stuff, but one of my professor friends made the final, "obvious" observation. The constant does go to 4/pi^{2} and here is how it works:
Say we have a Maclaurin series of the form ##\sum\limits_{k=0}^{\infty}{c_{k}x^{2k+1}}=c_{0}x+c_{1}x^{3}+c_{2}x^{5}+...##. Then applying Horner's Method:
##c_{0}x+c_{1}x^{3}+c_{2}x^{5}+c_{3}x^{7}+c_{4}x^{9}+...##
##=c_{0}x(1+\frac{c_{1}}{c_{0}}x^{2}+\frac{c_{2}}{c_{0}}x^{4}+\frac{c_{3}}{c_{0}}x^{6}+\frac{c_{4}}{c_{0}}x^{8}+...##
##=c_{0}x(1+\frac{c_{1}}{c_{0}}x^{2}(1+\frac{c_{2}c_{0}}{c_{0}c_{1}}x^{2}+\frac{c_{3}c_{0}}{c_{0}c_{1}}x^{4}+\frac{c_{4}c_{0}}{c_{0}c_{1}}x^{6}+...##
##=c_{0}x(1+\frac{c_{1}}{c_{0}}x^{2}(1+\frac{c_{2}}{c_{1}}x^{2}+\frac{c_{3}}{c_{1}}x^{4}+\frac{c_{4}}{c_{1}}x^{6}+...##
...
##=c_{0}x(1+\frac{c_{1}}{c_{0}}x^{2}(1+\frac{c_{2}}{c_{1}}x^{2}(1+\frac{c_{3}}{c_{2}}x^{2}(1+\frac{c_{4}}{c_{3}}x^{2}(1+...##
Here is the last piece of the puzzle:
Me: So I know this converges because I am basically using the ratio test of the Maclaurin Series for tangent, which I already know converges.
Dr. B: Well Zeda, you know that the limit of the ratio test is inversely related to the radius of convergence, and you know that the radius of converges for tangent about 0 is pi/2 before it gets wonky.
So if we were finding the radius of convergence, we would want to get r in the following:
##\lim\limits_{n\rightarrow\infty}{\frac{c_{n+1}(xr)^{2n+3}}{c_{n}(xr)^{2n+1}}}=1##
However, we know r, and we know it is centered about x=0 (because we are using a Maclaurin series), so we get:
##\lim\limits_{n\rightarrow\infty}{\frac{c_{n+1}(0r)^{2}}{c_{n}}}=1##
##\lim\limits_{n\rightarrow\infty}{r^{2}\frac{c_{n+1}}{c_{n}}}=1##
##\lim\limits_{n\rightarrow\infty}{\frac{c_{n+1}}{c_{n}}}=\frac{1}{r^{2}}##
Also, it happens to be that all the coefficients in the Maclaurin series of tan(x) are positive, so getting rid of the absolute value:
##\lim\limits_{n\rightarrow\infty}{\frac{c_{n+1}}{c_{n}}}=\frac{4}{\pi^{2}}##
Yay :D For extra fun, instead of using the shorthand c_{n} to refer to the coefficients, let's plug in actual values:
##\lim\limits_{n\rightarrow\infty}{\frac{B_{2n+4}4^{n+2}(4^{n+2}1)(2n+2)!}{B_{2n+2}4^{n+1}(4^{n+1}1)(2n+4)!}}=\frac{4}{\pi^{2}}##
##\lim\limits_{n\rightarrow\infty}{\frac{B_{2n+4}4(4^{n+2}1)(2n+2)!}{B_{2n+2}(4^{n+1}1)(2n+4)!}}=\frac{4}{\pi^{2}}##
##\lim\limits_{n\rightarrow\infty}{4\frac{B_{2n+4}(4^{n+2}1)}{B_{2n+2}(4^{n+1}1)(2n+3)(2n+4)}}=\frac{4}{\pi^{2}}##
##\lim\limits_{n\rightarrow\infty}{4\frac{B_{2n+4}4}{B_{2n+2}(2n+3)(2n+4)}}=\frac{4}{\pi^{2}}##
##\lim\limits_{n\rightarrow\infty}{16\frac{B_{2n+4}}{B_{2n+2}(2n+3)(2n+4)}}=\frac{4}{\pi^{2}}##
##\lim\limits_{n\rightarrow\infty}{\frac{B_{2n+4}}{B_{2n+2}(2n+3)(2n+4)}}=\frac{1}{4\pi^{2}}##
##\lim\limits_{n\rightarrow\infty}{\frac{B_{2n+2}(2n+3)(2n+4)}{B_{2n+4}}}=4\pi^{2}##
##\lim\limits_{n\rightarrow\infty}{\frac{B_{2n}(2n+1)(2n+2)}{B_{2n+2}}}=4\pi^{2}##
From this you can get that ##B_{2n+2}\approx \frac{B_{2n}(2n+1)(2n+2)}{4\pi^{2}}##. And further:
##B_{2n+4}\approx \frac{B_{2n+2}(2n+3)(2n+4)}{(2\pi)^{2}}\approx \frac{B_{2n}(2n+1)(2n+2)(2n+3)(2n+4)}{(2\pi)^{4}}##
And in general:
##B_{2n+2m}\approx \frac{B_{2n}(2n+2m)!(1)^{m}}{(2\pi)^{2m}(2n)!}##
So let's put that into practice. Say we know B_{20}=174611/330 and we want to estimate B_{30}:
##B_{20+2\cdot 5}\approx \frac{B_{20}(30)!(1)^{5}}{(2\pi)^{10}(20)!}##
##= \frac{174611/330(30)!}{(2\pi)^{10}(20)!}##
##= \frac{55016531182.1502685546875}{\pi^{10}}##
##B_{30}\approx 601581447.225687068979178849893094649083608778571...## according to Wolfram Alpha. Meanwhile the actual value:
##B_{30}\approx 601580873.900642368384303868174835916771400642368...##
This is a useful identity since the Bernoulli numbers are used in many applications, but they are tough to compute efficiently. I am pretty sure I have seen this identity before, and I've actually hypothesized it earlier in the year (in one of my previous notebooks), but it is awesome to actually prove it!