# COME ON CODE ON

A blog about programming and more programming.

## 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 $a^{p-1}\equiv1\pmod{p}.$
2> If p is a prime and $x^2\equiv1\pmod{p}$ or $\left(x-1\right)\left(x+1\right)\equiv0\pmod{p}$ then $x\equiv1\pmod{p}$ or $x\equiv-1\pmod{p}.$
3> If n is an odd prime then n-1 is an even number and can be written as $2^{s}.d$. By Fermat’s Little Theorem either $a^{d}\equiv1\pmod{n}$ or $a^{2^r\cdot d}\equiv -1\pmod{n}$ 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 $a^{d}\not\equiv1\pmod{n}$ and $a^{2^r\cdot d}\not\equiv -1\pmod{n}$ 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 $4^{-k}.$

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 $2^{s}.d$
For each iteration
Pick a random a in [1,N-1]
x = $a^{d}$ mod n
if x =1 or x = n-1
Next iteration
for r = 1 to s-1
x  = $x^{2}$ 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

Written by fR0DDY

September 18, 2010 at 12:23 PM

## 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;
while (T>0)
{
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 , , , , , , , , ,

## Programming Multiplicative Functions

In number theory, a multiplicative function is an arithmetic function f(n) of the positive integer n with the property that f(1) = 1 and whenever a and b are coprime, then f(ab) = f(a)f(b).

This type of functions can be programmed in quick time using the multiplicative formmula. Examples of multiplicative functions include Euler’s totient function, Möbius function and divisor function.

A multiplicative function is completely determined by its values at the powers of prime numbers, a consequence of the fundamental theorem of arithmetic. Thus, if n is a product of powers of distinct primes, say n = pa qb …, then f(n) = f(pa) f(qb) …

Lets take an example. If we need to calculate the divisor of 12 then div(12) = div(2231) or div(22)div(31). If we are calculating it for a range then it becomes even more easier. We only need to calculte the div(22)part because the div(31) part will already be present in the array. So lets look at an sample code for calculating the divisor of all numbers from 1 to 10000000.

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int A[10000000]={0};
int main()
{

int isprime[3163],d,n,e,p,k,i;

for (n=2;n<3163;n++)
isprime[n]=1;

//Sieve for Eratosthenes for Prime
//Storing the smallest prime which divides n.
//If A[n]=0 it means it is prime number.
for(d=2;d<3163;d++)
{
if(isprime[d])
{
for (n=d*d;n<3163;n+=d)
{
isprime[n]=0;
A[n]=d;
}
for (;n<=10000000;n+=d)
A[n]=d;
}
}

//Applying the formula
//Divisor(N)=Divisors(N/p^f(N,p))*(f(N,p)+1)
A[1]=1;
for(n=2;n<=10000000;n++)
{
if (A[n]==0)
A[n]=2;
else
{
p=A[n],k=n/p,e=2;
while (k%p==0)
k/=p,e++;
A[n]=A[k]*e;
}
}
printf("time=%.3lf sec.\n",
(double) (clock())/CLOCKS_PER_SEC);
while (scanf("%d",&i),i)
{
printf("%d\n",A[i]);
}
return 0;
}

If you run the above program you will be able to see how fast have we made this program. Similar techniques can be used to calculate the divisor function or the sum of squares of divisors which is also called sigma2. The general formula for calculating sigma function is :

sigma_x

or

-fR0DDY

Written by fR0DDY

April 13, 2009 at 6:31 AM

Posted in Maths, Programming

## 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 p2 = 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

Written by fR0DDY

March 10, 2009 at 3:34 PM

Posted in Maths, Programming

Tagged with , , , ,