## Prime number

Moderator: Board moderators

temper_3243
Experienced poster
Posts: 105
Joined: Wed May 25, 2005 7:23 am

### Prime number

I am trying to determine whether a number is prime or not. The numbers range from 1 to 2^64 -1. I am using data type unsigned long long in C.

I am trying it with Miller rabin test but it is failing .
Well i am changing the code. How do i prevent overflow in multiplication ?

Code: Select all

``````
#include<stdio.h>
int TRUE = 1;
int FALSE = 0;
int COMPOSITE = 1;
int PRIME = 0;
int NO_TIMES = 30;
int table[1500] = { 2, 3, 5, 7 }, count = 4;
int a[10000] = { 0 };
unsigned long long
mycomb (unsigned long long a, unsigned long long b, unsigned long long p)
{
int k;
a = a % p;
if (b == 0)
return 1;
if (b == 1)
return a;

k = mycomb ((a * a) % p, b / 2, p) % p;
k = k % p;
if (b % 2 == 1)
k = (k * a) % p;

return k % p;

}

int
fermat (unsigned long long m, unsigned long long u)
{
unsigned long long k;

k = mycomb (m, u - 1, u) % u;
if (k != 1)
return COMPOSITE;
return PRIME;
}

int
witness (unsigned long long a, unsigned long long n)
{
unsigned long long z = 0, u, j, prevx = 0, x;
u = n - 1;
while (u % 2 == 0)
{
u = u / 2;
z++;
}
x = mycomb (a, u, n);
x = x % n;

for (j = 1; j <= z; j++)
{
prevx = x;
x = (x * x) % n;
if (x == 1 && prevx != 1 && prevx != n - 1)
return TRUE;
}

if (x != 1)
return TRUE;

return FALSE;

}

int
miller_rabin (unsigned long long z)
{
int i;
unsigned long long l;
for (i = 0; i < count; i++)
{
l = table[i] % (z - 1);
if (witness (l, z) == TRUE)
return COMPOSITE;
}
return PRIME;
}

unsigned long long
primeornot (unsigned long long u)
{
unsigned long long z, k;

if (u == 2 || u == 3)
return PRIME;

if (u % 2 == 0)
return COMPOSITE;
if (u % 3 == 0)
return COMPOSITE;

if (u % 6 != 1 && u % 6 != 5)
return COMPOSITE;

if (u < 10000)
if (a[u] == 0)
return PRIME;

k = fermat (2, u);

if (k == COMPOSITE)
return COMPOSITE;
else
k = miller_rabin (u);

if (k == COMPOSITE)
return COMPOSITE;
else
return PRIME;

}

int
primes (void)
{
int i = 2, j = 2, k;
a[0] = 1;
a[1] = 1;
for (i = 2; i * j <= 10000; i++)
a[i * j] = 1;

for (i = 3; i <= 100; i = i + 2)
{
j = i;
if (a[i] == 1)
continue;
for (k = i; k * j <= 10000; k = k + 2)
a[k * j] = 1;

}
for (i = 0; i < 10000; i++)
if (a[i] == 0)
{
table[count++] = i;
}
return 0;

}

int
main ()
{
unsigned long long x, z, i;
primes ();
scanf (" %lld", &z);
for (i = 0; i < z; i++)
{
scanf (" %lld", &x);
if (primeornot (x) == PRIME)
printf ("YES\n");
else
printf ("NO\n");
}

return 0;
}

``````
Last edited by temper_3243 on Wed May 25, 2005 5:07 pm, edited 2 times in total.

mf
Guru
Posts: 1244
Joined: Mon Feb 28, 2005 4:51 am
Location: Zürich, Switzerland
Contact:

Code: Select all

``````for(k=0;k<count;k++){
z=(mycomb(table[k],u-1,u)%u);
if(z!=1 && z!= u-1)
return 0;
}``````
IMHO the code above implements (incorrectly) Fermat's test, not the Miller-Rabin's.

Also, there can be an overflow in this line:

Code: Select all

``k=mycomb((a*a)%p,b/2,p)%p;``
When a >= 2^32, the expression `(a*a)%p' will actually compute the quantity ((a * a) % 2^64) % p.

temper_3243
Experienced poster
Posts: 105
Joined: Wed May 25, 2005 7:23 am
Hi,
Thanks for the suggestion.I found my mistake.
I reimplemented the algorithm
Can you tell me for which inputs the above code breaks ?

How do i avoid overflow ?

frk_styc
New poster
Posts: 10
Joined: Sun Apr 17, 2005 5:00 pm
To avoid overflow you may use the following code:

Code: Select all

``````I64 mod_mul(I64 x,I64 y,I64 n)
{
if(!x) return 0;
return (((x&255)*y)%n+(mod_mul(x>>8,y,n)<<8)%n)%n;
}
``````
The function shown above computes (x * y) % n.

But in my practice I prefer implementing full-precision 64-bit multiplication and module operations because this latter way is much more efficient than what's shown above for it requires much fewer modulo operations of 64-bit integers when some floating-point operations are used.