Omnimaga
General Discussion => Other Discussions => Math and Science => Topic started by: Xeda112358 on October 28, 2013, 09:38:46 am
-
EDIT: To understand the method, see this document. (https://docs.google.com/document/d/1HVwAIPnpg4LS74rzoRZhHxTTAketUHWz-m-IyQQfByg/edit?usp=sharing)
For simple formulas:
arctan(x) : (9x+2x3)/(9+5x2) (at least 10 bits on [-1,1] use atan(x)=pi/2-atan(1/x) to extend to the whole real line)
ln(x+1) : x(8+x)/(8+5x) (at least 9 bits on [0,1])
A while ago, I mentioned a nice formula that I derived for estimating arctangent to within 10 bits of accuracy on [-1,1] (so within 1/1024 of the actual arctan function). I really like it because it is so simple and short, and requires few math operations for such accuracy:
(9x+2x3)/(9+5x2)
You can get away with using the following operations:
y=x*x
a=(y<<1)+9
b=(y<<2)+y+9
return a*x/b
multiply : 2
Divide : 1
add : 3
sub : 0
shift : 3
-
Wow that's some nice optimization you pulled off here Xeda. :thumbsup:
-
Indeed, pretty cool :)
(http://i.imgur.com/Zn80rTg.png)
Came from series expansion ?
(I'm too lazy to reverse-engineer the formula ... :D)
-
Actually, series expansions like Taylor series are really not good for arctangent (and ln()). Instead, I used a few terms for the continued fraction of arctangent (http://en.wikipedia.org/wiki/Inverse_trigonometric_functions#Continued_fractions_for_arctangent) and then I adjusted values to account for the error in the rest of the terms not included. So basically:
x/(1+x2/(3+4x2))
(next term in the continued fraction divides by 5+?, so I use x2 instead of 4x2)
x/(1+x2/(3+x2))
=x/((3+2x2)/(3+x2))
=x(3+x2)/(3+2x2)
=(9x+3x3)/(9+6x2)=f(x)
From there, I leave the constants the same (9 and 9) because those should be the same in the numerator and denominator (as x→0, this makes f(x)→x which is best for series expansion of arctangent around 0). Then I just adjusted the other two coefficients. It is possible that there are even better approximations using different numbers, too.
-
Ah, I see, nice :)
-
xeda, you and your dark magic..../me hides
-
That's okay, I plan to *try* to extend the precision much further than on [-1,1] and as well I want to try to make closer approximations. Currently, this is better than my CORDIC implementations for 24-bit floats.
-
Using atan(x)=pi/2-atan(1/x), it extends very nicely.
For example,
/ (9x+3x^3)/(9+6x^2) [if x <= 1]
\ pi/2 - (9x+3a^3)/(9+6a^2) (where a = 1/x) [otherwise]The biggest error is about 0.014, at x=1. See here (http://www.wolframalpha.com/input/?i=plot%7BAbs%5BArcTan%5Bx%5D+-+Piecewise%5B%7B%7B%289x%2B3x%5E3%29%2F%289%2B6x%5E2%29%2C+x+%3C%3D+1%7D%2C+%7Bpi%2F2-%289%281%2Fx%29%2B3%281%2Fx%29%5E3%29%2F%289%2B6%281%2Fx%29%5E2%29%2C+x+%3E+1%7D%7D%5D%5D%7D%2C+x+%3D+0..3)
I don't have any clever ideas to fix that, though.
-
Well, remember that (9x+3x^3)/(9+6x^2) is not the end of the modification (see the first post for the best fit dealing with 9 as the first coefficients).
To arrive to the estimation of (9x+2x3)/(9+52), you can observe that atan(1)=pi/4. Since 22/7 is a good approximation of pi, 11/14 is a good approximation of atan(1).
From here, you can guess the coefficients so that at x=1, (9+a)/(9+b)=11/24 and you get a=2, b=5
:) Using this, the biggest error is around .00065.
I am trying to get a formula of the form (ax+bx3)/(a+cx2+dx4) using the ratio of 353/452 to approximate atan(1).
-
So this has been brought up in IRC but no-one said it here yet: shouldn't a=(y<<1)+x+(x<<3) be a=(y<<1)*x+x+(x<<3) ?
-
Ah, yes, good catch :P I modified it a little more and it is actually in a better situation :)
-
Between yesterday and this morning, I have formulated a more sound process for coming up with these formulas. For more information, I made a document here. (https://docs.google.com/document/d/1HVwAIPnpg4LS74rzoRZhHxTTAketUHWz-m-IyQQfByg/edit?usp=sharing) (I am still modifying it and checking it for accuracy, otherwise, I would upload it, too).
As a result of yesterday's ventures, I managed to avoid the guess 'n check method and derive a formula for 13-bits of accuracy. It only requires one more multiplication, too:
(240x+115x3)/(240+195x2+17x4)
However, this is less suitable when using 8.8 fixed point math than the previous (then again, if you are using 8.8 fixed point, you probably won't need 13-bits of accuracy).
-
Applying the same method to ln(1+x), there is a similar formula to that of atan that is also simple and requires one less multiplication and one less shift for 9 bits of accuracy on [1,2]
EDIT: See the first post for the function.