Poder al cuadrar para exponentes negativos

No estoy seguro si el poder al cuadrar cuida el exponente negativo. Implementé el siguiente código que funciona solo para números positivos.

#include  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); } 

Ver 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 el exponente fraccional podemos hacer a continuación (método de Spektre):

  1. Supongamos que tiene x ^ 0.5 y luego calcula fácilmente la raíz cuadrada con 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 binario .

  2. Entonces, en caso de que tengas x ^ (1/3) debes reemplazar if mid*mid <= n if mid*mid*mid <= n y obtendrás la raíz cúbica de x. Lo que 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 de nuevo reducir el problema de encontrar la 5ta raíz de x.

  3. Sin embargo, en este momento se habrá dado cuenta de que este método solo funciona en el caso en que puede convertir la raíz a formato 1 / x. Entonces, ¿estamos atascados si no podemos convertir? No, todavía podemos seguir adelante ya que tenemos la voluntad.

  4. Convierta su punto flotante a punto fijo y luego calcule pow (a, b). Suponga que el número es 0.6, convirtiéndolo en (24, 8) rendimiento en coma flotante Piso (0.6 * 1 << 8) = 153 (10011001). Como sabes, 153 representa la parte fraccional, así que en el punto fijo esto (10011001) representa (2 ^ -1, 0, 0, 2 ^ -3, 2 ^ -4, 0, 0, 2 ^ 7). Así que podemos volver calcule el pow (a, 0.6) calculando la raíz 2,3, 4 y 7 de x en el punto fijo. Después de calcular, de nuevo necesitamos obtener el resultado en coma flotante dividiendo 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 registro :

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

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

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

    Generalmente se evalúa así:

    • Cómo funciona Math.Pow (etc.)

    entonces el exponente fraccionario puede ser evaluado: pow(x,y) = exp2(y*log2(x)) . Esto también se puede hacer en punto fijo:

    • punto fijo bignum pow
  2. 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 
  3. entero pow(a,b)=a^b donde b>=0

    Solo 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; } 
  4. entero pow(a,b)=a^b

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

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

    Puedes hacer algo como esto a^(1/bb) donde bb es un número entero. En realidad, esto está enraizando para que pueda 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 de c de MSB a LSB y evalúe si pow(c,bb)<=a luego abandone el bit como está más claro. Este es el ejemplo 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<>=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 ... así que b es el exponente fraccional que busca y bb es 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<>=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 firmar. Los 16 bits altos son enteros y los 16 bits bajos son partes fraccionarias.

  • esto es conversión fp -> fx: DWORD(float(x)*float(_fx32_one))
  • esto es conversión fp <- fx: float(DWORD(x))/float(_fx32_one))
  • fx32_mul(x,y) es x*y utiliza el ensamblador de la architecture 80386+ de 32 bits (puede volver a escribirlo en karatsuba o en 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 bit fraccional para la multiplicación: (a<<16)*(b<<16)=(a*b<<32) necesita retroceder en >>16 para obtener el resultado (a*b<<16) . Además, el resultado puede desbordar 32 bits, por lo tanto, utilizo 64 bits para el ensamblaje.

[edit2] Ejemplo de C ++ con punto fijo de 32 bits con signo fijo

Cuando pongas todos los pasos anteriores juntos 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<>=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 puntos flotantes
  • x,y son las representaciones más cercanas de punto fijo de a,b
  • c0 es resultado matemático pow
  • c1 es resultado fx32_pow
  • d es la diferencia

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

PD Mire esto:

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