Potencia al cuadrado para exponentes negativos
No estoy seguro de si la potencia 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 el exponente fraccionario podemos hacer lo siguiente (método Spektre):
Supongamos que tiene x^0.5, entonces puede calcular fácilmente la raíz cuadrada con este método: comience desde 0 hasta x/2 y siga comprobando que x^2 sea igual al resultado o no en el método de búsqueda binaria .
Entonces, en caso de que tenga x^(1/3), debe reemplazarlo
if mid*mid <= nyif mid*mid*mid <= nobtendrá la raíz cúbica de x. Lo mismo se aplica para x^(1/4), x^(1/5), etc. . 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, a estas alturas ya te habrás dado cuenta de que este método sólo funciona en el caso de que puedas convertir la raíz al formato 1/x. Entonces, ¿estamos estancados si no podemos convertirnos? No, todavía podemos seguir adelante porque tenemos la voluntad.
Convierta su punto flotante en punto fijo y luego calcule pow(a, b). Supongamos que el número es 0,6, al convertirlo a (24, 8) punto flotante se obtiene Floor(0,6*1<<8) = 153(10011001). Como sabes, 153 representa la parte fraccionaria, por lo que en punto fijo esto (10011001) representa (2^-1, 0, 0, 2^-3, 2^-4, 0, 0, 2^7). Así que podemos nuevamente 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 punto flotante dividiendo por 1<<8.
El código para el método anterior se puede encontrar en la respuesta aceptada.
También existe un método basado en registros :
x^y = exp2(y*log2(x))
Los ejemplos de números enteros son para intaritmética de 32 bits, DWORDson 32 bitsunsigned int
- flotante
pow(x,y)=x^y
Generalmente se evalúa así:
- Cómo funciona realmente Math.Pow (y demás)
entonces se puede evaluar el exponente fraccionario: pow(x,y) = exp2(y*log2(x)). Esto también se puede hacer en punto fijo:
- punto fijo bignum pow
- entero
pow(a,b)=a^bdondea>=0 , b>=0
Esto es fácil (ya lo tienes) hecho elevando al cuadrado:
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^bdondeb>=0
Solo agrega algunas ifs para manejar lo 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<0significa que 1/powiu(a,-b), como puede ver, el resultado no es un número entero en absoluto, ignore este caso o devuelva un valor flotante o agregue una variable multiplicadora (para que pueda evaluar PIecuaciones en aritmética pura de enteros). 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^bdondebes fraccionario
Puedes hacer algo como esto a^(1/bb)donde bbes un número entero. En realidad, esto es enraizamiento, por lo que puedes usar la búsqueda binaria para evaluar:
a^(1/2)essquare root(a)a^(1/bb)esbb_root(a)
así que haga una búsqueda binaria cde MSB a LSB y evalúe si pow(c,bb)<=aluego deja el bitestado como está, de lo contrario, aclare. Este es sqrtun ejemplo:
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 cambia el if (a*a>x)con if (pow(a,bb)>x)donde bb=1/b... entonces bel exponente fraccionario que buscas bbes un número entero. También mes el número de bits del resultado, así que cambie m=(bits(x)>>1);am=(bits(x)/bb);
[editar1] 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 bits altos 16son números enteros y 16los bits bajos son partes fraccionarias.
- esto es fp -> conversión fx:
DWORD(float(x)*float(_fx32_one)) - esta es la conversión fp <- fx:
float(DWORD(x))/float(_fx32_one)) fx32_mul(x,y)¿x*yUtiliza un ensamblador de arquitectura 80386+ de 32 bits (puede reescribirlo en karatsuba o cualquier otra cosa para que sea independiente de la plataforma)?fx32_sqrt(x)essqrt(x)
En punto fijo, debe tener en cuenta el desplazamiento de bits fraccionarios para la multiplicación: (a<<16)*(b<<16)=(a*b<<32)debe retroceder >>16para obtener el resultado (a*b<<16). Además, el resultado puede desbordar 32el bit, por lo que utilizo 64el resultado del bit en el ensamblaje.
[editar2] Ejemplo de C++ de potencia de punto fijo con signo de 32 bits
Cuando juntes 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,bson punto flotantex,yson las representaciones de punto fijo más cercanas dea,bc0es el resultado matemático powc1es el resultado fx32_powdes 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 Echa un vistazo a esto:
- ¿Cómo obtener una raíz cuadrada para una entrada de 32 bits en un solo ciclo de reloj?
- punto fijo log2,exp2
- implementación de entero C++ pow(x,1/y)
- implementación de entero C++ pow(a,b),log2(a),logb(a)
- poder de dominio real basado en matemáticas de dominio complejas
