实现基于整数的幂函数pow(int, int)的最有效方法

用C语言求一个整数的幂的最有效方法是什么?

// 2^3
pow(2,3) == 8


// 5^5
pow(5,5) == 3125
257427 次浏览
小开

一种非常特殊的情况是，当你需要2^(-x ^ y)时，其中x当然是负的y太大了，不能对int型进行移位。你仍然可以用浮点数在常数时间内完成2^x。

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);
}
使用double作为基类型，可以获得更多的2的幂。 (非常感谢评论者帮助整理这篇文章)。

也有可能在学习更多关于IEEE浮点数的知识后，会出现其他幂的特殊情况。

小开
最佳答案

平方求幂。

int ipow(int base, int exp)
{
int result = 1;
for (;;)
{
if (exp & 1)
result *= base;
exp >>= 1;
if (!exp)
break;
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);
}
小开

注意，平方求幂不是最优方法。这可能是一种适用于所有指数值的通用方法，但对于特定的指数值，可能有更好的序列，需要更少的乘法。

例如，如果你想计算x^15，用平方求幂的方法会给你:

x^15 = (x^7)*(x^7)*x
x^7 = (x^3)*(x^3)*x
x^3 = x*x*x

这一共有6次乘法。

事实证明，这可以通过addition-chain求幂使用“仅仅”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²)²，这也需要三次乘法)。

小开

这是对平方求幂效率的后续讨论。

这种方法的优点是它在log(n)时间内运行。例如，如果你要计算一个巨大的数，比如x^1048575(2^20 - 1)，你只需要循环20次，而不是使用朴素方法的100万+次。

此外，在代码复杂性方面，它比试图找到最优的乘法序列更简单，这是la Pramod的建议。

编辑:

我想我应该在有人指责我可能会溢出之前澄清一下。这种方法假设您有某种巨大的int库。

小开

如果要取2的a次方。最快的方法是按幂位移位。

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

如果你想得到一个整数的2的幂，最好使用shift选项:

pow(2,5)可以用1<<5代替

这样效率更高。

小开

另一个实现(在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);
}
小开

我用递归，如果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);
}
}
小开

迟到的人:

下面是一个解决方案，也是处理y < 0最好的，因为它可以。

  1. 它使用intmax_t的结果作为最大范围。没有提供不适合intmax_t的答案。
  2. powjii(0, 0) --> 1在这个例子中是一个共同的结果。

  3. pow(0,negative)，另一个未定义的结果，返回INTMAX_MAX

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

This code uses a forever loop for(;;) to avoid the final base *= base common in other looped solutions. That multiplication is 1) not needed and 2) could be int*int overflow which is UB.

小开

我已经实现了记忆所有计算权力的算法，然后在需要时使用它们。比如x^13等于(x^2)^2^2 * x^2 * x其中x^2^2是从表中取出来的而不是再计算一次。这基本上是@Pramod answer的实现(但在c#中)。 需要的乘法次数是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);
}
小开

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(exp)

power()函数为负exp和浮点基数工作。

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(exp)

小开

我的情况有点不同，我试图用一种力量创造一个面具，但我想无论如何我都要分享我找到的解决方案。

显然，它只适用于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;
}
小开

除了Elias的答案，当使用有符号整数实现时，会导致未定义行为，当使用无符号整数实现时，会导致高输入的不正确值，

下面是平方求幂的修改版本，它也适用于有符号整数类型，并且不会给出错误的值:

#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，但任何宽度(带符号或无符号)都可以使用，只需要稍加修改。然而，如果你需要使用一个非固定宽度的整数类型，你将需要通过(int)sqrt(INT_MAX)(在使用int的情况下)或类似的东西来改变SQRT_INT64_MAX，这应该是优化的，但它更丑陋，而不是一个C常量表达式。同样，将sqrt()的结果强制转换为int也不是很好，因为在完全平方的情况下，浮点精度很高，但我不知道有任何实现，其中INT_MAX -或任何类型的最大值-是完全平方，你可以接受。

小开

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


}
小开 
int pow(int const x, unsigned const e) noexcept
{
return !e ? 1 : 1 == e ? x : (e % 2 ? x : 1) * pow(x * x, e / 2);
//return !e ? 1 : 1 == e ? x : (((x ^ 1) & -(e % 2)) ^ 1) * pow(x * x, e / 2);
}

是的，它是递归的，但是一个好的优化编译器会优化递归。

小开

下面是一个O(1)算法用于计算x ** y，灵感来自这样的评论。它适用于32位signed int。

对于y的小值，它使用平方求幂。对于较大的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;
}
}
}
小开

我注意到gnu-GMP的标准指数平方算法有一些奇怪的地方:

我实现了两个几乎相同的函数——一个是幂模函数，使用最普通的二进制指数平方算法，

  • 标签______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进制的范例中。

我是否在代码/方法中遗漏了一些关键或明显的东西，可能会以令人困惑的方式歪曲结果?谢谢。


提交答案