Search code examples
cprimesgmp

How to improve finding Twin Primes


So I'm using the GMP library on the C language in order to find Twin primes above a certain value. While I'm confident that my strategy would work, the problem becomes the fact that it takes an immense amount of time (I understand that it becomes difficult in finding primes the higher you get.) Is there a way to optimize the search? Here's a snippet of the code that I have:

    mpz_ui_pow_ui(a, base, exponent);
    mpz_nextprime(b, a); // b is the next prime number after a.
                         // c and d will be prime + 2 and
                         // prime - 2.

    /* Fortunate of fortunalities, mpz_nextprime gives the next
       prime greater than what one adds in! */
    /* We need to test if numbers are prime too. */
    while (al == false) {
        mpz_add_ui (c, b, 2);
        mpz_add_ui (d, b, -2);
        if ((mpz_probab_prime_p(c, 15) == 2) ||
            (mpz_probab_prime_p(d, 15) == 2)) { // Returns 2
                                                // if c/d are
                                                // definitely
                                                // prime.
            mpz_set(firstprime,b); 
            al == true;
            break;
        }
        {
            mpz_nextprime(b, b); // b is the next prime number
                                 // after a. c and d will be
                                 // prime + 2 and prime - 2.
        }
    }
    printf("first twin is: ");
    mpz_out_str(stdout, 10, firstprime);
    printf("\n");
    printf("second twin is: ");
    if (mpz_probab_prime_p(c, 15) == 2) {
        mpz_out_str(stdout, 10, c);
    } else {
        mpz_out_str(stdout, 10, d);
    }
    printf ("\n");

Solution

  • Your algorithm is a little odd. You don't test whether b itself is a prime, but test one or both of b - 2 and b + 2. Then, if either of those is definitely a prime, you declare b is one of the twin primes.

    mpz_nextprime may return a non-prime, since it is using a probabilistic algorithm.

    @chqrlie is correct in pointing out that b - 2 had already been processed by mpz_nextprime. The only edge case is if the first call to mpz_nextprime resulted in a number only one or two away from a.

    Since you are willing to accept that b is only probably a prime, you should be happy if both are probably primes. So:

    /* a won't be prime */
    mpz_ui_pow_ui(a, base, exponent);
    
    if (exponent == 0) {
        mpz_nextprime(firstprime, a);
    } else {
        /* Handle the edge case of a - 1 and a + 1 being twins */
        mpz_sub_ui(b, a, 2);
        mpz_nextprime(firstprime, b);
    }
    
    for (;;) {
        mpz_add_ui(c, firstprime, 2);
        if (mpz_probab_prime_p(c, 15) > 0) {
            break;
        }
        /* Optimize out an mpz_set call, thanks @chqrlie */
        mpz_nextprime(firstprime, c);
    }
    

    Which would find probably twin primes. If you want at least one to be definitely prime, you can implement your own prime test, or add a mpz_probab_prime_p call for firstprime.