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

文章目录

scratch lenet(10): C语言计算log
- 1. 目的
- 2. 原理: $x = m * 2^p$
- 3. 公式展开
- - 3.1 公式分解
  - 3.2 获取 $m$
  - 3.3 获取 $p$
  - 3.4 Remez 算法
  - - Remez 算法用于 sin(x) 的快速计算
    - Remez 算法用于 exp 的近似
    - Remez 用于自然对数 ln(x) 的计算
- 4. 结果
- - 4.1 代码
  - 4.2 运行结果
- 5. References

1. 目的

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

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

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

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

2. 原理: $x = m * 2^p$

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

$x = m * 2^p$

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

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

其中:

$\in [1, 2]$
exp 减去 127 是规定，适用于 $x > 2^{-126}$
$p = e x p - 127$

3. 公式展开

3.1 公式分解

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

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

其中：

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

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

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

3.2 获取 $m$

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

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

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

网友 Crouching Kitten 提供了 Remez 算法的具体使用 ² ³.

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 >= w
        o=o+1;
        z = inp(y);
        [s,r4] = comp(z);
        if r4 < w
            break;
        end
        y = subs(y,s,z);
    end
    o
    fprintf('%f*x^2+%f*x+%f\n',z(4),z(3),z(2))
    r4
end

function [ x ] = inp( y )
    %UNTITLED 此处显示有关此函数的摘要
    %   此处显示详细说明
    a = zeros(4,4);
    b = zeros(4,1);
    for i =1:4
        a(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));
    end
    x = pinv(a)*b;
end

function [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);
    else
        r3 = abs(r1);
    end
    s = [s1,s2];
    r4 = r3 - abs(q(1));
end

function [y] = subs(y,s,q)     %替代，选取新点组，单一交换法
    p = length(s);
    for i = 1:p
        if s(i)>y(1)||s(i)<y(4)
        for j = 2:4
            if 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 == c2
                    y(j-1) = s(i);
                else
                    y(j) = s(i);
                end
            end
        end
        else 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 == c2
                    y(1) = s(i);
                else
                    y(4) = s(i);
                end
        else if s(i)> y(4) && s(i)<1
                c1 = 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 == c2
                    y(4) = s(i);
                else
                    y(1) = s(i);
                end
            end
            end
        end
    end
end

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

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

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

也可以使用 boost 中提供的 remez 函数 ¹ ⁷, 不过新版 boost 已经找不到 remez.hpp 头文件，需要使用 $b oos t <= 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 p
    unsigned 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 = 838607
    unsigned 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 result
    const 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

https://gist.github.com/LingDong-/7e4c4cae5cbbc44400a05fba65f06f23 ↩︎ ↩︎
Efficient implementation of natural logarithm (ln) and exponentiation - Crouching Kitten’s Answer ↩︎ ↩︎
https://en.wikipedia.org/wiki/Remez_algorithm ↩︎
【硬核】从零开始自己编写FOC 算法篇3：快速正弦算法 ↩︎
https://github.com/samhocevar/lolremez ↩︎ ↩︎
Remez算法的MATLAB实现 ↩︎
How to Find a Fast Floating-Point atan2 Approximation ↩︎