factorialtail.cpp

HugeCalc V5.1.0.1 ----> 这是一套绿色软件

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

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

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

#ifndef _UNICODE
	#pragma comment( lib, "../../HugeCalc_Dll_Import/Lib/HugeCalc.lib" )
#else
	#pragma comment( lib, "../../HugeCalc_Dll_Import/Lib/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 );

CHugeInt& ModPow10( CHugeInt &HugeInt )
{
	return HugeInt.ModPowTen( m_u32L );
}

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

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

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

	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] );
				}
			}
		}
	}

	for(i=0; m_u32L!=i; ++i)
	{
		ModPow10( dst[i] );
	}
}

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 )
		{
			ModPow10( *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) ));
				ModPow10( (*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] ));
			ModPow10( (*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 ); //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 );//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 ); //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 ); //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 );
		PolyMul( p1=ptmp4, p, p-m_u32L );
		for(j=0; m_u32L!=j; ++j,++p, ++p1)
//			F(i,j).Swap( tmp4[j] );
			ModPow5( (*p).Swap( *p1 ));
    }

	vHugeFun.swap( s_vHugeFun );
}

void Get_m5L( CHugeInt &m5L )
{
	static UINT32 s_u32L( 0 );
	static CHUGEINT_VECTOR vHugeFun;
	if ( s_u32L != m_u32L )
	{
		Calc_Fun( vHugeFun );
		s_u32L = m_u32L;
	}

	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;

	// 注意：如果未注册，所得到的结果 CarryParam 有可能被随机干扰，从而导致结果不正确！
	const CCarryParam CarryParam( 5, m_HugeN );
	(( m_HugeZ = m_HugeN ) -= (UINT32)CarryParam.SumOfDigits() ) /= 4;

	const U32_VECTOR * const pVector = CarryParam.GetLPVector();
	const U32_VECTOR::const_iterator p_begin = pVector->begin();
	const U32_VECTOR::const_iterator p_end = pVector->end();
	U32_VECTOR::const_iterator p = p_begin;

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

	U32_VECTOR::const_iterator pHead = ( p = p_begin ) + m_u32L;
	if ( p_end < pHead )
	{
		pHead = p_end;
	}

	const U32_VECTOR v32First( p_begin, pHead );
	const CCarryParam CarryParam0( 5, SIGN_POS, &v32First );
	CHugeInt t1, t2( CarryParam0 ), tmp;

	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;
	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 )
			{
				ModPow10( t1.Swap( tmp.Mul( t1 += *(--pV), t2 )));
			}
			ModPow10( m5L.Swap( tmp.Mul( t1 += *(--pV), m5L )));
		}

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

		++p;
		++u32Step;
	}
//	ModPow5( m5L );
}

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

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

//	m_HugeZ = u32Count5;
	ModPow10( m_HugeT.DecRShift( u32Count5 ));
}

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;

	// (twomL = 2^(-HugeExp) mod (5^HugeExp))
	const CHugeInt HugePhi(( CHugeInt( m_HugePow5 ) /= 5 ) *= 4 );	// Euler
	HugeExp %= HugePhi;

	CHugeInt twomL;
	if ( !HugeExp )
	{
		twomL = 1;
	}
	else if ( HugePhi.GetDigits() - HugeExp.GetDigits() <= ( m_HugePow5.GetDigits() / 2 ) )
	{
		twomL.PowMod( 2, HugeExp.Negate() += HugePhi, m_HugePow10 );
	}
	else
	{
		twomL.PowMod( ( ++CHugeInt( m_HugePow5 )) /= 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
}

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

	printf( "\nPlease input number of digits to get (no more than 255): " );
	scanf( "%d", &m_u32L );
	if ( m_u32L < 4 )
	{
		m_u32L = 4;
	}
	else if ( m_u32L > 255 )
	{
		m_u32L = 18;
	}

	m_u32Pow2 = ( 32 > m_u32L ) ? ( 1UL << m_u32L ) : 0;
	( m_HugePow2 = 1 )<<= m_u32L;
	( m_HugePow5 = 5 ).Pow( m_u32L );
	( m_HugePow10 = 1 ).DecLShift( m_u32L );
	m_HugeN = m_HugePow5;
	Get_m5L( m_HugeZ );

	printf("\nPlease input large number N, so that we could process for N!\nIf N is 0 then exit ...\n");

	for ( ; ; )
	{
		printf("\nN = ");
		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("Calcuate last %d non-zero digits of \"%s\":\n\t", m_u32L, m_HugeN.ShowStr( FS_BAND ));

		HugeCalc::EnableTimer( TRUE );
		HugeCalc::ResetTimer( 0 );

		Calc();

		HugeCalc::EnableTimer( FALSE );

		printf( m_HugeT.ShowStr( FS_NORMAL ) );

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

		m_HugeN = 0;
		m_HugeT = 0;
		m_HugeZ = 0;
	}

	system("pause");
	return 0;
}