不使用浮点型在C中建立对数函数 [英] Building a logarithm function in C without using float type
问题描述
我需要重写日志函数（基本 2 或基本 10 无关紧要）而无需使用 float 类型，但是我需要获得小数点后几位小数的精度。 （例如 float * 100 在整数类型中获取 2 个小数，例如：如果 1.4352 是结果，我的函数应该返回类似 143 （ int 类型）的内容，我知道最后两个数字是小数。
I need to rewrite the log function (base 2 or base 10 doesn't matter which) without using float type, but I need to get the precision of few decimal digits after the decimal point. ( like a float * 100 to get 2 decimals inside integer type eg: if the 1.4352 would be the result, my function should return something like 143 (int type) and I know that the last 2 numbers are the decimals.
我在stackoverflow上发现了一些方法，例如：
I found over the stackoverflow some methods like:
- How can I compute a base 2 logarithm without using the built-in math functions in C#?
但是它们全部返回 int 精度（避免小数点）。
but all of them returns int precision ( avoiding the decimals ).
我不知道如何解决这个问题的方法是：
I have no idea how to approach this so the question is:
如何编码（和或更改）整数 log 实现以支持十进制结果？
How to encode (and or change) integer log implementation to support decimal result?
推荐答案
您需要使用不动点精度/算术CS /数学。这意味着您使用整数类型变量，但某些位在小数点后。
You need to use fixed point precision/arithmetics/math for this. It means you use integer type variables but some of the bits are after the decimal point.
例如，假设8个小数位，则操作如下：
for example let assume 8 decimal bits so operations are done like this:
a = number1*256
b = number2*256
c=a+b // +
c=a-b // -
c=(a*b)>>8 // *
c=(a/b)<<8 // /
以下是通过 C ++ 进行二进制搜索的简单定点 log2 示例：
Here simple fixed point log2 example via binary search in 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;
}
//---------------------------------------------------------------------------
这里简单的测试（使用浮点数仅用于加载和打印，您也可以处理整数值的展位，或者由编译器评估常量）：
Here simple test (using floats just for loading and printing you can handle booth on integers too, or by compiler evaluated constants):
float(fx32_log2(float(125.67*float(_fx32_one)))) / float(_fx32_one)
其计算结果为： log2（125.67）= 6.98828125 我的获胜计算返回 6.97349648 非常接近。更精确的结果是您需要使用更多的小数位。整数和编译时间评估浮点示例：
This evaluates: log2(125.67) = 6.98828125 my win calc returns 6.97349648 which is pretty close. More precise result you need more fractional bits you need to use. Int and compile time evaluation float example:
(100*fx32_log2(125.67*_fx32_one))>>_fx32_fract_bits
返回 698 表示 6.98 乘以 100 。您还可以编写自己的加载和打印功能，以直接在定点和字符串之间进行转换。
returns 698 which means 6.98 as we multiplied by 100. You can also write your own load and print function to convert between fixed point and string directly.
要更改精度只需使用 _fx32_fract_bits 常数即可。无论如何，如果您的 C ++ 不知道 DWORD ，则它只是32位的 unsigned int 。如果您使用其他类型（例如 16 或 64 位），则只需相应地更改常量即可。
To change precision just play with _fx32_fract_bits constant. Anyway if your C++ does not know DWORD it is just 32bit unsigned int. If you are using different type (like 16 or 64 bit) then just change the constants accordingly.
有关更多信息，请查看：
For more info take a look at:
- Power by squaring for negative exponents
[Edit2] fx32_mul 基于32位算术，没有 asm 基准 2 ^ 16 O（n ^ 2）
fx32_mul on 32bit arithmetics without asm base 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 将此添加到代码开头
In case you do not have WORD,DWORD add this to start of code
typedef unsigned __int32 DWORD;
typedef unsigned __int16 WORD;
或以下方式：
typedef uint32_t DWORD;
typedef uint16_t WORD;
[Edit3] fx32_mul调试信息
对此进行调用并跟踪/断点（15个小数位）：
let call and trace/breakpoint this (15 fractional bits):
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
}
这篇关于不使用浮点型在C中建立对数函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
