simplificar - Potencia cuadrando para exponentes negativos

recursividad potencia de 2 (1)

Ejemplos enteros son para aritmética int de 32 bits, DWORD es unsigned int 32 bits

pow(x,y)=x^y flotante pow(x,y)=x^y

Generalmente se evalúa así:
- Cómo funciona Math.Pow (y así sucesivamente)
entonces el exponente fraccional se puede evaluar: pow(x,y) = exp2(y*log2(x)) . Esto se puede hacer también en un punto fijo:
- punto fijo bignum pow
entero pow(a,b)=a^b donde a>=0 , b>=0

Esto es fácil (ya lo tienes) hecho al cuadrar:

DWORD powuu(DWORD a,DWORD b) { int i,bits=32; DWORD d=1; for (i=0;i<bits;i++) { d*=d; if (DWORD(b&0x80000000)) d*=a; b<<=1; } return d; }
entero pow(a,b)=a^b donde b>=0

Simplemente agregue algunos if s para manejar el negativo a

int powiu(int a,DWORD b) { int sig=0,c; if ((a<0)&&(DWORD(b&1)) { sig=1; a=-a; } // negative output only if a<0 and b is odd c=powuu(a,b); if (sig) c=-c; return c; }
entero pow(a,b)=a^b

Entonces, si b<0 entonces significa 1/powiu(a,-b) Como puede ver, el resultado no es entero, por lo tanto, ignore este caso o devuelva un valor flotante o agregue una variable multiplicadora (para que pueda evaluar las ecuaciones PI en aritmética entera pura). Este es el resultado flotante:

float powfii(int a,int b) { if (b<0) return 1.0/float(powiu(a,-b)); else return powiu(a,b); }
entero pow(a,b)=a^b donde b es fraccional

Puede hacer algo como esto a^(1/bb) donde bb es entero. En realidad, esto es enraizamiento, por lo que puede usar la búsqueda binaria para evaluar:
- a^(1/2) es square root(a)
- a^(1/bb) es bb_root(a)
así que haga una búsqueda binaria para c de MSB a LSB y evalúe si pow(c,bb)<=a luego deje el bit como está claro. Este es un ejemplo de sqrt :

int bits(DWORD p) // count how many bits is p { DWORD m=0x80000000; int b=32; for (;m;m>>=1,b--) if (p>=m) break; return b; } DWORD sqrt(const DWORD &x) { DWORD m,a; m=(bits(x)>>1); if (m) m=1<<m; else m=1; for (a=0;m;m>>=1) { a|=m; if (a*a>x) a^=m; } return a; }
así que ahora solo cambie el if (a*a>x) con if (pow(a,bb)>x) donde bb=1/b ... entonces b es el exponente fraccionario que está buscando y bb es un número entero. También m es el número de bits del resultado, así que cambie m=(bits(x)>>1); a m=(bits(x)/bb);

[edit1] ejemplo de punto fijo sqrt

//--------------------------------------------------------------------------- const int _fx32_fract=16; // fractional bits count const int _fx32_one =1<<_fx32_fract; DWORD fx32_mul(const DWORD &x,const DWORD &y) // unsigned fixed point mul { 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; } 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) m-=_fx32_fract; 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; } //---------------------------------------------------------------------------

entonces este es un punto fijo sin signo. Los 16 bits altos son enteros y los 16 bits bajos son parte fraccionaria.

esta es la conversión de fp -> fx: DWORD(float(x)*float(_fx32_one))
esta es la conversión de fp <- fx: float(DWORD(x))/float(_fx32_one))
fx32_mul(x,y) es x*y usa ensamblador de arquitectura 80386+ 32bit (puede reescribirlo en karatsuba o cualquier otra cosa para ser independiente de la plataforma)
fx32_sqrt(x) es sqrt(x)

En el punto fijo, debe tener en cuenta el desplazamiento de bits fraccional para la multiplicación: (a<<16)*(b<<16)=(a*b<<32) necesita retroceder >>16 para obtener el resultado (a*b<<16) . Además, el resultado puede desbordar 32 bits, por lo tanto, utilizo el resultado de 64 bits en el ensamblaje.

[edit2] 32bit firmado punto fijo pow C ++ ejemplo

Cuando juntas todos los pasos anteriores, deberías tener algo como esto:

//--------------------------------------------------------------------------- //--- 32bit signed fixed point format (2os complement) //--------------------------------------------------------------------------- // |MSB LSB| // |integer|.|fractional| //--------------------------------------------------------------------------- const int _fx32_bits=32; // all bits count const int _fx32_fract_bits=16; // fractional bits count const int _fx32_integ_bits=_fx32_bits-_fx32_fract_bits; // integer bits count //--------------------------------------------------------------------------- const int _fx32_one =1<<_fx32_fract_bits; // constant=1.0 (fixed point) const float _fx32_onef =_fx32_one; // constant=1.0 (floating point) const int _fx32_fract_mask=_fx32_one-1; // fractional bits mask const int _fx32_integ_mask=0xFFFFFFFF-_fx32_fract_mask; // integer bits mask const int _fx32_sMSB_mask =1<<(_fx32_bits-1); // max signed bit mask const int _fx32_uMSB_mask =1<<(_fx32_bits-2); // max unsigned bit mask //--------------------------------------------------------------------------- float fx32_get(int x) { return float(x)/_fx32_onef; } int fx32_set(float x) { return int(float(x*_fx32_onef)); } //--------------------------------------------------------------------------- int fx32_mul(const int &x,const int &y) // x*y { int a=x,b=y; // asm has access only to local variables asm { // compute (a*b)>>_fx32_fract mov eax,a mov ebx,b mul eax,ebx // (edx,eax)=a*b mov ebx,_fx32_one div ebx // eax=(a*b)>>_fx32_fract mov a,eax; } return a; } //--------------------------------------------------------------------------- int fx32_div(const int &x,const int &y) // x/y { int a=x,b=y; // asm has access only to local variables asm { // compute (a*b)>>_fx32_fract mov eax,a mov ebx,_fx32_one mul eax,ebx // (edx,eax)=a<<_fx32_fract mov ebx,b div ebx // eax=(a<<_fx32_fract)/b mov a,eax; } return a; } //--------------------------------------------------------------------------- int fx32_abs_sqrt(int x) // |x|^(0.5) { int m,a; if (!x) return 0; if (x<0) x=-x; m=bits(x); // integer bits for (a=x,m=0;a;a>>=1,m++); // count all bits m-=_fx32_fract_bits; // compute result integer bits (half of x integer bits) if (m<0) m=0; m>>=1; 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; } //--------------------------------------------------------------------------- int fx32_pow(int x,int y) // x^y { // handle special cases if (!y) return _fx32_one; // x^0 = 1 if (!x) return 0; // 0^y = 0 if y!=0 if (y==-_fx32_one) return fx32_div(_fx32_one,x); // x^-1 = 1/x if (y==+_fx32_one) return x; // x^+1 = x int m,a,b,_y; int sx,sy; // handle the signs sx=0; if (x<0) { sx=1; x=-x; } sy=0; if (y<0) { sy=1; y=-y; } _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_uMSB_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 (int(y&m)) a=fx32_mul(a,x); } } // powering by rooting x^_y if (_y) { for (b=x,m=_fx32_one>>1;m;m>>=1) // use only fractional part { b=fx32_abs_sqrt(b); if (int(_y&m)) a=fx32_mul(a,b); } } // handle signs if (sy) { if (a) a=fx32_div(_fx32_one,a); else a=0; /*Error*/ } // underflow if (sx) { if (_y) a=0; /*Error*/ else if(int(y&_fx32_one)) a=-a; } // negative number ^ non integer exponent, here could add test if 1/_y is integer instead return a; } //---------------------------------------------------------------------------

Lo he probado así:

float a,b,c0,c1,d; int x,y; for (a=0.0,x=fx32_set(a);a<=10.0;a+=0.1,x=fx32_set(a)) for (b=-2.5,y=fx32_set(b);b<=2.5;b+=0.1,y=fx32_set(b)) { if (!x) continue; // math pow has problems with this if (!y) continue; // math pow has problems with this c0=pow(a,b); c1=fx32_get(fx32_pow(x,y)); d=0.0; if (fabs(c1)<1e-3) d=c1-c0; else d=(c0/c1)-1.0; if (fabs(d)>0.1) d=d; // here add breakpoint to check inconsistencies with math pow }

a,b son coma flotante
x,y son las representaciones de punto fijo más cercanas de a,b
c0 es el resultado matemático
c1 es el resultado de fx32_pow
d es la diferencia

La esperanza no olvidó algo trivial, pero parece que funciona correctamente. No olvide que el punto fijo tiene una precisión muy limitada, por lo que los resultados diferirán un poco ...

PD Mira esto:

¿Cómo obtener una raíz cuadrada para una entrada de 32 bits solo en un ciclo de reloj?
punto fijo log2, exp2
implementación entera de C ++ pow (x, 1 / y)

No estoy seguro de si el poder al cuadrado se ocupa del exponente negativo. Implementé el siguiente código que funciona solo para números positivos.

#include <stdio.h> int powe(int x, int exp) { if (x == 0) return 1; if (x == 1) return x; if (x&1) return powe(x*x, exp/2); else return x*powe(x*x, (exp-1)/2); }

Mirar https://en.wikipedia.org/wiki/Exponentiation_by_squaring no ayuda, ya que el siguiente código parece incorrecto.

Function exp-by-squaring(x, n ) if n < 0 then return exp-by-squaring(1 / x, - n ); else if n = 0 then return 1; else if n = 1 then return x ; else if n is even then return exp-by-squaring(x * x, n / 2); else if n is odd then return x * exp-by-squaring(x * x, (n - 1) / 2).

Editar: Gracias a amit, esta solución funciona tanto para números negativos como positivos:

float powe(float x, int exp) { if (exp < 0) return powe(1/x, -exp); if (exp == 0) return 1; if (exp == 1) return x; if (((int)exp)%2==0) return powe(x*x, exp/2); else return x*powe(x*x, (exp-1)/2); }

Para exponente fraccionario podemos hacer a continuación (método Spektre):

Suponga que tiene x ^ 0.5 y luego calcula fácilmente la raíz cuadrada mediante este método: comience de 0 a x / 2 y siga comprobando que x ^ 2 es igual al resultado o no en el método de búsqueda binaria .
Entonces, en caso de que tenga x ^ (1/3), debe reemplazar if mid*mid <= n a if mid*mid*mid <= n y obtendrá la raíz cúbica de x. Lo mismo se aplica para x ^ (1/4), x ^ (1/5) y así sucesivamente. En el caso de x ^ (2/5) podemos hacer (x ^ (1/5)) ^ 2 y nuevamente reducir el problema de encontrar la quinta raíz de x.
Sin embargo, en este momento se habrá dado cuenta de que este método funciona solo en el caso en que pueda convertir la raíz al formato 1 / x. Entonces, ¿estamos atrapados si no podemos convertir? No, todavía podemos seguir adelante, ya que tenemos la voluntad.
Convierta su punto flotante en punto fijo y luego calcule pow (a, b). Suponga que el número es 0.6, convirtiéndolo a (24, 8) el punto flotante produce Piso (0.6 * 1 << 8) = 153 (10011001). Como saben 153 representa la parte fraccionaria, así que en punto fijo esto (10011001) representa (2 ^ -1, 0, 0, 2 ^ -3, 2 ^ -4, 0, 0, 2 ^ 7). calcule el pow (a, 0.6) calculando la raíz 2,3, 4 y 7 de x en punto fijo. Después de calcular nuevamente necesitamos obtener el resultado en coma flotante dividiéndolo con 1 << 8.

El código para el método anterior se puede encontrar en la respuesta aceptada.

También hay un método basado en el registro :

x ^ y = exp2 (y * log2 (x))