求exp和log的针对单精度fpu实现的C函数。

achild

搞个小算法，用的stm32f7，带的有fpu。
里面的pow自己写简单，写出来比math库里面的快了好多倍。
exp和log速度那是相当慢。有没有指点一下这个要怎么写。查了一圈，没有一点思路。

takashiki

exp：麦克劳林公式、查表，速度应该能上去吧，精度不保证，因为每一步计算都有误差而且会累积
exp查表如下：e的1次方、2次方、4次方、...、64次方，也没几个，小数部分就是麦克劳林级数了，收敛很快，然后把指数拆了相乘就是了

ln暂时只会麦克劳林公式，分区间算

achild

takashiki 发表于 2020-12-15 17:58
exp：麦克劳林公式、查表，速度应该能上去吧，精度不保证，因为每一步计算都有误差而且会累积
exp查表如下 ...

谢了，有思路了

wye11083

takashiki 发表于 2020-12-15 17:58
exp：麦克劳林公式、查表，速度应该能上去吧，精度不保证，因为每一步计算都有误差而且会累积
exp查表如下 ...

你这算法弱爆了。如果以2为基数，那么可以化简为查表法，精度取决于表容量。。
log2=整数部分bitscan+除去最高位的查表的和，exp2=整数部分移位*小数部分查表的积。
这样可以用十几条指令快速计算出来。

对了，fpga可以2个周期出快速结果。64字节rom表可以做到5%精度。

takashiki

wye11083 发表于 2020-12-15 20:54
你这算法弱爆了。如果以2为基数，那么可以化简为查表法，精度取决于表容量。。
log2=整数部分bitscan+除 ...

ln(x)我是真不知道，也许牛顿法都比麦克劳林公式快
exp(x)可以做成两部分，指数是负数的先取绝对值最后取下倒数，整数部分可以做个大表，exp(0) ... exp(88)都放进去，再大就溢出了，小数部分麦克劳林公式收敛很快，迭代4次(1 + x + x^2/2 + x^3/6 + x^4/24)误差就只有千分之五了。至于你说的以2为基数我是真看不懂，比如要求exp(69.0003)应该怎么求，比我弱爆了的算法强多少？
库函数：       exp(69.0003f) = 925653556458082370000000000000.000000
我的算法：myexp(69.0003f) = 925655832077537080000000000000.000000
看起来差距还是很明显的，其实误差只有0.00024583921693310643843984592%

wye11083

takashiki 发表于 2020-12-16 06:32
ln(x)我是真不知道，也许牛顿法都比麦克劳林公式快
exp(x)可以做成两部分，指数是负数的先取绝对值最后取 ...

我指的是速度不是精度。

takashiki

wye11083 发表于 2020-12-16 07:08
我指的是速度不是精度。

那请教下暂不考虑精度的情况下，2为基数的速度模式是如何计算的？就简单点的吧，比如我说的exp(69.0003)

wye11083

takashiki 发表于 2020-12-16 07:15
那请教下暂不考虑精度的情况下，2为基数的速度模式是如何计算的？就简单点的吧，比如我说的exp(69.0003) ...

溢出了。exp4.35的话，先1左移4位，再乘以2^0.35，就是近似解。定点小数的话可以拿尾数移位然后查表。表在pc上做好。

takashiki

wye11083 发表于 2020-12-16 07:38
溢出了。exp4.35的话，先1左移4位，再乘以2^0.35，就是近似解。定点小数的话可以拿尾数移位然后查表。表 ...

exp(4.35) = 77.47846292526083
(1 << 4) * (2 ^ 0.35) = 20.392970037108196

你把这叫近似解？？？？指数函数、对数函数默认都是自然对数的底e，你把它换成2和楼主的原意根本不符好吧。对数可以用换底公式，指数怎么搞？

No.5

就爱看神仙打架的贴子，留个记号

wye11083

takashiki 发表于 2020-12-16 07:45
exp(4.35) = 77.47846292526083
(1

你算错了。我说了2为底。2^4.35不到32大于16。

wye11083

takashiki 发表于 2020-12-16 07:45
exp(4.35) = 77.47846292526083
(1

还有，这是近似计算，本身就是近似的，主要目的用于比较而不是出结果。2为底的定点指数和对数都可以硬件2个周期搞定。

takashiki

wye11083 发表于 2020-12-16 08:45
你算错了。我说了2为底。2^4.35不到32大于16。

到底是我算错了还是你语文表达和理解能力不行？难道我们处于不同的次元？
咱针对的都是LZ位的问题，fpu实现的C函数：pow、exp和log。fpu是浮点处理单元，而且LZ也说了是单精度浮点，你跟我扯定点？？？？
另外，exp就是以e为底的指数，你变换成2为底没问题，e^x = 2 ^ (log2e * x) = 2 ^ (1.4426950408889634 * x) 这都没问题，但根本不是你所说的exp(4.35) 就等价于 2^4.35的，指数函数的底不是你定义的，是有标准规范的。
总体来说，你这个定点算法对于楼主的浮点运算来说，一点价值都没有。

另外，Cortex M从M3开始就内置CLZ指令，他只有一个指令周期，计算定点的以2为底的对数也很简单，不会比你的FPGA慢多少，而且不用查表

achild

看了一圈，两个大佬真猛。taka的解释更贴近我的应用。

t3486784401

takashiki 发表于 2020-12-16 07:45
exp(4.35) = 77.47846292526083
(1

指数一样换底啊，本来和对数就是同一类东西：

exp(N) = 2^(1/ln2)^N = 2^(N/ln2)

然后 exp(4.35) = 2^(4.35/ln2) = 2^6.276

XA144F

把《中学数学用表》敲进去如何？

hexenzhou

exp我用坛子里介绍的Quake III游戏里用的，速度飞快。

achild

hexenzhou 发表于 2020-12-16 16:37
exp我用坛子里介绍的Quake III游戏里用的，速度飞快。

有帖子链接吗

achild

用了麦克劳林公式+查表法，级数到了6级，试了0.999999的小数，float结果只在最后两个有效位差了一点点。速度快了将近一百倍

achild

takashiki 发表于 2020-12-15 17:58
exp：麦克劳林公式、查表，速度应该能上去吧，精度不保证，因为每一步计算都有误差而且会累积
exp查表如下 ...

这个ln太难了啊，极限情况下，麦克劳林级数已经到了x^32次项了，刚刚收敛到有效数据第4位

takashiki

achild 发表于 2020-12-16 22:43
这个ln太难了啊，极限情况下，麦克劳林级数已经到了x^32次项了，刚刚收敛到有效数据第4位 ...

我知道ln的麦克劳林级数收敛很慢，但也不至于那么离谱吧。exp除的是阶乘，这个就是自然数顺序，收敛确实慢。因为尾数段在[1, 2)范围内，所以以1.5分段计算，因此猜测1.5附近时收敛最慢，此时迭代9次到x^10时有效位数能达到4位，其他的除了非常接近1的以外我测试迭代8次基本能到6位有效数字。

1. float myln(float f) {
   
2.     if(f < 0)                          // 没有意义的东西
   
3.         return -0.0 / 0.0;
   
4.     else if(f == 0)                    // 负无穷大
   
5.         return -1.0 / 0.0;
   
6.    
   
7.     // 分段计算：指数段(E)、尾数段(M)。符号位(S)已经排除掉了
   
8.     // float的内存排布：SEEEEEEE EMMMMMMM MMMMMMMM MMMMMMMM
   
9.     int pow = (*(int*)&f >> 23) - 127;                // 指数段
   
10.     *(int*)&f = *(int*)&f & 0x7FFFFF | 0x3F800000;    // 尾数段，就在半开半闭区间[1, 2)范围内，以1.5为分界分段
    
11.    
    
12.     int region = f >= 1.5 ? 1 : 0;
    
13.     if(region)
    
14.         f /= 2;
    
15. 
    
16.     // 麦克劳林级数，计算的是ln(1+x) = x - x^2/2 + x^3/3 - x^4/4 ... ，所以f先减去1
    
17.     f -= 1;
    
18.     int i;                            // 级数计算次数
    
19.     float ff = f, y = 0;
    
20.     for(i = 1; i < 8; i++, ff *= -f)
    
21.         y += ff / i;
    
22.         
    
23.     if(region)
    
24.         y += 0.6931471805599453;
    
25.    
    
26.     return pow * 0.6931471805599453 + y;
    
27. }

复制代码

如果你的编译器不支持除以浮点0请将最开始几行去掉，但在IEEE754中是有定义的。

在网上查到有用牛顿法的，迭代次数也很多，感觉也快不到哪儿去。泰勒级数ln(x + a)类型的还有很多展开式，其中有只有奇次项的，收敛似乎很快，但是我写了代码发现误差很感人，根本没法用

takashiki

achild 发表于 2020-12-16 22:43
这个ln太难了啊，极限情况下，麦克劳林级数已经到了x^32次项了，刚刚收敛到有效数据第4位 ...

再改一下，这个可以控制误差，有效位数轻松上到6~7位。每个数字的迭代次数不同，简单的也许一次就够了，复杂的二三次也很可能。

1. float myln(float f) {
   
2.     if(f < 0)                          // 没有意义的东西
   
3.         return -0.0 / 0.0;
   
4.     else if(f == 0)                    // 负无穷大
   
5.         return -1.0 / 0.0;
   
6.    
   
7.     // 分段计算：指数段(E)、尾数段(M)。符号位(S)已经排除掉了
   
8.     // float的内存排布：SEEEEEEE EMMMMMMM MMMMMMMM MMMMMMMM
   
9.     int pow = (*(int*)&f >> 23) - 127;                // 指数段
   
10.     *(int*)&f = *(int*)&f & 0x7FFFFF | 0x3F800000;    // 尾数段，就在半开半闭区间[1, 2)范围内
    
11.    
    
12.     // 麦克劳林级数，计算的是ln(1+x) = x - x^2/2 + x^3/3 - x^4/4 ... ，所以f先减去1
    
13.     int region = f >= 1.5 ? 1 : 0;
    
14.     if(region)
    
15.         f /= 2;
    
16. 
    
17.     f -= 1;
    
18.     int i = 1;                       // 级数计算次数
    
19.     float ff = f, y = 0;
    
20.     do {
    
21.         y += ff / i;
    
22.         ff *= -f;
    
23.         i ++;
    
24.         printf("  ... i=%d, ff=%f\n", i, ff);            // PC机上调试看迭代的次数
    
25.     } while((ff > 1e-5) || (ff < -1e-5));               // 误差控制
    
26.         
    
27.     if(region)
    
28.         y += 0.6931471805599453;
    
29.    
    
30.     return pow * 0.6931471805599453 + y;
    
31. }

复制代码

achild

takashiki 发表于 2020-12-17 07:57
再改一下，这个可以控制误差，有效位数轻松上到6~7位。每个数字的迭代次数不同，简单的也许一次就够了， ...

看明白了，我是按成十进制算的，好麻烦，效率还低。0.9的收敛速度比0.5的收敛速度确实慢太多了。！！