一些C实现的数学函数实现（估算）

最新推荐文章于 2021-07-13 18:02:09 发布

weixin_34162695

最新推荐文章于 2021-07-13 18:02:09 发布

阅读量209

点赞数

文章标签：嵌入式 python

原文链接：https://my.oschina.net/chaosannals/blog/668479

版权

2019独角兽企业重金招聘Python工程师标准>>>

之前搞过一段时间的单片机，当时用的芯片只有16M的空间用于写代码。用到了asin、acos、atan这些反三角函数，但是引用C标准库的时候发现，每多用一个这样的标准函数，空间多占用1M左右，也就是用到这3个函数的话就要占用3M，如果用到sin、cos这些，那么没有什么空间写其他代码了。所以就需要自己实现这些函数。

标准库的函数八成是用的查表法，而且表可能过大，我没有去看标准库的实现代码。那么我自己实现的函数和标准库的比起来有什么差别呢？其实就是容量小，但是误差大些，不过在我的应用里，就小数点好几位后的误差基本可以忽略。

0.实现Sqrt（平方根）和RSqrt函数（平方根倒数）

#include <stdint.h>
typedef union{
    float f;
    int32_t i;
}FIUnion;

float RSqrt(float value){
    size_t i;
    FIUnion buffer;
    buffer.f = value;
    buffer.i = 0x5f3759df - (buffer.i >> 1);
    for(i = 0;i != 4;++i)
        buffer.f *= (3.0f - value*buffer.f*buffer.f)*0.5f;
    return buffer.f;
}

float Sqrt(float value){
    return 1.0f / RSqrt(value);
}

这个算法在《算法心得》(Hacker's Delight)里“17.4 估算平方根”有提及，网上关于卡马克的传说里也有，大体上可以理解为从牛顿法演变过来的算法。这个例子做了4次迭代，精度还行，过多的迭代，单片机受不了。

1.通过CORDIC（坐标旋转数字计算方法）求反三角函数

#include <stddef.h>
#define MPI       3141593//π*1000000取整 
#define MUL       0.000001f//
#define NUM       21

float ATan2(float y,float x){
    static const int radin[NUM] = {
        785398,  463648,  244979,  124355,
        62419,  31240,  15624,  7812,
        3906,  1953,  977,  488,
        244,  122,  61,   31,
        15,   8,   4,   2,  1};
    size_t i;
    int tx,ty;
    int ix = x*1024.0f,iy = y*1024.0f;
    int result = 0;
    if(ix < 0){
        result = iy >= 0 ? MPI : -MPI;
        ix = -ix;
        iy = -iy;
    }
    for(i = 0;i != NUM;++i){
        if(iy == 0)break;
        if(iy > 0){
            tx = ix + (iy >> i);
            ty = iy - (ix >> i);
            result += radin[i];
        }else{
            tx = ix - (iy >> i);
            ty = iy + (ix >> i);
            result -= radin[i];
        }
        ix = tx;
        iy = ty;
    }
    return result * MUL;
}
float ASin(float value){
    float y = value;
    float x = 1.0f / RSqrt(1.0f - y*y);
    return ATan2(y,x);
}

float ACos(float value){
    float x = value;
    float y = 1.0f / RSqrt(1.0f - x*x);
    return ATan2(y,x);
}

float ATan(float value){
    return ATan2(value,1.0f);
}

这里用到了上面例子里的RSqrt，用于ASin和ACos。首先说说ATan2吧，其实就是传入“对边”和“邻边”的长，这样在一个坐标系下我们就可以认为这个角为向量(“对边”，“邻边”)与X正半轴的夹角。求解的过程也是一步步逼近，最后得到一个近似的估值。

基本的思路：首先是旋转这条边到一四象限，并累加旋转的角度。之后就是不停地旋转使其靠近X正半轴并累加旋转的角度。

radin数组是tan值为1,1/2,1/4....的角（这些值事先按计算器求得：) ）一直二分下去对应的弧度值的1000000倍（由于只有1000000倍，所以分到第21次时，是1了分不下去，要更高精度的话请用更大倍数）。用整数可以减少浮点计算，故意用了1000000倍来容纳小数点后的几位，因为用到的弧度值基本小于1。最后把结果乘以0.000001了。

至于 tx = ix + (iy >> i); ty = iy - (ix >> i);这种是坐标旋转的变体（x'/cos = x + y*tan; y'/cos = y - x*tan），无视向量的长度，除以cos来减少计算。而旋转的角特地选tan为1,1/2,1/4的，所以就相当于除2的i次方，而这种除2的次方又可以用右移代替。

最后的ATan是由于tan = y / x（“对”除“临”），所以当X = 1时，Y = tan 了，所以直接就这么用。写的时候就这思路，至于弊端的话，我发现误差会比ATan2大些。而且还是调用的ATan2，所以性能也差些。