/*
	File:	nf.java
	Approximate n factorial with Stirling's formula.
	cf. CLRS p. 55 Ch. 3

	Stirling approximation as a double overflows at 171!

	Use a Theta( lg( n ) ) algorithm for integer exponentiation.

	yanushka
 */

import java.math.*;

class Test
{
    public static void main( String [] args )
    {
	int	n = Entry.getInt( "Enter positive integer for factorial." );
	Factorial	f;

/*
	for ( int i = 6; i <= n; i++ ) {
	    f = new Factorial( i );
	    f.display();
	}
 */
	while ( n > 0 ) {
	    f = new Factorial( n );
	    f.display();
	    n = Entry.getInt( "Enter positive integer for factorial." );
	}
    }
}

/*
	The class Factorial computes n! exactly inside BigInteger
	and with Stirling's approximation of sqrt( 2 Pi n ) ( n / e )^n.
	It uses an efficient integer power algorithm.

	Its attributes are an int n, a BigInteger object named nf for the 
	exact value of n factorial, a double named snf and a BigDecimal 
	named bdSnf for Stirling's approximations.

	Its methods are a constructor, compute(), computeStirling(), 
	power() and display().
 */
class Factorial
{
    private final static int 		OVERFLOW = 170;
    private final static BigDecimal	HUNDRED = new BigDecimal( "100" );

    private int		n;
    private BigInteger	nf;
    private double	snf;
    private BigDecimal	bdSnf;

    public Factorial( int nn )
    {
	n = nn;
	compute();
    }

/*
	Compute n! with a loop using methods of the BigInteger class.
 */
    public void compute()
    {
	long	start = System.currentTimeMillis();

	nf = new BigInteger( "1" );
	for ( int i = 2; i <= n; i++ )
	    nf = nf.multiply( new BigInteger( new Integer( i ).toString() ) );

	long	finis = System.currentTimeMillis();
	System.out.println( "\t" + n + "! with n BigInteger products:    " 
				+ ( finis - start ) + " ms." );

	start = System.currentTimeMillis();
	if ( n < OVERFLOW )
	    computeStirling();
	else
	    findStirling();
	finis = System.currentTimeMillis();
	System.out.println( "\t" + n + "! with Stirling's approximation: " 
				+ ( finis - start ) + " ms." );
    }

/*
	Compute n! as a double with Stirling's approximation.
 */
    private void computeStirling()
    {
	double	dn = ( double )n;

	snf = Math.sqrt( 2 * Math.PI * dn ) * power( dn / Math.E );
    }

/*
	Compute n! as a BigDecimal with Stirling's approximation.
 */
    private void findStirling()
    {
	double		dn = ( double )n;
	BigDecimal	rtBD = new BigDecimal( 
					Math.sqrt( 2.0 * Math.PI * dn ) );

	bdSnf = powerBD( dn / Math.E ).multiply( rtBD );
    }

/*
	Compute x^n with a Theta( lg( n ) ) algorithm as a double.

	Initialize ans to 1 and pwr2 to x.
	while n > 0
	    if n is odd then multiply ans by pwr2.
	    Halve n.
	    Square pwr2.
 */
    private double power( double x )
    {
	double	pwr2 = x,
		ans = 1.0;
	int 	nn = n;

/*
	Invariant:  x^n = ans * pwr2^nn and nn >= 0.
 */
	while ( nn > 0 ) {
	    if ( nn % 2 == 1 )
		ans *= pwr2;
	    nn /= 2;
	    pwr2 *= pwr2;
	}

	return ans;
    }

/*
	Compute x^n with a Theta( lg( n ) ) algorithm as a BigDecimal.
 */
    private BigDecimal powerBD( double x )
    {
	BigDecimal	pwr2 = new BigDecimal( x ),
			ans = new BigDecimal( 1.0 );
	int 		nn = n;

/*
	Invariant:  x^n = ans * pwr2^nn and nn >= 0.
 */
	while ( nn > 0 ) {
	    if ( nn % 2 == 1 )
		ans = ans.multiply( pwr2 );
	    nn /= 2;
	    if ( nn > 0 )
	        pwr2 = pwr2.multiply( pwr2 );
	}

	return ans;
    }

/*
	Display n, snf exp( 1 / 12( n + 1 ) ), n!, snf exp( 1 / (12 n) ),
	and % error of snf / n!.
 */
    public void display()
    {
	if ( n < OVERFLOW ) {
	    System.out.println( n + "  " + 
			snf * Math.exp( 1 / ( 12.0 * ( n + 1 ) ) ) + "  " +
			nf.toString() + "  " + 
			snf * Math.exp( 1 / ( 12.0 * n ) ) + "  " +
			snf );

	    double	dnf = nf.doubleValue();
	    System.out.println( "\t\t" + 
			100.0 * Math.abs( snf - dnf ) / dnf + "%" );
	} else {
	    BigInteger	biSnf = bdSnf.toBigInteger();
	    BigDecimal 	loNf = bdSnf.multiply( new BigDecimal(
				Math.exp( 1 / ( 12.0 * ( n + 1 ) ) ) ) ),
	    		hiNf = bdSnf.multiply( new BigDecimal( 
					Math.exp( 1 / ( 12.0 * n ) ) ) );

	    System.out.println( n + "\n" + 
			loNf.toBigInteger().toString() + "\n\n" 
			+ nf.toString() + "\n\n" + 
			hiNf.toBigInteger().toString() + "\n\n" 
			+ biSnf.toString() );

	    BigDecimal	bdNf = new BigDecimal( nf ),
			absErr = hiNf.subtract( bdSnf ),
			pcErr = absErr.multiply( HUNDRED ).divide( bdNf, 0 );
	    System.out.println( "digits in absolute error: " + 
			absErr.toBigInteger().toString().length() + "  " +
				"digits in n!: " +
				nf.toString().length() +
				"\n\n  percent error " +
				pcErr.doubleValue() + "%" );

	    String 	hiS = hiNf.toBigInteger().toString(),
			nfS = nf.toString();

	    int		same = 0;

	    while ( hiS.charAt( same ) == nfS.charAt( same ) )
		same++;

	    System.out.println( "\n\nsame digits at hi end: " + same );
	}
    }
}

/*
	Some results follow.

	n	time for exact n! (ms)	time with Stirling (ms)	digits 	%error
	100	  15			    1			  158   0.083
	150	  13			    0			  263   0.055
	170	  14			   80			  307	0.049
	200	  29			   98			  375   0.0416
	500	  51			  835			 1135	0.0166
	1000	  60			 2629			 2568	0.0083
	5000	1743			53831			16326	0.0016

	Note for n >= 170 arithmetic occurs in the BigDecimal class.
 */