C语言 自然对数(ln)和幂运算的高效实现

ktca8awb  于 2022-12-17  发布在  其他
关注(0)|答案(8)|浏览(320)

我正在寻找C库<math.h>中提供的log()exp()函数的实现。我正在使用8位微控制器(OKI 411和431)。我需要计算Mean Kinetic Temperature。要求是我们应该能够尽可能快地计算MKT,并且使用尽可能少的代码内存。编译器在<math.h>中提供了log()exp()函数,但是调用其中一个函数并与库链接会导致代码大小增加5 KB,这将不适合我们使用的微处理器(OKI 411),因为我们的代码已经消耗了约12K的可用代码内存(约15K)。
我正在寻找的实现不应该使用任何其他的C库函数(如pow(),sqrt()等)。这是因为所有的库函数都打包在一个库中,即使调用一个函数,链接器也会将整个5K库带到代码内存中。
编辑
算法应正确,小数点后3位。

wqsoz72f

wqsoz72f1#

使用泰勒级数不是最简单也不是最快的方法。大多数专业的实现都使用近似多项式。我将展示如何在Maple(这是一个计算机代数程序)中使用Remez algorithm生成一个近似多项式。
对于3位精度,在Maple中执行以下命令:

with(numapprox):
Digits := 8
minimax(ln(x), x = 1 .. 2, 4, 1, 'maxerror')
maxerror

其响应为以下多项式:

-1.7417939 + (2.8212026 + (-1.4699568 + (0.44717955 - 0.056570851 * x) * x) * x) * x

最大误差为:0.000061011436

我们生成了一个近似ln(x)的多项式,但只在[1..2]区间内。增加区间是不明智的,因为这将增加最大误差甚至更多。

因此,首先求出2的最高次幂,它仍然小于数字(参见:What is the fastest/most efficient way to find the highest set bit (msb) in an integer in C?)。这个数字实际上是以2为底的对数。除以这个值,结果就进入了1..2区间。最后我们必须加上n*ln(2)才能得到最终结果。
数字〉= 1的示例实现:

float ln(float y) {
    int log2;
    float divisor, x, result;

    log2 = msb((int)y); // See: https://stackoverflow.com/a/4970859/6630230
    divisor = (float)(1 << log2);
    x = y / divisor;    // normalized value between [1.0, 2.0]

    result = -1.7417939 + (2.8212026 + (-1.4699568 + (0.44717955 - 0.056570851 * x) * x) * x) * x;
    result += ((float)log2) * 0.69314718; // ln(2) = 0.69314718

    return result;
}

不过如果你只打算在[1.0,2.0]区间内使用它,那么函数就像这样:

float ln(float x) {
    return -1.7417939 + (2.8212026 + (-1.4699568 + (0.44717955 - 0.056570851 * x) * x) * x) * x;
}
ss2ws0br

ss2ws0br2#

e^x的泰勒级数收敛非常快,你可以调整你的实现以达到你需要的精度。
对数的泰勒级数不太好...

kuuvgm7e

kuuvgm7e3#

如果你不需要浮点运算,你可以很容易地计算出一个近似的以2为底的分数对数,首先把你的值左移到32768或更高,然后把你这样做的次数存储在count中,然后重复一定的次数(取决于你想要的比例因子):

n = (mult(n,n) + 32768u) >> 16; // If a function is available for 16x16->32 multiply
count<<=1;
if (n < 32768) n*=2; else count+=1;

如果上述循环重复8次,则该数字以2为底的对数为count/256。如果重复10次,则为count/1024。如果重复11次,则为count/2048。实际上,此函数通过计算n**的整数2次幂对数(2^reps)来工作,但中间值经过缩放以避免溢出。

mbzjlibv

mbzjlibv4#

基本表与插值之间的值接近工作?如果值的范围是有限的(这可能是你的情况-我怀疑温度读数有很大的范围)和高精度是不需要的,它可能会工作。应该很容易在正常的机器上测试。
下面是关于函数的表表示的众多主题之一:Calculating vs. lookup tables for sine value performance?

xxe27gdn

xxe27gdn5#

招魂术。
我必须在有理数上实现对数。
我是这样做的:
根据维基百科,有哈雷-牛顿近似法

其可以用于非常高的精度。
使用牛顿法,迭代简化为(implementation),它具有对ln(x)的三次收敛,这比泰勒级数提供的收敛要好得多。

// Using Newton's method, the iteration simplifies to (implementation) 
// which has cubic convergence to ln(x).
public static double ln(double x, double epsilon)
{
    double yn = x - 1.0d; // using the first term of the taylor series as initial-value
    double yn1 = yn;

    do
    {
        yn = yn1;
        yn1 = yn + 2 * (x - System.Math.Exp(yn)) / (x + System.Math.Exp(yn));
    } while (System.Math.Abs(yn - yn1) > epsilon);

    return yn1;
}

这不是C,而是C#,但我相信任何能用C编程的人都能从中推导出C代码。
此外,由于
对数n(x)= ln(x)/ln(n)。
因此,您也实现了logN。

public static double log(double x, double n, double epsilon)
{
    return ln(x, epsilon) / ln(n, epsilon);
}

其中epsilon(误差)是最小精度。
至于速度,你可能更好地使用内嵌硬件,但正如我所说的,我用它作为一个基础,以实现对数的有理数类工作与任意精度。
在某些情况下,任意的精确度可能比速度更重要。
然后,对有理数使用对数恒等式:
对数B(x/y)=对数B(x)-对数B(y)

nvbavucw

nvbavucw6#

除了Crouching Kitten的答案给了我灵感之外,您还可以构建一个伪递归(最多1次自调用)对数,以避免使用多项式。

ln(x) :=
If (x <= 0)
    return NaN
Else if (!(1 <= x < 2))
    return LN2 * b + ln(a)
Else
    return taylor_expansion(x - 1)

这是相当有效和精确的,因为在[1; 2)泰勒级数收敛得快得多,如果我们的输入是正的但不在这个范围内,我们在第一次调用ln时得到这样一个数1〈= a〈2。
你可以从浮点数x中的数据找到'B'作为你的无偏指数,从浮点数x的尾数找到'a'(a和x是完全相同的浮点数,但是现在是exponent biased_0而不是exponent biased_b)。LN 2应该以十六进制浮点表示法IMO作为一个宏保存。你也可以使用http://man7.org/linux/man-pages/man3/frexp.3.html
还有诀窍

unsigned long tmp = *(ulong*)(&d);

对于“内存强制转换”,在处理浮点内存方式时,了解double到unsignedlong而不是“值强制转换”是非常有用的,因为位操作符将根据编译器而导致警告或错误。

ekqde3dh

ekqde3dh7#

在没有<math.h>的情况下,在C中可能计算ln(x)和expo(x):

static double expo(double n) {
    int a = 0, b = n > 0;
    double c = 1, d = 1, e = 1;
    for (b || (n = -n); e + .00001 < (e += (d *= n) / (c *= ++a)););
    // approximately 15 iterations
    return b ? e : 1 / e;
}

static double native_log_computation(const double n) {
    // Basic logarithm computation.
    static const double euler = 2.7182818284590452354 ;
    unsigned a = 0, d;
    double b, c, e, f;
    if (n > 0) {
        for (c = n < 1 ? 1 / n : n; (c /= euler) > 1; ++a);
        c = 1 / (c * euler - 1), c = c + c + 1, f = c * c, b = 0;
        for (d = 1, c /= 2; e = b, b += 1 / (d * c), b - e/* > 0.0000001 */;)
            d += 2, c *= f;
    } else b = (n == 0) / 0.;
    return n < 1 ? -(a + b) : a + b;
}

static inline double native_ln(const double n) {
    //  Returns the natural logarithm (base e) of N.
    return native_log_computation(n) ;
}

static inline double native_log_base(const double n, const double base) {
    //  Returns the logarithm (base b) of N.
    return native_log_computation(n) / native_log_computation(base) ;
}

在线试用

u5i3ibmn

u5i3ibmn8#

根据@Crouching Kitten上面的自然对数回答,如果你需要它在输入小于1时保持精确,你可以添加一个简单的比例因子。下面是一个我在微控制器中使用的C++示例。它的比例因子为256,输入精确到1/256 = ~0.04。最大为2^32/256 = 16777215(由于uint 32变量溢出)。
有趣的是,即使在没有FPU的STMF 103 Arm M3上,下面的浮点实现也比libfixmath中的16位定点实现快得多(例如3倍或更快)(也就是说,此浮点实现仍需要数千个周期,因此仍不快)

#include <float.h>

float TempSensor::Ln(float y) 
{
    // Algo from: https://stackoverflow.com/a/18454010
    // Accurate between (1 / scaling factor) < y < (2^32  / scaling factor). Read comments below for more info on how to extend this range

    float divisor, x, result;
    const float LN_2 = 0.69314718; //pre calculated constant used in calculations
    uint32_t log2 = 0;

    //handle if input is less than zero
    if (y <= 0)
    {
        return -FLT_MAX;
    }

    //scaling factor. The polynomial below is accurate when the input y>1, therefore using a scaling factor of 256 (aka 2^8) extends this to 1/256 or ~0.04. Given use of uint32_t, the input y must stay below 2^24 or 16777216 (aka 2^(32-8)), otherwise uint_y used below will overflow. Increasing the scaing factor will reduce the lower accuracy bound and also reduce the upper overflow bound. If you need the range to be wider, consider changing uint_y to a uint64_t
    const uint32_t SCALING_FACTOR = 256;
    const float LN_SCALING_FACTOR = 5.545177444; //this is the natural log of the scaling factor and needs to be precalculated

    y = y * SCALING_FACTOR; 
    
    uint32_t uint_y = (uint32_t)y;
    while (uint_y >>= 1) // Convert the number to an integer and then find the location of the MSB. This is the integer portion of Log2(y). See: https://stackoverflow.com/a/4970859/6630230
    {
        log2++;
    }

    divisor = (float)(1 << log2);
    x = y / divisor;    // FInd the remainder value between [1.0, 2.0] then calculate the natural log of this remainder using a polynomial approximation
    result = -1.7417939 + (2.8212026 + (-1.4699568 + (0.44717955 - 0.056570851 * x) * x) * x) * x; //This polynomial approximates ln(x) between [1,2]

    result = result + ((float)log2) * LN_2 - LN_SCALING_FACTOR; // Using the log product rule Log(A) + Log(B) = Log(AB) and the log base change rule log_x(A) = log_y(A)/Log_y(x), calculate all the components in base e and then sum them: = Ln(x_remainder) + (log_2(x_integer) * ln(2)) - ln(SCALING_FACTOR)

    return result;
}

相关问题