抱歉不是 python 编码器,但这里是通用方法(不受库或语言限制):
-
定义
所以你有 2 个(或 N)浮点数 a,b 并希望有 2 个整数 aa,bb 这样:
a/b == aa/bb
-
接近
浮点数只是整数尾数,以 2 为底的整数指数向左(如果为负指数,则向右)移动,所以:
a = sign(a)*mantisa(a)*2^exponent(a) = sign(a)*(mantisa(a)<<exponent(a))
b = sign(b)*mantisa(b)*2^exponent(b) = sign(b)*(mantisa(b)<<exponent(b))
所以如果我们移动两个 a,b 数字,那么较大数量的尾数的 msb(最高有效位)将转到某个整数变量的 msb将a,b 转换为整数而不改变它们的比率(除非由于目标变量数据类型的位宽较小而减少了尾数的某些位)。这就像将数字与相同的常数相乘。
-
从a,b中提取指数
这可以通过直接将指数位提取为整数并从中减去偏差以使其有符号或使用log2()数学函数来简单地完成。
-
计算shift
我们需要将a,b 的尾数位移shift 位或将a,b 乘以2^shift,这样更大的幅度数将是最大的,它仍然适合整数变量。因此,如果我假设 32 位有符号整数,我们希望较大数量的 msb 将是 30(位从 0 编号,我们希望保留最后一位,所以我们仍然可以申请标志)。
计算很简单:
shift=max( exponent(a), exponent(b) );
shift=30-shift;
// shift-=_f32_man_bits; // this is just in case of bit-shifting
-
位移或乘以a,b并构造结果
所以只需将a,b 转换为整数,如上一个项目符号中所述。之后,您可以将结果整数除以它们的 GCD 或将它们右移直到 a 或 b 的 lsb 非零(删除尾随零)。
这里是二进制的小例子:
exponent(b)=2 exponent(a)=-3
| |
| 0.0010101110b <- a
100.01101b <- b
--------------------------------------------------------------------------
_f32_man_bits = 23 // 32 bit float has 24 bit mantisa but first one is implicit
shift = 30 - max(exponent(b),exponent(a)) = 30 - 2 = 28
--------------------------------------------------------------------------
????????????????????????????????.0000000000b <- 32 bit integer variable
00000010101110000000000000000000.0000000000b <- a * (1 << shift) = mantissa(a)|(1<<_f32_man_bits) << (shift + exponent(a) - _f32_man_bits)
01000110100000000000000000000000.0000000000b <- b * (1 << shift) = mantissa(b)|(1<<_f32_man_bits) << (shift + exponent(b) - _f32_man_bits)
|
msb is zero so sign can still be applied ...
删除尾随零可以这样完成:
// remove trailing zeros
for (;((aa|bb)&1)==0;)
{
aa>>=1;
bb>>=1;
}
上面的例子会变成:
0000001010111b
0100011010000b
GCD 的除法可以这样完成(去除尾随零之后):
// divide by GCD
for (int d=3;(d<=aa)&&(d<=bb);d+=2)
while ((aa%d)+(bb%d)==0)
{ aa/=d; bb/=d; }
最后应用符号。
这里是 C++ 浮动示例(乘法):
void f32_ratio0(int &aa,int &bb,float a,float b) // aa/bb = a/b
{
// IEEE 754 constants
const DWORD _f32_man_bits=23; // mantisa bits (without implicit one)
// variables
int shift,d;
int expa,siga;
int expb,sigb;
// extract parts of a,b
siga=(a<0.0); a=fabs(a); sigb=(b<0.0); b=fabs(b);
expa=floor(log(a)/log(2.0)); expb=floor(log(b)/log(2.0));
// compute shift
shift=expa; if (shift<expb) shift=expb; // max(expa,expb)
shift=30-shift; // shift msb of bigger mantisa to 30th bit of integer
// construct result
aa=float(a*pow(2.0,shift));
bb=float(b*pow(2.0,shift));
// remove trailing zeros
for (;((aa|bb)&1)==0;)
{
aa>>=1;
bb>>=1;
}
// divide by GCD
for (d=3;(d<=aa)&&(d<=bb);d+=2)
while ((aa%d)+(bb%d)==0)
{ aa/=d; bb/=d; }
// sign
if (siga) aa=-aa;
if (sigb) bb=-bb;
}
这里是 C++ 整数示例(移位):
void f32_ratio1(int &aa,int &bb,float a,float b) // aa/bb = a/b
{
// IEEE 754 constants
const DWORD _f32_sig =0x80000000; // sign
const DWORD _f32_exp =0x7F800000; // exponent
const DWORD _f32_exp_sig=0x40000000; // exponent sign
const DWORD _f32_exp_bia=0x3F800000; // exponent bias
const DWORD _f32_exp_lsb=0x00800000; // exponent LSB
const DWORD _f32_man =0x007FFFFF; // mantisa
const DWORD _f32_man_msb=0x00400000; // mantisa MSB
const DWORD _f32_man_bits=23; // mantisa bits (without implicit one)
const DWORD _f32_exp_bias=127; // exponent bias
// float bits access
union
{
float f; // 32bit floating point
DWORD u; // 32 bit uint
} y;
// variables
int shift,d;
int mana,expa,siga;
int manb,expb,sigb;
// extract parts of a
y.f=a;
mana=(y.u&_f32_man)|_f32_exp_lsb;
expa=((y.u&_f32_exp)>>_f32_man_bits)-_f32_exp_bias;
siga=(y.u&_f32_sig);
// extract parts of b
y.f=b;
manb=(y.u&_f32_man)|_f32_exp_lsb;
expb=((y.u&_f32_exp)>>_f32_man_bits)-_f32_exp_bias;
sigb=(y.u&_f32_sig);
// compute shift
shift=expa; if (shift<expb) shift=expb; // max(expa,expb)
shift=(30-_f32_man_bits)-shift; // shift msb of bigger mantisa to 30th bit of integer
// construct result
d=shift+expa; aa=mana; if (d<0) aa>>=-d; else if (d>0) aa<<=d;
d=shift+expb; bb=manb; if (d<0) bb>>=-d; else if (d>0) bb<<=d;
// remove trailing zeros
for (;((aa|bb)&1)==0;)
{
aa>>=1;
bb>>=1;
}
// divide by GCD
for (d=3;(d<=aa)&&(d<=bb);d+=2)
while ((aa%d)+(bb%d)==0)
{ aa/=d; bb/=d; }
// sign
if (siga) aa=-aa;
if (sigb) bb=-bb;
}
其中DWORD 是任何无符号的 32 位数据类型,例如:
typedef unsigned __int32 DWORD;
double 精度将以相同的方式完成,只有常量发生变化,并且需要 64bit 或 2x32bit 变量来存储整数尾数和结果...
精度取决于指数的相对距离。如果数字相差太大,则生成的数字将不适合目标整数,从而导致较小的幅度数字在以下情况下转换为零:
abs( exponent(a) - exponent(b) ) >= 31
同样,如果整数使用更大的位宽,则 31 将相应地改变 ...
现在你的例子:
// a b a/b
0.50000 / 1.00000 = 0.500000 // floats
// aa bb aa/bb
1 / 2 = 0.500000 // ratio0
1 / 2 = 0.500000 // ratio1
// a b a/b
0.50000 / 0.60000 = 0.833333 // floats
// aa bb aa/bb
4194304 / 5033165 = 0.833333 // ratio0
4194304 / 5033165 = 0.833333 // ratio1
请注意,0.6 不能完全用浮点数表示,因此 aa,bb 的值很大!!!要解决这个问题,您需要添加舍入,但为此您需要知道阈值告诉你数字的哪一部分要四舍五入...在不知道浮点数的目标范围或精度的情况下,我无法安全地实现它...
如果您想保留更多浮点数之间的比率,而不是简单地将它们添加到函数中。
这里是 3 个变量的浮动 C++ 示例:
void f32_ratio0(int &aa,int &bb,int &cc,float a,float b,float c) // aa/bb/cc = a/b/c
{
// IEEE 754 constants
const DWORD _f32_man_bits=23; // mantisa bits (without implicit one)
// variables
int shift,d;
int expa,siga;
int expb,sigb;
int expc,sigc;
// extract parts of a,b
siga=(a<0.0); a=fabs(a); sigb=(b<0.0); b=fabs(b); sigc=(c<0.0); c=fabs(c);
expa=floor(log(a)/log(2.0)); expb=floor(log(b)/log(2.0)); expc=floor(log(c)/log(2.0));
// compute shift
shift=expa; // max(expa,expb)
if (shift<expb) shift=expb;
if (shift<expc) shift=expc;
shift=30-shift; // shift msb of bigger mantisa to 30th bit of integer
// construct result
aa=float(a*pow(2.0,shift));
bb=float(b*pow(2.0,shift));
cc=float(c*pow(2.0,shift));
// remove trailing zeros
for (;((aa|bb|cc)&1)==0;)
{
aa>>=1;
bb>>=1;
cc>>=1;
}
// divide by GCD
for (d=3;(d<=aa)&&(d<=bb)&&(d<=cc);d+=2)
while ((aa%d)+(bb%d)+(cc%d)==0)
{ aa/=d; bb/=d; cc/=d; }
// sign
if (siga) aa=-aa;
if (sigb) bb=-bb;
if (sigc) cc=-cc;
}
和您的示例结果:
// a b c
0.50000 / 0.60000 / 1.00000
// aa bb cc
4194304 / 5033165 / 8388608
[Edit1]N案例算法
-
提取N的部分浮点数O(N)
所以我们有浮点数 a0,a1,a2,...,a(N-1) 并需要整数指数 e0,e1,... 和尾数 m0,m1,... 和符号 s0,s1,...。对于 32 位浮点数(使用上述示例中的 // IEEE 754 常量):
int i,m[N],e[N],s[N];
float a[N]={ ... your numbers here ... };
unsigned __int32 *u=(unsigned __int32*)a,i;
for (i=0;i<N;i++)
{
m[i]=(u[i]&_f32_man)|_f32_exp_lsb;
a[i]=((u[i]&_f32_exp)>>_f32_man_bits)-_f32_exp_bias;
s[i]=(u[i]&_f32_sig);
}
-
计算 shift 其 O(N)
所以首先计算 e[i] O(N) 的最大值,然后是 shift 本身 O(1)
// shift = max(e[0...N-1])
int shift;
for (shift=e[0],i=1;i<N;i++)
if (shift<e[i])
shift=e[i];
// shift
shift=30-shift;
-
应用移位并构造结果O(N)
for (i=0;i<N;i++)
{
int d=shift+e[i]-_f32_man_bits;
if (d<0) m[i]>>=-d;
else if (d>0) m[i]<<= d;
if (s[i]) m[i]=-m[i];
}
结果在m[]。