COME ON CODE ON

A blog about programming and more programming.

Posts Tagged ‘factorial

Combination

with 16 comments

In mathematics a combination is a way of selecting several things out of a larger group, where (unlike permutations) order does not matter. More formally a k-combination of a set S is a subset of k distinct elements of S. If the set has n elements the number of k-combinations is equal to the binomial coefficient. In this post we will see different methods to calculate the binomial.

1. Using Factorials
We can calculate nCr directly using the factorials.
nCr = n! / (r! * (n-r)!)

#include<iostream>
using namespace std;

long long C(int n, int r)
{
    long long f[n + 1];
    f[0]=1;
    for (int i=1;i<=n;i++)
        f[i]=i*f[i-1];
    return f[n]/f[r]/f[n-r];
}

int main()
{
	int n,r,m;
	while (~scanf("%d%d",&n,&r))
	{
		printf("%lld\n",C(n, min(r,n-r)));
	}
}

But this will work for only factorial below 20 in C++. For larger factorials you can either write big factorial library or use a language like Python. The time complexity is O(n).
If we have to calcuate nCr mod p(where p is a prime), we can calculate factorial mod p and then use modular inverse to find nCr mod p. If we have to find nCr mod m(where m is not prime), we can factorize m into primes and then use Chinese Remainder Theorem(CRT) to find nCr mod m.

#include<iostream>
using namespace std;
#include<vector>

/* This function calculates (a^b)%MOD */
long long pow(int a, int b, int MOD)
{
	long long x=1,y=a; 
	while(b > 0)
	{
		if(b%2 == 1)
		{
			x=(x*y);
			if(x>MOD) x%=MOD;
		}
		y = (y*y);
		if(y>MOD) y%=MOD; 
		b /= 2;
	}
	return x;
}

/* 	Modular Multiplicative Inverse
	Using Euler's Theorem
	a^(phi(m)) = 1 (mod m)
	a^(-1) = a^(m-2) (mod m) */
long long InverseEuler(int n, int MOD)
{
	return pow(n,MOD-2,MOD);
}

long long C(int n, int r, int MOD)
{
	vector<long long> f(n + 1,1);
	for (int i=2; i<=n;i++)
		f[i]= (f[i-1]*i) % MOD;
	return (f[n]*((InverseEuler(f[r], MOD) * InverseEuler(f[n-r], MOD)) % MOD)) % MOD;
}

int main()
{    
	int n,r,p;
	while (~scanf("%d%d%d",&n,&r,&p))
	{
		printf("%lld\n",C(n,r,p));
	}
}

2. Using Recurrence Relation for nCr
The recurrence relation for nCr is C(i,k) = C(i-1,k-1) + C(i-1,k). Thus we can calculate nCr in time complexity O(n*r) and space complexity O(n*r).

#include<iostream>
using namespace std;
#include<vector>

/*
	C(n,r) mod m
	Using recurrence:
	C(i,k) = C(i-1,k-1) + C(i-1,k)
	Time Complexity: O(n*r)
	Space Complexity: O(n*r)
*/

long long C(int n, int r, int MOD)
{
	vector< vector<long long> > C(n+1,vector<long long> (r+1,0));

	for (int i=0; i<=n; i++)
	{
		for (int k=0; k<=r && k<=i; k++)
			if (k==0 || k==i)
				C[i][k] = 1;
			else
				C[i][k] = (C[i-1][k-1] + C[i-1][k])%MOD;
	}
	return C[n][r];
}
int main()
{
	int n,r,m;
	while (~scanf("%d%d%d",&n,&r,&m))
	{
		printf("%lld\n",C(n, r, m));
	}
}

We can easily reduce the space complexity of the above solution by just keeping track of the previous row as we don’t need the rest rows.

#include<iostream>
using namespace std;
#include<vector>

/*
	Time Complexity: O(n*r)
	Space Complexity: O(r)
*/
long long C(int n, int r, int MOD)
{
	vector< vector<long long> > C(2,vector<long long> (r+1,0));

	for (int i=0; i<=n; i++)
	{
		for (int k=0; k<=r && k<=i; k++)
			if (k==0 || k==i)
				C[i&1][k] = 1;
			else
				C[i&1][k] = (C[(i-1)&1][k-1] + C[(i-1)&1][k])%MOD;
	}
	return C[n&1][r];
}

int main()
{
	int n,r,m,i,k;
	while (~scanf("%d%d%d",&n,&r,&m))
	{
		printf("%lld\n",C(n, r, m));
	}
}

3. Using expansion of nCr
Since

C(n,k) = n!/((n-k)!k!)
         [n(n-1)...(n-k+1)][(n-k)...(1)]
       = -------------------------------
           [(n-k)...(1)][k(k-1)...(1)]


We can cancel the terms: [(n-k)…(1)] as they appear both on top and bottom, leaving:

    n (n-1)     (n-k+1)
    - ----- ... -------
    k (k-1)       (1)

which we might write as:

C(n,k) = 1,                 if k = 0
       = (n/k)*C(n-1, k-1), otherwise

#include<iostream>
using namespace std;

long long C(int n, int r)
{
	if (r==0) return 1;
    else return C(n-1,r-1) * n / r;
}

int main()
{
	int n,r,m;
	while (~scanf("%d%d",&n,&r))
	{
		printf("%lld\n",C(n, min(r,n-r)));
	}
}

4. Using Matrix Multiplication
In the last post we learned how to use Fast Matrix Multiplication to calculate functions having linear equations in logarithmic time. Here we have the recurrence relation C(i,k) = C(i-1,k-1) + C(i-1,k).
If we take k=3 we can write,
C(i-1,1) + C(i-1,0) = C(i,1)
C(i-1,2) + C(i-1,1) = C(i,2)
C(i-1,3) + C(i-1,2) = C(i,3)
Now on the left side we have four variables C(i-1,0), C(i-1,1), C(i-1,2) and C(i-1,3).
On the right side we have three variables C(i,1), C(i,2) and C(i,3).
We need those two sets to be the same, except that the right side index numbers should be one higher than the left side index numbers. So we add C(i,0) on the right side. NOw let’s get our all important Matrix.

(.   .   .   .)  ( C(i-1,0) )   ( C(i,0) )
(.   .   .   .)  ( C(i-1,1) ) = ( C(i,1) )
(.   .   .   .)  ( C(i-1,2) )   ( C(i,2) )
(.   .   .   .)  ( C(i-1,3) )   ( C(i,3) )

The last three rows are trivial and can be filled from the recurrence equations above.

(.   .   .   .)  ( C(i-1,0) )   ( C(i,0) )
(1   1   0   0)  ( C(i-1,1) ) = ( C(i,1) )
(0   1   1   0)  ( C(i-1,2) )   ( C(i,2) )
(0   0   1   1)  ( C(i-1,3) )   ( C(i,3) )

The first row, for C(i,0), depends on what is supposed to happen when k = 0. We know that C(i,0) = 1 for all i when k=0. So the matrix reduces to

(.   .   .   .)  ( C(i-1,0) )   ( C(i,0) )
(1   1   0   0)  ( C(i-1,1) ) = ( C(i,1) )
(0   1   1   0)  ( C(i-1,2) )   ( C(i,2) )
(0   0   1   1)  ( C(i-1,3) )   ( C(i,3) )

And this then leads to the general form:

               i
(.   .   .   .)  ( C(0,0) )   ( C(i,0) )
(1   1   0   0)  ( C(0,1) ) = ( C(i,1) )
(0   1   1   0)  ( C(0,2) )   ( C(i,2) )
(0   0   1   1)  ( C(0,3) )   ( C(i,3) )

For example if we wan’t C(4,3) we just raise the above matrix to the 4th power.

              4
(1   0   0   0)  ( 1 )   ( 1 )
(1   1   0   0)  ( 0 ) = ( 4 )
(0   1   1   0)  ( 0 )   ( 6 )
(0   0   1   1)  ( 0 )   ( 4 )

Here’s a C++ code.

#include<iostream>
using namespace std;

/*    
	C(n,r) mod m
	Using Matrix Exponentiation
	Time Complexity: O((r^3)*log(n))
	Space Complexity: O(r*r)
*/

long long MOD;

template< class T >
class Matrix
{
	public:
		int m,n;
		T *data;
		
		Matrix( int m, int n );
		Matrix( const Matrix< T > &matrix );
		
		const Matrix< T > &operator=( const Matrix< T > &A );
		const Matrix< T > operator*( const Matrix< T > &A );
		const Matrix< T > operator^( int P );
		
		~Matrix();
};

template< class T >
Matrix< T >::Matrix( int m, int n )
{
	this->m = m;
	this->n = n;
	data = new T[m*n];
}

template< class T >
Matrix< T >::Matrix( const Matrix< T > &A )
{
	this->m = A.m;
	this->n = A.n;
	data = new T[m*n];
	for( int i = 0; i < m * n; i++ )
		data[i] = A.data[i];
}

template< class T >
Matrix< T >::~Matrix()
{
	delete [] data;
}

template< class T >
const Matrix< T > &Matrix< T >::operator=( const Matrix< T > &A )
{
	if( &A != this )
	{
		delete [] data;
		m = A.m;
		n = A.n;
		data = new T[m*n];
		for( int i = 0; i < m * n; i++ )
			data[i] = A.data[i];
	}
	return *this;
}

template< class T >
const Matrix< T > Matrix< T >::operator*( const Matrix< T > &A )
{
	Matrix C( m, A.n );
	for( int i = 0; i < m; ++i )
		for( int j = 0; j < A.n; ++j )
		{
			C.data[i*C.n+j]=0;
			for( int k = 0; k < n; ++k )
				C.data[i*C.n+j] = (C.data[i*C.n+j] + (data[i*n+k]*A.data[k*A.n+j])%MOD)%MOD;		
		}
	return C;
}

template< class T >
const Matrix< T > Matrix< T >::operator^( int P )
{
	if( P == 1 ) return (*this);
	if( P & 1 ) return (*this) * ((*this) ^ (P-1));
	Matrix B = (*this) ^ (P/2);
	return B*B;
}

long long C(int n, int r)
{
	Matrix<long long> M(r+1,r+1);
	for (int i=0;i<(r+1)*(r+1);i++)
		M.data[i]=0;
	M.data[0]=1;
	for (int i=1;i<r+1;i++)
	{
		M.data[i*(r+1)+i-1]=1;
		M.data[i*(r+1)+i]=1;
	}
	return (M^n).data[r*(r+1)];
}

int main()
{
	int n,r;
	while (~scanf("%d%d%lld",&n,&r,&MOD))
	{
		printf("%lld\n",C(n, r));
	}
}

5. Using the power of prime p in n factorial
The power of prime p in n factorial is given by
εp = ⌊n/p⌋ + ⌊n/p2⌋ + ⌊n/p3⌋…
If we call the power of p in n factorial, the power of p in nCr is given by
e = countFact(n,i) – countFact(r,i) – countFact(n-r,i)
To get the result we multiply p^e for all p less than n.

#include<iostream>
using namespace std;
#include<vector>

/* This function calculates power of p in n! */
int countFact(int n, int p)
{
	int k=0;
	while (n>0)
	{
		k+=n/p;
		n/=p;
	}
	return k;
}

/* This function calculates (a^b)%MOD */
long long pow(int a, int b, int MOD)
{
	long long x=1,y=a; 
	while(b > 0)
	{
		if(b%2 == 1)
		{
			x=(x*y);
			if(x>MOD) x%=MOD;
		}
		y = (y*y);
		if(y>MOD) y%=MOD; 
		b /= 2;
	}
	return x;
}

long long C(int n, int r, int MOD)
{
	long long res = 1;
	vector<bool> isPrime(n+1,1);
	for (int i=2; i<=n; i++)
		if (isPrime[i])
		{
			for (int j=2*i; j<=n; j+=i)
				isPrime[j]=0;
			int k = countFact(n,i) - countFact(r,i) - countFact(n-r,i);  
			res = (res * pow(i, k, MOD)) % MOD;
		}
	return res;
}

int main()
{
	int n,r,m;
	while (scanf("%d%d%d",&n,&r,&m))
	{
		printf("%lld\n",C(n,r,m));
	}
}

6. Using Lucas Theorem
For non-negative integers m and n and a prime p, the following congruence relation holds:
\binom{m}{n}\equiv\prod_{i=0}^k\binom{m_i}{n_i}\pmod p,
where
m=m_kp^k+m_{k-1}p^{k-1}+\cdots +m_1p+m_0,
and
n=n_kp^k+n_{k-1}p^{k-1}+\cdots +n_1p+n_0
are the base p expansions of m and n respectively.
We only need to calculate nCr only for small numbers (less than equal to p) using any of the above methods.

#include<iostream>
using namespace std;
#include<vector>

long long SmallC(int n, int r, int MOD)
{
	vector< vector<long long> > C(2,vector<long long> (r+1,0));

	for (int i=0; i<=n; i++)
	{
		for (int k=0; k<=r && k<=i; k++)
			if (k==0 || k==i)
				C[i&1][k] = 1;
			else
				C[i&1][k] = (C[(i-1)&1][k-1] + C[(i-1)&1][k])%MOD;
	}
	return C[n&1][r];
}

long long Lucas(int n, int m, int p)
{
	if (n==0 && m==0) return 1;
	int ni = n % p;
	int mi = m % p;
	if (mi>ni) return 0;
	return Lucas(n/p, m/p, p) * SmallC(ni, mi, p);
}

long long C(int n, int r, int MOD)
{
	return Lucas(n, r, MOD);
}

int main()
{
    
    int n,r,p;
	while (~scanf("%d%d%d",&n,&r,&p))
	{
          printf("%lld\n",C(n,r,p));
    }
}

7. Using special n! mod p
We will calculate n factorial mod p and similarly inverse of r! mod p and (n-r)! mod p and multiply to find the result. But while calculating factorial mod p we remove all the multiples of p and write
n!* mod p = 1 * 2 * … * (p-1) * 1 * 2 * … * (p-1) * 2 * 1 * 2 * … * n.
We took the usual factorial, but excluded all factors of p (1 instead of p, 2 instead of 2p, and so on). Lets call this strange factorial.
So strange factorial is really several blocks of construction:
1 * 2 * 3 * … * (p-1) * i
where i is a 1-indexed index of block taken again without factors p.
The last block could be not full. More precisely, there will be floor(n/p) full blocks and some tail (its result we can compute easily, in O(P)).
The result in each block is multiplication 1 * 2 * … * (p-1), which is common to all blocks, and multiplication of all strange indices i from 1 to floor(n/p).
But multiplication of all strange indices is really a strange factorial again, so we can compute it recursively. Note, that in recursive calls n reduces exponentially, so this is rather fast algorithm.
Here’s the algorithm to calculate strange factorial.

int factMOD(int n, int MOD)
{
	long long res = 1;
	while (n > 1)
	{
		long long cur = 1;
		for (int i=2; i<MOD; ++i)
			cur = (cur * i) % MOD;
		res = (res * powmod (cur, n/MOD, MOD)) % MOD;
		for (int i=2; i<=n%MOD; ++i)
			res = (res * i) % MOD;
		n /= MOD;
	}
	return int (res % MOD);
}

But we can still reduce our complexity.
By Wilson’s Theorem, we know (n-1)!\ \equiv\ -1 \pmod n for all primes n. SO our method reduces to:

long long factMOD(int n, int MOD)
{
	long long res = 1;
	while (n > 1)
	{
		res = (res * pow(MOD - 1, n/MOD, MOD)) % MOD;
		for (int i=2, j=n%MOD; i<=j; i++)
			res = (res * i) % MOD;
		n/=MOD;
	}
	return res;
}

Now in the above code we are calculating (-1)^(n/p). If (n/p) is even what we are multiplying by 1, so we can skip that. We only need to consider the case when (n/p) is odd, in which case we are multiplying result by (-1)%MOD, which ultimately is equal to MOD-res. SO our method again reduces to:

long long factMOD(int n, int MOD)
{
	long long res = 1; 
	while (n > 0)
	{
		for (int i=2, m=n%MOD; i<=m; i++)
			res = (res * i) % MOD;
		if ((n/=MOD)%2 > 0) 
			res = MOD - res;
	}
	return res;
}

Finally the complete code here:

#include<iostream>
using namespace std;
#include<vector>

/* This function calculates power of p in n! */
int countFact(int n, int p)
{
	int k=0;
	while (n>=p)
	{
		k+=n/p;
		n/=p;
	}
	return k;
}

/* This function calculates (a^b)%MOD */
long long pow(int a, int b, int MOD)
{
	long long x=1,y=a; 
	while(b > 0)
	{
		if(b%2 == 1)
		{
			x=(x*y);
			if(x>MOD) x%=MOD;
		}
		y = (y*y);
		if(y>MOD) y%=MOD; 
		b /= 2;
	}
	return x;
}

/* 	Modular Multiplicative Inverse
	Using Euler's Theorem
	a^(phi(m)) = 1 (mod m)
	a^(-1) = a^(m-2) (mod m) */
long long InverseEuler(int n, int MOD)
{
	return pow(n,MOD-2,MOD);
}

long long factMOD(int n, int MOD)
{
	long long res = 1; 
	while (n > 0)
	{
		for (int i=2, m=n%MOD; i<=m; i++)
			res = (res * i) % MOD;
		if ((n/=MOD)%2 > 0) 
			res = MOD - res;
	}
	return res;
}

long long C(int n, int r, int MOD)
{
	if (countFact(n, MOD) > countFact(r, MOD) + countFact(n-r, MOD))
		return 0;

	return (factMOD(n, MOD) *
			((InverseEuler(factMOD(r, MOD), MOD) * 
			InverseEuler(factMOD(n-r, MOD), MOD)) % MOD)) % MOD;
}

int main()
{    
	int n,r,p;
	while (~scanf("%d%d%d",&n,&r,&p))
	{
		printf("%lld\n",C(n,r,p));
	}
}

-fR0D

Written by fR0DDY

July 31, 2011 at 5:30 PM

Number of zeores and digits in N Factorial in Base B

with 2 comments

Question : What is the number of trailing zeroes and the number of digits in N factorial in Base B.

For example,
20! can be written as 2432902008176640000 in decimal number system while it is equal to “207033167620255000000” in octal number system and “21C3677C82B40000” in hexadecimal number system. That means that 10 factorial has 4 trailing zeroes in Base 10 while it has 6 trailing zeroes in Base 8 and 4 trailing zeroes in Base 16. Also 10 factorial has 19 digits in Base 10 while it has 21 digits in Base 8 and 16 digits in Base 16. Now the question remains how to find it?

Now we can break the Base B as a product of primes :
B = ap1 * bp2 * cp3 * …
Then the number of trailing zeroes in N factorial in Base B is given by the formulae
min{1/p1(n/a + n/(a*a) + ….), 1/p2(n/b + n/(b*b) + ..), ….}.

And the number of digits in N factorial is :
floor (ln(n!)/ln(B) + 1)

Here’s a sample C++ code :

#include<iostream>
using namespace std;
#include<math.h>

int main()
{
    int N,B,i,j,p,c,noz,k;
    while (scanf("%d%d",&N,&B)!=EOF)
    {
          noz=N;
          j=B;
          for (i=2;i<=B;i++)
          {
              if (j%i==0)
              {   
                 p=0;
                 while (j%i==0)
                 {
                       p++;
                       j/=i;
                 }
                 c=0;
                 k=N;
                 while (k/i>0)
                 {
                       c+=k/i;
                       k/=i;
                 }
                 noz=min(noz,c/p);
              }
          }
          double ans=0;
          for (i=1;i<=N;i++)
          {
              ans+=log(i);
          }
          ans/=log(B);ans+=1.0;
          ans=floor(ans);
          printf("%d %.0lf\n",noz,ans);
    }
}

NJOY!
-fR0DDY

Written by fR0DDY

February 17, 2010 at 5:47 PM

Posted in Maths, Programming

Tagged with , , , , , ,

Last Non-zero Digit of Factorial

with 26 comments

If you are not familiar with factorial read this.

Question: Find the last non-zero digit of n!

Approach  1
If you think you can calculate the factorial first and then divide out all the zeroes, you can but only to a certain extent, i mean to the point you can calculate n factorial. What if you had to find the last non-zero digit of a very large number whose factorial you cannot calculate.

Approach  2
This is the most popular technique and quite handy too. We know that  n! = 1 x 2 x 3 …x (n-1) x n
or n! = (2^a2)(3^a3)(5^a5)(7^a7)…
or n! = (2^a2)(5^a5)(3^a3)(7^a7)…

Now we know that the number of zeroes in n1 is equal to a5. So n! divided by 10^a5 will be equal to factorial without trailing digits. Lets call this n! without trailing zeroes as (n!)’ .So,
(n!)’ = (2^(a2-a5))(3^a3)(7^a7)…

Let L(k) denote the least significant non-zero decimal digit of the integer k. So from above we can infer that
L(n!) = L((n!)’) = L[(2^(a2-a5))(3^a3)(7^a7)…]
= L[(3^a3)(7^a7)…]*L[(2^(a2-a5))]

This logic is implemented in the C/C++ function below.

int last_digit_factorial(int N)
{
    int i,j,ans=1,a2=0,a5=0,a;

    for (i = 1; i <= N; i++)
    {
        j = i;
        //divide i by 2 and 5
        while (j%2==0)
        {
              j /= 2;
              a2++;
        }
        while (j%5==0)
        {
              j/=5;
              a5++;
        }

        ans = (ans*(j%10))%10;
    }

    a=a2-a5;

    for (i = 1; i <= a; i++)
        ans=(ans * 2) %10;

    return ans;
}

Approach 3
Fact :The powers of 2 3 7 and 1 have cyclic last digit.
So we divide all the primes into one of these groups and then take out the power accordingly.

power1[4]={1, 1, 1, 1}
power2[4]={6, 2, 4, 8}
power3[4]={1, 3, 9, 7}
power7[4]={1, 7, 9, 3}
power9[4]={1, 9, 1, 9}

So for example, if we have to find last digit in 6!
6! = 1 * 2 * 3 * 4 * 5 * 6
6! = 2^3 * (3^2) * 10
6!’ = 2^3 * (3^2)

We know the power of 3 has cyclic last digit, so 3^2 last digit is 9 or power3[2 mod 4]
L(6!) = L(6!’) = power2[3 mod 4] * power3[2 mod 4] mod 10
= 32 mod 10 = 2

But why code this method when we have a simpler approach.

Approach 4

Theory : Since we know that only a factor 5 brings in zero and we always have more powers of 2 than of 5, so the value of L(n!) is easily computed recursively as L(L(n)L((n-1)!)) unless n is a multiple of 5, in which case we need more information. As a result, the values of L(n!) come in fixed strings of five, as shown below :

1 5 10 15 20 25
L(n!) 1264 22428 88682 88682 44846 44846

Thus if we know any value of 5n we can get the value of 5n+j. For example if we know value of L(15!) we know the value of L(16!), L(17!) up-to L(19!). But how to get the value of 5n. If we map the starting values of the L(5n!) we get like 2 8 8 4 4 8 2 2 6. They come in another fixed strings of five. So we need to know L(25n!) to know L((25n+5j)!) for j = 0,1,2,3,4. Continuing in this way if we tabulate the values of L((5^t n)!) we get :

L(n!) 1264 22428 88682 88682 44846 44846
L(5n!) 2884 48226 24668 48226 48226 86442
L(25n!) 4244 82622 82622 28488 46866 64244
L(125n!) 8824 68824 26648 68824 42286 26648
L(625n!) 1264 22428 88682 88682 44846 44846

So we see that the pattern for L(625n!) is same as L(n!). Hence we can conclude that the pattern for L((5^j n)!) is the same as for L((5^(j+4) n)!). So we can use a modulo 4 function to get the next results. Also we can see below that each of the four row have a starting digit as 0 , 2 , 4, 6 and 8.

Result :

k mod 4
0 06264 22428 44846 66264 88682
1 02884 24668 48226 62884 86442
2 04244 28488 46866 64244 82622
3 08824 26648 42286 68824 84462

This is all we need to get the last digit of n! or L(n!). First we convert the given number n into base 5. So we have
n = d_0 + d_1*5 + d_2*5^2 + … + d_h*5^h
where d_0 is the least significant digit.

Now we enter the above table at the row h (mod 4) in the block whose first digit is 0 (because the coefficient of 5^(h+1) is zero), and determine the digit in the (d_h)th position of this block. Let this digit be denoted by s_h. Then we enter the table at row h-1 (mod 4) in the block that begins with s_h, and determine the digit in the (d_(h-1))th position of this block. Let this digit be denote by s_(h-1). We continue in this way down to s_0, which is the least significant non-zero digit of n!.

Example :
To illustrate, consider the case of the decimal number n=1592. In the base 5 this is n=22332. Now we enter the above table at row k=4=0 (mod 4) in the block beginning with 0, which is 06264. The leading digit of n (in the base 5) is 2, so we check the digit in position 2 of this block to give L((2*5^4)!) = 2. Then we enter the table at row k=3 (mod 4) in the block beginning with 2, which is 26648, to find L((2*5^4 + 2*5^3)!) = 6.

Then in the row k=2 (mod 4), the block beginning with 6 is 64244, and we find L((2*5^4 + 2*5^3 + 3*5^2)!) = 4. From this we know we’re in the block 48226 in row k=1 (mod 4), so we have L((2*5^4 + 2*5^3 + 3*5^2 + 3*5)!) = 2. Finally, we enter the row k=0 (mod 4) in block 22428 to find the result

L(1592!) = L((2*5^4 + 2*5^3 + 3*5^2 + 3*5 + 2)!) = 4

Thus if there are k digits in the base 5 representation of n then we only need k lookups in the table to get the last non-zero digit. Finally a C++ program for the algorithm :

#include
using namespace std;

int A[4][5][5] = {
{0,6,2,6,4,2,2,4,2,8,4,4,8,4,6,6,6,2,6,4,8,8,6,8,2},
{0,2,8,8,4,2,4,6,6,8,4,8,2,2,6,6,2,8,8,4,8,6,4,4,2},
{0,4,2,4,4,2,8,4,8,8,4,6,8,6,6,6,4,2,4,4,8,2,6,2,2},
{0,8,8,2,4,2,6,6,4,8,4,2,2,8,6,6,8,8,2,4,8,4,4,6,2}
};

char num[200];
char* dto5(int n)
{
     int b=5;
     int j,l;
     register int i=0;
     do
     {
           j=n%b;
           num[i++]=(j<10) ? (j+'0'): ('A'+j-10);
     }while((n/=b)!=0);

     num[i]='';
     l=strlen(num);
     reverse(num,num+l);
     return num;
}

int last_digit_factorial(int N)
{
    char *num=dto5(N);
    int l = strlen(num),s=0,i;
    for (i=0;i<l;i++)
    {
        s=A[(l-1-i)%4][s/2][num[i]-'0'];
    }
    return s;
}

int main()
{
  int N;

  while (scanf("%d", &N) == 1)
  {
        printf("%d\n", last_digit_factorial(N));
  }
  return 0;
}

Can we make it even simpler. Yes.

Let’s write n as 5k + r

So,
(n)! = (5k + r)!
= 1·2·3…k.(k+1)…(5k-1)·(5k)·(5k+1)…(5k+r)
(Seperating out multiples of 5)
= {5·10·15…(5k)} · {1·2·3·4·6·7…(5k-4)·(5k-3)·(5k-2)·(5k-1)·(5k+1)…(5k+r)}
(Expanding multiples of 5)
= {5·(2·5)·(3·5)…(k·5)} · {1·2·3·4·6·7…(5k-4)·(5k-3)·(5k-2)·(5k-1)·(5k+1)…(5k+r)}
= 5^k {1·2·3…k} · {1·2·3·4·6·7…(5k-4)·(5k-3)·(5k-2)·(5k-1)·(5k+1)…(5k+r)}
= 5^k · k! · {1·2·3·4·6·7…(5k-4)·(5k-3)·(5k-2)·(5k-1)·(5k+1)…(5k+r)}
= 5^k · k! · product{(5·i+1)·(5·i+2)·(5·i+3)·(5·i+4), i = 0 to k – 1} · (5k+1)…(5k+r)
(Multiplying and dividing by 2^k)
= 5^k · 2^k · k! · product((5i+1)·(5i+2)·(5i+3)·(5i+4)/2, i = 0 to k – 1) · (5k+1)…(5k+r)
= 10^k · k! · product((5i+1)·(5i+2)·(5i+3)·(5i+4)/2, i = 0 to k – 1) · (5k+1)…(5k+r)

Now to find the last digit of n! or L(n), we remove the power of 10 first and move to next type.

So,
L(n) = L(5k+r)
L(n) = L(k! · product((5i+1)·(5i+2)·(5i+3)·(5i+4)/2, i = 0 to k – 1))
L(n) = [L(k) · {product((5i+1)·(5i+2)·(5i+3)·(5i+4)/2, i = 0 to k – 1) mod 10} · {(5k+1)…(5k+r) mod 10 }] mod 10

Now let us first find P = (5i+1)·(5i+2)·(5i+3)·(5i+4)/2 mod 10

To calculate this we will have to take (mod 2) value and (mod 5) value and then using chinese remainder take out (mod 10) value.

Clearly P mod 2 = 0.

Modular inverse of 2 (mod 5) is 3.

So,
P mod 5 = 1·2·3·4·3 mod 5
= 72 mod 5
= 2 mod 5

Using the chinese remainder theorem we get,

P mod 10 = 2 mod 10

Hence,

L(n) = [L(k) · {product((5i+1)·(5i+2)·(5i+3)·(5i+4)/2, i = 0 to k – 1) mod 10} · {(5k+1)…(5k+r) mod 10 }] mod 10
= [L(k) · {product(2, i = 0 to k – 1) mod 10} · L(r)] mod 10
= [L(k) · L(r) · {2^k mod 10}] mod 10
= 2^k · L(k) · L(r) mod 10

Also,
2^k mod 10 = 1 when k = 6 when k = 4i
= 2 when k = 4i+1
= 4 when k = 4i+2
= 8 when k = 4i+3

Hence,

L(0) = L(1) = 1
L(2) = 2
L(3) = 6
L(4) = 4

L(n) = L(5k+r) = 2^k·L(k)·L(r) (mod 10)

And here’s a simpler code.

#include<iostream>
using namespace std;
 
int P(int K)
{
        int A[]={6,2,4,8};
        if (K<1) return 1;
        return A[K%4];
}
 
int L(int N)
{
        int A[]={1,1,2,6,4};
        if (N<5) return A[N];
        return (P(N/5)*L(N/5)*L(N%5))%10;
}
 
int main()
{
        int N;
        while (scanf("%d", &N) == 1)
                printf("%d\n", L(N));
        return 0;
}

NJOY!!!
fR0D

Written by fR0DDY

June 20, 2009 at 6:12 PM

Posted in Programming

Tagged with , , , , , ,

Factorial

with 6 comments

In mathematics, the factorial of a non-negative integer n, denoted by n!, is the product of all positive integers less than or equal to n. For example,  4! = 1 x 2 x 3 x 4 = 24

In programming, factorial comes to use at many places. Sometimes you simply have to calculate factorials but in C/C++ you are limited by the range of integers. Even long long can store only up-to 20! . 

int main()
{
    long long i,j=1;
    for (i=1;i<21;i++)
    {
        j*=i;
        printf("%lld! = %lld\n",i,j);
    }
}&#91;/sourcecode&#93;

So you have to find factorials using strings and multiply the way you used to do in std I or II. If you haven't done so do now, you will learn a lot. Or if you know some language like Python or Haskell, which allow larger numbers you need not worry. For example,  a code to accept t number of test cases and then print the factorial of next t numbers(less than 100) in python would be :

&#91;sourcecode language='python'&#93;

f=&#91;1&#93;*101
for i in range(2,101):
    f&#91;i&#93;=f&#91;i-1&#93;*i

n=input()
for i in range(0,n):
    x=input()
    print f&#91;x&#93;&#91;/sourcecode&#93;

Sometimes the question can be to find n! as a power of primes. To know the power of prime in a factorial you can use the formula:
<em>                  ε</em><sub>p</sub> = ⌊n/p⌋ + ⌊n/<em>p</em><sup>2</sup>⌋ + ⌊n/<em>p</em><sup>3</sup>⌋...

where <em>ε</em><sub>p</sub> is the power of prime p in factorial n.

<strong>Number of trailing zeroes in n Factorial</strong>

Finding number of trailing zeroes is very easy. Because you can get a zero only through a multiple of 5. So the number of trailing zeroes is equal to the multiples of five which is :
<em>                         ε</em><sub>5</sub> = ⌊n/5⌋ + ⌊n/<em>5</em><sup>2</sup>⌋ + ⌊n/<em>5</em><sup>3</sup>⌋...

A simple c-sharp program to get t test cases and then output the number of trailing zeroes in the factorial of next t numbers would be :

using System;
 

class Program
{
    public static void Main()
    {
        int T, N, sum, x;
        T = int.Parse(Console.ReadLine());
        while (T>0)
        {
            N = int.Parse(Console.ReadLine());
            x = 5;
            sum = 0;
            while (N / x > 0)
            {
                sum += N / x;
                x *= 5;
            }
            Console.WriteLine(sum);
            T--;
        }
    }
}

I will be writing soon on how to find the last non-zero digit and if possible also last k non zero digits.

ciao
fR0D

Written by fR0DDY

June 20, 2009 at 12:39 PM

Posted in Programming

Tagged with , , , , , , , , ,

%d bloggers like this: