Наиболее эффективный способ реализации целочисленной степенной функции pow(int, int)
Каков наиболее эффективный способ поднять целое число до степени другого целого числа в C?
// 2^3
pow(2,3) == 8
// 5^5
pow(5,5) == 3125
22 ответа
Возведение в степень путем возведения в квадрат.
int ipow(int base, int exp)
{
int result = 1;
for (;;)
{
if (exp & 1)
result *= base;
exp >>= 1;
if (!exp)
break;
base *= base;
}
return result;
}
Это стандартный метод модульного возведения в степень для огромных чисел в асимметричной криптографии.
Обратите внимание, что возведение в квадрат по квадратам - не самый оптимальный метод. Вероятно, это лучший способ, который вы можете использовать в качестве общего метода, который работает для всех значений показателя, но для конкретного значения показателя может быть лучшая последовательность, которая требует меньше умножений.
Например, если вы хотите вычислить x^15, метод возведения в степень путем возведения в квадрат даст вам:
x^15 = (x^7)*(x^7)*x
x^7 = (x^3)*(x^3)*x
x^3 = x*x*x
Это всего 6 умножений.
Оказывается, это можно сделать, используя "всего лишь" 5 умножений посредством возведения в степень сложения.
n*n = n^2
n^2*n = n^3
n^3*n^3 = n^6
n^6*n^6 = n^12
n^12*n^3 = n^15
Не существует эффективных алгоритмов для нахождения этой оптимальной последовательности умножений. Из Википедии:
Проблема нахождения кратчайшей цепочки сложений не может быть решена динамическим программированием, поскольку она не удовлетворяет предположению об оптимальной подструктуре. То есть недостаточно разложить мощность на меньшие степени, каждая из которых вычисляется минимально, поскольку цепочки сложения для меньших степеней могут быть связаны (для совместного использования вычислений). Например, в кратчайшей цепочке сложения для a¹⁵, указанной выше, подзадача для a⁶ должна быть вычислена как (a³)², поскольку a³ используется повторно (в отличие, скажем, от a⁶ = a²(a²)², что также требует трех умножений).
Если вам нужно поднять 2 до степени. Самый быстрый способ сделать это - сдвинуться с места силой.
2 ** 3 == 1 << 3 == 8
2 ** 30 == 1 << 30 == 1073741824 (A Gigabyte)
Вот метод в Java
private int ipow(int base, int exp)
{
int result = 1;
while (exp != 0)
{
if ((exp & 1) == 1)
result *= base;
exp >>= 1;
base *= base;
}
return result;
}
int pow( int base, int exponent)
{ // Does not work for negative exponents. (But that would be leaving the range of int)
if (exponent == 0) return 1; // base case;
int temp = pow(base, exponent/2);
if (exponent % 2 == 0)
return temp * temp;
else
return (base * temp * temp);
}
power() функция для работы только для целых
int power(int base, unsigned int exp){
if (exp == 0)
return 1;
int temp = power(base, exp/2);
if (exp%2 == 0)
return temp*temp;
else
return base*temp*temp;
}
Сложность = O(log(exp))
power() функция работать для отрицательного опыта и плавающей базы.
float power(float base, int exp) {
if( exp == 0)
return 1;
float temp = power(base, exp/2);
if (exp%2 == 0)
return temp*temp;
else {
if(exp > 0)
return base*temp*temp;
else
return (temp*temp)/base; //negative exponent computation
}
}
Сложность = O(log(exp))
Чрезвычайно специализированный случай, когда вам нужно сказать 2^(-x к y), где x, конечно, отрицателен, а y слишком велик, чтобы делать смещение на int. Вы все еще можете делать 2 ^ х в постоянное время, прикручивая поплавком.
struct IeeeFloat
{
unsigned int base : 23;
unsigned int exponent : 8;
unsigned int signBit : 1;
};
union IeeeFloatUnion
{
IeeeFloat brokenOut;
float f;
};
inline float twoToThe(char exponent)
{
// notice how the range checking is already done on the exponent var
static IeeeFloatUnion u;
u.f = 2.0;
// Change the exponent part of the float
u.brokenOut.exponent += (exponent - 1);
return (u.f);
}
Вы можете получить больше степеней 2, используя удвоение в качестве базового типа. (Большое спасибо комментаторам за помощь в выравнивании этого поста).
Также существует вероятность того, что, узнав больше о плавающих элементах IEEE, могут проявиться и другие особые случаи возведения в степень.
Если вы хотите получить значение целого числа для 2, возведенное в степень чего-либо, всегда лучше использовать параметр shift:
pow(2,5) можно заменить на 1<<5
Это гораздо эффективнее.
Так же, как продолжение комментариев об эффективности возведения в степень путем возведения в квадрат.
Преимущество такого подхода заключается в том, что он выполняется за время log(n). Например, если вы собирались вычислить что-то огромное, например, x^1048575 (2^20 - 1), вам нужно всего лишь пройти цикл 20 раз, а не 1 миллион +, используя наивный подход.
Кроме того, с точки зрения сложности кода это проще, чем пытаться найти наиболее оптимальную последовательность умножений, как предложено Прамодом.
Редактировать:
Я думаю, мне следует уточнить, прежде чем кто-то пометит меня на предмет потенциального переполнения. Этот подход предполагает, что у вас есть какая-то огромная библиотека.
Поздно на вечеринку:
Ниже приведено решение, которое также касается y < 0 как можно лучше.
- Он использует результат
intmax_tдля максимальной дальности. Там нет положения для ответов, которые не вписываются вintmax_t, powjii(0, 0) --> 1что является общим результатом для этого случая.pow(0,negative)другой неопределенный результат, возвращаетINTMAX_MAXintmax_t powjii(int x, int y) { if (y < 0) { switch (x) { case 0: return INTMAX_MAX; case 1: return 1; case -1: return y % 2 ? -1 : 1; } return 0; } intmax_t z = 1; intmax_t base = x; for (;;) { if (y % 2) { z *= base; } y /= 2; if (y == 0) { break; } base *= base; } return z; }
Этот код использует цикл навсегда for(;;) чтобы избежать финала base *= base распространено в других зацикленных решениях. Это умножение 1) не нужно и 2) может быть int*int переполнение, которое является UB.
int pow(int const x, unsigned const e) noexcept
{
return !e ? 1 : 1 == e ? x : pow(x * x, e / 2) * (e % 2 ? x : 1);
}
Да, это рекурсивно, но хороший оптимизирующий компилятор оптимизирует рекурсию.
В дополнение к ответу Элиаса, который приводит к неопределенному поведению при реализации с целыми числами со знаком и к неправильным значениям для высокого ввода при реализации с целыми числами без знака,
Вот модифицированная версия Exponentiation by Squaring, которая также работает со знаковыми целочисленными типами и не дает неправильных значений:
#include <stdint.h>
#define SQRT_INT64_MAX (INT64_C(0xB504F333))
int64_t alx_pow_s64 (int64_t base, uint8_t exp)
{
int_fast64_t base_;
int_fast64_t result;
base_ = base;
if (base_ == 1)
return 1;
if (!exp)
return 1;
if (!base_)
return 0;
result = 1;
if (exp & 1)
result *= base_;
exp >>= 1;
while (exp) {
if (base_ > SQRT_INT64_MAX)
return 0;
base_ *= base_;
if (exp & 1)
result *= base_;
exp >>= 1;
}
return result;
}
Соображения для этой функции:
(1 ** N) == 1
(N ** 0) == 1
(0 ** 0) == 1
(0 ** N) == 0
Если произойдет какое-либо переполнение или упаковка, return 0;
я использовал int64_t, но любая ширина (со знаком или без знака) может быть использована с небольшими изменениями. Однако если вам нужно использовать целочисленный тип без фиксированной ширины, вам нужно будет изменить SQRT_INT64_MAX от (int)sqrt(INT_MAX) (в случае использования int) или что-то подобное, что должно быть оптимизировано, но это уродливее, а не константное выражение C. Также кастинг результат sqrt() для int не очень хорошо из-за точности с плавающей запятой в случае идеального квадрата, но, как я не знаю ни одной реализации, где INT_MAX - или максимум любого типа - это идеальный квадрат, вы можете жить с этим.
Еще одна реализация (на Java). Может быть не самое эффективное решение, но число итераций такое же, как у экспоненциального решения.
public static long pow(long base, long exp){
if(exp ==0){
return 1;
}
if(exp ==1){
return base;
}
if(exp % 2 == 0){
long half = pow(base, exp/2);
return half * half;
}else{
long half = pow(base, (exp -1)/2);
return base * half * half;
}
}
Более общее решение с учетом отрицательной степени
private static int pow(int base, int exponent) {
int result = 1;
if (exponent == 0)
return result; // base case;
if (exponent < 0)
return 1 / pow(base, -exponent);
int temp = pow(base, exponent / 2);
if (exponent % 2 == 0)
return temp * temp;
else
return (base * temp * temp);
}
Решение O (log N) в Swift ...
// Time complexity is O(log N)
func power(_ base: Int, _ exp: Int) -> Int {
// 1. If the exponent is 1 then return the number (e.g a^1 == a)
//Time complexity O(1)
if exp == 1 {
return base
}
// 2. Calculate the value of the number raised to half of the exponent. This will be used to calculate the final answer by squaring the result (e.g a^2n == (a^n)^2 == a^n * a^n). The idea is that we can do half the amount of work by obtaining a^n and multiplying the result by itself to get a^2n
//Time complexity O(log N)
let tempVal = power(base, exp/2)
// 3. If the exponent was odd then decompose the result in such a way that it allows you to divide the exponent in two (e.g. a^(2n+1) == a^1 * a^2n == a^1 * a^n * a^n). If the eponent is even then the result must be the base raised to half the exponent squared (e.g. a^2n == a^n * a^n = (a^n)^2).
//Time complexity O(1)
return (exp % 2 == 1 ? base : 1) * tempVal * tempVal
}
Я реализовал алгоритм, который запоминает все вычисленные мощности, а затем использует их при необходимости. Так, например, x^13 равно (x^2)^2^2 * x^2^2 * x, где x ^ 2 ^ 2 это взято из таблицы вместо того, чтобы вычислять это снова. Необходимое количество умножения: Ceil(Log n)
public static int Power(int base, int exp)
{
int tab[] = new int[exp + 1];
tab[0] = 1;
tab[1] = base;
return Power(base, exp, tab);
}
public static int Power(int base, int exp, int tab[])
{
if(exp == 0) return 1;
if(exp == 1) return base;
int i = 1;
while(i < exp/2)
{
if(tab[2 * i] <= 0)
tab[2 * i] = tab[i] * tab[i];
i = i << 1;
}
if(exp <= i)
return tab[i];
else return tab[i] * Power(base, exp - i, tab);
}
Я использую рекурсивный, если exp четный,5^10 =25^5.
int pow(float base,float exp){
if (exp==0)return 1;
else if(exp>0&&exp%2==0){
return pow(base*base,exp/2);
}else if (exp>0&&exp%2!=0){
return base*pow(base,exp-1);
}
}
Вот алгоритм O(1) для вычисленияx ** y, вдохновленный этим комментарием . Он работает для 32-битных подписанныхint.
Для малых значений используется возведение в степень путем возведения в квадрат. Для больших значенийy, есть только несколько значенийxгде результат не переполняется. В этой реализации используется таблица поиска для чтения результата без вычислений.
При переполнении стандарт C разрешает любое поведение, включая сбой. Тем не менее, я решил выполнить граничную проверку индексов LUT, чтобы предотвратить нарушение доступа к памяти, что может быть неожиданным и нежелательным.
Псевдокод:
If `x` is between -2 and 2, use special-case formulas.
Otherwise, if `y` is between 0 and 8, use special-case formulas.
Otherwise:
Set x = abs(x); remember if x was negative
If x <= 10 and y <= 19:
Load precomputed result from a lookup table
Otherwise:
Set result to 0 (overflow)
If x was negative and y is odd, negate the result
C-код:
#define POW9(x) x * x * x * x * x * x * x * x * x
#define POW10(x) POW9(x) * x
#define POW11(x) POW10(x) * x
#define POW12(x) POW11(x) * x
#define POW13(x) POW12(x) * x
#define POW14(x) POW13(x) * x
#define POW15(x) POW14(x) * x
#define POW16(x) POW15(x) * x
#define POW17(x) POW16(x) * x
#define POW18(x) POW17(x) * x
#define POW19(x) POW18(x) * x
int mypow(int x, unsigned y)
{
static int table[8][11] = {
{POW9(3), POW10(3), POW11(3), POW12(3), POW13(3), POW14(3), POW15(3), POW16(3), POW17(3), POW18(3), POW19(3)},
{POW9(4), POW10(4), POW11(4), POW12(4), POW13(4), POW14(4), POW15(4), 0, 0, 0, 0},
{POW9(5), POW10(5), POW11(5), POW12(5), POW13(5), 0, 0, 0, 0, 0, 0},
{POW9(6), POW10(6), POW11(6), 0, 0, 0, 0, 0, 0, 0, 0},
{POW9(7), POW10(7), POW11(7), 0, 0, 0, 0, 0, 0, 0, 0},
{POW9(8), POW10(8), 0, 0, 0, 0, 0, 0, 0, 0, 0},
{POW9(9), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{POW9(10), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}
};
int is_neg;
int r;
switch (x)
{
case 0:
return y == 0 ? 1 : 0;
case 1:
return 1;
case -1:
return y % 2 == 0 ? 1 : -1;
case 2:
return 1 << y;
case -2:
return (y % 2 == 0 ? 1 : -1) << y;
default:
switch (y)
{
case 0:
return 1;
case 1:
return x;
case 2:
return x * x;
case 3:
return x * x * x;
case 4:
r = x * x;
return r * r;
case 5:
r = x * x;
return r * r * x;
case 6:
r = x * x;
return r * r * r;
case 7:
r = x * x;
return r * r * r * x;
case 8:
r = x * x;
r = r * r;
return r * r;
default:
is_neg = x < 0;
if (is_neg)
x = -x;
if (x <= 10 && y <= 19)
r = table[x - 3][y - 9];
else
r = 0;
if (is_neg && y % 2 == 1)
r = -r;
return r;
}
}
}
Мой случай немного отличается, я пытаюсь создать маску из власти, но я решил поделиться решением, которое нашел в любом случае.
Очевидно, это работает только для степеней 2.
Mask1 = 1 << (Exponent - 1);
Mask2 = Mask1 - 1;
return Mask1 + Mask2;
Если вы знаете экспоненту (и это целое число) во время компиляции, вы можете использовать шаблоны, чтобы развернуть цикл. Это можно сделать более эффективным, но я хотел продемонстрировать здесь основной принцип:
#include <iostream>
template<unsigned long N>
unsigned long inline exp_unroll(unsigned base) {
return base * exp_unroll<N-1>(base);
}
Мы прекращаем рекурсию, используя специализацию шаблона:
template<>
unsigned long inline exp_unroll<1>(unsigned base) {
return base;
}
Экспонента должна быть известна во время выполнения,
int main(int argc, char * argv[]) {
std::cout << argv[1] <<"**5= " << exp_unroll<5>(atoi(argv[1])) << ;std::endl;
}
Я заметил кое-что странное в стандартном алгоритме экспоненциального возведения в квадрат сgnu-GMP:
Я реализовал 2 почти идентичные функции - функцию степени по модулю, использующую самый ванильный алгоритм бинарного экспоненциального возведения в квадрат,
- помечен
______2()
затем еще один, в основном та же концепция, но переназначенная на деление на 10 в каждом раунде вместо деления на 2,
- помечен
______10()
.
( time ( jot - 1456 9999999999 6671 | pvE0 |
gawk -Mbe '
function ______10(_, __, ___, ____, _____, _______) {
__ = +__
____ = (____+=_____=____^= \
(_ %=___=+___)<_)+____++^____—
while (__) {
if (_______= __%____) {
if (__==_______) {
return (_^__ *_____) %___
}
__-=_______
_____ = (_^_______*_____) %___
}
__/=____
_ = _^____%___
}
}
function ______2(_, __, ___, ____, _____) {
__=+__
____+=____=_____^=(_%=___=+___)<_
while (__) {
if (__ %____) {
if (__<____) {
return (_*_____) %___
}
_____ = (_____*_) %___
--__
}
__/=____
_= (_*_) %___
}
}
BEGIN {
OFMT = CONVFMT = "%.250g"
__ = (___=_^= FS=OFS= "=")(_<_)
_____ = __^(_=3)^--_ * ++_-(_+_)^_
______ = _^(_+_)-_ + _^!_
_______ = int(______*_____)
________ = 10 ^ 5 + 1
_________ = 8 ^ 4 * 2 - 1
}
-
GNU Awk 5.1.1, API: 3.1 (GNU MPFR 4.1.0, GNU MP 6.2.1)
.
($++NF = ______10(_=$___, NR %________ +_________,_______*(_-11))) ^!___'
out9: 48.4MiB 0:00:08 [6.02MiB/s] [6.02MiB/s] [ <=> ]
in0: 15.6MiB 0:00:08 [1.95MiB/s] [1.95MiB/s] [ <=> ]
( jot - 1456 9999999999 6671 | pvE 0.1 in0 | gawk -Mbe ; )
8.31s user 0.06s system 103% cpu 8.058 total
ffa16aa937b7beca66a173ccbf8e1e12 stdin
($++NF = ______2(_=$___, NR %________ +_________,_______*(_-11))) ^!___'
out9: 48.4MiB 0:00:12 [3.78MiB/s] [3.78MiB/s] [<=> ]
in0: 15.6MiB 0:00:12 [1.22MiB/s] [1.22MiB/s] [ <=> ]
( jot - 1456 9999999999 6671 | pvE 0.1 in0 | gawk -Mbe ; )
13.05s user 0.07s system 102% cpu 12.821 total
ffa16aa937b7beca66a173ccbf8e1e12 stdin
По причинам, крайне нелогичным и неизвестным мне, для самых разных входных данных, которые я использовал, вариант div-10 почти всегда быстрее. Это совпадение хэшей между двумя сделало это действительно сбивающим с толку, несмотря на то, что компьютеры явно не были встроены и для парадигмы с основанием 10.
Я упускаю что-то важное или очевидное в коде/подходе, что может привести к искажению результатов? Спасибо.
Игнорируя особый случай 2, возведенный в степень, наиболее эффективным будет простая итерация.
int pow(int base, int pow) {
int res = 1;
for(int i=pow; i<pow; i++)
res *= base;
return res;
}
РЕДАКТИРОВАТЬ: Как было указано, это не самый эффективный способ... при условии, что вы определяете эффективность как циклы процессора, что, я думаю, достаточно справедливо.