fft.cpp

支持无限位大数的加减乘除平方运算。运用了FFT

///////////////////////////////////////////////////////////
/// 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 < NumberOne.NumberSize && j < NumberTwo.NumberSize){
		Result.NumberArray[count+1] = GetCarry( (NumberOne.NumberArray[i] + NumberTwo.NumberArray[j] + Result.NumberArray[count]) );
		Result.NumberArray[count++] = (long)(NumberOne.NumberArray[i++] + NumberTwo.NumberArray[j++] + Result.NumberArray[count]+100)%10;
	}

	while( i < NumberOne.NumberSize ){
		Result.NumberArray[count+1] = GetCarry( (NumberOne.NumberArray[i] + Result.NumberArray[count]) );
		Result.NumberArray[count++] = (long)(NumberOne.NumberArray[i++] + Result.NumberArray[count]+100)%10;
	}
	while( j < NumberTwo.NumberSize ){
		Result.NumberArray[count+1] = GetCarry( (NumberTwo.NumberArray[j] + Result.NumberArray[count]) );
		Result.NumberArray[count++] = (long)(NumberTwo.NumberArray[j++] + Result.NumberArray[count]+100)%10;
	}

	if( Result.NumberArray[count] >0 )
		count++;
	Result.DotPosition = DotPosition;
	Result.NumberSize = count;
	return Result;
}

////////////////////////////////////// Subtraction ////////////////////////
Number Subtraction( Number NumberOne, Number NumberTwo )
{
	long i;
	Number Result;
	Result.NumberArray = new double[ max( NumberOne.NumberSize, NumberTwo.NumberSize) + 10];//

	if( Campare( NumberOne, NumberTwo) == 1 ){
		for( i = 0; i < NumberTwo.NumberSize; i++ )
			NumberTwo.NumberArray[i] = - NumberTwo.NumberArray[i];
		return Add( NumberOne, NumberTwo );
	}
	else {
		for( i = 0; i < NumberTwo.NumberSize; i++ )
			NumberOne.NumberArray[i] = - NumberOne.NumberArray[i];
		Result = Add( NumberTwo, NumberOne );
		Result.NumberArray[Result.NumberSize-1] = -Result.NumberArray[Result.NumberSize-1];
		return Result;
	}
}
/////////////////////////////////////////////////// Multiple /////////////////////////
Number Multiple( Number &NumberOne, Number &NumberTwo )
{
	Number Result;
	Result.NumberArray = new double[ NumberOne.NumberSize + NumberTwo.NumberSize + 10];
	long MaxLength = 8192;

	while (MaxLength < 2*(NumberOne.NumberSize + NumberTwo.NumberSize) )
	  MaxLength *= 2;
	if( MaxLength < 128 )
		MaxLength = 128;

	InitializeFFT( MaxLength);
	
	MulWithFFT( NumberOne.NumberArray, NumberOne.NumberSize, NumberTwo.NumberArray, NumberTwo.NumberSize, Result.NumberArray );
	Result.NumberSize = NumberOne.NumberSize + NumberTwo.NumberSize - 1;
	Result.DotPosition = NumberOne.DotPosition + NumberTwo.DotPosition;
	ResultDeal( Result );
	return Result;
}
////////////////////////////////////////// Square //////////////////////////////////////
Number Square( Number NumberOne ) 
{
	return Multiple( NumberOne, NumberOne );
}

//////////////////////////////////////// Divid ///////////////////////////////////
Number Divid( Number NumberOne, Number NumberTwo )
{
	long i,j;
	Number Result;
	Result.NumberArray = new double[Precision*3];//

	memset( Result.NumberArray, 0, Precision*3*sizeof(double));
	Result.NumberArray[0] = (long)(10/NumberTwo.NumberArray[ NumberTwo.NumberSize - 1]);
	Result.DotPosition = NumberTwo.NumberSize ;
	Result.NumberSize = NumberTwo.NumberSize + 1;
	
	for( i = 0; ; i++ )
	{
		Result = Subtraction( Add( Result, Result) , Multiple( Square( Result ), NumberTwo ) );
		if( 2*Result.NumberSize + NumberTwo.NumberSize > Precision*3 )
			break;
	}
	return Multiple( NumberOne, Result);
}

//////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////// main //////////////////////////////////////
int main()
{
	int count;
	string str1, str, str2 ;
	Number NumberOne, NumberTwo, Result;

	printf("\n\t\t\tFFT Based BigNumber Calculator \n\nInput Form : [ A Operator B ]; Operatior Include +, -, *, /\n\t     Or [A 2] : calculator the square of A\n\n");
	while( 1 )
	{
		cout << "Input : ";
		cin >> str1 >> str;
		if( str == "2")
		{
			InitialInput( NumberOne, str1 );
			Result = Square( NumberOne );
			OutPut( Result, Result.NumberSize );
			continue;
		}
		cin >> str2;
		if( str1 == "exit" || str1 == "EXIT" )
			break;
		
		////////////////// input error ////////
		if( InitialInput( NumberOne, str1 ) == -1 || InitialInput( NumberTwo, str2 )==-1 )
		{
			cout << "Error : Invalid Number\n\n";
			continue;
		}
		if( str.length() != 1 ){
			cout << "Error : Invalid Operator\n\n";
			continue;
		}
		if( str[0] == '/' && str2 == "0" )
		{
			cout << "Error : Divid By Zero\n\n";
			continue;
		}
		/////////////////////////////
		
		switch( str[0] )
		{
		case '*':
			Result = Multiple( NumberOne, NumberTwo );
			OutPut( Result, Result.NumberSize );
			break;
		case '+':
			Result = Add( NumberOne, NumberTwo );
			OutPut( Result,Result.NumberSize );
			break;
		case '-':
			Result = Subtraction( NumberOne, NumberTwo );
			OutPut( Result,Result.NumberSize );
			break;
		case '/':
			Result = Divid( NumberOne, NumberTwo );
			OutPut( Result , Precision );
			break;
		default:
			cout << "Invalid Operator"<<endl;
		}
	}
}