一道编程面试题引出的，用泰勒级数展开式自行实现C语言的atan反正切函数
0x00 缘起
今天别人问我一道编程题：已知向量a=x+yi,编写程序实现a转换为a=A∠a形式，要求：不能用库函数，纯代码实现。
毕业良久，匮乏的数学知识都随岁月而逝了，于是一通网上搜索。
发现这里涉及到几个知识点：
1. 向量的表示方法，实际上最后还是使用了复数的幅角表示法。实际上原理是一样的，只是表述方式不同。
2. math.h 数学函数库中的几个函数实现。sqrt ，pow，atan
0x01 向量
这题的核心实际上是数学的公式的实现，于是一通搜索。
找到了这个：
还有这个：
基本上可以明确 向量转几何意义，再转成角度计算和向量长度的计算。对应原题中，长度就是A 角度就是α
那么列出程序如下：
typedef struct { double amp; // 幅值 double angle; // 弧度 } stRet_T; stRet_T convert(double x, double y) { stRet_T ret; ret.amp = sqrt(x*x + y*y); ret.angle = ( atan(y/x) / M_PI ) * 180; return ret; }
0x02 库函数实现
题目中有要求不能用库函数，因此，需要手动实现库函数，已知需要实现的函数 sqrt 和 atan
开方函数的实现很容易，
float msqrt(float a) { double x,y; x=0.0; y=a/2; while(x!=y) { x=y; y=(x+a/x)/2; } return x; }
反正切函数怎么实现呢？
找了一通网上的资源，有说用cordic算法的，有说用查表法计算的。最后测试了代码都不对。
返璞归真，就抓一下arctan的数学形式，于是找到了arctan转成泰勒级数展开式，泰勒级数就简单了，求和处理。
于是，atan实现如下：
float calc(float x, int n) { float a = mpow(-1, n); float c = 2*n+1; float b = mpow(x, c); // printf("n: %d, a: %f, b: %f, c: %f\n", n, a, b, c); return (a*b/c); } float matan(float x, int level) { float ret = 0; int i = 0; for (i=0; i<level; i++) { ret += calc(x, i); // printf("lev[%d]: %f\n", i, calc(x,i)); } return ret; }
完成后，测试发现：级数越少误差越高。需要级数到10000才能和系统库的值大体一致。测试如下：
printf("atanx %f, mmatanx %f, level: %d\n", atan(1), mmatan(1,10), 10); printf("atanx %f, mmatanx %f, level: %d\n", atan(1), mmatan(1,100), 100); printf("atanx %f, mmatanx %f, level: %d\n", atan(1), mmatan(1,1000), 1000); printf("atanx %f, mmatanx %f, level: %d\n", atan(1), mmatan(1,10000), 10000); // atanx 0.785398, mmatanx 0.760460, level: 10 // atanx 0.785398, mmatanx 0.782898, level: 100 // atanx 0.785398, mmatanx 0.785148, level: 1000 // atanx 0.785398, mmatanx 0.785375, level: 10000
0x03 总结
数学简直是一门艺术，把复杂的变换轻松化解成为了加减乘除。果然加减乘除才是科学密码。
数学YYDS，泰勒级数YYDS
完整代码如下：
1 #include <stdio.h> 2 #include <math.h> 3 4 typedef struct { 5 double amp; // 幅值 6 double angle; // 弧度 7 } stRet_T; 8 /** 9 * 10 * @param x 实部 11 * @param y 虚部 12 * @return 13 * amp 幅值 √(x²+y²) 14 * angle 幅度转角度 arctan(y/x) 15 */ 16 17 stRet_T convert(double x, double y) 18 { 19 stRet_T ret; 20 21 ret.amp = sqrt(x*x + y*y); 22 ret.angle = ( atan(y/x) / M_PI ) * 180; 23 24 return ret; 25 } 26 27 float msqrt(float a) 28 { 29 double x,y; 30 x=0.0; 31 y=a/2; 32 while(x!=y) 33 { 34 x=y; 35 y=(x+a/x)/2; 36 } 37 return x; 38 } 39 40 41 float mpow(float a,int b) 42 { 43 float r=a; 44 if(b>0) 45 { 46 while(--b) 47 r*=a; 48 49 } 50 else if(b<0) 51 { 52 while(++b) r*=a; 53 r=1.0/r; 54 } 55 else r=1; 56 return r; 57 } 58 59 float calc(float x, int n) 60 { 61 float a = mpow(-1, n); 62 float c = 2*n+1; 63 float b = mpow(x, c); 64 // printf("n: %d, a: %f, b: %f, c: %f\n", n, a, b, c); 65 66 return (a*b/c); 67 } 68 69 float mmatan(float x, int level) 70 { 71 float ret = 0; 72 int i = 0; 73 for (i=0; i<level; i++) { 74 ret += calc(x, i); 75 // printf("lev[%d]: %f\n", i, calc(x,i)); 76 } 77 return ret; 78 } 79 80 81 stRet_T mconvert(double x, double y) 82 { 83 stRet_T ret; 84 85 ret.amp = msqrt(x*x + y*y); 86 ret.angle = ( mmatan(y/x, 1000) / M_PI ) * 180; 87 88 return ret; 89 } 90 91 int main() { 92 93 printf("atanx %f, mmatanx %f, level: %d\n", atan(1), mmatan(1,10), 10); 94 printf("atanx %f, mmatanx %f, level: %d\n", atan(1), mmatan(1,100), 100); 95 printf("atanx %f, mmatanx %f, level: %d\n", atan(1), mmatan(1,1000), 1000); 96 printf("atanx %f, mmatanx %f, level: %d\n", atan(1), mmatan(1,10000), 10000); 97 98 double x = 110* msqrt(3); 99 double y = 110; 100 101 stRet_T r = convert(x,y); 102 printf("lib: amp: %f, rad: %f\n", r.amp, r.angle); 103 104 stRet_T r2 = mconvert(x,y); 105 printf("mine: amp: %f, rad: %f\n", r2.amp, r2.angle); 106 107 return 0; 108 109 110 // 220∠30°=220cos30°+j220sin30°=110√3+110j 111 112 return 0; 113 }
测试效果如下：
atanx 0.785398, mmatanx 0.760460, level: 10 atanx 0.785398, mmatanx 0.782898, level: 100 atanx 0.785398, mmatanx 0.785148, level: 1000 atanx 0.785398, mmatanx 0.785375, level: 10000 lib: amp: 220.000000, rad: 30.000000 mine: amp: 220.000000, rad: 30.000004
参考文献
[1]. 怎样把相量形式转换代数形式？_百度知道
