factorialtail.cpp

RSA的应用需要用到很多大素数。如何生成大素数

// FactorialTail.cpp : Defines the entry point for the console application.
//

//#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>

#include "../../HugeCalc_Dll_Import/HugeCalc.h"	// 公共接口
#include "../../HugeCalc_Dll_Import/HugeInt.h"	// 10进制系统
#include "../../HugeCalc_Dll_Import/HugeIntX.h"	// 16进制系统

#ifndef _UNICODE
	#pragma comment( lib, "../../HugeCalc_Dll_Import/HugeCalc.lib" )
#else
	#pragma comment( lib, "../../HugeCalc_Dll_Import/HugeCalcU.lib" )
#endif

//	Project -> Setting -> C/C++ -> Code Generation --> Use run-time library:
//		Win32 Debug:	Debug Multithreaded DLL
//		Win32 Release:	Multithreaded DLL
CHugeInt m_HugeN;
CHugeInt m_HugeT;
CHugeInt m_HugeZ;
UINT32 m_u32L;

UINT32 m_u32Pow2( 0 );
CHugeInt m_HugePow2( 1 );
CHugeInt m_HugePow5( 5 );
CHugeInt m_HugePow10( 1 );

const UINT32 GetVal( const CHugeInt &HugeInt )
{
	SINT64 s64Num = -1;
	HugeInt.CanConvertSINT64( s64Num );
	return ( s64Num <= 0x7FFFFFFF) ? (UINT32)s64Num : -1;
}

CHugeInt &ModPow5( CHugeInt &HugeInt )
{
	switch( CompareAbs( HugeInt %= m_HugePow10, m_HugePow5 ))
	{
		case 0: HugeInt = 0;
		case -1: return HugeInt;
		default: break;
	}

	if ( !m_u32Pow2 )
	{
		return HugeInt -= (( CHugeInt( HugeInt ) *= m_HugePow2 ).RShift( m_u32L ) *= m_HugePow5 );
	}
	else
	{
		return HugeInt -= (( CHugeInt( HugeInt ) *= m_u32Pow2 ).RShift( m_u32L ) *= m_HugePow5 );
	}
}

void PolyMul( const CHUGEINT_VECTOR::iterator dst, const CHUGEINT_VECTOR::iterator src1, const CHUGEINT_VECTOR::const_iterator src2, const BOOL bMod5 = FALSE )
{
	UINT32 i,j;
	for( i=0; m_u32L!=i; ++i )
	{
		dst[i] = 0;
		src1[i] %= m_HugePow10;
	}

	CHugeInt tmp;

	for(i=0; m_u32L!=i; ++i)
	{
		if ( !(!src1[i]) )
		{
			for(j=0; m_u32L-i!=j; ++j)
			{
				if ( !(!src2[j]) )
				{
					dst[i+j] += tmp.Mul( src1[i], src2[j] );
				}
			}
		}
	}

	if ( bMod5 )
	{
		for(i=0; m_u32L!=i; ++i)
		{
			ModPow5( dst[i] );
		}
	}
	else
	{
		for(i=0; m_u32L!=i; ++i)
		{
			dst[i] %= m_HugePow10;
		}
	}
}

VOID Calc_Fun( CHUGEINT_VECTOR &vHugeFun )
{
//	#define TRIANGLE(x,y)  triangle[(x)*s_u32L+(y)]
//	#define F(x,y)         s_vHugeFun[(x)*s_u32L+(y)]
//	#define F_T(x,y)       f_t[(x)*s_u32L+(y)]

	static UINT32 s_u32L( 0 );
	static CHUGEINT_VECTOR s_vHugeFun;

	if ( s_vHugeFun.size() < vHugeFun.size() )
	{
		s_vHugeFun.swap( vHugeFun );
	}
	
	vHugeFun.clear();

	UINT32 i,j,k;
	if ( s_u32L == m_u32L )
	{
		vHugeFun.swap( s_vHugeFun );
		return;
	}
	else if ( s_u32L > m_u32L )
	{
		vHugeFun.resize( m_u32L * m_u32L, CHugeInt() );
		CHUGEINT_VECTOR::iterator pV = vHugeFun.begin();
		CHUGEINT_VECTOR::const_iterator p1 = NULL;
		CHUGEINT_VECTOR::const_iterator p2 = NULL;

		for ( i = 0, p1 = s_vHugeFun.begin(); m_u32L != i; ++i, p1 += s_u32L )
		{
			for ( j = 0, p2 = p1; m_u32L != j; ++j, ++pV, ++p2 )
			{
// 				ModPow5( *pV = F( i, j ));
				ModPow5( *pV = *p2 );
			}
		}
		return;
	}
	s_u32L = m_u32L;

	const UINT32 L2 = m_u32L * m_u32L;

	s_vHugeFun.clear();
	s_vHugeFun.resize( L2, CHugeInt() );

	CHUGEINT_VECTOR triangle( L2, CHugeInt() );
	CHUGEINT_VECTOR f_t( L2, CHugeInt() );
	CHUGEINT_VECTOR tmp1( m_u32L, CHugeInt() );
	CHUGEINT_VECTOR tmp2( m_u32L, CHugeInt() );
	CHUGEINT_VECTOR tmp3( m_u32L, CHugeInt() );
	CHUGEINT_VECTOR tmp4( m_u32L, CHugeInt() );
	CHUGEINT_VECTOR tmp5( m_u32L, CHugeInt() );
	CHugeInt tmp,tmpp,tmppp;

	const CHUGEINT_VECTOR::iterator pfun = s_vHugeFun.begin();
	const CHUGEINT_VECTOR::iterator ptriangle = triangle.begin();
	const CHUGEINT_VECTOR::iterator pf_t = f_t.begin();
	const CHUGEINT_VECTOR::iterator ptmp1 = tmp1.begin();
	const CHUGEINT_VECTOR::iterator ptmp2 = tmp2.begin();
	const CHUGEINT_VECTOR::iterator ptmp3 = tmp3.begin();
	const CHUGEINT_VECTOR::iterator ptmp4 = tmp4.begin();
	const CHUGEINT_VECTOR::iterator ptmp5 = tmp5.begin();

	CHUGEINT_VECTOR::iterator p = NULL;
	CHUGEINT_VECTOR::iterator p1 = NULL;
	CHUGEINT_VECTOR::iterator p2 = NULL;

	//Initialize Pascal Triangle
//	TRIANGLE(0,0) = 1;
//	TRIANGLE(1,0) = 1;
//	TRIANGLE(1,1) = 1;
	*( p = ptriangle ) = 1;
	*( p += m_u32L ) = 1;
	*( ++p ) = 1;
	for(i=2; m_u32L!=i; ++i)
	{
//		TRIANGLE(i,0) = 1;
//		TRIANGLE(i,i) = 1;
//		for(j=1; i!=j; ++j)
//			( TRIANGLE(i,j) =TRIANGLE(i-1,j-1) ) += TRIANGLE(i-1,j);
		p = ptriangle + i * m_u32L;	// TRIANGLE(i,0)
		p1 = p + i;						// TRIANGLE(i,i)
		p2 = p - m_u32L;				// TRIANGLE(i-1,0)
		*p = 1;
		while ( p1 != ++p )
		{
			*p = *p2;
			*p += *(++p2);
		}
		*p = 1;
	}

//	F(0,0)=1;
//	F(0,1)=1;
//	//Initialize F(1,x)
//	F(1,0)=24;
//	F(1,1)=50*5;
//	F(1,2)=35*5*5;
//	F(1,3)=10*5*5*5;
//	F(1,4)=1*5*5*5*5;
	*( p = pfun ) = 1;
	*( p + 1 ) = 1;
	*( p += m_u32L ) = 24;
	*( ++p ) = 50*5;
	*( ++p ) = 35*5*5;
	*( ++p ) = 10*5*5*5;
	*( ++p ) = 1*5*5*5*5;

	//set tmp1[i] to be 5^i
	tmp1[0]=1;
	for(i=1; m_u32L!=i; ++i) ( tmp1[i] = tmp1[i-1] ) *= 5;

	for(j=0; m_u32L!=j; ++j)
	{
		for(k=0, p=ptriangle+j*m_u32L, p1=ptmp1; \
			k<=j; \
			++k, ++p, ++p1 )
		{
			ModPow5( *p *= *p1 );
		}
	}

	for(i=2; m_u32L!=i; ++i)
	{
		//Set F(i,x) now
		for(j=0, p1=pfun+(i-1)*m_u32L; \
			m_u32L!=j; \
			++j, ++p1)
		{
			for(k=0, p=pf_t+j*m_u32L, p2=ptriangle+j*m_u32L; \
				k<=j; \
				++k, ++p, ++p2)
			{				
//				ModPow5( F_T(j, k ).Mul( F(i-1,j), TRIANGLE(j,k) ));
				ModPow5( p->Mul( *p1, *p2 ));
			}
		}

		for(j=0, p=ptmp4, p1=pfun+(i-1)*m_u32L, p2=ptmp1; \
			m_u32L!=j; \
			++j, ++p, ++p1, ++p2)
		{ //tmp4 is F(i-1,j)*5^j mod 5^m_u32L
// 			ModPow5( tmp4[j].Mul( F(i-1,j), tmp1[j] ));
			ModPow5( p->Mul( *p1, *p2 ));
		}

		for(j=0; m_u32L!=j; ++j) tmp3[j]=0;
		for(j=0; m_u32L!=j; ++j)
		{//Adding F(i-1,j)*(5x+1)^j mod 5^m_u32L into tmp3
			for(k=0, p=ptmp3, p1=pf_t+j*m_u32L; \
				k<=j; \
				++k, ++p, ++p1)
			{
//				tmp3[k] += F_T(j, k );
 				*p += *p1;
			}
		}
		//Next get Multiplication of polynomial tmp4, tmp3 into tmp5
		PolyMul( ptmp5, ptmp3, ptmp4, FALSE ); //tmp5=tmp3*tmp4;

		tmp2[0]=1; //set tmp2 to 2^j
		for(j=1; m_u32L!=j; ++j) ( tmp2[j] = tmp2[j-1] ) *= 2;

		for(j=0; m_u32L!=j; ++j) tmp3[j]=0;
		for(j=0; m_u32L!=j; ++j)
		{ //Adding F(i-1,j)*(5x+2)^j mod 5^m_u32L into tmp3
			for(k=0, p=ptmp3, p1=pf_t+j*m_u32L, p2=ptmp2+j; \
				k<=j; \
				++k, ++p, ++p1,--p2)
			{
// 				tmp3[k] += tmp.Mul( F_T(j, k ), tmp2[j-k] );//*2^(j-k)
				*p += tmp.Mul( *p1, *p2 );//*2^(j-k)
			}
		}

		PolyMul( ptmp4, ptmp3, ptmp5, FALSE );//tmp4=tmp3*tmp5

		tmp2[0]=1; //set tmp2 to 3^j
		for(j=1; m_u32L!=j; ++j) ( tmp2[j] = tmp2[j-1] ) *= 3;

		for(j=0; m_u32L!=j; ++j) tmp3[j]=0;
		for(j=0; m_u32L!=j; ++j)
		{ //Adding F(i-1,j)*(5x+3)^j mod 5^m_u32L into tmp3
			for(k=0, p=ptmp3, p1=pf_t+j*m_u32L, p2=ptmp2+j; \
				k<=j; \
				++k, ++p, ++p1,--p2)
			{
// 				tmp3[k] += tmp.Mul( F_T(j, k ), tmp2[j-k] );//*3^(j-k)
				*p += tmp.Mul( *p1, *p2 );//*3^(j-k)
			}
		}

		PolyMul( ptmp5, ptmp3, ptmp4, FALSE ); //tmp5=tmp3*tmp4;

		tmp2[0] = 1; //set tmp2 to 4^j
		for(j=1; m_u32L!=j; ++j) (tmp2[j]=tmp2[j-1]) *= 4;

		for(j=0; m_u32L!=j; ++j) tmp3[j]=0;
		for(j=0; m_u32L!=j; ++j)
		{ //Adding F(i-1,j)*(5x+4)^j mod 5^m_u32L into tmp3
			for(k=0, p=ptmp3, p1=pf_t+j*m_u32L, p2=ptmp2+j; \
				k<=j; \
				++k, ++p, ++p1,--p2)
			{
// 				tmp3[k] += tmp.Mul( F_T(j, k ), tmp2[j-k] );//*4^(j-k)
				*p += tmp.Mul( *p1, *p2 );//*4^(j-k)
			}
		}
		PolyMul( pfun+i*m_u32L, ptmp3, ptmp5, TRUE ); //tmp4=tmp3*tmp5;
	}

	for(i=1,p=pfun+i*m_u32L; m_u32L!=i; ++i)
	{
// 		PolyMul( ptmp4, pfun+i*m_u32L, pfun+(i-1)*m_u32L, TRUE );
		PolyMul( p1=ptmp4, p, p-m_u32L, TRUE );
		for(j=0; m_u32L!=j; ++j,++p, ++p1)
// 			F(i,j).Swap( tmp4[j] );
			p->Swap( *p1 );
    }

	triangle.clear();
	f_t.clear();
	tmp1.clear();
	tmp2.clear();
	tmp3.clear();
	tmp4.clear();
	tmp5.clear();

	vHugeFun.swap( s_vHugeFun );

#if 0
	FILE * fptr = fopen("F36.txt","w");
	for( i=1,p=pfun+i*m_u32L; m_u32L!=i; ++i)
	{
		fprintf(fptr,"\n\t");
		for(j=0; m_u32L!=j; ++j,++p)
		{
//			fprintf(fptr, F(i,j).ConvertToStr(FS_NORMAL) );
			fprintf(fptr, p->ConvertToStr(FS_NORMAL) );
			fprintf(fptr,", ");
		}
	}
	fflush(fptr);
	fclose(fptr);
#endif
}

void Get_m5L( CHugeInt &m5L )
{
	CARRY_PARAM CarryParam;
	U32_VECTOR &vU32Num = CarryParam.vU32Num;
	CarryParam.u32Carry = 5;
	// 注意：如果未注册，所得到的结果 CarryParam 有可能被随机干扰，从而导致结果不正确！
	m_HugeN.ConvertToCarry( CarryParam );
	const U32_VECTOR::const_iterator p_begin = vU32Num.begin();
	const U32_VECTOR::const_iterator p_end = vU32Num.end();

	U32_VECTOR::const_iterator p = p_begin;

	UINT32 u32Sum = 0;
	p = p_begin - 1;
	while ( p_end != ++p )
	{
		u32Sum += *p;
	}
	(( m_HugeZ = m_HugeN ) -= u32Sum ) /= 4;

	BOOL bOdd = FALSE;
	p = p_begin + m_u32L;
	while ( p_end > p )
	{
		if ( (*p&1) )
		{
			bOdd = !bOdd;
		}
		p += 2;
	}
	if ( bOdd )
	{
		--( m5L = m_HugePow5 );
	}
	else
	{
		m5L = 1;
	}

	CHUGEINT_VECTOR vHuge5( 5, CHugeInt());
	vHuge5[ 0 ] = 0;
	( vHuge5[ 1 ] = 5 ).Pow( m_u32L - 1 );
	( vHuge5[ 2 ] = vHuge5[ 1 ] ) *= 2;
	( vHuge5[ 3 ] = vHuge5[ 2 ] ) += vHuge5[ 1 ];
	( vHuge5[ 4 ] = vHuge5[ 2 ] ) *= 2;
	
	CHugeInt t1, t2, tmp;
// 	ModPow5( t2 = m_HugeN );
	U32_VECTOR::const_iterator pHead = p_begin + m_u32L;
	{
		if ( p_end < pHead )
		{
			pHead = p_end;
		}
		CARRY_PARAM CarryParamL;
		CarryParamL.nSign = 1;
		CarryParamL.u32Carry = 5;
		CarryParamL.vU32Num.swap( U32_VECTOR( p_begin, pHead ));
		t2 = CarryParamL;
	}

	static UINT32 s_u32L( 0 );
	static CHUGEINT_VECTOR vHugeFun;
	if ( s_u32L != m_u32L )
	{
		Calc_Fun( vHugeFun );
		s_u32L = m_u32L;
	}

	const CHUGEINT_VECTOR::const_iterator pV_end = vHugeFun.end();
	CHUGEINT_VECTOR::const_iterator pV_mem = vHugeFun.begin();
	CHUGEINT_VECTOR::const_iterator pV = pV_mem;

	UINT32 i, j;
	UINT32 u32Step = 0;
	pHead = (--( p = p_begin )) + m_u32L;
	while( p_end != ++p )
	{
		if ( pV_end != pV_mem )
		{
			pV_mem += m_u32L;
		}

		for( i=0; (*p)!=i; ++i )
		{
			--t2;
			t1 = 0;
			pV = pV_mem;
			for( j=m_u32L-1; 0!=j; --j )
			{
				t1.Swap( tmp.Mul( t1 += *(--pV), t2 ) %= m_HugePow10 );
			}
			m5L.Swap( tmp.Mul( t1 += *(--pV), m5L )) %= m_HugePow10;
		}

		t2 /= 5;
		if ( p_end > ++pHead )
		{
			t2 += vHuge5[ *pHead ];
		}

		++u32Step;
	}
	ModPow5( m5L );
}

void Calc_Short( VOID ){
	UINT32 u32N = GetVal( m_HugeN );
	m_HugeT.Factorial( u32N );

	UINT32 u32Count5 = 0;
	while( 0 != u32N )
	{
		u32Count5 += ( u32N /= 5 );
	}

//	m_HugeZ = u32Count5;
	m_HugeT.RShift( u32Count5 ) %= m_HugePow10;
	printf( m_HugeT.ConvertToStr() );
	printf("\n");
}

void Calc(){
	if( m_HugeN < m_u32L*4 )
	{
		Calc_Short();
		return;
	}

	CHugeInt m5L;
	Get_m5L( m5L );  //m5L = final_result % (5^m_u32L)

	CHugeInt HugeExp( m_HugeZ );
	HugeExp += m_u32L;

	CHugeInt twomL( 2 );

	//	mpz_powm(twomL,twomL,HugeExp,m_HugePow5);// (twomL = 2^(-HugeExp) mod (5^HugeExp))
	if ( m_HugePow5 < 0x10000000 )
	{
		const UINT32 u32Mod = GetVal( m_HugePow5 );
		const UINT32 u32Phi = u32Mod / 5 * 4;
		twomL = HugeCalc::PowMod( 2, u32Phi - (UINT32)(HugeExp % u32Phi), u32Mod );
	}
	else
	{
		CHugeInt HugePhi( m_HugePow5 );
		( HugePhi /= 5 ) *= 4;	// Euler
		HugeExp %= HugePhi;
		if ( !HugeExp )
		{
			twomL = 1;
		}
		else
		{
			HugeExp.Negate() += HugePhi;
			twomL.PowMod( 2, HugeExp, m_HugePow10 );
		}
	}

	ModPow5( twomL *= m5L );			//twomL = (m5L*2^(-m_u32L)) % (5^m_u32L)
 	m_HugeT.Swap( twomL ) <<= m_u32L;	//2^m_u32L * (m5L*2^(-m_u32L))%(5^m_u32L), the final m_HugeT
	printf( m_HugeT.ConvertToStr() );
	printf("\n");
}

int
main(void)
{
	printf("Call %s", HugeCalc::GetVersion());
	if ( HC_LICENSE_ALL == HugeCalc::GetLicenseLevel() )
	{
		printf("\n");
	}
	else
	{
		printf(" (Unregistered)\n");
	}

	if ( HC_LICENSE_NONE == HugeCalc::GetLicenseLevel() )
	{
		printf( "\r\n您未通过 HugeCalc 的认证许可！\r\n\r\n" \
				"解决方案可选下列方案之一：\r\n" \
				"    一、请移动和[或]修改文件名为：/CopyrightByGuoXianqiang/.../HugeCalc.exe；\r\n" \
				"    二、请至 HugeCalc.chm 相关页面进行注册（一劳永逸）。\r\n\r\n" );

		system("pause");
		return -1;
	}

	while (TRUE)
	{
		printf("\nPlease input large number N, so that we could process for N!\n(if N is 0 then exit) N=");
		char buffer[1032] = { 0 };
		scanf("%s",buffer);
		if ( 0 > (m_HugeN = buffer).GetSign())
		{
			m_HugeN.Negate();
		}

		if ( !m_HugeN )
			return 0;

		if( m_HugeN.GetDigits()>1024 )
		{
			printf("Current only process number whose length no more than 1024\n");
			return -1;
		}

		printf("Please input number of digits to get (no more than 256):");
		scanf("%d",&m_u32L);
		if(m_u32L>256){
			printf("Input Out of range\n");
			return -1;
		}
		if(m_u32L<4)m_u32L=4;
		m_u32Pow2 = ( 32 > m_u32L ) ? ( 1UL << m_u32L ) : 0;

		( m_HugePow2 = 1 )<<= m_u32L;
		( m_HugePow5 = 5 ).Pow( m_u32L );
		( m_HugePow10 = 1 ).LShift( m_u32L );

		printf("Calcuate last %d non-zero digits of %s\n", m_u32L, m_HugeN.ConvertToStr());
		
		HugeCalc::EnableTimer();
		HugeCalc::ResetTimer();

		Calc();
		
		#if 1
			printf("timer: %s s\n", HugeCalc::ShowTimer( TRUE, FALSE ) );
		#else
			printf("timer: %s\n", HugeCalc::ShowTimer( TRUE, TRUE ) );
		#endif

		m_HugeN = 0;
		m_HugeT = 0;
		m_HugeZ = 0;
		m_HugePow2 = 1;
		m_HugePow5 = 5;
		m_HugePow10 = 1;
	}

	system("pause");
	return 0;
}