#include <iostream>
#include <cassert>

using namespace std;


// sieve_fact
// Finds the prime powers of n!
//
// Inputs: an array of size n+1
// Outputs: the prime powers of n! in the array, or 0 for non-prime powers
// Invariant: product of i*factor[i]

void sieve_fact(int *factor, int n)
{
    // Set an array with each item

    for(int i=0; i<=n; ++i)
    {
        factor[i] = 1;
    }

    // A simple number sieve
    for(int p=2; p<=n; ++p)
    {
        if(factor[i])   // i is prime
        {
            for(int q=p*2; q<=n; q+=p)
            {
                if(factor[q])
                {
                    // We must factorize q
                    // q consists only of factors >= p

                    int r;
                    for(r=q; r%p==0; r/=p)
                        factor[p] += factor[q];

                    if(r>p)
                    {
                        assert(factor[r]!=0);
                        factor[r] += factor[q];
                    }
                    else
                        assert(r==1);
     
                    // We "tick off" q, setting it to zero
                    factor[q]=0;
                }
            }
        }
    }
}



int multiplications = 0;

int multiply(int a, int b)
{
    if(a==1) return b;
    if(b==1) return a;

    multiplications++;
    return (a*b); // % 10099;
}


int mypow(int a, int b)
{
    if(b==0) return 1;
    else if(b==1) return a;
    else if(b==2) return multiply(a,a);
    else if(b==3) return multiply(multiply(a,a),a);

    int p = 1;
    int m = a;

    for(; b; b=b>>1)
    {
        if(b&1)
            p = multiply(p, m);
        if(b>1)
            m = multiply(m, m);
    }

    return p;
}


// mult_factors1
// Multiplies an array of factors
//
// Input: an array of size n+1
// Output: product of i * factor[i]

int mult_factors1(int *factor, int n)
{
    int m=1;
    for(int i=2; i<=n; ++i)
    {
        if(factor[i])
        {
            m = multiply(m, mypow(i, factor[i]));
        }
    }
    return m;
}

// mult_factors2
// Multiplies an array of factors
//
// Input: an array of size n+1
// Output: product of i * factor[i]

int mult_factors2(int *factor, int n)
{
    if(n<2) return 1;   // Does not work for these numbers

    int fact=1;
    int m=2;
    int lastFactor=factor[2];

    for(int i=3; i<=n; ++i)
    {
        if(factor[i])
        {
            if(factor[i] != lastFactor)
            {
                fact = multiply(fact, mypow(m, lastFactor - factor[i]));
                lastFactor = factor[i];
            }
            m = multiply(m, i);
        }
    }

    assert(lastFactor==1);

    return multiply(fact, m);
}

template<class N>
N naivefac(N n) {
    N x = 1;
    while (n > 1) x*=n, n--;
    return x;
}


int fac(int n) {
    int x, y;

    x = (n & 1) ? n-- : 1;
    y = n;
    while (n) {
        x = multiply(x,y);
        y += (n -= 2);
    }
    return x;
}

template<class N>
N fac1(N n1, N n0=2)
{
    N m=1;

    for(; n0<n1; ++n0)
        m *= n0;

    return m;
}

template<class N>
N mult_range(N n0, N n1)
{
    N n=1;
    for(; n0<n1; ++n0)
        n *= n0;

    return n;
}

int apascal(int row, int col)
{
    if(col==0) return 1;
    if(col>row) return 0;

    return apascal(row-1, col) - apascal(row-1, col-1);
}

template<class N>
void derive_coefficients(N *c, int len)
{
    N *m = new N[len+1];
    int i, j;
    N n;

    for(i=0, n=1; i<=len; ++i, n+=len)
    {
        m[i] = mult_range(n, n+N(len));

        c[i] = 0;
        for(j=0; j<=i; ++j)
            c[i] += apascal(i,i-j)*m[j];
    }
    
    delete [] m;
}

// sfac
// Performs factorial by multiplying in steps of "step"
// The array of coefficients, "coeff", is overwritten
template<class N>
N sfac(N n, int step, N*coeff)
{
    N n0=n-n%step;    
    if(n0==0) n0=1;

    N f = 1;

    while(n>n0) f*=n, n--;    

    for(N n0=step; n0<=n; n0+=step)
    {
        f *= coeff[0];

        for(int i=0; i<step; ++i)
        {
            coeff[i] += coeff[i+1];
        }
    }

    return f;
}


template<class N>
N quickfac(N n)
{
    N coeffs[11] = { 3628800, 
        670438944000, 
        107686468915200, 
        2750919796800000, 
        25623563263200000,
        116882616060000000,
        297274516560000000,
        443961000000000000,
        387555840000000000,
        183254400000000000,
        36288000000000000 };

    return sfac(n, 10, coeffs);
}


void report(int n)
{
    int *factor = new int[n+1];
    sieve_fact(factor, n);
 
    cout << n << "! = ";

    if(n<2) cout << "1";

    for(int i=2; i<=n; ++i)
    {
        if(factor[i])
        {
            if(i>2) cout <<" * ";
            cout << i;
            if(factor[i]>1) cout << "^" << factor[i];
        }
    }

    multiplications = 0;
    cout << "\nNaive:       " << n << "! = " << naivefac(n);
    cout << "  multiplications = " << multiplications;

    multiplications=0;
    cout << "\nEricson:     " << n << "! = " << fac(n);
    cout << "  multiplications = " << multiplications;

    multiplications=0;
    cout << "\nAlgorithm 1: " << n << "! = " << mult_factors1(factor, n);
    cout << "  multiplications = " << multiplications;

    multiplications=0;
    cout << "\nAlgorithm 2: " << n << "! = " << mult_factors2(factor, n);
    cout << "  multiplications = " << multiplications;

    delete [] factor;
}


int mults = 0;
int adds = 0;

// TestNum is a class that performs modular arithmetic, and
// measures the number of additions and multiplications
// This is useful for testing algorithms that might otherwise overflow

template<class I=long, I M=10099>
class TestNum
{
    I n;
public:

    operator I() const { return n; }
    // TestNum(I m) : n(m%M) { }

    template<class N>
    TestNum(const N &m) : n(m%M) { }

    TestNum() { }
    TestNum<I,M> operator*(TestNum m)
    {
        mults++;
        return (n*m)%M;
    }

    TestNum<I,M> operator+(TestNum m)
    {
        adds++;
        return (n+m)%M;
    }

    TestNum<I,M> &operator+=(TestNum m)
    {
        adds++;
        n = (n + m)%M;
        return *this;
    }

    TestNum<I,M> &operator=(TestNum m)
    {
        n=m%M;
        return *this;
    }

    TestNum<I,M> &operator++()
    {
        adds++;
        n = (n+1)%M;
        return *this;
    }

    TestNum<I,M> &operator--()
    {
        adds++;
        n--;
        return *this;
    }

    TestNum<I,M> &operator*=(TestNum m)
    {
        mults++;
        n = (n*m)%M;
        return *this;
    }
};

typedef TestNum<__int64, __int64(1)<<60> Int64;
typedef TestNum<int, 10099> Mod32;

int main(int argc, char* argv[])
{
    if(argc!=2)
    {
        cout << "Usage: fact [n]\n";
        return 1;
    }

    int x2[3];
    derive_coefficients(x2, 2);
    for(int i=0; i<3; ++i)
        cout << "Coeff: " << x2[i] << endl;

    int x3[4];
    derive_coefficients(x3, 3);
    for(int i=0; i<4; ++i)
        cout << "Coeff: " << x3[i] << endl;

    cout << endl;

    const int N=10;
    Int64 xn[N+1];
    derive_coefficients(xn, N);
    for(int i=0; i<=N; ++i)
        cout << "Coeff: " << xn[i] << endl;

    for(int i=0; i<13; ++i)
    {
        Int64 c[N+1];
        memcpy(c, xn, sizeof(c));
        Int64 f = sfac<Int64>(i, N, c);
        cout << i << ": " << naivefac(Int64(i)) << "=" << f << "=" << quickfac(Int64(i)) << endl;
    }

    cout << "Quickfac: ";
    mults=0;
    adds=0;
    cout << "10000! = " << quickfac(Mod32(10000)) << endl;
    cout << "Adds = " << adds << ", mults=" << mults << endl;

    cout << "Naivefac: ";
    mults=0;
    adds=0;
    cout << "10000! = " << naivefac(Mod32(10000)) << endl;
    cout << "Adds = " << adds << ", mults=" << mults << endl;

return 0;



    report(atoi(argv[1]));
    return 0;
}