资讯

历史

科技

环境与自然

成长

游戏

财经

文学与艺术

美食

健康

家居

文化

情感

汽车

三农

军事

旅行

运动

教育

生活

星座命理

一种快速的幂运算方法（底数是自然数e，指数是浮点数）

创作时间:

作者:

@小白创作中心

一种快速的幂运算方法（底数是自然数e，指数是浮点数）

引用

CSDN

https://blog.csdn.net/Yemiekai/article/details/108431606

本文介绍了一种快速计算自然数e的浮点数次幂的方法，这种方法在深度学习等需要大量幂运算的场景中具有重要的应用价值。文章详细介绍了方法的原理、实现代码以及应用场景，对于从事相关领域工作的技术人员具有较高的参考价值。

问题

【给一个浮点数y，现在需要求出e的y次幂的值】。

对于这个问题，最直接的方法是使用库函数，例如在C++中<math.h>头文件提供了exp()函数，Python里通过import math使用math.exp()。这些方法精度较高，但是速度相当慢。

在深度学习（DeepLearning）中经常需要花费大量时间进行幂运算，典型场景是使用激活函数和计算概率分布的时候。例如在SoftMax层通常需要进行底数是e，指数是浮点数的幂运算。

提高幂运算的速度，能有效提高实际应用的速度。据说一些消费级的NVIDIA显卡都是把双精度给砍了，有时候你甚至是用半精度在训练。对于大多数神经网络计算而言，近似精度是完全足够的，并且可以节省很多时间。

方法

有一些其它的快速幂运算方法，如查找表，使用线性插值等。

这里参考文章《A Fast, Compact Approximation of the Exponential Function》的方法，能够以较少的精度损失换取明显的速度提升。

经测试，这种方法比库函数快几倍到几十倍，具体看指数有多复杂。据说在某些特定的值上误差范围有±10%，这个需要根据实际需求权衡精度与速度。

假设目标机器是大端字节序，double类型为64位，float类型为32位，int类型为32位，short类型为32位。

double版本

(version1)

inline double fast_exp(double y){
    union{
        double d;
        int x[2];
    }data = {y};
    
    data.x[0] = 0;
    data.x[1] = (int)(y * 1512775 + 1072632447);
    
    return data.d;
}

(version2)

inline double fast_exp(double y){
    double d;
    *(reinterpret_cast<int*>(&d) + 0) = 0;
    *(reinterpret_cast<int*>(&d) + 1) = static_cast<int>(y * 1512775  + 1072632447);
    return d;
}

以上两段代码意思是一样的，只是实现方式不一样。version1用联合体是为了能分别拿到一块64位数据的高32位和低32位，version2通过修改指针的类型，也是为了拿到高位和低位。

float版本

(version1)

inline float fast_exp(float y) {
    float d;
    *(reinterpret_cast<short*>(&d) + 0) = 0;
    *(reinterpret_cast<short*>(&d) + 1) = static_cast<short>(184 * y + (16256-7));
    return d;
}

因为在Cortex-A7的Neon Intrinsics中没有双精度浮点数的类型，只能用到float，所以参考别人的文章写了一个float版本的实现，以便使用Neon加速计算。

(version2)参考自https://www.itread01.com/content/1550634858.html，尚未验证

union
{
    uint32_t i;
    float f;
}v;
v.i=(1<<23)*(1.4426950409*val+126.94201519f);
return v.f;

原理

根据IEEE754-1985标准（IEEE Standard for Binary Floating-Point Arithmetic），一个浮点数可以通常用以下形式表示：
[
(-1)^s \cdot(1+m)\cdot2^{x-x_0}
]
公式1

其中s是符号位，m是尾数（一串内存里的二进制的数字），x是指数项，x0是偏置（bias）。

对于一个64位的浮点数，尾数m占52位，指数项x占11位，偏置x0=1023，在内存空间占8个字节：

(图1) 双精度浮点数的内存排列

以上是浮点数的表示法及其数据存放特点。先记住。

然后… 先来看看2的y次幂怎么求：

现在输入一个浮点数y，需要计算2的y次幂。

y是一个浮点数，它的表达式为y=(-1)^s \cdot(1+m)\cdot2^{x-x_0}。观察发现，y的表达式里面就含有2次幂 “2^{x-x_0}”，我们正好需要计算2次幂，把x-x0换成y不就行了？妙啊。

上面讲了思路，具体怎么操作呢？看回图1，指数项x在内存的[53~63]位，把y放到对应的位上面，就完成了替换。

(图2) 和图1一样的意思

首先把y当成int数，然后加上x0，也就是y+1023（根据规范，双精度浮点数的偏置项bias是x0=1023）（可能加上x0是为了消掉x-x0中的x0，把2^{x-x_0}变成2^y）。

然后把结果左移20位（乘以2^{20}），就对应到指数项所在的坑（图2中绿色格子），由此把指数项换成了y。

结合上述步骤，求2的y次幂的方法就是：取出y的高32位（图1中的i），让它等于2^{20}\cdot(y+1023)即可。得到的结果就是(-1)^s \cdot(1+m)\cdot2^y，这里还有m，后面再讲怎么处理。

通用表达式是：令y的高32位i = ay + (b-c)

求e的y次幂的时候，式中a = 2^{20} / ln(2)，b = 1023\cdot2^{20}，c的经验值是60801，c是用于减少误差的。

为什么a是这个值，不是2^{20}吗。因为这是在求e的y次幂。前面原理讲是针对2的y次幂讲的，求e的y次幂的时候需要变一下，看下面的推导：

2^a = e^{ln2^a}=e^{a{\cdot}ln2}

令y = a{\cdot}ln2，则a = y {\cdot} {\frac{1}{ln2}}，上面的式子变成：

2^{y {\cdot} {\frac{1}{ln2}} }=e^y

1/ln2是一个常数，约为1.442695....，因此e的y次幂可以通过求2的y次幂得到，过程是一样的，变换一下y，把输入的y乘上1/ln2即可。

c为什么是68243，这有点复杂，请看原文作者的推导。

所以：
a = 2^{20} / ln(2) = 1512775，
(b-c)=1023\cdot2^{20}-60801=1072632447

就和代码里的数值对应上了(见double版本的version1)。

单精度的计算方法类似，根据单精度浮点数的存储方式改一下参数就可以了。

注意：
这种快速幂运算的方法对输入数据y是有要求的，对于double版本而言，输入y大概要在[-700,700]的区间，超出范围算法失效。对于float版本而言，在[-10,10]之间是没问题的。

关于m：为什么代码里把低32位的数据置零，因为这一步是为了把公式1中的m置零，保证只有指数项。再看图2，只是把低32位的m置零了，高32位还有20个m不用管？确实没有管，原文作者说保留这部分的m不仅没什么影响，反而有助于提高精度。

总结

上面的原理只是大概近似的理解，并不是很深刻。原文只讲述了做法和过程，给了一条公式，没有详细解释原因，可能需要更深入的研究。根据这种方法修改到float类型上也能work，看来原理是没问题的。有兴趣的可以再看看原文《A fast, compact approximation of the exponential function》。另外需要结合浮点数的原理，参考IEEE754规范《754-1985 - IEEE Standard for Binary Floating-Point Arithmetic》。