快速幂算法
博客园 快速幂算法
CSDN 快速幂算法详解
算法导论：Strassen算法
整数运算
基本原理：
(a + b) % p = (a % p + b % p) % p （1）
(a - b) % p = (a % p - b % p ) % p （2）
**(a * b) % p = (a % p * b % p) % p （3）**
最后一个mod p 是为了保证结果不超过p，多个因子连续的乘积取模的结果等于每个因子取模后的乘积再取模的结果。提前取模以避免结果过大的爆表。
a*b
快速乘法的基本思想 ，是二进制和乘法分配律的结合（浮点数不满足结合律），比如说，13 ==(1101)2 ，4 *13等于4 *(1101)2 ，用分配律展开得到4 *13 == 4 *(1000+100+1)2，我们不难观察出，快速幂可以通过判断当前的位（bit）是1还是0，推断出是否需要做求和操作，每次移动到下一位（bit）时，就对 ans进行 *2操作 ，等待是否求和。由于除以2和位移操作是等效的，因此这也可以看作是二分思想的应用，这种算法将b进行二分从而减少了不必要的运算，时间复杂度是log（n）。
a^b
快速幂其实可以看作是快速乘法的特例，在快速幂中，我们不再对ans进行 *2操作，因为在a ^b中b的意义已经从乘数变成了指数，但是我们可以仍然把b写成二进制，举例说明：此时，我们将4 *13改为4 ^13，13=（1101）2 ，二进制13写开我们得到（1000+100+1），注意，这里的所有二进制是指数，指数的相加意味着底数相乘，因此有4 ^13 == (4 ^8) * (4 ^4) * (4 ^1)。再注意到指数之间的2倍关系，我们就可以用很少的几个变量，完成这一算法。这样，我们就将原本用循环需要O(n)的算法，改进为O（logN）的算法。
- 快速乘法
//快速乘法
int qmul(int a,int b)
{ // 根据数据范围可选择long long
int ans=0;
while(b)
{
if( b&1)ans+=a; //按位与完成位数为1的判断
b>>=1;a<<=1; //位运算代替/2和*2
}
return ans;
}
如果涉及到快速乘法取模，则需要进行一些微小改动
//快速乘法取模
int qmul_mod(int a,int b,int mod)
{
int ans=0;
while(b)
{
if((b%=mod)&1)ans+=a%=mod;//这里需要b%=mod 以及a%=mod
b>>=1;a<<=1;
}
return ans%mod; //ans也需要对mod取模
}
- 快速幂
//快速幂 a^b
int qpow(int a,int b)
{
if(a==0) return 0;//attention
int ans=1;
while(b)
{
if(b&1)ans*=a;//和快速乘法的区别
b>>=1;a*=a;//区别，同上
}
return ans;
}
以及含有取模的快速幂：
int qpow_mod(int a,int b,int mod){
if(a==0)return 0;
int ans=1;
while(b){
if(b&1)ans=(ans%mod)*(a%mod);//如果确定数据不会爆的话，可写成 ans*=a%=mod;
b>>=1;a*=a%=mod;//等价于a=(a%mod)*(a%mod)，且将一个模运算通过赋值代替，提高了效率
}
return ans%mod;//数据不会爆的话，这里的%运算会等价于第5中不断重复的 ans%mod
}
如果我们对于性能还有更进一步的要求，那么也就是减少取模运算了，那么我们需要确定数据范围不会爆掉
在这样的前提下，我们可以只用原先1/4的取模运算量完成快速幂
int qpow_mod(int a,int b,int mod){
if(!a)return 0;
int ans=1;
while(b){
if(b&1)ans*=a%=mod;//这里的模运算只有一个
b>>=1;a*=a;//这里的模运算没有了
}
return ans%mod;
}
矩阵幂
依然可以按上述方法实现，底数改为矩阵，如 lc.509 斐波那契 那样做，注意到，我们可以构造函数，同样可以重载字符来完成。
当指数足够大，(N>300)，我们考虑Strassen算法。
矩阵原本相乘需要8次乘法：
r = a * e + b * g ;
s = a * f + b * h ;
t = c * e + d * g;
u = c * f + d * h;
运用strassen提出的Strassen算法，可以将8次乘法降为7次乘法，虽然只是一次乘法，但是其实一次乘法耗时要比加减法多很多。
处理的方法是写成：
p1 = a * ( f - h )
p2 = ( a + b ) * h
p3 = ( c +d ) * e
p4 = d * ( g - e )
p5 = ( a + d ) * ( e + h )
p6 = ( b - d ) * ( g + h )
p7 = ( a - c ) * ( e + f )
r = p5 + p4 + p6 - p2
s = p1 + p2
t = p3 + p4
u = p5 + p1 - p3 - p7
复杂度：O( n^lg7 ) ~= O( n^2.81 )
/*
Strassen Algorithm Implementation in C++
Coded By: Seyyed Hossein Hasan Pour MatiKolaee in May 5 2010 .
Mazandaran University of Science and Technology,Babol,Mazandaran,Iran
--------------------------------------------
Email : Master.huricane@gmail.com
YM : Deathmaster_nemessis@yahoo.com
Updated may 09 2010.
*/
#include <iostream>
#include <cstdlib>
#include <iomanip>
#include <ctime>
#include <windows.h>
using namespace std;
int Strassen(int n, int** MatrixA, int ** MatrixB, int ** MatrixC);//Multiplies Two Matrices recrusively.
int ADD(int** MatrixA, int** MatrixB, int** MatrixResult, int length );//Adds two Matrices, and places the result in another Matrix
int SUB(int** MatrixA, int** MatrixB, int** MatrixResult, int length );//subtracts two Matrices , and places the result in another Matrix
int MUL(int** MatrixA, int** MatrixB, int** MatrixResult, int length );//Multiplies two matrices in conventional way.
void FillMatrix( int** matrix1, int** matrix2, int length);//Fills Matrices with random numbers.
void PrintMatrix( int **MatrixA, int MatrixSize );//prints the Matrix content.
int main()
{
int MatrixSize = 0;
int** MatrixA;
int** MatrixB;
int** MatrixC;
clock_t startTime_For_Normal_Multipilication ;
clock_t endTime_For_Normal_Multipilication ;
clock_t startTime_For_Strassen ;
clock_t endTime_For_Strassen ;
time_t start,end;
srand(time(0));
cout<<setw(45)<<"In the name of GOD";
cout<<endl<<setw(60)<<"Strassen Algorithm Implementation in C++ "
<<endl<<endl<<setw(50)<<"By Seyyed Hossein Hasan Pour"
<<endl<<setw(60)<<"Mazandaran University of Science and Technology"
<<endl<<setw(40)<<"May 9 2010";
cout<<"\nPlease Enter your Matrix Size(must be in a power of two(eg:32,64,512,..): ";
cin>>MatrixSize;
int N = MatrixSize;//for readiblity.
MatrixA = new int *[MatrixSize];
MatrixB = new int *[MatrixSize];
MatrixC = new int *[MatrixSize];
for (int i = 0; i < MatrixSize; i++)
{
MatrixA[i] = new int [MatrixSize];
MatrixB[i] = new int [MatrixSize];
MatrixC[i] = new int [MatrixSize];
}
FillMatrix(MatrixA,MatrixB,MatrixSize);
//*******************conventional multiplication test
cout<<"Phase I started: "<< (startTime_For_Normal_Multipilication = clock());
MUL(MatrixA,MatrixB,MatrixC,MatrixSize);
cout<<"\nPhase I ended: "<< (endTime_For_Normal_Multipilication = clock());
cout<<"\nMatrix Result... \n";
PrintMatrix(MatrixC,MatrixSize);
//*******************Strassen multiplication test
cout<<"\nMultiplication started: "<< (startTime_For_Strassen = clock());
Strassen( N, MatrixA, MatrixB, MatrixC );
cout<<"\nMultiplication: "<<(endTime_For_Strassen = clock());
cout<<"\nMatrix Result... \n";
PrintMatrix(MatrixC,MatrixSize);
cout<<"Matrix size "<<MatrixSize;
cout<<"\nNormal mode "<<(endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication)<<" Clocks.."<<(endTime_For_Normal_Multipilication - startTime_For_Normal_Multipilication)/CLOCKS_PER_SEC<<" Sec";
cout<<"\nStrassen mode "<<(endTime_For_Strassen - startTime_For_Strassen)<<" Clocks.."<<(endTime_For_Strassen - startTime_For_Strassen)/CLOCKS_PER_SEC<<" Sec\n";
system("Pause");
return 0;
}
/*
in order to be able to create a matrix without any limitaion in c++,
one way is to create it using pointers.
as you see by using a pointer to pointer strategy we can make a multi-
dimensional Matrix of any size . The notation also makes us capable of
creating a matrix with VARIABLE size at runtime ,meaning we can resize
the size of our matrix at runtime , shrink it or increase it , your choice.
what we do is simple , first we make a pointer of pointer variable , this
means that our first pointer will point to another pointer which again
this pointer ,points to sth else(we can make it point to an array) .
int **A;
will declare the variable , we now need to expand it .
now make a pointer based array and allocate the memory dynamicly
A = new int *[desired_array_row];
this gives us a one diminsional pointer based array,now you want a 2D array?
big deal,lets make one.
we use for() to achieve this goal , remember when i said we are going to make
a variable which is a pointer of pointer ? which meant any location pointed to somewhere else
, we made a pointer based array , a one diminsional one , just up there ,
and you know this fatct that an array is consits of individual blocks right?
and the fact that each block can be used just like a solo variable.
so simply if we could write
A = new int *[any_size];
cant we do it to all of our indiviual array blocks which are just like the solo variable ?
so this means that if we could do it with A, and get an array , we can use the same method
to make different arrays for different block of the array we made in first place.
we use for() to iterate through all of the blocks of the previously made array, and
then for each block we create a single array .
for ( int i = 0; i < desired_array_row; i++)
A[i] = new int [desired_column_size];
after this for , we can enjoy our 2D array wich can be access like any ordinary array we know.
just use the conventional notation for accessing array blocks for either reading or writing.( A[i][j])
and remember to free the space we allocated for our 2D array at the end of the program .
we do such a thing this way:
for ( int i = 0; i < your_array_row; i++)
{
delete [] A[i];
}
delete[] A;
.using this method you can make any N-diminsional array, you just need to use for with right iteration.
*/
int Strassen(int N, int **MatrixA, int **MatrixB, int **MatrixC)
{
int HalfSize = N/2;
int newSize = N/2;
if ( N <= 64 )//choosing the threshhold is extremely important, try N<=2 to see the result
{
MUL(MatrixA,MatrixB,MatrixC,N);
}
else
{
int** A11;
int** A12;
int** A21;
int** A22;
int** B11;
int** B12;
int** B21;
int** B22;
int** C11;
int** C12;
int** C21;
int** C22;
int** M1;
int** M2;
int** M3;
int** M4;
int** M5;
int** M6;
int** M7;
int** AResult;
int** BResult;
//making a 1 diminsional pointer based array.
A11 = new int *[newSize];
A12 = new int *[newSize];
A21 = new int *[newSize];
A22 = new int *[newSize];
B11 = new int *[newSize];
B12 = new int *[newSize];
B21 = new int *[newSize];
B22 = new int *[newSize];
C11 = new int *[newSize];
C12 = new int *[newSize];
C21 = new int *[newSize];
C22 = new int *[newSize];
M1 = new int *[newSize];
M2 = new int *[newSize];
M3 = new int *[newSize];
M4 = new int *[newSize];
M5 = new int *[newSize];
M6 = new int *[newSize];
M7 = new int *[newSize];
AResult = new int *[newSize];
BResult = new int *[newSize];
int newLength = newSize;
//making that 1 diminsional pointer based array , a 2D pointer based array
for ( int i = 0; i < newSize; i++)
{
A11[i] = new int[newLength];
A12[i] = new int[newLength];
A21[i] = new int[newLength];
A22[i] = new int[newLength];
B11[i] = new int[newLength];
B12[i] = new int[newLength];
B21[i] = new int[newLength];
B22[i] = new int[newLength];
C11[i] = new int[newLength];
C12[i] = new int[newLength];
C21[i] = new int[newLength];
C22[i] = new int[newLength];
M1[i] = new int[newLength];
M2[i] = new int[newLength];
M3[i] = new int[newLength];
M4[i] = new int[newLength];
M5[i] = new int[newLength];
M6[i] = new int[newLength];
M7[i] = new int[newLength];
AResult[i] = new int[newLength];
BResult[i] = new int[newLength];
}
//splitting input Matrixes, into 4 submatrices each.
for (int i = 0; i < N / 2; i++)
{
for (int j = 0; j < N / 2; j++)
{
A11[i][j] = MatrixA[i][j];
A12[i][j] = MatrixA[i][j + N / 2];
A21[i][j] = MatrixA[i + N / 2][j];
A22[i][j] = MatrixA[i + N / 2][j + N / 2];
B11[i][j] = MatrixB[i][j];
B12[i][j] = MatrixB[i][j + N / 2];
B21[i][j] = MatrixB[i + N / 2][j];
B22[i][j] = MatrixB[i + N / 2][j + N / 2];
}
}
//here we calculate M1..M7 matrices .
//M1[][]
ADD( A11,A22,AResult, HalfSize);
ADD( B11,B22,BResult, HalfSize);
Strassen( HalfSize, AResult, BResult, M1 ); //now that we need to multiply this , we use the strassen itself .
//M2[][]
ADD( A21,A22,AResult, HalfSize); //M2=(A21+A22)B11
Strassen(HalfSize, AResult, B11, M2); //Mul(AResult,B11,M2);
//M3[][]
SUB( B12,B22,BResult, HalfSize); //M3=A11(B12-B22)
Strassen(HalfSize, A11, BResult, M3); //Mul(A11,BResult,M3);
//M4[][]
SUB( B21, B11, BResult, HalfSize); //M4=A22(B21-B11)
Strassen(HalfSize, A22, BResult, M4); //Mul(A22,BResult,M4);
//M5[][]
ADD( A11, A12, AResult, HalfSize); //M5=(A11+A12)B22
Strassen(HalfSize, AResult, B22, M5); //Mul(AResult,B22,M5);
//M6[][]
SUB( A21, A11, AResult, HalfSize);
ADD( B11, B12, BResult, HalfSize); //M6=(A21-A11)(B11+B12)
Strassen( HalfSize, AResult, BResult, M6); //Mul(AResult,BResult,M6);
//M7[][]
SUB(A12, A22, AResult, HalfSize);
ADD(B21, B22, BResult, HalfSize); //M7=(A12-A22)(B21+B22)
Strassen(HalfSize, AResult, BResult, M7); //Mul(AResult,BResult,M7);
//C11 = M1 + M4 - M5 + M7;
ADD( M1, M4, AResult, HalfSize);
SUB( M7, M5, BResult, HalfSize);
ADD( AResult, BResult, C11, HalfSize);
//C12 = M3 + M5;
ADD( M3, M5, C12, HalfSize);
//C21 = M2 + M4;
ADD( M2, M4, C21, HalfSize);
//C22 = M1 + M3 - M2 + M6;
ADD( M1, M3, AResult, HalfSize);
SUB( M6, M2, BResult, HalfSize);
ADD( AResult, BResult, C22, HalfSize);
//at this point , we have calculated the c11..c22 matrices, and now we are going to
//put them together and make a unit matrix which would describe our resulting Matrix.
for (int i = 0; i < N/2 ; i++)
{
for (int j = 0 ; j < N/2 ; j++)
{
MatrixC[i][j] = C11[i][j];
MatrixC[i][j + N / 2] = C12[i][j];
MatrixC[i + N / 2][j] = C21[i][j];
MatrixC[i + N / 2][j + N / 2] = C22[i][j];
}
}
// dont forget to free the space we alocated for matrices,
for (int i = 0; i < newLength; i++)
{
delete[] A11[i];delete[] A12[i];delete[] A21[i];
delete[] A22[i];
delete[] B11[i];delete[] B12[i];delete[] B21[i];
delete[] B22[i];
delete[] C11[i];delete[] C12[i];delete[] C21[i];
delete[] C22[i];
delete[] M1[i];delete[] M2[i];delete[] M3[i];delete[] M4[i];
delete[] M5[i];delete[] M6[i];delete[] M7[i];
delete[] AResult[i];delete[] BResult[i] ;
}
delete[] A11;delete[] A12;delete[] A21;delete[] A22;
delete[] B11;delete[] B12;delete[] B21;delete[] B22;
delete[] C11;delete[] C12;delete[] C21;delete[] C22;
delete[] M1;delete[] M2;delete[] M3;delete[] M4;delete[] M5;
delete[] M6;delete[] M7;
delete[] AResult;
delete[] BResult ;
}//end of else
return 0;
}
int ADD(int** MatrixA, int** MatrixB, int** MatrixResult, int MatrixSize )
{
for ( int i = 0; i < MatrixSize; i++)
{
for ( int j = 0; j < MatrixSize; j++)
{
MatrixResult[i][j] = MatrixA[i][j] + MatrixB[i][j];
}
}
return 0;
}
int SUB(int** MatrixA, int** MatrixB, int** MatrixResult, int MatrixSize )
{
for ( int i = 0; i < MatrixSize; i++)
{
for ( int j = 0; j < MatrixSize; j++)
{
MatrixResult[i][j] = MatrixA[i][j] - MatrixB[i][j];
}
}
return 0;
}
int MUL( int** MatrixA, int** MatrixB, int** MatrixResult, int MatrixSize )
{
for (int i=0;i<MatrixSize ;i++)
{
for (int j=0;j<MatrixSize ;j++)
{
MatrixResult[i][j]=0;
for (int k=0;k<MatrixSize ;k++)
{
MatrixResult[i][j]=MatrixResult[i][j]+MatrixA[i][k]*MatrixB[k][j];
}
}
}
return 0;
}
void FillMatrix( int** MatrixA, int** MatrixB, int length)
{
for(int row = 0; row<length; row++)
{
for(int column = 0; column<length; column++)
{
MatrixB[row][column] = (MatrixA[row][column] = rand() %5);
//matrix2[row][column] = rand() % 2;//ba hazfe in khat 50% afzayeshe soorat khahim dasht
}
}
}
void PrintMatrix(int **MatrixA,int MatrixSize)
{
cout<<endl;
for(int row = 0; row<MatrixSize; row++)
{
for(int column = 0; column<MatrixSize; column++)
{
cout<<MatrixA[row][column]<<"\t";
if ((column+1)%((MatrixSize)) == 0)
cout<<endl;
}
}
cout<<endl;
}
876
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
