MATLAB: Next prime number using While loops

eratosthenesprime numbersoft-lockwhile loop

We have to write a function called next_prime that takes a scalar positive integer input n. Use a while-loop to find and return k, the smallest prime number that is greater than n. we may use isprime function. plz help..

Best Answer

There have now been multiple answers to this homework, all of which are fairly equivalent. (The response by Shubham Sangle has a lot of good points.) Let me offer a few ideas.

A good code to solve for the next prime would try to avoid calling isprime as much as possible, because calls to isprime are expensive. Of course, this depends on how large a number you are looking at

isprime(10000000000037)
ans =
  logical
   1
   
timeit(@() isprime(10000000000037))
ans =
            0.009439529724

10000000000037 is close to as large a number you can use isprime on, working in double precision. If a single test for primality costs roughly 0.01 seconds (at least on my computer, your mileage may vary) then how long will it take to find the next prime after some number?

As it turns out, 10000000000037 is the next prime that exceeds 1e13. After that, we see primes at 10000000000051, 10000000000099, 10000000000129, 10000000000183, 10000000000259...

So out as far as 1e13, the spaces between primes seem to be running around 30 to 70, sometimes more, sometims less. That means, if you wanted to find the next prime after one of those primes, you might need to test perhaps as many as 70 numbers, or more. (However, you might get unlucky, if you find yourself at the beginning of a rather large hole between primes.)

If a single test by isprime takes 0.01 seconds out there, then your code for nextprime might take as much as one second of CPU time. That is not too bad, as long as one second is acceptable to you. Even so, there are some tricks you might employ.

The idea I will suggest here relies on the concept of k-rough numbers. Thus a rough number is one that has no small prime divisors. A k-rough number is one that has no prime divisors less than k. Some definitions will allow k as the smallest divisor, whereas I won't lose any sleep over the distinction.

https://en.wikipedia.org/wiki/Rough_number

A test in advance for small divisors is not a bad thing. For example, you can rule out all even numbers in your search. And you can rule out all multiples of 3, 5, 7, etc. That is, suppose I picked an arbitary large integer. What is the probability that it happens to be divisible by some small primes? (I call this pre-test a partial sieve.) Thus, 50% of the time, some large number will be divisible by 2. But 2/3 of the time (67%), it will be divisible by EITHER 2, OR 3. That means we could exclude 67% of the numbers, merely by knowing in advance if they were divisuble by EITHER 2, or by 3.

100*(1 - prod(1 - 1./[2 3]))
ans =
          66.6666666666667
100*(1 - prod(1 - 1./[2 3 5 7 11 13 17 19]))
ans =
          82.8975977582789
100*(1 - prod(1 - 1./primes(100)))
ans =
          87.9682709525065
100*(1 - prod(1 - 1./primes(1000)))
ans =
          91.9034736493157

So, if we excluded all numbers that are divisible by any small prime below 1000, then we could reduce the number of calls to isprime by almost 92%. We can carry this pretty far of course. So, if we pretest for divisibility by all primes below 1e8, the exclusion rate goes down to only roughly 97%.

100*(1 - prod(1 - 1./primes(1e8)))
ans =
          96.9520278389414

That is important when a single call to isprime costs 30 minutes of CPU time. (This only happens if you are looking at numbers with many thousands of decimal digits though, using perhaps the sym version of isprime.)

As such, a decent version of nextprime might look like this:

function nextp = nextprime(N)
  % solve for the next prime that comes immediately after N
  
  % insure that N is itself an integer.
  if N == ceil(N)
    % then N is an integer already. We need to start looking at the NEXT
    % prime number though, so we need to start our search at N+1.
    N = N + 1;
  else
    % N was not an integer at all. So ceil will suffice to increase N
    N = ceil(N);
  end
  
  % is N even? If so, then since we need ony test odd numbers, since
  % even numers are not prime, then we can increment it to an odd number.
  % however, just in case N was so small that the next prime would have
  % been 2, special case that.
  if N < 2
    nextp = 2;
    return
  elseif rem(N,2) == 0
    % we only need to look at odd numbers after this.
    N = N + 1;
  end
  
  % Set up a partial sieve. Here, I'll set the sieve to check for divisibility
  % by primes that do not exceed 1000. We don't need to check for divisibility
  % by 2 though, since we will just step around all even numbers anyway, and
  % we have already excluded the case where N was less than 2. So we never
  % need to worry about the only even prime at this point.
  partialsieve = primes(1000);
  partialsieve(1) = [];
  
  % finally, this is just a while loop. The test in the while becomes a no-op,
  % because we break out when a prime gets hit. Use a while instead of a for
  % because we have no idea how long the search will go. It does take a very long
  % time before the intervening distance between two primes gets very large, so
  % the while loop will terminate before long.
  while true
    if ~all(rem(N,partialsieve))
      % no need to test this candidate
      N = N + 2;
    else
      if isprime(N)
        nextp = N;
        % just break out of the loop. No need to even reset flag.
        break
      else
        N = N + 2;
      end
    end
  end
end

This code is actually pretty efficient, allowing us to exclude tests with isprime for 92% of the numbers we might otherwise try. In fact, it looks like it had to make only one call to isprime in this example run:

timeit(@() nextprime(1e13))
ans =
            0.009934246724

I can make that claim, since that is roughly the time that MATLAB took for ONE call to isprime out there. And in this next test, there should have been only 5 calls to isprime.

timeit(@() nextprime(10000000000183))
ans =
            0.047855878724

The nextprime code above is efficient, since it is MUCH faster to test for divisibility by even a rather lengthy list of integers using the rem test, then it it to make one call to isprime out there.

How long will the search for the next prime require? Luckily, the primes are scattered randomly, but it takes a seriously long time before the inter-prime stride becomes terribly large. At the same time, it can be shown this stride can be arbitrarily large if you go out far enough.

As you can see below, it generally will not take that many tests before the next prime is found, but it will help if you can avoid some of those tests with careful code.

>> max(diff(primes(1e6)))
ans =
   114
>> max(diff(primes(1e7)))
ans =
   154
>> max(diff(primes(1e8)))
ans =
   220
>> max(diff(primes(1e9)))
ans =
   282

Will this all be significant for large numbers? Yes, very much so. If you plot the maximum interprime stride, as a function of the size of the numbers we are looking at, you will find that on a plot with logged x axis, the curve should be roughly linear.

semilogx(10.^[6 7 8 9],[114 154 220 282],'-o')
grid on

I'll just fit a linear polynomial to it, as:

interprimepoly = polyfit([6 7 8 9],[114 154 220 282],1)
interprimepoly =
                        57                      -235

This will be adequate for an extrapolation, since I don't care if it is only wildly approximate.

polyval(interprimepoly,[13, 100000])
ans =
                       506                   5699765

As such, out in the vicinity of 1e13, we may possibly be pushed to test several hundred numbers to see if they are prime. If we can intelligently exclude perhaps 90% of those tests, then we can cut the time required by up to 90%.

Further out, looking at numbers with tens or hundreds of thousands of decimal digits, the importance of avoiding millions of calls to isprime will be quite significant.

Related Solutions

MATLAB: Function to find the next prime number…

Method 1 This functionality does not run in MATLAB and can be used only via MuPAD Notebook Interface.

To create an MuPAD interface use the following code

 mupad
nb = allMuPADNotebooks

Then a interface screen shall pop up where you can proceed by using the nextprime(num) function.

Method 2

You could also create a user defined function to compute the next prime number. This function takes only a non-negative integers as an argument

 function p = nextprime(n)
if (isprime(n))
    p=n;
else
    while(~isprime(n))
        n=n+1;
    end
    p=n;
end
end

MATLAB: I wrote a code to produce twin primes, however it is pretty strenuous for the computer and I need some high numbers to check a conjecture. Is there a better way to get twin primes. Note: these twin primes are represented by 6m’s where 6m+-1 are pri

Ok, it looks like you are still interested in this, so here is an approach targetted more at what I hear from your last comment.

Simple would be, if you just want to find all twin primes less than some value, just use primes to get all primes, then find all primes that are exactly 2 larger than the one before.

N = 1e7;
P = primes(N);
ind = find(diff(P) == 2);
TP = [P(ind)',P(ind+1)'];
size(TP)
ans =
       58980           2
TP(1:10,:)
ans =
     3     5
     5     7
    11    13
    17    19
    29    31
    41    43
    59    61
    71    73
   101   103
   107   109

But there is no formula to find all such twin primes, nor to predict when a twin prime pair will occur, except that it must be a pair of the form 6*m+[-1 1].

Yes, you could use isprime to test each pair of numbers of that form, but that would be extremely inefficient, far less good than just calling primes directly.

So, is there something slightly better? For example, could we write a variation of a sieve, targetted specifically at twin primes? For example, suppose I wanted to find all twin primes less than 1e9? primes(1e9) is probably too big for me to do on my computer. But I know that twin primes must be of the form 6*m +/-1.

tic
% find all twin prime pairs where both are less than N
N = 1e9;
Plist = primes(sqrt(N));
% Keep the small twin primes from Plist
twin_primes = Plist(diff(Plist) == 2);
twin_primes = [twin_primes;twin_primes+2]';
% We don't need to check 6*m-1 for divisibility by 2 or 3.
Plist(1:2) = [];
 
% m will correspond to the index of a potential twin prime, 6*m-1.
% At the end of this comutation, if m(k) is true, then
% 6*m(k) +/- 1 will be a twin prime. Note that all of this will miss the
% very first twin prime pair, [3 5]. But all other twin prime pairs are of
% the form 6*m +/- 1.
m = true(1,floor((N-2)/6));
 
for p_i = Plist
  % which elements of m correspond to 6*m-1 that are divisible by p_i?
  % Find them using the modular inverse of 6, modulo p_i.
  p_i_inv = minv(6,p_i);
  
  m1 = p_i_inv;
  if p_i == (6*m1 - 1)
    m1 = m1 + p_i;
  end
  m(m1:p_i:end) = false;
  
  m2 = mod((p_i - 1)*p_i_inv,p_i);
  
  if p_i == (6*m2 - 1)
    m2 = m2 + p_i;
  end
  m(m2:p_i:end) = false;
  
end
m = find(m);
 
twin_primes =[twin_primes;[6*m-1;6*m+1]'];
toc

So reasonably fast.

Elapsed time is 9.773307 seconds.

This ran pretty quickly on my computer, since the only RAM it required was a vector of size N/6 bytes. So for N=1e9, m only required 0.17 gigabytes. And the loop was also pretty efficient.

Note that the above code uses my minv tool, which is found in my vpi toolbox, although when used on doubles, it returns double precision output.

whos m twin_primes
  Name                   Size                    Bytes  Class     Attributes
  m                      1x3424019            27392152  double              
  twin_primes      3424506x2                  54792096  double

So it found every one of the 3.4 million twin prime pairs for which neither of them exceed 1e9. But since it did not need to explicitly test any of them explicitly for primality, nor did it need to generate the entire list of primes below 1e9, the computation was actually pretty fast. While my computer is moderately old, I'd guess I could easily push N out further out than 1e10. (I just tried it for N=1e10, and it took 110 seconds to find all 27412679 twin prime pairs less than 1e10. Not bad for a 10 year old computer.)

I'm not sure if this will be of any real value to you, but it will indeed generate all such twin primes less than N.

Best Answer

Related Solutions

MATLAB: Function to find the next prime number…

MATLAB: I wrote a code to produce twin primes, however it is pretty strenuous for the computer and I need some high numbers to check a conjecture. Is there a better way to get twin primes. Note: these twin primes are represented by 6m’s where 6m+-1 are pri

Related Question