Miller Rabin Primality Test

Miller Rabin Primality Test is a probabilistic test to check whether a number is a prime or not. It relies on an equality or set of equalities that hold true for prime values, then checks whether or not they hold for a number that we want to test for primality.

Theory

1> Fermat’s little theorem states that if p is a prime and 1 ≤ a < p then
2> If p is a prime and or then or
3> If n is an odd prime then n-1 is an even number and can be written as . By Fermat’s Little Theorem either or for some 0 ≤ r ≤  s-1.
4> The Miller–Rabin primality test is based on the contrapositive of the above claim. That is, if we can find an a such that and for all 0 ≤ r ≤  s-1 then a is witness of compositeness of n and we can say n is not a prime. Otherwise, n may be a prime.
5> We test our number N for some random a and either declare that N is definitely a composite or probably a prime. The probably that a composite number is returned as prime after k itereations is

Algorithm

Input :A number N to be tested and a variable iteration-the number 
of 'a' for which algorithm will test N.
Output :0 if N is definitely a composite and 1 if N is probably a prime.

Write N as 
For each iteration
	Pick a random a in [1,N-1]
	x =  mod n
	if x =1 or x = n-1 
		Next iteration
	for r = 1 to s-1
		x  =  mod n
		if x = 1
			return false
		if x = N-1
			Next iteration
	return false
return true

Here’s a python implementation :

import random
def modulo(a,b,c):
        x = 1
        y = a
        while b>0:
                if b%2==1:
                        x = (x*y)%c
                y = (y*y)%c
                b = b/2
        return x%c
        
def millerRabin(N,iteration):
        if N<2:
                return False
        if N!=2 and N%2==0:
                return False
        
        d=N-1
        while d%2==0:
                d = d/2
        
        for i in range(iteration):
                a = random.randint(1, N-1)
                temp = d
                x = modulo(a,temp,N)
                while (temp!=N-1 and x!=1 and x!=N-1):
                        x = (x*x)%N
                        temp = temp*2
                
                if (x!=N-1 and temp%2==0):
                        return False
        
        return True

-fR0DDY

Factorial

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

Prime Numbers

Suppose that you have to find all prime numbers till say some value N or let’s say sum of all prime numbers till N. How would you proceed. Let’s play around a bit. A grade 5 student would tell you that a prime has only two divisors 1 and itself and would code it like this :

void basic1(int N)
{
     int i,j,k,c;
     long long s=0;
     for (i=1;i<=N;i++)
     {
         c=0;
         for (j=1;j<=i && c<3;j++)
             if (i%j==0)
                c++;
         if (c==2)
            s+=i;
     }
     cout<<s<<endl;
}

A little smarter kid would tell you that a prime number has no divisor between 2 and N/2 and would code it like this :

void basic2(int N)
{
     int i,j,k,F;
     long long s=0;
     for (i=2;i<=N;i++)
     {
         F=1;
         for (j=2;j<=i/2 && F==1;j++)
             if (i%j==0)
                F=0;
         if (F)
            s+=i;
     }
     cout<<s<<endl;
}

An even smarter kid will tell you that a prime number has no divisor between 2 and its root and would code it like this :

void moderate1(int N)
{
     vector<int> p(N/2,0);
     int i,j,k,F,c=0;
     long long s=0;
     for (i=2;i<=N;i++)
     {
         F=1;
         for (j=2;j*j<=i && F==1;j++)
             if (i%j==0)
                F=0;
         if (F)
         {
            p[c++]=i;
            s+=i;
         }
     }
     cout<<s<<endl;
}

A good maths student can tell you that a prime number has no prime divisors between 2 and its root.

void moderate2(int N)
{
     vector<int> p(N/2,0);
     p[0]=2;
     int i,j,F,c=1;
     long long s=2;
     for (i=3;i<=N;i+=2)
     {
         F=1;
         for (j=0; p[j]*p[j]<=i && F==1; j++)
             if (i%p[j] == 0)
                F=0;
         if (F)
         {
            p[c++]=i;
            s+=i;
         }
     }
     cout<<s<<endl;
}

A good programmer will tell you that Sieve of Eratosthenes is the best way to find the list of prime numbers.
Its algorithm is as follows :

1. Create a contiguous list of numbers from two to some highest number n.
2. Strike out from the list all multiples of two (4, 6, 8 etc.).
3. The list’s next number that has not been struck out is a prime number.
4. Strike out from the list all multiples of the number you identified in the previous step.
5. Repeat steps 3 and 4 until you reach a number that is greater than the square root of n (the highest number in the list).
6. All the remaining numbers in the list are prime.

And the C++ implementation would be :

void sieve(int N)
{
     int x=sqrt(N),i,j;
     vector<bool> S(N+1,0);
     for (i=4;i<=N;i+=2)
         S[i]=1;
     for (i=3;i<=x;i+=2)
         if (!S[i])
            for (j=i*i;j<=N;j+=2*i)
                S[j]=1;

     long long s=0;
     for (i=2;i<=N;i++)
         if (!S[i])
            s+=i;

     printf("%lld\n",s);
}

But then we can optimize even the Sieve of Eratosthenes. If you look closer you will realise that apart from 2 all the other even numbers are of no use and just waste our time and memory. So we should get rid of them. We can do that by sieving only odd numbers. Let an element with index i correspond to number 2*i+1. So in the sieve we only need to go upto (N-1)/2. Also if p = 2*i+1 p² = 4*i^2 + 4*i + 1 which is represented by index 2*i*(i+1). Also the step will be of order 2*i+1. So the code would look something like this:

void improved_sieve(int N)
{
     int M=(N-1)/2;
     int x=(floor(sqrt(N))-1)/2,i,j;
     vector<bool> S(M+1,0);
     for (i=1;i<=x;i++)
         if (!S[i])
            for (j=2*i*(i+1);j<=M;j+=(2*i+1))
                S[j]=1;

     long long s=2;
     for (i=1;i<=M;i++)
         if (!S[i])
            s+=(2*i+1);

     printf("%lld\n",s);
}

and since i am learning Python. Here’s the python code as well :

import math
def improved_sieve(N):
    M=(N-1)/2
    x=(int(math.sqrt(N))-1)/2
    S=[]
    for i in range(M+1):
        S.append(False);

    for i in range (1,x+1):
        if (S[i]==False):
            for j in range (2*i*(i+1),M+1,2*i+1):
                S[j]=True

    s=2;
    for i in range (1,M+1):
        if (S[i]==False):
            s+=(2*i+1)

    print s

N=input()
improved_sieve(N)

Lets compare all these codes on runtime :

Limit10^N where N=	basic1	basic2	moderate1	moderate2	sieve	improved_sieve
4	0.079	0.026	0.003	0.007	0.002	0.001
5	7.562	2.611	0.050	0.117	0.020	0.011
6	–	–	1.160	2.218	0.208	0.121
7	–	–	26.350	44.946	2.183	1.269

All time are in seconds. If you notice that though algorithm 4 seems to be better than 3 it takes more time because of the array refrences that are made several times. Choose the best.

Happy Coding!
-fR0D