【问题标题】:Building a logarithm function in C without using float type在不使用浮点类型的情况下在 C 中构建对数函数
【发布时间】:2017-06-25 18:01:18
【问题描述】:

我需要在不使用float类型的情况下重写日志函数(基数2或基数10无所谓),但我需要获得小数点后几位小数的精度。 (就像float * 100 得到整数类型内的2 小数例如:如果1.4352 是结果,我的函数应该返回类似143int 类型)的东西,我知道最后两个数字是小数。

我在 stackoverflow 上发现了一些方法,例如:

但它们都返回 int 精度(避免小数)。

我不知道如何解决这个问题,所以问题是:

如何编码(和或更改)整数log 实现以支持十进制结果?

【问题讨论】:

  • 这是作业吗?
  • 不,不是。我想在个人项目中使用它。我做了一个函数女巫检查最高位的位置,但使用它我只返回结果的整数部分,忽略 zecimals,我需要至少 2 个......它将帮助我获得指数函数。 P.S:我不必使用任何库
  • 嗯 log2(1024^100) 给你 log2(1024)*100 但这不会有太大帮助(因为你不能在 C 中有效地计算它)......
  • 你需要问一个具体的问题。是什么阻碍了您编写代码?展示我们的尝试,解释其中的问题并提出一个具体问题,以帮助您改进代码。
  • 我敢打赌,你是在没有 FPU 的 MCU 上。无论如何都要使用固定精度和二进制搜索

标签: c algorithm logarithm


【解决方案1】:

您需要为此使用定点精度/算术/数学。这意味着您使用整数类型变量,但有些位在小数点之后。

例如,假设有 8 个十进制位,因此操作如下:

a = number1*256
b = number2*256
c=a+b // +
c=a-b // -
c=(a*b)>>8 // *
c=(a/b)<<8 // /

这里是简单的定点log2示例,通过C++中的二分查找:

//---------------------------------------------------------------------------
const DWORD _fx32_bits      =32;                            // all bits count
const DWORD _fx32_fract_bits= 8;                            // fractional bits count
const DWORD _fx32_integ_bits=_fx32_bits-_fx32_fract_bits;   // integer bits count
//---------------------------------------------------------------------------
const DWORD _fx32_one       =1<<_fx32_fract_bits;           // constant=1.0 (fixed point)
const DWORD _fx32_fract_mask=_fx32_one-1;                   // fractional bits mask
const DWORD _fx32_integ_mask=0xFFFFFFFF-_fx32_fract_mask;   // integer bits mask
const DWORD _fx32_MSB_mask=1<<(_fx32_bits-1);               // max unsigned bit mask
//---------------------------------------------------------------------------
DWORD bits(DWORD p)             // count how many bits is p
    {
    DWORD m=0x80000000; DWORD b=32;
    for (;m;m>>=1,b--)
     if (p>=m) break;
    return b;
    }
//---------------------------------------------------------------------------
DWORD fx32_mul(DWORD x,DWORD y)
    {
    // this should be done in asm with 64 bit result !!!
    DWORD a=x,b=y;              // asm has access only to local variables
    asm {                       // compute (a*b)>>_fx32_fract
        mov eax,a               // eax=a
        mov ebx,b               // ebx=b
        mul eax,ebx             // (edx,eax)=eax*ebx
        mov ebx,_fx32_one
        div ebx                 // eax=(edx,eax)>>_fx32_fract
        mov a,eax;
        }
    return a;
    // you can also do this instead but unless done on 64bit variable will overflow
    return (x*y)>>_fx32_fract_bits;
    }
//---------------------------------------------------------------------------
DWORD fx32_sqrt(const DWORD &x) // unsigned fixed point sqrt
    {
    DWORD m,a;
    if (!x) return 0;
    m=bits(x);                  // integer bits
    if (m>_fx32_fract_bits) m-=_fx32_fract_bits; else m=0;
    m>>=1;                      // sqrt integer result is half of x integer bits
    m=_fx32_one<<m;             // MSB of result mask
    for (a=0;m;m>>=1)           // test bits from MSB to 0
        {
        a|=m;                   // bit set
        if (fx32_mul(a,a)>x)    // if result is too big
         a^=m;                  // bit clear
        }
    return a;
    }
//---------------------------------------------------------------------------
DWORD fx32_exp2(DWORD y)       // 2^y
    {
    // handle special cases
    if (!y) return _fx32_one;                    // 2^0 = 1
    if (y==_fx32_one) return 2;                  // 2^1 = 2
    DWORD m,a,b,_y;
    // handle the signs
    _y=y&_fx32_fract_mask;      // _y fractional part of exponent
     y=y&_fx32_integ_mask;      //  y integer part of exponent
    a=_fx32_one;                // ini result
    // powering by squaring x^y
    if (y)
        {
        for (m=_fx32_MSB_mask;(m>_fx32_one)&&(m>y);m>>=1);     // find mask of highest bit of exponent
        for (;m>=_fx32_one;m>>=1)
            {
            a=fx32_mul(a,a);
            if (DWORD(y&m)) a<<=1;  // a*=2
            }
        }
    // powering by rooting x^_y
    if (_y)
        {
        for (b=2<<_fx32_fract_bits,m=_fx32_one>>1;m;m>>=1)      // use only fractional part
            {
            b=fx32_sqrt(b);
            if (DWORD(_y&m)) a=fx32_mul(a,b);
            }
        }
    return a;
    }
//---------------------------------------------------------------------------
DWORD fx32_log2(DWORD x)    // = log2(x)
    {
    DWORD y,m;
    // binary search from highest possible integer power of 2 to avoid overflows (log2(integer bits)-1)
    for (y=0,m=_fx32_one<<(bits(_fx32_integ_bits)-1);m;m>>=1)
        {
        y|=m;   // set bit
        if (fx32_exp2(y)>x) y^=m; // clear bit if result too big
        }
    return y;
    }
//---------------------------------------------------------------------------

这里是简单的测试(使用浮点数仅用于加载和打印,您也可以处理整数的展位,或通过编译器评估的常量):

float(fx32_log2(float(125.67*float(_fx32_one)))) / float(_fx32_one)

这会评估:log2(125.67) = 6.98828125 我的获胜计算返回 6.97349648,这非常接近。更精确的结果,您需要更多需要使用的小数位。 Int 和编译时评估浮点示例:

(100*fx32_log2(125.67*_fx32_one))>>_fx32_fract_bits

返回698,这意味着6.98乘以100。您也可以编写自己的加载和打印函数,直接在fixed pointstring之间进行转换。

要更改精度,只需使用_fx32_fract_bits 常量即可。无论如何,如果您的 C++ 不知道 DWORD 它只是 32 位 unsigned int。如果您使用不同的类型(如1664 位),则只需相应地更改常量即可。

欲了解更多信息,请查看:

[Edit2] fx32_mul 在没有asm 基础2^16 O(n^2) 的 32 位算术上2^16 O(n^2)

DWORD fx32_mul(DWORD x,DWORD y)
    {
    const int _h=1; // this is MSW,LSW order platform dependent So swap 0,1 if your platform is different
    const int _l=0;
    union _u
        {
        DWORD u32;
        WORD u16[2];
        }u;
    DWORD al,ah,bl,bh;
    DWORD c0,c1,c2,c3;
    // separate 2^16 base digits
    u.u32=x; al=u.u16[_l]; ah=u.u16[_h];
    u.u32=y; bl=u.u16[_l]; bh=u.u16[_h];
    // multiplication (al+ah<<1)*(bl+bh<<1) = al*bl + al*bh<<1 + ah*bl<<1 + ah*bh<<2
    c0=(al*bl);
    c1=(al*bh)+(ah*bl);
    c2=(ah*bh);
    c3= 0;
    // propagate 2^16 overflows (backward to avoid overflow)
    c3+=c2>>16; c2&=0x0000FFFF;
    c2+=c1>>16; c1&=0x0000FFFF;
    c1+=c0>>16; c0&=0x0000FFFF;
    // propagate 2^16 overflows (normaly to recover from secondary overflow)
    c2+=c1>>16; c1&=0x0000FFFF;
    c3+=c2>>16; c2&=0x0000FFFF;
    // (c3,c2,c1,c0) >> _fx32_fract_bits
    u.u16[_l]=c0; u.u16[_h]=c1; c0=u.u32;
    u.u16[_l]=c2; u.u16[_h]=c3; c1=u.u32;
    c0 =(c0&_fx32_integ_mask)>>_fx32_fract_bits;
    c0|=(c1&_fx32_fract_mask)<<_fx32_integ_bits;
    return c0;
    }

如果您没有WORD,DWORD,请将其添加到代码开头

typedef unsigned __int32 DWORD;
typedef unsigned __int16 WORD;

或者这个:

typedef uint32_t DWORD;
typedef uint16_t WORD;

[Edit3] fx32_mul 调试信息

调用并跟踪/断点这个(15 个小数位):

fx32_mul(0x00123400,0x00230056);

这是:

0x00123400/32768 * 0x00230056/32768 =
36 * 70.00262451171875 = 2520.094482421875

所以:

DWORD fx32_mul(DWORD x,DWORD y) // x=0x00123400 y=0x00230056
    {
    const int _h=1;
    const int _l=0;
    union _u
        {
        DWORD u32;
        WORD u16[2];
        }u;
    DWORD al,ah,bl,bh;
    DWORD c0,c1,c2,c3;
    // separate 2^16 base digits
    u.u32=x; al=u.u16[_l]; ah=u.u16[_h]; // al=0x3400 ah=0x0012
    u.u32=y; bl=u.u16[_l]; bh=u.u16[_h]; // bl=0x0056 bh=0x0023
    // multiplication (al+ah<<1)*(bl+bh<<1) = al*bl + al*bh<<1 + ah*bl<<1 + ah*bh<<2
    c0=(al*bl);        // c0=0x00117800
    c1=(al*bh)+(ah*bl);// c1=0x0007220C
    c2=(ah*bh);        // c2=0x00000276
    c3= 0;             // c3=0x00000000
    // propagate 2^16 overflows (backward to avoid overflow)
    c3+=c2>>16; c2&=0x0000FFFF; // c3=0x00000000 c2=0x00000276
    c2+=c1>>16; c1&=0x0000FFFF; // c2=0x0000027D c1=0x0000220C
    c1+=c0>>16; c0&=0x0000FFFF; // c1=0x0000221D c0=0x00007800
    // propagate 2^16 overflows (normaly to recover from secondary overflow)
    c2+=c1>>16; c1&=0x0000FFFF; // c2=0x0000027D c1=0x0000221D
    c3+=c2>>16; c2&=0x0000FFFF; // c3=0x00000000 c2=0x0000027D
    // (c3,c2,c1,c0) >> _fx32_fract_bits
    u.u16[_l]=c0; u.u16[_h]=c1; c0=u.u32; // c0=0x221D7800
    u.u16[_l]=c2; u.u16[_h]=c3; c1=u.u32; // c1=0x0000027D
    c0 =(c0&_fx32_integ_mask)>>_fx32_fract_bits; // c0=0x0000443A
    c0|=(c1&_fx32_fract_mask)<<_fx32_integ_bits; // c0=0x04FA443A
    return c0; // 0x04FA443A -> 83510330/32768 = 2548.53302001953125
    }

【讨论】:

  • 评论不用于扩展讨论;这个对话是moved to chat
  • 在 log2 期间,无需运行完整的幂运算:您已经有一堆先前计算过的位 - 在每一步中您所需要的只是一个试验定点乘以表中的常数。事实上,您既不需要求幂也不需要平方根。只是相乘、相加和比较。但是此时您可能会发现,如果使用模拟浮点的典型方式,您将获得更少的乘数,因此从浮点软件库中获取逐字代码并降低其准确性会产生更快的代码!
猜你喜欢
  • 2015-07-19
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2011-03-11
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多