Skip to content

Instantly share code, notes, and snippets.

@pmslavin

pmslavin/atan.py

Created September 5, 2017 20:28
Show Gist options
  • Save pmslavin/fc29d0a9169666155de514e69bc509f0 to your computer and use it in GitHub Desktop.
Save pmslavin/fc29d0a9169666155de514e69bc509f0 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
atan.py
def atan_ts(z):
""" ___ oo (-1)^2n+1 * x^(2n+1)
atan(z) = \ --------------------
/__ n=0 2n+1
"""
a = 0
for n in range(17):
q = 2.0*n+1
a += (-1)**n * z**q / q
return a
def atan_id(z):
"""
atan(x) = atan(c) + atan((x-c)/(1+xc)) where c is small
Use series for arguments near centre of convergence interval
"""
if z<= 0.2:
return atan_ts(z)
ac = atan_ts(0.2)
return ac + atan_ts((z-0.2)/(1.0+0.2*z))
def octantify(x,y):
"""
Rotates coordinates to lie within first octant
"""
o = 0
t = 0
if y<0:
pass ; x = -x ; y = -y ; o += 4
if x<=0:
t = x ; x = y ; y = -t ; o += 2
if x<y:
t = y-x ; x = x+y ; t = t ; o += 1
return x,y,o
PI = 3.1415926536
BRAD_PI = 1<<14
ATAN_ONE = 0x1000
LUT_STRIDE = ATAN_ONE/128
lut = [ 0x0000,0x0146,0x028C,0x03D2,0x0517,0x065D,0x07A2,0x08E7,
0x0A2C,0x0B71,0x0CB5,0x0DF9,0x0F3C,0x107F,0x11C1,0x1303,
0x1444,0x1585,0x16C5,0x1804,0x1943,0x1A80,0x1BBD,0x1CFA,
0x1E35,0x1F6F,0x20A9,0x21E1,0x2319,0x2450,0x2585,0x26BA,
0x27ED,0x291F,0x2A50,0x2B80,0x2CAF,0x2DDC,0x2F08,0x3033,
0x315D,0x3285,0x33AC,0x34D2,0x35F6,0x3719,0x383A,0x395A,
0x3A78,0x3B95,0x3CB1,0x3DCB,0x3EE4,0x3FFB,0x4110,0x4224,
0x4336,0x4447,0x4556,0x4664,0x4770,0x487A,0x4983,0x4A8B,
# 64
0x4B90,0x4C94,0x4D96,0x4E97,0x4F96,0x5093,0x518F,0x5289,
0x5382,0x5478,0x556E,0x5661,0x5753,0x5843,0x5932,0x5A1E,
0x5B0A,0x5BF3,0x5CDB,0x5DC1,0x5EA6,0x5F89,0x606A,0x614A,
0x6228,0x6305,0x63E0,0x64B9,0x6591,0x6667,0x673B,0x680E,
0x68E0,0x69B0,0x6A7E,0x6B4B,0x6C16,0x6CDF,0x6DA8,0x6E6E,
0x6F33,0x6FF7,0x70B9,0x717A,0x7239,0x72F6,0x73B3,0x746D,
0x7527,0x75DF,0x7695,0x774A,0x77FE,0x78B0,0x7961,0x7A10,
0x7ABF,0x7B6B,0x7C17,0x7CC1,0x7D6A,0x7E11,0x7EB7,0x7F5C,
# 128
0x8000,0x80A2
]
def atan2_lut(y,x):
"""
Direct lookup of values from table
"""
x,y,phi = octantify(x,y)
phi *= BRAD_PI/4;
t = (y << 12)/x
return phi + lut[t/LUT_STRIDE]/8 / float(BRAD_PI) * PI
def atan2_luli(y,x):
"""
Lookup with linear interpolation within interval
"""
x,y,phi = octantify(x,y)
phi *= BRAD_PI/4;
t = (y << 12)/x
h = t%LUT_STRIDE
fa = lut[t/LUT_STRIDE]
fb = lut[t/LUT_STRIDE+1]
luli = phi + (fa + ((fb-fa)*h >> 5))/8
return PI * luli/BRAD_PI
if __name__ == "__main__":
print("atan_ts(0.1) = {0}".format(atan_ts(0.1)))
print("atan_ts(0.2) = {0}".format(atan_ts(0.2)))
print("atan_ts(0.5) = {0}".format(atan_ts(0.5)))
print("atan_ts(1.0) = {0}".format(atan_ts(1.0)))
print("")
print("atan_id(0.1) = {0}".format(atan_id(0.1)))
print("atan_id(0.2) = {0}".format(atan_id(0.2)))
print("atan_id(0.5) = {0}".format(atan_id(0.5)))
print("atan_id(1.0) = {0}".format(atan_id(1.0)))
print("")
print("atan2_lut(0.1) = {0}".format(atan2_lut(1,10)))
print("atan2_lut(0.2) = {0}".format(atan2_lut(1,5)))
print("atan2_lut(0.5) = {0}".format(atan2_lut(1,2)))
print("atan2_lut(1.0) = {0}".format(atan2_lut(1,1)))
print("")
print("atan2_luli(0.1) = {0}".format(atan2_luli(1,10)))
print("atan2_luli(0.2) = {0}".format(atan2_luli(1,5)))
print("atan2_luli(0.5) = {0}".format(atan2_luli(1,2)))
print("atan2_luli(1.0) = {0}".format(atan2_luli(1,1)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment