извините, не кодер 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 числа, чтобы мсб (старший значащий бит) мантиссы с большим значением величины перешел в мсб некоторой целочисленной переменной, которую вы преобразовали a,b в целые числа безизменение их отношения (если только некоторые биты мантиссы не обрезаны из-за меньшей ширины в битах целевого типа данных переменной).Это как умножение чисел на одну и ту же константу.
извлечение показателей из a,b
, которые можно просто сделать путем непосредственного извлечения экспонентыбиты как целое число и вычитание из него, чтобы сделать его подписанным или с помощью математической функции log2().
вычислить shift
нам нужно сдвинуть мантиссы на a,b на shift бит или умножить a,b на 2^shift, чтобы большее значение величины было наибольшим, и оно все еще вписывается в целочисленную переменную.Таким образом, если я предположу, что 32 битовое целое число со знаком, мы хотим, чтобы мсб большего значения было равно 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 в целое число, как описано в предыдущем пункте.После этого вы можете разделить целые числа Resultign на их GCD или сдвинуть их вправо до тех пор, пока lsb a или b не станет ненулевым (удалить конечные нули).
Здесь небольшой пример в двоичном:
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 !!! Чтобы решить, что вам нужно добавить округление, но для этого вам нужно знать порог, который сообщает вам, какую часть числа округлить ... Не зная целевого диапазона значений с плавающей запятой или точности, я не могу реализоватьэто безопасно ...
Если вы хотите сохранить соотношение между большим числом операций с плавающей точкой, чем просто добавить их в функцию.
Вот пример плавающего C ++ для 3 переменных:
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[].