目前我计算日志如下：
#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];
}
如您所见，我只对 range 感兴趣0.1 - 2.0。但是，我需要更精确的0. 我怎样才能实现非线性计算？还有什么方法可以在这个函数中使用一些插值来获得更好的精度？
#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 附近获得更精确的值，您将需要：
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 时更精确，还是希望结果在给定范围内更精确（就像现在一样）？
将您的函数“原点”从零或 0.1 转移到 1.0。
for (i=-478;i<523;i++) {
double j = 1.0 + (double)i / 523.0;
myloglut[i+478] = log(j);
}
此函数仅选择两个点：1.0 和 2.0 作为 1.0 + (523.0 / 523.0) == 2.0。
那么第一个值是：
myloglut[0] = log(0.0860420650095602);
更“自然”的大小是 973，这将使除数为 512（确切地说），第一个条目将是 51/512 = ~0.099609375。
这需要多准确和多快？你可以用分段切比雪夫近似做一些很好的事情。您可以使用切比雪夫近似的阶数来控制精度（高阶 = 更慢但更准确）。我还建议通过将你的 double 分解为尾数（1 到 2 之间）和指数（2 的幂，其对数就是指数 times log(2)，你可以预先计算）来减少参数。
我不认为您可以在 [0.1, 2] 上提出任何非常准确的信息，而无需在需要日志时进行更多算术运算或使用巨大的表并引发所有缓存和不可预测的内存访问问题。但是，如果您有足够的时间，请考虑进行分段切比雪夫近似。（如果你想让我向你展示使用切比雪夫近似的代码，请在评论中告诉我，我会更新这篇文章。）
编辑：使用切比雪夫近似的对数代码。精确到 1e-5。
double mylog(double x) {
static double logtwo = log(2);
static double tbls[4][5] = {{0,0,0,0,0},
{-.9525934480e-2,-.87402539075,-1.135790603,1.1519051721,-1.7063112037},
{.53892330786e-3,-1.0117355213,-.4085197450,-.6242237228,0},
{.60393290013e-6,-1.0001523639,-.4940510719,-.4058961978,0}};
if (x>1) return -mylog(1/x);
int expo,e2;
x = 1-frexp(x, &expo);
double y = frexp(x, &e2);
if (e2 < -3) e2 = -3;
double *tbl = tbls[-e2];
return expo*logtwo + tbl[0]+x*(tbl[1]+x*(tbl[2]+x*(tbl[3]+x*tbl[4])));
}
我使用 Maple 计算了 Chebyshev 近似值，并将它们扩展为传统的多项式以提高速度。
如果您想要非常接近 1 的良好精度，您可以更改if (e2 < -3) e2 = -3该行并将其添加到对应于泰勒近似值{0,-1,-.5,-1/3.,-.25}的末尾。tbls如果您希望它更快，请在 1/2 和 3/4 之间计算一个更好的 log(x) 三次近似值，并将其存储在tbls.