Algorithm to find all primes from 2 to 1000 not working

Question

Here's a piece of code to compute all primes from 2 to 1000 using the statement, that a number n is a prime number iff:

In the first version I think that I implemented the algorithm correctly:

public class Giuga {

    public static void main(String[] args){
        int n = 2;

        while(n<=1000){
            int k = 1;
            long sum = 0;
            while(k<=n-1){
                sum = sum+(long)Math.pow((double)k,(double)n-1);
                k++;
            }
            if(sum%n==n-1){
                System.out.println(n + " is a prime.");
            }
            n++;
        }
    }
}

But, since the variable sum grows rapidly, an overflow happens and after the prime number 17 there will be no output anymore.

To prevent that I have to use this:

Well, I did that and here is my 2. version:

public class Giuga {

    public static void main(String[] args){
        int n = 2;

        while(n<=1000){
            int k = 1;
            long sum = 0;
            while(k<=n-1){
                sum = sum+((long)Math.pow((double)k%n,(double)n-1))%n; //Here are the changes
                k++;
            }
            if(sum%n==n-1){
                System.out.println(n + " is a prime.");
            }
            n++;
        }
    }
}

I think I did it correctly, but now the output stops after the prime number 13.

I'm trying to find my mistake for quite some time now. What am I doing wrong? There must be 168 primes from 2 to 1000.

Answer 1

As has been pointed out, doubles, which only have about 16 digits of precision, aren't precise enough to maintain the correct remainders for calculations on high enough numbers.

You can switch to longs and perform you own modular exponentiation.

int k = 1;
long sum = 0;
while(k<=n-1){
    long pow = 1;
    for (int i = 0; i < n - 1; i++)
        pow = (pow * k) % n;
    sum = (sum + pow)%n;
    k++;
}

This algorithm could be improved by changing this straightforward modular exponentiation to using modular exponentiation by repeated squaring, and it's not the most efficient prime finding algorithm, but it is now correct.

2 is a prime.
3 is a prime.
5 is a prime.
7 is a prime.
11 is a prime.
13 is a prime.
17 is a prime.
19 is a prime.
23 is a prime.
29 is a prime.
31 is a prime.

(snip)

977 is a prime.
983 is a prime.
991 is a prime.
997 is a prime.

To make it modular exponentiation by repeated squaring, replace

long pow = 1;
for (int i = 0; i < n - 1; i++)
    pow = (pow * k) % n;

with

long pow = 1;
long square = k;
int exp = n - 1;
while (exp > 0)
{
    if ((exp & 1) == 1)
    {
        pow = (pow * square) % n;
    }
    square = (square * square) % n;
    exp >>= 1;
}

It tests each bit of the exponent in succession, and multiplies the current square in to pow if it's set.

Answer 2

Java doubles (i.e., the outputs of pow) do not represent large numbers precisely enough to yield a correct remainder. You should switch to modular exponentiation.

Answer 3

You can always use the BigInteger class for large number computations

private static boolean isPrime(int n) {
    BigInteger N = new BigInteger(String.valueOf(n));
    BigInteger N_MINUS_1 = N.subtract(BigInteger.ONE); 

    BigInteger sum = BigInteger.ZERO;
    for (int k = 1;  k < n; k++)
        sum = sum.add(new BigInteger(String.valueOf(k)).modPow(N_MINUS_1,N)).mod(N);
    return sum.equals(N_MINUS_1);
}

It's interesting though that this is a variation of Fermat's little theorem, for each k in the sum Σ

k^(n-1)%n should be 1, else this number is not prime! So, if we find a k that k^(n-1)%n != 1, we can stop calculations. The above algorithm can be rewritten as:

private static boolean isPrimeFermat(int n) {
    BigInteger N = new BigInteger(String.valueOf(n));
    BigInteger N_MINUS_1 = N.subtract(BigInteger.ONE); 

    for (int k = 1;  k < n; k++){
        if (new BigInteger(String.valueOf(k)).modPow(N_MINUS_1, N).equals(BigInteger.ONE) == false)
            return false;
    }       
    return true;
}

Voilà!

Answer 4

Your powers are too big as well (guess what 999^999 is). So you have to calculate the power as well using stepwise (ab) mod n = ((a mod n)(b mod n)) mod n (modular exponentiation):

public class Giuga {

    public static void main(String[] args) {
        int count = 0;
        for (int i = 2; i < 1000; i++) {
            if (isPrime(i)) {
                count++;
                System.out.printf("%d is prime\n", i);
            }
        }
        System.out.printf("Found %d primes.\n", count);
    }

    // checks if number is a prime
    private static boolean isPrime(int number) {
        int bigSum = 0;
        for (int k = 1; k <= number - 1; k++) {
            bigSum = (bigSum + calcPowMod(k, number - 1, number)) % number;
        }
        return bigSum % number == number - 1;
    }

    // calculates (a^b)%mod, making sure that intermediate results are always < max(a^2,mod)
    private static int calcPowMod(int a, int b, int mod) {
        int pow = 1;
        for (int k = 1; k <= b; k++) {
            pow = (pow * a) % mod;
        }
        return pow;
    }
}