在不使用浮点类型的情况下在 C 中构建对数函数 [英] Building a logarithm function in C without using float type
问题描述
我需要在不使用 float 类型的情况下重写日志函数(基本 2 或基本 10 无关紧要)，但我需要得到小数点后几位小数的精度.(就像一个 float * 100 在整数类型中获取 2 小数例如:如果 1.4352 将是结果，我的函数应该返回类似143(int 类型)，我知道最后 2 个数字是小数.
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:
但它们都返回 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?
推荐答案
为此您需要使用定点精度/算术/数学.这意味着您使用整数类型变量，但有些位在小数点之后.
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 我的 win calc 返回 6.97349648 非常接近.更精确的结果需要使用更多的小数位.int 和编译时评估浮点示例:
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:
[Edit2] fx32_mul 在没有 asm 基础的 32 位算术上 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屋！
