Poder por quadratura para expoentes negativos

Não tenho certeza se o poder por quadratura cuida do expoente negativo. Eu implementei o código a seguir, que funciona apenas 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); } 

Olhando para https://en.wikipedia.org/wiki/Exponentiation_by_squaring não ajuda, pois o código a seguir parece errado.

  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). 

Edit: Graças a amit esta solução funciona para números negativos e 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 o expoente fracionário , podemos fazer abaixo (método Spektre):

  1. Suponha que você tenha x ^ 0,5, então você calcula facilmente a raiz quadrada por este método: comece de 0 a x / 2 e continue verificando x ^ 2 é igual ao resultado ou não no método de busca binária .

  2. Então, caso você tenha x ^ (1/3) você tem que replace if mid*mid <= n para if mid*mid*mid <= n e você irá obter a raiz cúbica de x.Same se aplica para x ^ (1/4), x ^ (1/5) e assim por diante. No caso de x ^ (2/5), podemos fazer (x ^ (1/5)) ^ 2 e novamente reduzir o problema de encontrar a quinta raiz de x.

  3. No entanto, a essa altura, você teria percebido que esse método funciona apenas quando você pode converter a raiz para o formato 1 / x. Então estamos presos se não podemos converter? Não, ainda podemos ir em frente, pois temos a vontade.

  4. Converta seu ponto flutuante em ponto fixo e depois calcule pow (a, b). Suponha que o número seja 0,6, convertendo isso para (24, 8) rendimentos de ponto flutuante Floor (0,6 * 1 << 8) = 153 (10011001). Como você sabe, 153 representa a parte fracionária de modo que em ponto fixo, isto (10011001) representa (2 ^ -1, 0, 0, 2 ^ -3, 2 ^ -4, 0, 0, 2 ^ 7). calcule o pow (a, 0.6) calculando as raízes 2,3, 4 e 7 de x em ponto fixo. Após o cálculo, novamente precisamos obter o resultado em ponto flutuante, dividindo com 1 << 8.

O código para o método acima pode ser encontrado na resposta aceita.

Existe também um método baseado em log :

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

Exemplos inteiros são para aritmética int 32 bits, DWORD é 32bit unsigned int

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

    Geralmente é avaliado assim:

    • Como o Math.Pow (e assim por diante) funciona

    então o expoente fracionário pode ser avaliado: pow(x,y) = exp2(y*log2(x)) . Isso também pode ser feito em um ponto fixo:

    • ponto fixo bignum pow
  2. inteiro pow(a,b)=a^b onde a>=0 , b>=0

    Isso é fácil (você já tem isso) feito por quadratura:

      DWORD powuu(DWORD a,DWORD b) { int i,bits=32; DWORD d=1; for (i=0;i 
  3. inteiro pow(a,b)=a^b onde b>=0

    Basta adicionar alguns if s para lidar com o negativo

      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. inteiro pow(a,b)=a^b

    Então, se b<0 , significa 1/powiu(a,-b) Como você pode ver, o resultado não é inteiro, então ignore este caso ou retorne um valor flutuante ou adicione uma variável multiplicadora (assim você pode avaliar as equações PI em Aritmética Integral Pura). Este é o resultado float:

      float powfii(int a,int b) { if (b<0) return 1.0/float(powiu(a,-b)); else return powiu(a,b); } 
  5. inteiro pow(a,b)=a^b onde b é fracional

    Você pode fazer algo assim a^(1/bb) onde bb é inteiro. Na realidade, isso é o enraizamento, então você pode usar a pesquisa binária para avaliar:

    • a^(1/2) é square root(a)
    • a^(1/bb) é bb_root(a)

    então faça uma busca binária por c de MSB para LSB e avalie se pow(c,bb)<=a então deixe o bit como está claro. Este é o exemplo 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; } 

    então agora apenas mude o if (a*a>x) com if (pow(a,bb)>x) onde bb=1/b ... então b é o expoente fracional que você está procurando e bb é inteiro. Também m é o número de bits do resultado, então mude m=(bits(x)>>1); para m=(bits(x)/bb);

[edit1] exemplo de sqrt de ponto fixo

 //--------------------------------------------------------------------------- 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; } //--------------------------------------------------------------------------- 

então este é o ponto fixo sem sinal. Altos 16 bits são inteiros e baixos 16 bits são parte fracionária.

  • isso é fp -> conversão fx: DWORD(float(x)*float(_fx32_one))
  • isto é fp <- conversão fx: float(DWORD(x))/float(_fx32_one))
  • fx32_mul(x,y) é x*y ele usa o assembler de arquitetura 80386+ 32bit (você pode reescrevê-lo para karatsuba ou qualquer outra coisa que seja independente de plataforma)
  • fx32_sqrt(x) é sqrt(x)

    No ponto fixo, você deve estar ciente do deslocamento de bits fracionários para multiplicação: (a<<16)*(b<<16)=(a*b<<32) você precisa voltar para >>16 para obter o resultado (a*b<<16) . Além disso, o resultado pode ultrapassar 32 bits, portanto, eu uso o resultado de 64 bits na assembly.

[editar2] 32bit assinado ponto fixo pow C ++ example

Quando você coloca todos os passos anteriores juntos, você deve ter algo assim:

 //--------------------------------------------------------------------------- //--- 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; } //--------------------------------------------------------------------------- 

Eu testei assim:

 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 são ponto flutuante
  • x,y são representações de ponto fixo mais próximas de a,b
  • c0 é resultado de matemática
  • c1 é o resultado fx32_pow
  • d é diferença

    a esperança não esqueceu algo trivial, mas parece que funciona corretamente. Não se esqueça que o ponto fixo tem uma precisão muito limitada, pelo que os resultados diferem um pouco ...

PS Dê uma olhada nisso:

  • Como obter uma raiz quadrada para input de 32 bits em apenas um ciclo de clock?
  • ponto fixo log2, exp2
  • Integer C ++ pow (x, 1 / y) implementação