///////////////////////////////////////////////////////////
/// FFT Based Big Number Calculator
// Author :05hhb 王美宏.
// PS : 调整Peicision 的大小可以改变除法输出结果的有效位数。
/////////////////////////////////////////////////////////
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <string>
#include <cctype>
#include <iostream>
using namespace std;
const double PI = acos( -1 );
typedef double Real;
struct Complex {
Real R,I;
};
struct Number
{
double *NumberArray;
long DotPosition;
long NumberSize;
};
const Precision = 500;
long FFTLengthMax;
Complex * OmegaFFT;
Complex * ArrayFFT0, * ArrayFFT1;
Complex * ComplexCoef;
double FFTSquareWorstError;
long AllocatedMemory;
long max( long a, long b)
{
return a > b ? a : b;
}
//////////////////////////////////////////// initial FFt //////////////////////////
void InitializeFFT(long MaxLength)
{
long i;
Real Step;
FFTLengthMax = MaxLength;
OmegaFFT = (Complex *) malloc(MaxLength/2*sizeof(Complex));
ArrayFFT0 = (Complex *) malloc(MaxLength*sizeof(Complex));
ArrayFFT1 = (Complex *) malloc(MaxLength*sizeof(Complex));
ComplexCoef = (Complex *) malloc(MaxLength*sizeof(Complex));
Step = 2.0*PI/(double)MaxLength;
for ( i = 0; 2*i < MaxLength; i++ ) {
OmegaFFT[i].R = cos(Step*(double)i);
OmegaFFT[i].I = sin(Step*(double)i);
}
FFTSquareWorstError=0.;
AllocatedMemory = 7*MaxLength*sizeof(Complex)/2;
}
/////////////////////////////////////////// recursively FFT /////////////////////
void RecursiveFFT(Complex * Coef, Complex * FFT, long Length, long Step, long Sign)
{
long i, OmegaStep;
Complex * FFT0, * FFT1, * Omega;
Real tmpR, tmpI;
if (Length==2) {
FFT[0].R = Coef[0].R + Coef[Step].R;
FFT[0].I = Coef[0].I + Coef[Step].I;
FFT[1].R = Coef[0].R - Coef[Step].R;
FFT[1].I = Coef[0].I - Coef[Step].I;
return;
}
FFT0 = FFT;
RecursiveFFT(Coef ,FFT0,Length/2,Step*2,Sign);
FFT1 = FFT+Length/2;
RecursiveFFT(Coef+Step,FFT1,Length/2,Step*2,Sign);
Omega = OmegaFFT;
OmegaStep = FFTLengthMax/Length;
for ( i = 0; 2*i < Length; i++, Omega += OmegaStep) {
tmpR = Omega[0].R*FFT1[i].R-Sign*Omega[0].I*FFT1[i].I;
tmpI = Omega[0].R*FFT1[i].I+Sign*Omega[0].I*FFT1[i].R;
FFT1[i].R = FFT0[i].R - tmpR;
FFT1[i].I = FFT0[i].I - tmpI;
FFT0[i].R = FFT0[i].R + tmpR;
FFT0[i].I = FFT0[i].I + tmpI;
}
}
//////////////////////////////////////////// FFT ////////////////////////////
void FFT(Real * Coef, long Length, Complex * FFT, long NFFT)
{
long i;
for (i=0; i<Length; i++) {
ComplexCoef[i].R = Coef[i];
ComplexCoef[i].I = 0.;
}
for (; i<NFFT; i++)
ComplexCoef[i].R = ComplexCoef[i].I = 0.;
RecursiveFFT(ComplexCoef,FFT,NFFT,1,1);
}
////////////////////////////////////// Inverse FFT /////////////////////////////
void InverseFFT(Complex * FFT, long NFFT, Real * Coef, long Length)
{
long i;
Real invNFFT = 1.0/NFFT, tmp;
RecursiveFFT(FFT, ComplexCoef, NFFT, 1, -1);
for (i=0; i<Length; i++) {
tmp = ComplexCoef[i].R/NFFT;
Coef[i] = floor(0.5+tmp);
if ((tmp-Coef[i])*(tmp-Coef[i])>FFTSquareWorstError)
FFTSquareWorstError = (tmp-Coef[i])*(tmp-Coef[i]);
}
}
/////////////////////////////////////////////// Convolution ///////////////////////////
void Convolution(Complex * A, Complex * B, long NFFT, Complex * C)
{
long i;
Real tmpR, tmpI;
for (i=0; i<NFFT; i++) {
tmpR = A[i].R*B[i].R-A[i].I*B[i].I;
tmpI = A[i].R*B[i].I+A[i].I*B[i].R;
C[i].R = tmpR;
C[i].I = tmpI;
}
}
//////////////////////////////////////////// Multiple with FFT ///////////////////
void MulWithFFT(Real * ACoef, long ASize,Real * BCoef, long BSize, Real * CCoef)
{
long NFFT = 2;
while (NFFT<ASize+BSize)
NFFT *= 2;
if (NFFT>FFTLengthMax) {
printf("Error, FFT Size is too big in MulWithFFT\n");
return;
}
FFT(ACoef, ASize, ArrayFFT0, NFFT);
FFT(BCoef, BSize, ArrayFFT1, NFFT);
Convolution(ArrayFFT0,ArrayFFT1,NFFT,ArrayFFT0);
InverseFFT(ArrayFFT0,NFFT,CCoef, ASize+BSize-1);
}
///////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////// initial the input ///////////////////
int InitialInput( Number &Num, string str )
{
long i, count;
Num.NumberArray = new double[ str.length()+10 ];
Num.DotPosition = str.find( '.' );
if( Num.DotPosition !=-1 )
{
str.erase( Num.DotPosition ,1);
Num.DotPosition = str.length() - Num.DotPosition ;
}
else Num.DotPosition = 0;
for( i = str.length() - 1, count=0; i >= 0; i-- )
{
if( (str[i] > '9' || str[i] < '0') && str[i]!='.')
return -1;
Num.NumberArray[count++] = str[i] - '0';
}
Num.NumberSize = count;
return 1;
}
//////////////////////////////////////////// deal wit the result /////////////////////////////
void ResultDeal( Number &Result )
{
long i, Carry;
double Temp;
for( i = 0, Carry = 0; i < Result.NumberSize ; i++ )
{
Temp = Result.NumberArray[i];
Result.NumberArray[i] = (long)(Result.NumberArray[i] + Carry ) % 10;
Carry = ( Temp + Carry ) / 10;
}
if( Carry )
{
Result.NumberArray[i] = Carry;
Result.NumberSize++;
}
/////////// delete the additional aero
while( Result.NumberSize > Result.DotPosition+1 && Result.NumberArray[Result.NumberSize-1] == 0 )
Result.NumberSize--;
////////////////////////////////
}
//////////////////////////////////////////////// OutPut ///////////////////////////////
void OutPut( Number Result , long count)
{
long i;
cout << "Result : ";
for( i = Result.NumberSize - 1; i >= Result.NumberSize - count; i-- )
{
if( i == Result.DotPosition - 1 )
cout << ".";
cout << Result.NumberArray[i];
}
cout << endl <<endl;
}
/////////////////////////////////////////////// calculate the carry /////////////////////
long GetCarry( long a )
{
if( a >= 0 )
return a/10;
if( a < 0 )
return a/10 - 1;
}
/////////////////////////////////////////////// campare two number ////////////////////
long Campare( Number NumberOne, Number NumberTwo)
{
if( NumberOne.NumberSize - NumberOne.DotPosition > NumberTwo.NumberSize - NumberTwo.DotPosition )
return 1;
else if( NumberOne.NumberSize - NumberOne.DotPosition < NumberTwo.NumberSize - NumberTwo.DotPosition )
return -1;
else {
long i = NumberOne.NumberSize - 1, j = NumberTwo.NumberSize - 1;
while( i >= 0 && j >= 0 && NumberOne.NumberArray[i] == NumberTwo.NumberArray[j] )
i--, j--;
if( j == -1 || NumberOne.NumberArray[i] > NumberTwo.NumberArray[j] )
return 1;
if( i == -1 || NumberOne.NumberArray[i] < NumberTwo.NumberArray[j] )
return -1;
}
}
/////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////// Add /////////////////////////////////////
Number Add( Number NumberOne, Number NumberTwo )
{
long i = 0, j = 0, count=0;
long DotPosition = max( NumberOne.DotPosition, NumberTwo.DotPosition);
long Size = max (NumberOne.NumberSize - NumberOne.DotPosition, NumberTwo.NumberSize - NumberTwo.DotPosition ) + DotPosition + 5;
Number Result;
Result.NumberArray = new double [ Size + 10 ];//
memset( Result.NumberArray, 0, Size*sizeof(double));
if( NumberOne.DotPosition >= NumberTwo.DotPosition ){
while( i < NumberOne.DotPosition - NumberTwo.DotPosition ){
Result.NumberArray[count+1] = GetCarry( NumberOne.NumberArray[i] + Result.NumberArray[count]);
Result.NumberArray[count++] = (long)(NumberOne.NumberArray[i++] + Result.NumberArray[count] + 100)%10;;
}
}
else {
while( j < NumberTwo.DotPosition - NumberOne.DotPosition ){
Result.NumberArray[count+1] = GetCarry( NumberTwo.NumberArray[j] + Result.NumberArray[count]);
Result.NumberArray[count++] = (long)( NumberTwo.NumberArray[j++] + Result.NumberArray[count] + 100)%10;
}
}
while( i <
评论0