opencv中的exp32f函式

Kkkzzzz發表於2021-10-17
OpenCV函式

exp32f

opencv的exp函式和cmath的exp函式在精度上存在一定差異，通過查詢原始碼，發現了這麼一段實現。程式碼如下：

點選檢視程式碼 
#define EXPTAB_SCALE 6
#define EXPTAB_MASK  ((1 << EXPTAB_SCALE) - 1)
 
#define EXPPOLY_32F_A0 .9670371139572337719125840413672004409288e-2
 
static const double expTab[EXPTAB_MASK + 1] = {
    1.0 * EXPPOLY_32F_A0,
    1.0108892860517004600204097905619 * EXPPOLY_32F_A0,
    1.0218971486541166782344801347833 * EXPPOLY_32F_A0,
    1.0330248790212284225001082839705 * EXPPOLY_32F_A0,
    1.0442737824274138403219664787399 * EXPPOLY_32F_A0,
    1.0556451783605571588083413251529 * EXPPOLY_32F_A0,
    1.0671404006768236181695211209928 * EXPPOLY_32F_A0,
    1.0787607977571197937406800374385 * EXPPOLY_32F_A0,
    1.0905077326652576592070106557607 * EXPPOLY_32F_A0,
    1.1023825833078409435564142094256 * EXPPOLY_32F_A0,
    1.1143867425958925363088129569196 * EXPPOLY_32F_A0,
    1.126521618608241899794798643787 * EXPPOLY_32F_A0,
    1.1387886347566916537038302838415 * EXPPOLY_32F_A0,
    1.151189229952982705817759635202 * EXPPOLY_32F_A0,
    1.1637248587775775138135735990922 * EXPPOLY_32F_A0,
    1.1763969916502812762846457284838 * EXPPOLY_32F_A0,
    1.1892071150027210667174999705605 * EXPPOLY_32F_A0,
    1.2021567314527031420963969574978 * EXPPOLY_32F_A0,
    1.2152473599804688781165202513388 * EXPPOLY_32F_A0,
    1.2284805361068700056940089577928 * EXPPOLY_32F_A0,
    1.2418578120734840485936774687266 * EXPPOLY_32F_A0,
    1.2553807570246910895793906574423 * EXPPOLY_32F_A0,
    1.2690509571917332225544190810323 * EXPPOLY_32F_A0,
    1.2828700160787782807266697810215 * EXPPOLY_32F_A0,
    1.2968395546510096659337541177925 * EXPPOLY_32F_A0,
    1.3109612115247643419229917863308 * EXPPOLY_32F_A0,
    1.3252366431597412946295370954987 * EXPPOLY_32F_A0,
    1.3396675240533030053600306697244 * EXPPOLY_32F_A0,
    1.3542555469368927282980147401407 * EXPPOLY_32F_A0,
    1.3690024229745906119296011329822 * EXPPOLY_32F_A0,
    1.3839098819638319548726595272652 * EXPPOLY_32F_A0,
    1.3989796725383111402095281367152 * EXPPOLY_32F_A0,
    1.4142135623730950488016887242097 * EXPPOLY_32F_A0,
    1.4296133383919700112350657782751 * EXPPOLY_32F_A0,
    1.4451808069770466200370062414717 * EXPPOLY_32F_A0,
    1.4609177941806469886513028903106 * EXPPOLY_32F_A0,
    1.476826145939499311386907480374 * EXPPOLY_32F_A0,
    1.4929077282912648492006435314867 * EXPPOLY_32F_A0,
    1.5091644275934227397660195510332 * EXPPOLY_32F_A0,
    1.5255981507445383068512536895169 * EXPPOLY_32F_A0,
    1.5422108254079408236122918620907 * EXPPOLY_32F_A0,
    1.5590044002378369670337280894749 * EXPPOLY_32F_A0,
    1.5759808451078864864552701601819 * EXPPOLY_32F_A0,
    1.5931421513422668979372486431191 * EXPPOLY_32F_A0,
    1.6104903319492543081795206673574 * EXPPOLY_32F_A0,
    1.628027421857347766848218522014 * EXPPOLY_32F_A0,
    1.6457554781539648445187567247258 * EXPPOLY_32F_A0,
    1.6636765803267364350463364569764 * EXPPOLY_32F_A0,
    1.6817928305074290860622509524664 * EXPPOLY_32F_A0,
    1.7001063537185234695013625734975 * EXPPOLY_32F_A0,
    1.7186192981224779156293443764563 * EXPPOLY_32F_A0,
    1.7373338352737062489942020818722 * EXPPOLY_32F_A0,
    1.7562521603732994831121606193753 * EXPPOLY_32F_A0,
    1.7753764925265212525505592001993 * EXPPOLY_32F_A0,
    1.7947090750031071864277032421278 * EXPPOLY_32F_A0,
    1.8142521755003987562498346003623 * EXPPOLY_32F_A0,
    1.8340080864093424634870831895883 * EXPPOLY_32F_A0,
    1.8539791250833855683924530703377 * EXPPOLY_32F_A0,
    1.8741676341102999013299989499544 * EXPPOLY_32F_A0,
    1.8945759815869656413402186534269 * EXPPOLY_32F_A0,
    1.9152065613971472938726112702958 * EXPPOLY_32F_A0,
    1.9360617934922944505980559045667 * EXPPOLY_32F_A0,
    1.9571441241754002690183222516269 * EXPPOLY_32F_A0,
    1.9784560263879509682582499181312 * EXPPOLY_32F_A0,
};
 
static const double exp_prescale = 1.4426950408889634073599246810019 * (1 << EXPTAB_SCALE);
static const double exp_postscale = 1./(1 << EXPTAB_SCALE);
static const double exp_max_val = 3000.*(1 << EXPTAB_SCALE); // log10(DBL_MAX) < 3000
 
const double* getExpTab64f()
{
    return expTab;
}
 
const float* getExpTab32f()
{
    static float expTab_f[EXPTAB_MASK+1];
    static std::atomic<bool> expTab_f_initialized(false);
    if (!expTab_f_initialized.load())
    {
        for( int j = 0; j <= EXPTAB_MASK; j++ )
            expTab_f[j] = (float)expTab[j];
        expTab_f_initialized = true;
    }
    return expTab_f;
}
 
void exp32f( const float *_x, float *y, int n )
{
    const float* const expTab_f = getExpTab32f();
 
    const float
    A4 = (float)(1.000000000000002438532970795181890933776 / EXPPOLY_32F_A0),
    A3 = (float)(.6931471805521448196800669615864773144641 / EXPPOLY_32F_A0),
    A2 = (float)(.2402265109513301490103372422686535526573 / EXPPOLY_32F_A0),
    A1 = (float)(.5550339366753125211915322047004666939128e-1 / EXPPOLY_32F_A0);
 
    int i = 0;
    const Cv32suf* x = (const Cv32suf*)_x;
    float minval = (float)(-exp_max_val/exp_prescale);
    float maxval = (float)(exp_max_val/exp_prescale);
    float postscale = (float)exp_postscale;
 
    for( ; i < n; i++ )
    {
        float x0 = x[i].f;
        x0 = std::min(std::max(x0, minval), maxval);
        x0 *= (float)exp_prescale;
        Cv32suf buf;
 
        int xi = cv::saturate_cast<int>(x0);
        x0 = (x0 - xi)*postscale;
 
        int t = (xi >> EXPTAB_SCALE) + 127;
        t = !(t & ~255) ? t : t < 0 ? 0 : 255;
        buf.i = t << 23;
 
        y[i] = buf.f * expTab_f[xi & EXPTAB_MASK] * ((((x0 + A1)*x0 + A2)*x0 + A3)*x0 + A4);
    }
}

這是從opencv4原始碼中摘出的一段計算\(e^x\)的函式（32位浮點）。初看起來可能有些迷糊，但仔細思考和推導之後，稍微明白了一些。想要理解這段程式碼，主要有三點需要清楚：第一點，從計算\(e^x\)轉換到計算\(2^x\)；第二點，使用四階泰勒展開逼近\(2^x\); 第三點，expTab_f的作用。

從\(e^x\)轉到\(2^x\)。這個比較簡單，對輸入\(x\)預先乘以\(log_2(e)\)就可以了，exp_prescale就是起這樣的作用（這裡預先左移了6位，相當於乘以64，這個先不管）。轉換之後對\(x\)進行了整數和小數部分的分離，整數部分為xi,小數部分為x0。對於整數部分，按照IEEE754浮點數的表示方式，可以直接移動到階碼部分。

int t = (xi >> EXPTAB_SCALE) + 127;
t = !(t & ~255) ? t : t < 0 ? 0 : 255;
buf.i = t << 23;

這個時候整數部分已經計算完成（即buf.f，這裡用了union的儲存方式，方便機器數級別的浮點和整數的轉換），剩下只要計算小數部分，即\(2^{x0}\)。

y[i] = buf.f * expTab_f[xi & EXPTAB_MASK] * ((((x0 + A1)x0 + A2)x0 + A3)*x0 + A4);

小數部分的計算不是很直接，這裡用到了級數近似的方式。

\(2^x\)的多項式展開。對\(2^x\)進行泰勒展開可以得到

\(2^x = 1 + xln2 + \frac{x^2ln^{2}\ 2}{2} + \frac{x^3ln^{3}\ 2}{3!} + \frac{x^4ln^{4}\ 2}{4!} + \cdots\)

\(A0 = \frac{ln^{4}\ 2}{4!},\)

\(A1 = \frac{ln^{3}\ 2}{3!A0},\)

\(A2 = \frac{ln^{2}\ 2}{2!A0},\)

\(A3 = \frac{ln2}{A0},\)

\(A4 = \frac{1}{A0}.\)

將以上代入((((x0 + A1)*x0 + A2)*x0 + A3)*x0 + A4)可以得到\(2^x\)的四階多項式展開，不過每項係數都多除以了一個A0。這時候就需要expTab_f了。

expTab_f。expTab_f有64項，其作用除了消去\(2^x\)多項式展開中的A0，還有就是乘上之前忽略的小數部分。一開始，輸入x乘以64，相當於小數點向右移動了6位，因而有6位小數移動到了整數部分（按照二進位制表示來考慮），不過在計算整數階碼時沒有考慮這部分，所以之後需要把這6位小數重新乘上去，根據整數部分低6位的值計算對應的小數部分，這就是expTab_f所表示的內容。比如整數的低6位為1時，表示小數\(2^{1/64}=1.0108892860517005\)，低6位為3時，表示小數\(2^{1/32+1/64}=1.0330248790212284\)，低6位為63時（所有位為1），表示小數\(2^{1-1/64}=1.978456026387951\)。

一些疑問

至此，這段程式碼基本就搞清楚了。對此，還有幾個問題：第一，為什麼需要預先左移6位，從運算角度考慮並沒有起到什麼實質作用；第二，為什麼\(2^x\)多項式展開的係數都要除以A0，如果不除會有什麼不同。對於第一個問題，我猜可能是精度考慮，移動到整數部分之後，這幾位的小數精度可以保證，而且移動的位數越多，近似的精度越高。對於第二個問題，我猜也可能跟精度有關，但並不太清楚是怎樣影響的。

相關文章