Sieve of Eratosthenes

The Sieve of Eratosthenes is simplest method to produce a table of prime numbers.  We start
 with a table of integers 2 to n, and starting with 2, we take the current number in the 
 table and remove all its higher multiples.  Once the currentt number in the table exceeds
 the square root of n, then  all the numbers are prime.

 The following function is a straightforward implementation of the algorithm. The
 variable P in the function starts as a list of numbers from 2 to n.  We construct a table 
 of the multiples of the  current prime in  the varianle x. These are removed from P
 using the Complement command.

SieveOfEratosthenes[n_]:= 
         Module[ {t,P,i=1,x},
         t=N[Sqrt[n]];
         P=Range[2,n];
         While[ P[[i]] <= t,
              x=Table[k P[[i]],
                        {k,P[[i]], Floor[n/P[[i]]]}];
              P=Complement[P,x];
              i++
              ];
          Return[P]]

SieveOfEratosthenes[100]

You can store the list of primes in a file  by using the >> redirection operator.

SieveOfEratosthenes[100]>> primetable

The previous command created a file called primetable that contained the primes up to
 100. These can be read later by the ReadList command so that the primes can be used
 in functions such as primality testing and factoring.

primes=Flatten[ReadList["primetable", Expression ]]

primes[[10]]

Warning:  You cannot use the Sieve directly for producing a  very large table of primes in
 Mathematica.  Since Mathematica is slow, it is best to do this incrementally. The above
 function can be used to produce a table of primes up to 10^5 in a few minutes. Then 
 modify the function to generate the list of primes between two number n and m, where we 
 we already know the primes up to square root of n. Then you can go up in steps of
 10^5 until a large table of primes is constructed.

Simple Test of Primality

The following function demonstrates a test of primality using the  the table of primes constructed 
 above. The variable primes containes the list of primes  read from the file primetable using
 the ReadList command. Since primetable contained a list of primes up to 100, the function
 is correct  for input values under 10,000.

(* The function Isprime is using a list of primes in
  primes as a global variable *) 
Isprime[n_]:= Module[ {i=1,t },
   t=N[Sqrt[n]];
   While[ primes[[i]] <=t,
    If[Mod[n,primes[[i]]]==0, Return[False]];
    i++];
   Return[True]]

Isprime[25]

Isprime[37]

isprime[49999]

isprime[99997]

Isprime[99991]

Isprime[ Prime[1000]Prime[3000]]//Timing

Isprime[Prime[4000] Prime[5000]]//Timing