scratch lenet(10): C语言计算log

1. 目的

用于复现 LeNet-5 时，损失函数的计算。损失函数中的惩罚项（正则项）出现了 $\log ()$ 函数， 而我的设定是复现过程不能用 C 标准库的数学库。

具体的实现依赖于选择的数学公式，了解到的有两种：

1. 泰勒级数展开。缺点是比较慢。
2. 结合 IEEE-754 二进制表示法 和 Remez 算法。速度快。

本文使用第二种方法，即 IEEE-754 二进制 + Remez 算法, 考察的数据类型是 float。

2. 原理: $x=m\ast 2^p$

对于 float 类型变量 x, 它可以表示为

$$x=m\ast 2^p$$

为啥呢？把 y 的二进制展开一下就可以理解了¹：

   sign  exp        frac
0b 0     [00000000] 00000000000000000000000
value = (-1)^s * m * 2 ^ (exp - 127)

其中:

• $m=1.0+0.frac\in [1,2]$
• exp 减去 127 是规定，适用于 $x>2^{-126}$
• $p=exp-127$

3. 公式展开

3.1 公式分解

也就是对于 2.1 小节公式，由于等号左右都大于0，可以分别计算对数，得到的公式

$$\text{ln}(x)=\text{ln}(m\ast 2^p)=\text{ln}(m)+\text{ln}(2^p)=\text{ln}(m)+\text{ln}(2)p$$

其中：

• $m\in [1,2]$, 可以从 $x$ 的二进制表示的 “frac” 部分获取到
• $\text{ln}(2)$ 是常数 0.69314718
• $p$ 可以从 $x$ 的二进制表示的 “exp” 部分获取到

对于 $\text{ln}(m)$, 由于 $m\in [1,2]$, 可用 Remez 算法近似算出。因此剩下的问题为：

• 获取 $m$
• 获取 $p$
• 应用 Remez 算法

3.2 获取 $m$

把 $x$ 的二进制表示中，exp 部分修改为 127， 使得 value 表达式等于：
$$value=(-1{)}^s\ast M\ast 2^{127-127}=(-1{)}^s\ast m$$
并且 x 的二进制的其他部分是不变的。这就得到了 m 的二进制表示。对应代码如下：

float m_ln(float x)
{unsigned int bx = *(unsigned int *) (&x);// exp:  00000000// frac: 0b11111111111111111111111 = 838607unsigned int bm = (127 << 23) | (bx & 8388607);float m = *(float *) (&bm);printf("m = %f\n", m); ...
}

3.3 获取 $p$

根据 $p=\text{exp}-127$, 要获取 $p$ 就需要先获取 exp. 在 C 语言中标准库已经用了 exp 这个名字，因此我们用 ex 来表示它：

float m_ln(float x)
{unsigned int bx = *(unsigned int *) (&x);unsigned int ex = bx >> 23;signed int p = (signed int)ex - (signed int)127;printf("p = %d\n", p);...
}

3.4 Remez 算法

Remez 算法被 Maple 等数值计算软件采用，精度基本够用，速度够快，解决了给定区间 $x\in [a,b]$ 上用多项式逼近给定函数 $f(x)$ 的 minimax() 问题的解。理论发展过程大概是： Chebyshev 证明了迭代求解过程的收敛性， Remez 则给出了具体的迭代求解过程。

网友 Crouching Kitten 提供了 Remez 算法的具体使用 ^{2 3}.

Remez 算法用于 sin(x) 的快速计算

B站视频 【硬核】从零开始自己编写FOC 算法篇3：快速正弦算法⁴ 则给出了一些理论公式，以及基于 Remez C++ 库 LolRemez⁵ 算出了一些数据，从而在 DSP 上得到了比 DSP 数学库还快还准的 sin(x) 的实现。公式：

Remez 算法用于 exp 的近似

博客⁶ 给出了 exp 的 Remez 的 Octave/MatLab 实现。

% https://loki-pup.github.io/posts/2019-9-30-Remez%E7%AE%97%E6%B3%95%E7%9A%84MATLAB%E5%AE%9E%E7%8E%B0-2019
function [] = remez(y)w = 10^(-6);r4 = 1;o = 0;while r4 >= wo=o+1;z = inp(y);[s,r4] = comp(z);if r4 < wbreak;endy = subs(y,s,z);endofprintf('%f*x^2+%f*x+%f\n',z(4),z(3),z(2))r4
endfunction [ x ] = inp( y )%UNTITLED 此处显示有关此函数的摘要%   此处显示详细说明a = zeros(4,4);b = zeros(4,1);for i =1:4a(i,1) = (-1)^(i-1);a(i,2) = 1;a(i,3) = y(i);a(i,4) = y(i)^(2);b(i,1) = exp(y(i));endx = pinv(a)*b;
endfunction [s,r4] = comp(q)f1 = @(x)exp(x)-q(2)-q(3)*x-q(4)*x^2;f2 = @(x)q(4)*x^2+q(2)+q(3)*x-exp(x);[s1,r1] = fminbnd(f1,-1,1);[s2,r2] = fminbnd(f2,-1,1);if abs(r1)<=abs(r2)r3 = abs(r2);elser3 = abs(r1);ends = [s1,s2];r4 = r3 - abs(q(1));
endfunction [y] = subs(y,s,q)     %替代，选取新点组，单一交换法p = length(s);for i = 1:pif s(i)>y(1)||s(i)<y(4)for j = 2:4if s(i)<y(j)c1 = sign(exp(y(j-1))-q(2)-q(3)*y(j-1)-q(4)*(y(j-1))^2);c2 = sign(exp(s(i))-q(2)-q(3)*s(i)-q(4)*(s(i))^2);if c1 == c2y(j-1) = s(i);elsey(j) = s(i);endendendelse if s(i)> -1 && s(i)<y(1)c1 = sign(exp(y(1))-q(2)-q(3)*y(1)-q(4)*y(1)^2);c2 = sign(exp(s(i))-q(2)-q(3)*s(i)-q(4)*(s(i))^2);if c1 == c2y(1) = s(i);elsey(4) = s(i);endelse if s(i)> y(4) && s(i)<1c1 = sign(exp(y(4))-q(2)-q(3)*y(4)-q(4)*y(4)^2);c2 = sign(exp(s(i))-q(2)-q(3)*s(i)-q(4)*(s(i))^2);if c1 == c2y(4) = s(i);elsey(1) = s(i);endendendendend
end

Remez 用于自然对数 ln(x) 的计算

对于 $ln(m),m\in [1,2]$ 而言， Remez 的4阶展开为：(来自²)
$$-1.7417939+(2.8212026+(-1.4699568+(0.44717955-0.056570851\ast x)\ast x)\ast x)\ast x;$$
也可以用3阶的展开，精度差一点，计算量少一点：
$$-1.49278+(2.11263+(-0.729104+0.10969\ast x)\ast x)\ast x;$$

尝试使用 Remez C++ 库 LolRemez⁵ 获取同样的结果， 不过编译失败了， 而 automake / P4 那一套我完全不熟悉，暂时放弃这个库。

也可以使用 boost 中提供的 remez 函数 ^{1 7}, 不过新版 boost 已经找不到 remez.hpp 头文件， 需要使用 $boost<=1.56$ 的版本，暂时没有尝试

#include <boost/math/tools/remez.hpp>boost::math::tools::remez_minimax<double> approx([](const double& x) { return log(x); },4, 0, 1, 2, false, 0, 0, 64
);

4. 结果

4.1 代码

// Author: Zhuo Zhang
// Homepage: https://github.com/zchrissirhcz
// Create Date: 2023.06.24// logarithm for natural number `e`
float m_log(float x)
{// x = m * 2^p, m \in [1, 2]// ln(x) = ln(m * 2^p) = ln(m) + ln(2) * p// determine punsigned int bx = *(unsigned int *) (&x);unsigned int ex = bx >> 23;signed int p = (signed int)ex - (signed int)127;// determine m// exp:  00000000// frac: 0b11111111111111111111111 = 838607unsigned int bm = (127 << 23) | (bx & 8388607);float m = *(float *) (&bm);// printf("m = %f\n", m); // determine ln(m) by Remez algorithm for m in [1, 2]float ln_m_approx_4th_order = -1.7417939 + (2.8212026 + (-1.4699568 + (0.44717955 - 0.056570851 * m) * m) * m) * m;float ln_m_approx_3rd_order = -1.49278 + (2.11263 +(-0.729104 + 0.10969 * m) * m) * m;float ln_m_approx = ln_m_approx_4th_order;// combine the resultconst float ln2 = 0.6931471806;float res = ln_m_approx + ln2 * p;return res;
}#include <stdio.h>
#include <stdbool.h>
#include <math.h>
int main()
{float x;while (true){scanf("%f", &x);float y1 = m_log(x);float y2 = log(x);printf("m_log(%f) = %f\n", x, y1);printf("log(%f) = %f\n", x, y2);}return 0;
}

4.2 运行结果

zz@Legion-R7000P% gcc log.c -lm
zz@Legion-R7000P% ./a.out 
0.0123
m_log(0.012300) = -4.398208
log(0.012300) = -4.398156
1.234
m_log(1.234000) = 0.210295
log(1.234000) = 0.210261
2.345
m_log(2.345000) = 0.852274
log(2.345000) = 0.852285
3.456
m_log(3.456000) = 1.240081
log(3.456000) = 1.240112

5. References

---

1. https://gist.github.com/LingDong-/7e4c4cae5cbbc44400a05fba65f06f23 ↩︎ ↩︎

2. Efficient implementation of natural logarithm (ln) and exponentiation - Crouching Kitten’s Answer ↩︎ ↩︎

3. https://en.wikipedia.org/wiki/Remez_algorithm ↩︎

4. 【硬核】从零开始自己编写FOC 算法篇3：快速正弦算法 ↩︎

5. https://github.com/samhocevar/lolremez ↩︎ ↩︎

6. Remez算法的MATLAB实现 ↩︎

7. How to Find a Fast Floating-Point atan2 Approximation ↩︎