c - 非线性 lut 对数精度

Question

目前我计算日志如下:

#define MAXLOG 1001
double myloglut[MAXLOG];
void MyLogCreate()
{
    int i;
    double exp, expinc;
    expinc = (2.0 - 0.1) / MAXLOG;
    for (i = 0, exp = 0.1; i <= MAXLOG; ++i, exp += expinc)
            myloglut[i] = log(exp);
    myloglut[478] = 0; // this one need to be precise
}

double MyLog(double v)
{
    int idx = (int)((MAXLOG*(v - 0.1)) / (2.0 - 0.1));
    return myloglut[idx];
}

如您所见，我只对 0.1 - 2.0 范围感兴趣。但是，我需要 0 附近更精确。我怎样才能实现非线性计算？还有什么方法可以在此函数中使用一些插值以获得更好的精度？

Answer 1

最新版本

#include <stdio.h>                  // for input/output.
#include <math.h>                   // for mathmatic functions (log, pow, etc.)

// Values
#define MAXELM 1000                 // Array size
#define MINVAL 0.1                  // Minimum x value
#define MAXVAL 1.9                  // Maximum x value
#define EXPVAR 1.4                  // Exponent which makes the variation non linear. If set to 1, the variation will be linear.
#define ACRTPT (MINVAL + MAXVAL)/2  // Accurate point. This value is used to know where to compute with maximum accuracy. Can be set to a fixed value.
// Behavior
#define STRICT 0                    // if TRUE: Return -1 instead of the floored (or closest if out of range) offset when (x) hasn't been calculated for this value.
#define PNTALL 0                    // if TRUE: Print all the calculated values.
#define ASKFOR 1                    // if TRUE: Ask for a x value then print the calculated ln value for it.

// Global vars
double results[MAXELM];             // Array to store computed values.

// Func: offset to var conversion
double getvar(int offset)
{
    double x = (double)MINVAL + ((double)MAXVAL - (double)MINVAL) * (double)offset / (double)MAXELM;

    if(x >= (double)ACRTPT)
        x = pow(x - (double)ACRTPT, (double)EXPVAR) + (double)ACRTPT;
    else
        x = -pow((double)ACRTPT - x, (double)EXPVAR) + (double)ACRTPT;
    // This ^ is the equation used when NONLIN = 1; to have a non linear repartition. Feel free to change it. The inverse equation is in `int getoffset(double)`.
    return x;
}

// Func: var to offset conversion
int getoffset(double var)
{
    double x = var;

    if(x >= (double)ACRTPT)
        x = pow(x - (double)ACRTPT, 1.0/(double)EXPVAR) + (double)ACRTPT;
    else
        x = -pow((double)ACRTPT - x, 1.0/(double)EXPVAR) + (double)ACRTPT;
    // This ^ is the equation used when NONLIN = 1; to calculate offset with a non linear repartition. Feel free to change it (but it must be the inverse of the one in
    // `double getvar(int)` for this to work.). These equations are tied, so you cannot modify one without modifying the other. They are here because
    // `pow(negative, non-integer)` always returns `-nan` instead of the correct value. This 'trick' uses the fact that (-x)^(1/3) == -(x^(1/3)) to cicumvent the
    // limitation.

    int offset = (x - (double)MINVAL) * (double)MAXELM / ((double)MAXVAL - (double)MINVAL);
#if STRICT
    if(getvar(offset) != var)
        return -1;
    return (offset < 0)?-1:(offset > (MAXELM - 1))?-1:offset;
#else
    return (offset < 0)?0:(offset > (MAXELM - 1))?MAXELM - 1:offset;
#endif
}

// Func: main.
int main(int argc, char* argv[])
{
    int offset;
    for(offset = 0; offset < MAXELM; offset++)
        results[offset] = log(getvar(offset));

#if PNTALL
    for(offset = 0; offset < MAXELM; offset++)
    {
        printf("[log(%lf) = %lf] ", getvar(offset), results[offset]);
        if(!((offset + 1) % 6))
            printf("\n");
    }
    printf("\n");
#endif

#if ASKFOR
    double x;
    printf("log(x) for x = ");
    scanf("%lf", &x);
    if((offset = getoffset(x)) < 0)
        printf("ERROR: Value for x = %lf hasn't been calculated\n", x);
    else
        printf("results[%d]: log(%lf) = %lf\n", offset, getvar(offset), results[offset]);
#endif

    return 0;
}

最新版本的特点:

• 使用固定大小的数组。
• 仅计算存储的值(不会计算一个数组单元格的多个值)。
• 使用函数从值获取偏移量以及从偏移量获取值，因此您不必存储已计算出 log 的值。

与上一个版本相比的优点:

• 不使用cbrt，而是使用pow。
• 允许在编译时指定微积分变量的增长。 (因此这些值或多或少地围绕准确点进行分组 (ACRTPT))

<小时/>

第三版

#include <stdio.h>                  // for input/output.
#include <math.h>                   // for mathmatic functions (log, pow, etc.)

// Values
#define MAXELM 1000                 // Array size
#define MINVAL 0.1                  // Minimum x value
#define MAXVAL 1.9                  // Maximum x value
#define ACRTPT (MINVAL + MAXVAL)/2  // Accurate point. This value is used to know where to compute with maximum accuracy. Can be set to a fixed value.
// Behavior
#define NONLIN 1                    // if TRUE: Calculate log values with a quadratic distribution instead of linear distribution.
#define STRICT 1                    // if TRUE: Return -1 instead of the floored (or closest if out of range) offset when (x) hasn't been calculated for this value.
#define PNTALL 0                    // if TRUE: Print all the calculated values.
#define ASKFOR 1                    // if TRUE: Ask for a x value then print the calculated ln value for it.

// Global vars
double results[MAXELM];             // Array to store computed values.

// Func: offset to var conversion
double getvar(int offset)
{
    double x = (double)MINVAL + ((double)MAXVAL - (double)MINVAL) * (double)offset / (double)MAXELM;
#if NONLIN
    x = pow((x - ACRTPT), 3) + ACRTPT;
    // This ^ is the equation used when NONLIN = 1; to have a non linear repartition. Feel free to change it. The inverse equation is in `int getoffset(double)`.
#endif
    return x;
}

// Func: var to offset conversion
int getoffset(double var)
{
#if NONLIN
    int offset = ((
        cbrt(var - ACRTPT) + ACRTPT
    // This ^ is the equation used when NONLIN = 1; to calculate offset with a non linear repartition. Feel free to change it (but it must be the inverse of the one in
    // `double getvar(int)` for this to work.)
                    ) - (double)MINVAL) * (double)MAXELM / ((double)MAXVAL - (double)MINVAL);
#else
    int offset = (var - (double)MINVAL) * (double)MAXELM / ((double)MAXVAL - (double)MINVAL);
#endif
#if STRICT
    if(getvar(offset) != var)
        return -1;
    return (offset < 0)?-1:(offset > (MAXELM - 1))?-1:offset;
#else
    return (offset < 0)?0:(offset > (MAXELM - 1))?MAXELM - 1:offset;
#endif
}

// Func: main.
int main(int argc, char* argv[])
{
    int offset;
    for(offset = 0; offset < MAXELM; offset++)
        results[offset] = log(getvar(offset));

#if PNTALL
    for(offset = 0; offset < MAXELM; offset++)
    {
        printf("[log(%lf) = %lf] ", getvar(offset), results[offset]);
        if(!((offset + 1) % 6))
            printf("\n");
    }
    printf("\n");
#endif

#if ASKFOR
    double x;
    printf("log(x) for x = ");
    scanf("%lf", &x);
    if((offset = getoffset(x)) < 0)
        printf("ERROR: Value for x = %lf hasn't been calculated\n", x);
    else
        printf("results[%d]: log(%lf) = %lf\n", offset, getvar(offset), results[offset]);
#endif

    return 0;
}

这个版本比以前的版本更干净、更容易维护。如果您还有什么需要，请留言。

您可以使用文件顶部的宏配置其行为。

特点:

• 使用固定大小的数组。
• 仅计算存储的值(不会计算一个数组单元格的多个值)。
• 使用函数从值获取偏移量以及从偏移量获取值，因此您不必存储已计算出 log 的值。

<小时/>

第二个版本

好吧，这是我的第二个解决方案。请参阅下面的原始评论。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MIN_INC 0.001    // This is the minimum increment. If its set to 0, when tmp will be equal to avg, it will never leave this state, since INC_MUL * (tmp - avg)^2 will be 0.
#define INC_MUL 0.2      // This is a number which influences the precision you will get. The smaller it is, the more precise you will be, and the greater will be your result array cardinality.

typedef struct {
    double offset;
    double value;    // value = log(offset). Since the results are not linarly widespread, this is pretty important.
} logCalc;

// Here, we need to use a pointer on a logCalc pointer, since we want to actually SET the address of the logCalc pointer, not the address of one of its copies.
int MyLogCreate(logCalc** arr, double min, double max)
{
    if((*arr) != NULL)
        return 0;
    unsigned int i = 0;
    double tmp, avg = (max + min) / 2.0;
    for( ; min < avg; min += (INC_MUL * ((avg - min) * (avg - min)) + MIN_INC))
    {
        (*arr) = (logCalc*)realloc((*arr), sizeof(logCalc) * (i + 1));
        (*arr)[i].offset  = min;
        (*arr)[i++].value = log(min);
    }
    for(tmp = avg ; tmp < max; tmp += (INC_MUL * ((tmp - avg) * (tmp - avg)) + MIN_INC))
    {
        (*arr) = (logCalc*)realloc((*arr), sizeof(logCalc) * (i + 1));
        (*arr)[i].offset  = tmp;
        (*arr)[i++].value = log(tmp);
    }
    return i;
}

int main(int argc, char** argv)
{
    logCalc *myloglut = NULL;
    unsigned int i,
        t = MyLogCreate(&myloglut, .1, 1.9);
    for(i = 0; i < (t-1); i++)
    {
        printf("[log(%lf) = %lf], ", myloglut[i].offset, myloglut[i].value);
        if(!((i+1)%6))         // Change 6 to what's best for your terminal $COLUMNS
            printf("\n");
    }
    printf("\n");
    free(myloglut);
    return 0;
}

<小时/>

原始评论

计算的线性度来自于您使用线性增量的事实。在 for 循环的每次迭代中，您将 exp 增加 (2.0 - 0.1)/MAXLOG。

要获得 0 附近更精确的值，您需要:

1. 定义更大的范围 - 更大的数组 -(能够存储更多 0 左右的值)
2. 使用非线性增量。此增量可能取决于 i (或取决于 exp，具体取决于您的操作方式)，因此您可以准确地知道您要尝试的数字的“偏移量”计算(以及您需要增加 exp 的量)。当然，你会在0左右计算出更多的结果。

这是我当前的实现:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define CALCULATE_UNTIL 2.0
#define PRECISE_UNTIL   1.0

typedef struct {
    double offset;
    double value;
} logCalc;

logCalc *myloglut = NULL;

int MyLogCreate()
{
    double exp = 0.1;
    int i;
    for (i = 0; exp <= CALCULATE_UNTIL; exp += (exp < PRECISE_UNTIL)?0.0001898:0.001898)
    {
        myloglut = realloc(myloglut, sizeof(logCalc) * (i + 1));
        myloglut[i].offset = exp;
        myloglut[i++].value = (i == 4780)?0:log(exp);
    }
    return i; // So you know how big the array is. Don't forget to free(myloglut); at the end of your code.
}

int main(int argc, char** argv)
{
    int i,
    t = MyLogCreate();
    for(i = 0; i < t; i++)
    {
        printf("[log(%lf) = %lf], ", myloglut[i].offset, myloglut[i].value);
        if(!(i%6))    // For formatting purposes.
            printf("\n");
    }
    printf("\n");
    free(myloglut);
    return 0;
}

我创建了一个新类型来存储 exp 的值，这对于了解结果的对数可能很有用。

更新:我不确定您想做什么。您想要精确到 log(x) = 0 左右还是 x = 0 左右？对于第一种情况，我可能必须再次重写代码才能使其按您的要求工作。另外，您是否希望结果在接近 0 时更加精确，或者您希望结果在给定范围内(就像现在一样)更加精确？