Сила в квадрате для отрицательных показателей

Question

Сила в квадрате для отрицательных показателей

Я не уверен, что сила путем возведения в квадрат учитывает отрицательный показатель. Я реализовал следующий код, который работает только для положительных чисел.

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

Просмотр https://en.wikipedia.org/wiki/Exponentiation_by_squaring не помогает, так как следующий код кажется неправильным.

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

Изменить: Благодаря этому это решение работает как для отрицательных, так и для положительных чисел:

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

Для дробного показателя мы можем сделать ниже (метод Спектре):

Предположим, что у вас есть x^0.5, тогда вы легко вычисляете квадратный корень с помощью этого метода: начните с 0 до x/2 и продолжайте проверять, что x^2 равно результату или нет в методе двоичного поиска.
Так что, если у вас есть х ^(1/3), вы должны заменить if mid*mid <= n в if mid*mid*mid <= n и вы получите кубический корень из x. То же самое относится к x^(1/4), x^(1/5) и так далее. В случае x^(2/5) мы можем сделать (x^(1/5))^2 и снова уменьшить проблему нахождения 5-го корня из x.
Однако к этому времени вы бы поняли, что этот метод работает только в том случае, когда вы можете конвертировать корень в формат 1/x. Так застряли ли мы, если мы не можем конвертировать? Нет, мы все еще можем идти вперед, поскольку у нас есть воля.
Преобразуйте вашу плавающую точку в фиксированную точку и затем вычислите pow (a, b). Предположим, что число равно 0,6, преобразовав его в (24, 8) с плавающей запятой и получит Floor(0.6*1<<8) = 153(10011001). Как вы знаете, 153 представляет дробную часть, поэтому в фиксированной точке это (10011001) представляет (2^-1, 0, 0, 2^-3, 2^-4, 0, 0, 2^7). Так что мы можем снова вычислите pow(a, 0.6), вычисляя корни 2,3, 4 и 7 в x в фиксированной точке. После вычисления нам снова нужно получить результат с плавающей точкой, разделив на 1<<8.

Код для вышеуказанного метода можно найти в принятом ответе.

Существует также метод на основе журнала:

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

6

c algorithm recursion math

Источник

user3903331 20 июн '15 в 21:16

1 ответ

Решение

Другие вопросы по тегам c algorithm recursion math

user2521214 21 июн '15 в 08:11 2015-06-21 08:11 · Accepted Answer · 2015-06-21 08:11

Целочисленные примеры для 32 бит int арифметика, DWORD 32-битный unsigned int

плавающий pow(x,y)=x^y
Обычно оценивается так:
- Как на самом деле работает Math.Pow (и так далее)
таким образом, дробный показатель можно оценить: pow(x,y) = exp2(y*log2(x)), Это можно сделать также в фиксированной точке:
- bignum pow с фиксированной точкой

целое число pow(a,b)=a^b где a>=0 , b>=0

Это легко (у вас это уже есть) сделать путем возведения в квадрат:

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

целое число pow(a,b)=a^b где b>=0

Просто добавь несколько if s для обработки негатива 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;
     }

целое число pow(a,b)=a^b
Так что если b<0 тогда это значит 1/powiu(a,-b) Как видите, результат не является целым числом, поэтому либо проигнорируйте этот случай, либо верните плавающее значение, либо добавьте переменную множителя (чтобы вы могли оценить PI уравнения на чистой целочисленной арифметике). Это результат с плавающей точкой:
```
    float powfii(int a,int b)
     {
     if (b<0) return 1.0/float(powiu(a,-b));
     else return powiu(a,b);
     }
```
целое число pow(a,b)=a^b где b дробный
Вы можете сделать что-то вроде этого a^(1/bb) где bb целое число На самом деле это рутинг, поэтому вы можете использовать бинарный поиск для оценки:
- a^(1/2) является square root(a)
- a^(1/bb) является bb_root(a)
так что бинарный поиск c от MSB до LSB и оценить, если pow(c,bb)<=a затем оставьте bit как это еще понятно. Это 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;
        }
```
так что теперь просто измените if (a*a>x) с if (pow(a,bb)>x) где bb=1/b... так b дробный показатель, который вы ищете, и bb целое число Также m количество бит результата, поэтому изменить m=(bits(x)>>1); в m=(bits(x)/bb);

[edit1] Пример с фиксированной точкой 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;
    }
//---------------------------------------------------------------------------

так что это беззнаковая фиксированная точка. Высоко 16 биты целые и младшие 16 биты являются дробной частью.

это конверсия fp -> fx: DWORD(float(x)*float(_fx32_one))
это fp <- fx преобразование: float(DWORD(x))/float(_fx32_one))
fx32_mul(x,y) является x*y он использует ассемблер 80386+ 32-битной архитектуры (вы можете переписать его на karatsuba или любой другой, чтобы быть независимым от платформы)
fx32_sqrt(x) является sqrt(x)
В фиксированной точке вы должны знать о дробном сдвиге битов для умножения: (a<<16)*(b<<16)=(a*b<<32) вам нужно вернуться назад >>16 чтобы получить результат (a*b<<16), Также результат может переполниться 32 немного поэтому я использую 64 Немного результата в сборке.

[edit2] 32-битный со знаком пример с фиксированной точкой Pow C++

Когда вы объедините все предыдущие шаги, у вас должно получиться что-то вроде этого:

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

Я проверил это так:

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 являются плавающей точкой
x,y являются ближайшими представлениями с фиксированной точкой a,b
c0 это математический результат
c1 такое результат fx32_pow
d разница
Надежда не забыла чего-то тривиального, но, похоже, все работает правильно. Не забывайте, что фиксированная точка имеет очень ограниченную точность, поэтому результаты будут немного отличаться...

PS Посмотрите на это: