Shift-and-add algorithms

计算对数

为了计算${\log }_{b}x,x>1$，我们可以预计算表格 (写在硬件上)：
$$A_k={\log }_b(1+\frac{1}{2^k}),k=0,1,\cdots ,n-1$$
然后可以做分解：
$${\log }_{b}x\approx {\log }_b\prod _{k=0}^{n-1}(1+\frac{1}{2^k}{)}^{a_k}=\sum _{k=0}^{n-1}{a}_k{A}_k$$

设置$x_0=1$，依次计算$x_i=x_{i-1}(1+\frac{1}{2^k})$数次，保证$x_i\le x$，以确定$a_k$的大小。直到$n$长的预计算表格用完。

绝对误差为：$O({\log }_b(1+\frac{1}{2^n}))$

一、预计算表格

void make_shiftadd_ltab(double b)
{
    double l1b = 1.0 / log(b);
    double s = 1.0;
    for (ulong k=0; k<ltab_n; ++k)
    {
        shiftadd_ltab[k] = log(1.0+s) * l1b;  // == log_b(1+1/2^k)
        s *= 0.5;
    }
}

二、计算对数${\log }_{b}x$

double shiftadd_log(double x, ulong n)
{
    if ( n>=ltab_n )  n = ltab_n;
    double t = 0.0;
    double e = 1.0;
    double v = 1.0;

    for (ulong k=1; k<n; ++k)
    {
        v *= 0.5;  // v == (1>>k)

        double u;
        bool d;
        while ( 1 )
        {
            u = e + e * v;  // u=e;  u+=(e>>k);
            d = ( u<=x );

            if ( d==false )  break;

            t += shiftadd_ltab[k];
            e = u;
        }
    }
    return  t;
}

计算指数

为了计算$b^x,b>1$，我们同样预计算表格 (写在硬件上)：
$$A_k={\log }_b(1+\frac{1}{2^k}),k=0,1,\cdots ,n-1$$
然后可以做分解：
$$b^x\approx \prod _{i=0}^{n-1}(1+\frac{1}{2^k}{)}^{a_k}=b^{{\sum }_{i=0}^{n-1}{a}_k{A}_k}$$

设置$x_0=0$，依次计算$x_i=x_{i-1}+(1+\frac{1}{2^k})$数次，保证$x_i\le x$，以确定$a_k$的大小。直到$n$长的预计算表格用完。

相对误差为：$O(1+\frac{1}{2^n})$

一、预计算表格

void make_shiftadd_ltab(double b)
{
    double l1b = 1.0 / log(b);
    double s = 1.0;
    for (ulong k=0; k<ltab_n; ++k)
    {
        shiftadd_ltab[k] = log(1.0+s) * l1b;  // == log_b(1+1/2^k)
        s *= 0.5;
    }
}

二、计算指数$b^x$

double shiftadd_exp(double x, ulong n)
{
    if ( n>=ltab_n )  n = ltab_n;
    double t = 0.0;
    double e = 1.0;
    double v = 1.0;

    for (ulong k=1; k<n; ++k)
    {
        v *= 0.5;  // v == (1>>k)

        double u;
        bool d;
        while ( 1 )
        {
            u = t + shiftadd_ltab[k];
            d = ( u<=x );

            if ( d==false )  break;

            t = u;
            e += e * v;  // e+=(e>>k);
        }
    }
    return  e;
}