Prime number

Let's talk about algorithms!

Moderator: Board moderators

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

Prime number

Post by temper_3243 » Wed May 25, 2005 7:28 am

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:

Post by mf » Wed May 25, 2005 9:30 am

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

Post by temper_3243 » Wed May 25, 2005 5:08 pm

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

Post by frk_styc » Sat May 28, 2005 4:24 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.

Post Reply

Return to “Algorithms”