Script doesn't seem to be execute properly

I was hoping to acquire some help on how to make my program work efficiently and not take a substantial amount of time to finish:

clear variables
a=...;
p=nextprime(a);
count=0;
limit=200000;
tic
while isprime((p-1)/2)~=1
    a=a+1;
p=nextprime(a);
 count=count + 1;
 if count>limit
        break
 end
end
toc

This program outputs a number p greater than a such that p is prime and (p-1)/2 is prime. However I've noticed that for any number a greater than approximately 15 digits, the program will take an absurd amount of time to finish, which isn't ideal since I need to test numbers of the order 10^50.

Answer 1

[Walter Roberson]

Beyond about 4E15 the distance between adjacent representable doubles becomes greater than 1. p becomes forced to be even (and so not a prime) and p-1 becomes the same as p .

You can do marginally better by switching to uint64, which gets you to about 1.8E19 . But you cannot get beyond that using ordinary numeric forms.

You need to switch to a variable precision toolbox, such as Symbolic Toolbox, or John D'Errico's File Exchange contribution for variable precision integers.

[Manuel Barros]

I have attempted to input a as a symbolic number but it still take a rather long time for the aforementioned scale of prime numbers I'm trying to calculate. I'm trying to see a way to make my code more efficient to the point that it shouldn't take more than a few minutes to find such a number, but no success so far.

[Walter Roberson]

prime tests for a series of numbers are sometimes more efficient as you test more values at one time . It depends on the implementation .

Is the question to find the first such number greater than a or to find some such number greater than a?

[Manuel Barros]

It is to find the first prime greater than a such that (p-1)/2 is also prime

[Walter Roberson]

Then it is going to depend upon the quality of implementation of isprime() or nextprime() . There is a possibility that it might be faster to test

test_vals = p : 2 : p + 10000;
candidate_mask = isprime(test_vals);
next_few_primes = test_vals(candidate_mask);

instead of looping doing nextprime().

But that is going to depend on how the isprime() and nextprime() are implemented in the symbolic package.

Answer 2

[Christopher Creutzig]

In your code, you spend a lot of time computing the same prime over and over again. Do not start the search at a+1 for the second search, but start after the prime you already found.

It might also be marginally faster to look for the next prime q starting at a/2 such that p=2*q+1 is also prime.

>> tic
>> a = sym('12345678901234567890');
>> q = nextprime(fix(a/2));
>> while ~isprime(2*q+1), q = nextprime(q+1); end
>> toc
Elapsed time is 4.304840 seconds.
>> [q, 2*q+1]
 
ans =
 
[ 6172839450617290091, 12345678901234580183]

[Stephen Cobeldick]

+1 impressively fast for symbolic math!

[Manuel Barros]

Yes thank you, I happened to notice this underlying issue a while after. :)