“ICSIlog”的编译时间表生成

Question

以下 C 代码用于在运行时生成查找表，以帮助实现“ICSI”日志algorithm（引用自https://github.com/mgbellemare/SkipCTS/blob/master/src/icsilog.cpp）：

    /*
This method fills a given array of floats with the information necessary to compute the icsi_log. This method has to be called before any call to icsi_log.
Parameters:
    n is the number of bits used from the mantissa (0<=n<=23). Higher n means higher accuracy but slower execution. We found that a good value for n is 14.
    lookup_table requires a float* pointing to a continuous (preallocated) memory array of 2^n*sizeof(float) bytes.
Return values: void
*/
void fill_icsi_log_table(const int n, float *lookup_table)
{
    float numlog;
    int incr,i,p;
    int *const exp_ptr = ((int*)&numlog);
    int x = *exp_ptr; /*x is the float treated as an integer*/
    x = 0x3F800000; /*set the exponent to 0 so numlog=1.0*/
        *exp_ptr = x;
    incr = 1 << (23-n); /*amount to increase the mantissa*/
    p = 1 << n;
    for(i=0;i<p;i++)
    {
        lookup_table[i] = (float) log2(numlog); /*save the log of the value*/
        x += incr;
        *exp_ptr = x; /*update the float value*/
    }
}


/* ICSIlog V 2.0 */
void fill_icsi_log_table2(const unsigned precision, float* const   pTable)
{
    /* step along table elements and x-axis positions
      (start with extra half increment, so the steps intersect at their midpoints.) */
    float oneToTwo = 1.0f + (1.0f / (float)( 1 <<(precision + 1) ));
    int i;
    for(i = 0;  i < (1 << precision);  ++i )
    {+
        // make y-axis value for table element
        pTable[i] = logf(oneToTwo) / 0.69314718055995f;

        oneToTwo += 1.0f / (float)( 1 << precision );
    }
}

有没有一种方法可以调整这些函数中的任何一个，以便在编译时使用模板和 C++11-amenable 单行返回 constexpr 函数来生成查找表，类似于以下结构？

/** Range generation,
 * from http://stackoverflow.com/questions/13313980/populate-an-array-using-constexpr-at-compile-time **/
template<unsigned... Is> struct seq{};

template<unsigned N, unsigned... Is>
struct gen_seq : gen_seq<N-1, N-1, Is...>{};

template<unsigned... Is>
struct gen_seq<0, Is...> : seq<Is...>{};

/** A table consisting of indexes and values,
 * which will all be computed at compile-time **/
template<unsigned N>
struct Table
{
    unsigned indexes[N];
    double  values[N];

    static constexpr unsigned length = N;
};


template< typename LambdaType, unsigned... Is>
constexpr Table< sizeof...(Is) > TableGenerator(seq<Is...>, LambdaType evalFunc)
{
    return {{ Is... }, { evalFunc(Is)... }};
}

template<unsigned N, typename LambdaType>
constexpr Table<N> TableGenerator( LambdaType evalFunc )
{
    return TableGenerator(gen_seq<N>(), evalFunc);
}



/** Function that computes a value for each index **/
constexpr double myFunc( unsigned idx )
{ 
    return sin(0.2 * idx) + cos(0.5*idx);
}

Answer 1

以this example 为起点和表格生成代码的“v2.0”变体：

  /* ICSIlog V 2.0 */
void fill_icsi_log_table2(const unsigned precision, float* const   pTable)
{
    /* step along table elements and x-axis positions
      (start with extra half increment, so the steps intersect at their midpoints.) */
    float oneToTwo = 1.0f + (1.0f / (float)( 1 <<(precision + 1) ));
    int i;
    for(i = 0;  i < (1 << precision);  ++i )
    {
        // make y-axis value for table element
        pTable[i] = logf(oneToTwo) / 0.69314718055995f;

        oneToTwo += 1.0f / (float)( 1 << precision );
    }
}

这个递归模板结构：

#include <math.h>

#define PRECISION (4)

constexpr float table_log(float oneToTwo)
{
    return logf(oneToTwo) / 0.69314718055995f;
}

template<size_t c, size_t precision,  float* const* pTable>
struct ForLoop {
    template<template <size_t, size_t,  float* const*> class Func>
    static void iterate(float oneToTwo) {
        ForLoop<c - 1, precision, pTable>::template 
        iterate<Func>(Func<c - 1, precision, pTable>()(oneToTwo));
    }
};

template<size_t precision,  float* const* pTable>
struct ForLoop<0, precision, pTable> {
    template<template <size_t, size_t,  float* const*> class Func>
    static void iterate(float oneToTwo) {
        Func<0, precision, pTable>()(oneToTwo);
    }
};

template <size_t index, size_t precision,  float* const *pTable>
struct LogTabe {
    float operator()(float oneToTwo) {
        float a = table_log(oneToTwo);
        (*pTable)[(1 << precision) - index] = a;
        return oneToTwo + 1.0f / (float)(1 << precision);
    }
};

static float *const table = new float[1 << PRECISION];
extern float *const table;

int main() {
    ForLoop<(1 << PRECISION) + 1, PRECISION, &table>::iterate<LogTabe>(1.0f + (1.0f / (float)( 1 << (PRECISION + 1))));
}

用gcc x86-64 8.1编译，-std=c++11 -O1，生成与原代码和asm输出一致的输出表：

  mov rax, QWORD PTR table[rip]
  mov DWORD PTR [rax], 0x3d35d69b
  mov DWORD PTR [rax+4], 0x3e0462c4
  mov DWORD PTR [rax+8], 0x3e567af2
  mov DWORD PTR [rax+12], 0x3e92203d
  mov DWORD PTR [rax+16], 0x3eb7110e
  mov DWORD PTR [rax+20], 0x3eda3f60
  mov DWORD PTR [rax+24], 0x3efbd42b
  mov DWORD PTR [rax+28], 0x3f0df989
  mov DWORD PTR [rax+32], 0x3f1d5da0
  mov DWORD PTR [rax+36], 0x3f2c2411
  mov DWORD PTR [rax+40], 0x3f3a58fe
  mov DWORD PTR [rax+44], 0x3f480731
  mov DWORD PTR [rax+48], 0x3f553848
  mov DWORD PTR [rax+52], 0x3f61f4e6
  mov DWORD PTR [rax+56], 0x3f6e44cd
  mov DWORD PTR [rax+60], 0x3f7a2f04
  mov DWORD PTR [rax+64], 0x3f88759c
  mov eax, 0
  ret
_GLOBAL__sub_I_main:
  sub rsp, 8
  mov edi, 64
  call operator new[](unsigned long)
  mov QWORD PTR table[rip], rax
  add rsp, 8
  ret

显示表值已在编译时成功预计算。然而，最近版本的 Clang 拒绝编译由 cmets 中 max66 给出的反对意见的代码，即“cmath”和“math.h”库函数不是严格的 constexpr（但因为无论如何它都是在编译时评估的，所以 Taylor作为 constexpr 函数实现的任意精度的级数扩展本身可能会很好地替代。）