add SIN and ARCTAN@ approximation algorithm

add SIN and ARCTAN@ approximation algorithm

ARCTAN/arctan_approx.c

@@ -1,54 +1,78 @@
-#include <stdio.h>
-#include <math.h>
+////////////////////////
+// 近似arctan2计算方法
+////////////////////////
 
-#define MAX2(a,b) (((a)>(b))?(a):(b))
+#define PI 3.1415926535897932384626433832795
+#define SIGN(x) (((x)>=0)?((x)>0):-1)
+#define ABS(x) ((x)>0?(x):(-(x)))
 
-// 几种近似tanh计算方法
-// 代码清单 2 - 12
-float tanh_1(float x)
+float algo1(float i, float q) { return i*q/(i*i+0.28125f*q*q); }
+float algo2(float i, float q) { return i*q/(i*i+0.28086f*q*q); }
+float algo3(float i, float q)
 {
-    if (x < -3) return -1;
-    if (x >  3) return  1;
-    return x*(27.0f+x*x)/(27.0f+9.0f*x*x);
+    float x=q/i;
+    return (float)(PI/4.0)*x+0.285f*x*(1.0f-ABS(x));
 }
-
-float tanh_2(float x)
+float algo4(float i, float q)
 {
-    if (x>(float) 3.4f) return  1;
-    if (x<(float)-3.4f) return -1;
-
-    x*= 1.0f/3.4f;
-    x*= fabsf(x)-2.0f;
-    return x*(fabsf(x)-2.0f);
+    float x=q/i;
+    return (float)(PI/4.0)*x+0.273f*x*(1.0f-ABS(x));
 }
-
-float tanh_3(float x)
+float algo5(float i, float q)
 {
-    if (x>(float) 3.44540506f) return  1;
-    if (x<(float)-3.44540506f) return -1;
-
-    x*= 1.0f/3.44540506f;
-    x*= fabsf(x)-2.00002276f;
-    return x*(fabsf(x)-2.0000286f);
+    float x=q/i;
+    return (float)(PI/4.0)*x+x*(0.186982f-0.191942f*x*x);
+}
+float algo6(float i, float q)
+{
+    float x=q/i;
+    return (float)(PI/4.0)*x-x*(ABS(x)-1.0f)*(0.2447f+0.0663f*ABS(x));
 }
 
-int main()
+typedef float (*algo_t)(float, float);
+
+// 计算ARCTAN近似，使用不同的算法
+// 输出数据范围是-PI~+PI
+float arctan_approx(float i, float q, algo_t algo)
 {
-    float e1=0, e2=0, e3=0, y0, y1, y2, y3;
+    if (ABS(i)>ABS(q))
+        return (i>0)? algo(i,q):SIGN(q)*3.141592654f+algo(i,q);
+    else
+        return (q>0)? (float)(PI/2.0)-algo(q,i):-(float)(PI/2.0)-algo(q,i);
+}
 
-    for (float x=-3.5f; x<3.5f; x+=1e-3f)
+////////////////////////
+// 单元测试
+////////////////////////
+
+#include <stdio.h>
+#include <math.h>
+#define N 10000
+#define MAX(a,b) (((a)>(b))?(a):(b))
+
+// 测试近似算法误差
+float test_algo(algo_t algo)
+{
+    float t,ta,err=0;
+    for (int n=0; n<N; n++) 
     {
-        y0=(float)tanh(x);
-        y1=tanh_1(x);
-        y2=tanh_2(x);
-        y3=tanh_3(x);
-
-        e1=MAX2(e1,fabsf(y0-y1));
-        e2=MAX2(e2,fabsf(y0-y2));
-        e3=MAX2(e3,fabsf(y0-y3));
+        t=(float)(n-N/2)/(float)(N+1)*(float)M_PI*2.0f;
+        ta=arctan_approx(cos(t),sin(t),algo);
+        err=MAX(err,fabs(t-ta));
     }
-    printf("[INF] tanh_1, difference: %e\n",e1);
-    printf("[INF] tanh_2, difference: %e\n",e2);
-    printf("[INF] tanh_3, difference: %e\n",e3);
+    return err;
+}
+
+#define RAD2DEG(x) ((x)*(float)(180.0/M_PI))
+
+int main()
+{
+    printf("algo1 err. max: %f(deg)\n", RAD2DEG(test_algo(algo1)));
+    printf("algo2 err. max: %f(deg)\n", RAD2DEG(test_algo(algo2)));
+    printf("algo3 err. max: %f(deg)\n", RAD2DEG(test_algo(algo3)));
+    printf("algo4 err. max: %f(deg)\n", RAD2DEG(test_algo(algo4)));
+    printf("algo5 err. max: %f(deg)\n", RAD2DEG(test_algo(algo5)));
+    printf("algo6 err. max: %f(deg)\n", RAD2DEG(test_algo(algo6)));
+
     return 0;
 }

ARCTAN/arctan_approx.py

@@ -31,21 +31,31 @@ def arctan2(i,q,algo):
 # 单元测试
 ####################
 if __name__ == '__main__':
-    t=np.linspace(-np.pi,np.pi,10001)[:-1]
+    N=10000
+    t=np.linspace(-np.pi,np.pi,N+1)[:-1]
+
     t1=np.array([arctan2(i,q,algo1) for i,q in zip(np.cos(t),np.sin(t))])
     t2=np.array([arctan2(i,q,algo2) for i,q in zip(np.cos(t),np.sin(t))])
     t3=np.array([arctan2(i,q,algo3) for i,q in zip(np.cos(t),np.sin(t))])
     t4=np.array([arctan2(i,q,algo4) for i,q in zip(np.cos(t),np.sin(t))])
     t5=np.array([arctan2(i,q,algo5) for i,q in zip(np.cos(t),np.sin(t))])
     t6=np.array([arctan2(i,q,algo6) for i,q in zip(np.cos(t),np.sin(t))])
-plt.clf()
-plt.plot(np.rad2deg(t),np.rad2deg(t-t1))
-plt.plot(np.rad2deg(t),np.rad2deg(t-t2))
-plt.plot(np.rad2deg(t),np.rad2deg(t-t3))
-plt.plot(np.rad2deg(t),np.rad2deg(t-t4))
-plt.plot(np.rad2deg(t),np.rad2deg(t-t5))
-plt.plot(np.rad2deg(t),np.rad2deg(t-t6))
 
-plt.legend(['err. algo. 1','err. algo. 2','err. algo. 3','err. algo. 4','err. algo.5','err. algo. 6'])
-plt.grid(True)
-plt.show()
+    print('algo1 err. %f (max: %f)'%(np.mean(np.abs(np.rad2deg(t-t1))),np.max(np.abs(np.rad2deg(t-t1)))))
+    print('algo2 err. %f (max: %f)'%(np.mean(np.abs(np.rad2deg(t-t2))),np.max(np.abs(np.rad2deg(t-t2)))))
+    print('algo3 err. %f (max: %f)'%(np.mean(np.abs(np.rad2deg(t-t3))),np.max(np.abs(np.rad2deg(t-t3)))))
+    print('algo4 err. %f (max: %f)'%(np.mean(np.abs(np.rad2deg(t-t4))),np.max(np.abs(np.rad2deg(t-t4)))))
+    print('algo5 err. %f (max: %f)'%(np.mean(np.abs(np.rad2deg(t-t5))),np.max(np.abs(np.rad2deg(t-t5)))))
+    print('algo6 err. %f (max: %f)'%(np.mean(np.abs(np.rad2deg(t-t6))),np.max(np.abs(np.rad2deg(t-t6)))))
+
+    plt.clf()
+    plt.plot(np.rad2deg(t),np.rad2deg(t-t1))
+    plt.plot(np.rad2deg(t),np.rad2deg(t-t2))
+    plt.plot(np.rad2deg(t),np.rad2deg(t-t3))
+    plt.plot(np.rad2deg(t),np.rad2deg(t-t4))
+    plt.plot(np.rad2deg(t),np.rad2deg(t-t5))
+    plt.plot(np.rad2deg(t),np.rad2deg(t-t6))
+
+    plt.legend(['err. algo. 1','err. algo. 2','err. algo. 3','err. algo. 4','err. algo. 5','err. algo. 6'])
+    plt.grid(True)
+    plt.show()

ARCTAN/build.py

@@ -1,3 +1,4 @@
 import os
-os.system('gcc -o ./artan_approx arctan_approx.c -std=c99')
-os.system('arctan_approx')
+os.system('gcc arctan_approx.c -Wall -lm -o test')
+os.system('test')
+

ARCTAN/readme.md

@@ -1,6 +1,6 @@
-# arctan的近似计算方法
+# arctan2的近似计算方法
 1. 分别用python和C实现
 2. C程序使用命令：
-gcc -o ./arctan_approx arctan_approx.c -std=c99 
-编译，编译结果是arctan_approx(*.exe)，运行它得到近似arctan计算方法误差
-(在build.py可也可以完成编译和测试)
+gcc arctan_approx.c -Wall -lm -o test
+编译，编译结果是test，运行它得到近似arctan计算方法误差
+(在build.py可也可以完成编译和测试)

SIN/build.py

@@ -0,0 +1,4 @@
+import os
+os.system('gcc sin_approx.c -Wall -lm -o test')
+os.system('test')
+

SIN/readme.md

@@ -0,0 +1,6 @@
+# sin的近似计算方法
+1. 分别用python和C实现
+2. C程序使用命令：
+gcc sin_approx.c -Wall -lm -o test
+编译，编译结果是test，运行它得到近似sin计算方法误差
+(在build.py可也可以完成编译和测试)

SIN/sin_approx.c

@@ -0,0 +1,72 @@
+////////////////////////
+// 近似sin计算方法
+////////////////////////
+
+#define PI 3.1415926535897932384626433832795
+#define SIGN(x) (((x)>=0)?((x)>0):-1)
+#define ABS(x) ((x)>0?(x):(-(x)))
+
+float algo1(float x)
+{
+    float y=x/(float)PI;
+    return 4.0f*y*(1.0f-y);
+}
+
+float algo2(float x)
+{
+//    return -0.417698f*x*x+1.312236f*x-0.050465f;
+    return (float)(60.0*(PI*PI-12.0)/(PI*PI*PI*PI*PI))*x*x-
+           (float)(60.0*(PI*PI-12.0)/(PI*PI*PI*PI))*x+
+           (float)(12.0*(PI*PI-10.0)/(PI*PI*PI));
+}
+
+float algo3(float x)
+{
+    return 16.0f*x*((float)PI-x)/((float)(5.0*PI*PI)-4.0f*x*((float)PI-x));
+}
+
+typedef float (*algo_t)(float);
+
+// 计算ARCTAN近似，使用不同的算法
+// 输入数据范围是-PI~+PI
+float sin_approx(float x, algo_t algo)
+{
+    if ((x>(float)(PI*2.0f)) || (x<-(float)(PI*2.0f)))
+    {
+        int k=x/(float)(PI*2.0f);
+        x-=(float)k*(float)(PI*2.0);
+    }
+    if (x>PI) x-=(float)(PI*2.0);
+
+    return SIGN(x)*algo(ABS(x));
+}
+
+////////////////////////
+// 单元测试
+////////////////////////
+
+#include <math.h>
+#include <stdio.h>
+#define N 10000
+#define MAX(a,b) (((a)>(b))?(a):(b))
+
+// 测试近似算法误差
+float test_algo(algo_t algo)
+{
+    float t,err=0;
+    for (int n=0; n<N; n++)
+    {
+        t=(float)(n-N/2)/(float)(N+1)*(float)M_PI*2.0f;
+        err=MAX(err,ABS(sin_approx(t,algo)-sin(t)));
+    }
+    return err;
+}
+
+int main()
+{
+    printf("algo1 err. max: %f\n", test_algo(algo1));
+    printf("algo2 err. max: %f\n", test_algo(algo2));
+    printf("algo3 err. max: %f\n", test_algo(algo3));
+
+    return 0;
+}

SIN/sin_approx.py

@@ -0,0 +1,51 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+########################
+# 近似sin计算方法
+########################
+
+def algo1(x): 
+    y=x/np.pi
+    return 4*y*(1-y)
+
+def algo2(x): 
+    if True:
+        return -0.417698*x*x+1.312236*x-0.050465
+    else:
+        return (60*(np.pi**2-12)/np.pi**5)*x**2-\
+               (60*(np.pi**2-12)/np.pi**4)*x+\
+               (12*(np.pi**2-10)/np.pi**3)\
+
+def algo3(x):
+    return 16*x*(np.pi-x)/(5*np.pi**2-4*x*(np.pi-x))
+
+def sin_approx(x,algo):
+    x%=np.pi*2.0
+    if x>np.pi:x-=np.pi*2.0
+    return np.sign(x)*algo(np.abs(x))
+
+
+####################
+# 单元测试
+####################
+if __name__ == '__main__':
+    N=10000
+    x=np.linspace(-np.pi,np.pi,N+1)[:-1]
+    s=np.sin(x)
+
+    s1=np.array([sin_approx(x_,algo1) for x_ in x])
+    s2=np.array([sin_approx(x_,algo2) for x_ in x])
+    s3=np.array([sin_approx(x_,algo3) for x_ in x])
+
+    print('algo1 err. %f (max: %f)'%(np.mean(np.abs(s-s1)),np.max(np.abs(s-s1))))
+    print('algo2 err. %f (max: %f)'%(np.mean(np.abs(s-s2)),np.max(np.abs(s-s2))))
+    print('algo3 err. %f (max: %f)'%(np.mean(np.abs(s-s3)),np.max(np.abs(s-s3))))
+
+    plt.clf()
+    plt.plot(s-s1)
+    plt.plot(s-s2)
+    plt.plot(s-s3)
+    plt.legend(['algo. 1 err.', 'algo. 2 err.','algo. 3 err.'])
+    plt.grid(True)
+    plt.show()