【log计算的标量实现】

CMATH

---

计算机是如何实现log计算的

一、总体思路

所有的对数函数计算核心都是利用多项式展开（泰勒级数）然后多项式求和计算结果。为了性能或者精度的要求可能会对展开后的求和式子做进一步优化，最终会封装一个 函数出来。其余的对数函数都是使用换底公式来套 函数做的最底层实现，随着大量图形运算的需求提升， 函数实现得好不好直接决定你电脑快不快。

二、使用步骤

1.换底公式

可以看出所有的对数函数计算都可以转换成两个ln相除，比如：

double log10(double x) {
    return log(x) / log(10); // log == ln
}

所以我们只需要解决Ln函数就可以

2. ln函数实现源码

#include "fdlibm.h"

#ifdef __STDC__
static const double
#else
static double
#endif
ln2_hi  =  6.93147180369123816490e-01,	/* 3fe62e42 fee00000 */
ln2_lo  =  1.90821492927058770002e-10,	/* 3dea39ef 35793c76 */
two54   =  1.80143985094819840000e+16,  /* 43500000 00000000 */
Lg1 = 6.666666666666735130e-01,  /* 3FE55555 55555593 */
Lg2 = 3.999999999940941908e-01,  /* 3FD99999 9997FA04 */
Lg3 = 2.857142874366239149e-01,  /* 3FD24924 94229359 */
Lg4 = 2.222219843214978396e-01,  /* 3FCC71C5 1D8E78AF */
Lg5 = 1.818357216161805012e-01,  /* 3FC74664 96CB03DE */
Lg6 = 1.531383769920937332e-01,  /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01;  /* 3FC2F112 DF3E5244 */

static double zero   =  0.0;

#ifdef __STDC__
	double __ieee754_log(double x)
#else
	double __ieee754_log(x)
	double x;
#endif
{
	double hfsq,f,s,z,R,w,t1,t2,dk;
	int k,hx,i,j;
	unsigned lx;

	hx = __HI(x);		/* high word of x */
	lx = __LO(x);		/* low  word of x */

	k=0;
	if (hx < 0x00100000) {			/* x < 2**-1022  */
	    if (((hx&0x7fffffff)|lx)==0) 
		return -two54/zero;		/* log(+-0)=-inf */
	    if (hx<0) return (x-x)/zero;	/* log(-#) = NaN */
	    k -= 54; x *= two54; /* subnormal number, scale up x */
	    hx = __HI(x);		/* high word of x */
	} 
	if (hx >= 0x7ff00000) return x+x;
	k += (hx>>20)-1023;
	hx &= 0x000fffff;
	i = (hx+0x95f64)&0x100000;
	__HI(x) = hx|(i^0x3ff00000);	/* normalize x or x/2 */
	k += (i>>20);
	f = x-1.0;
        //  到这，第一步标记，得出 k, f, s 的值了
	if((0x000fffff&(2+hx))<3) {	/* |f| < 2**-20 */
	    if(f==zero) if(k==0) return zero;  else {dk=(double)k;
				 return dk*ln2_hi+dk*ln2_lo;}
	    R = f*f*(0.5-0.33333333333333333*f);
	    if(k==0) return f-R; else {dk=(double)k;
	    	     return dk*ln2_hi-((R-dk*ln2_lo)-f);}
	}
 	s = f/(2.0+f); 
	dk = (double)k;
	z = s*s;
	i = hx-0x6147a;
	w = z*z;
	j = 0x6b851-hx;
	t1= w*(Lg2+w*(Lg4+w*Lg6)); 
	t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); 
	i |= j;
	R = t2+t1;
	if(i>0) {
	    hfsq=0.5*f*f;
	    if(k==0) return f-(hfsq-s*(hfsq+R)); else
		     return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
	} else {
	    if(k==0) return f-s*(f-R); else
		     return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
	}
}

下面挑重点讲解一下

1. 参数约化（Argument Reduction）

将输入的数字 x 重整化为这个形式 $$x=2^k\times (1+f)$$
其中$$\frac{\sqrt{2}}{2}<(1+f)<\sqrt{2}$$以确保 1+f 小到足够合适，要不然之后换元的泰勒展开精度会很差。
然后很容易发现规律 $$lnx=\ln \left[2^k\times (1+f)\right]=k\ln 2+\ln (1+f)$$
其中 k*ln2口算即可，我们进一步来关注ln(1 + f) 的计算。

2.计算1+f的对数值

将第一步的 1 + f 进行计算，其结果因范围控制得比较小，所以很容易近似到一个很精确的值。计算方法：
取中间代换变量 $$s=\frac{f}{2+f}$$
（为什么这么代换呢？因为泰勒对 lnx 函数的展开在 x 接近 1 的时候最好），则有泰勒展开：
$$\begin{array}{rl}\ln (1+f)& =\ln (1+s)-\ln (1-s)\\ & =2s+\frac{2}{3}{s}^3+\frac{2}{5}{s}^5+\cdots \\ & =2s+s\cdot R(z)\end{array}$$
我们用一个 Remez 算法来求解 R(z) ，逼近到 14 次方项，此时的精度可以到 $2^{-58.45}$ 这个级别，完全足够满足 64 位浮点运算要求了，于是上式的 R 就可以等于：
$$R(z)\approx Lg1\cdot s+Lg2\cdot s+Lg3\cdot s+Lg4\cdot s+Lg5\cdot s+Lg6\cdot s+Lg7\cdot s$$
其中每一项都要带对应指数，式子可满足精度 $|(Lg1\cdot s+\cdots +Lg7\cdot s)-R(z)|<2^{-58.45}$，且$Lgi，i=1,2,3\cdots 7$都是常数，从上面的代码可以看到这些常数的值。
然后根据 $s=\frac{f}{2+f}$，我们有进一步的代换关系：$2s=\frac{2f}{2+f}=f-sf$
于是有：$$\begin{array}{rl}\ln (1+f)& =2s+s\cdot R(z)\\ & =f-sf+s\cdot R(z)\\ & =f-s(f-R(z))\end{array}$$
则：$$\ln (x)=k\ln 2+\ln (1+f)=k\ln 2+f-s(f-R(z))$$
另外，如果要求更高精度的话，重复代换关系可得：$2s=f-sf=f-\frac{ff}{2}+\frac{sff}{2}$
于是有：
$$\begin{array}{rl}\ln (1+f)& =2s+s\cdot R(z)\\ & =f-\frac{ff}{2}+\frac{sff}{2}+s\cdot R(z)\\ & =f-\frac{ff}{2}+s\left(\frac{ff}{2}+R(z)\right)\end{array}$$
则：
$\ln (x)=k\ln 2+\ln (1+f)=k\ln 2+f-\frac{ff}{2}+s\left(\frac{ff}{2}+R(z)\right)$

总览

其中的 ln2 是个常数， R(z) 项也是通过 7 个常数和 s 的值计算得到的，所以这个方法依赖了 8 个常数和 3 个变量 k, s, f ，变量关系为 $x=2^k\times (1+f),s=\frac{f}{2+f}$