Construire une fonction logarithmique en C sans utiliser le type float

Je dois réécrire la fonction de journalisation (base 2 ou base importe peu) sans utiliser de type float , mais je dois obtenir la précision de quelques chiffres après la virgule décimale. (comme un float * 100 pour obtenir 2 décimales dans un type entier, par exemple: si le résultat est 1.4352 , ma fonction devrait renvoyer quelque chose comme 143 (type int ) et je sais que les 2 derniers chiffres sont les décimales.

J’ai trouvé sur le stackoverflow des méthodes comme:

  • Comment puis-je calculer un logarithme en base 2 sans utiliser les fonctions mathématiques intégrées en C #?

mais tous renvoient une précision int (en évitant les décimales).

Je ne sais pas comment aborder cela, alors la question est:

Comment encoder (et / ou changer) l’implémentation du log entier pour supporter le résultat décimal?

Pour cela, vous devez utiliser une précision de point fixe / arithmétique / maths. Cela signifie que vous utilisez des variables de type entier, mais que certains bits se trouvent après le point décimal.

Par exemple, supposons 8 bits décimaux pour que les opérations se fassent comme ceci:

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

Voici un exemple simple de log2 point fixe via une recherche binary en 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<>=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; } //--------------------------------------------------------------------------- 

Voici un test simple (en utilisant des flottants uniquement pour le chargement et l’impression, vous pouvez également gérer les kiosques sur des entiers ou avec des constantes évaluées par le compilateur):

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

Cela évalue: log2(125.67) = 6.98828125 mon gain de 6.97349648 qui est assez proche. Résultat plus précis, vous avez besoin de plus de bits fractionnaires que vous devez utiliser. Int et comstack time evaluation float example:

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

renvoie 698 ce qui signifie 6.98 , multiplié par 100 . Vous pouvez également écrire vos propres fonctions de chargement et d’impression pour convertir directement entre point fixe et chaîne .

Pour changer de précision, il suffit de jouer avec la constante _fx32_fract_bits . Quoi qu'il en soit, si votre C ++ ne connaît pas DWORD il ne s'agit que de 32 bits unsigned int . Si vous utilisez un type différent ( 16 ou 64 bits, par exemple), il vous suffit de modifier les constantes en conséquence.

Pour plus d'informations, consultez:

  • Puissance en quadrature pour les exposants négatifs

[Edit2] fx32_mul sur une arithmétique 32 bits sans base asm 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; } 

Si vous n'avez pas WORD,DWORD ajoute au début du code

 typedef unsigned __int32 DWORD; typedef unsigned __int16 WORD; 

ou ca:

 typedef uint32_t DWORD; typedef uint16_t WORD; 

[Edit3] Informations de débogage de fx32_mul

laissez call et trace / breakpoint ceci (15 bits fractionnaires):

 fx32_mul(0x00123400,0x00230056); 

Lequel est:

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

Alors:

 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 }