next prime number using While loops

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..

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/6200-tutorial-how-to-ask-a-question-on-answers-and-get-a-fast-answer

[Masum Mushfiq Ishti]

This worked for me:

function k=next_prime(n)
if (~isscalar(n)||(n<0))
    fprintf('error\n');
    return
end
n=n+1;
while true
    if isprime(n)
        k=n;
        break;
    else
        n=n+1;
    end
end

[Suleman]

what do you mean by break ? @Masum Mushfiq Ishti

[DGM]

See break

Answer 1

[John D'Errico]

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.

[Jan]

@John D'Errico: Thanks for this excellent analysis.

The code fails for detecting primes < 1000. nextprime(2) replies 1009. In the condition:

~all(rem(N,partialsieve))

3 is recognized as divisible by 3, but it is still a prime.

[Jan]

@John D'Errico: Beside the above mentioned problem, your method can be improved:

while true
if ~all(rem(N, partialsieve))
elseif isprime(N)
    break
end
N = N + 2;
end

If partialsieve are the primes(1000), using this as test is cheaper than calling isprime(N), where primes(sqrt(N)) is calculated and tested. But here we can exploit 2 two details:

1. Most non-primes are found with a shorter test vector already, e.g. primes(100). Splitting partialsieve in a low and high part and testing the low part at first is cheaper.
2. isprime is based on the same algorithm again and tests with primes(sqrt). Then primes(1000) are tested twice.

The solution is to join both parts and omit the vectorization, but test the values in a loop. Then the smaller prime factors are found earlier:

% Just as demonstration - !!! NOT WORKING !!!
p  = primes(sqrt(N));
np = numel(p);
while true
   isp = 1;          % Assume N is a prime at first
   for ip = 2:np     % Omit the leading p(1)==2, because we process odd values only
      if rem(N, p(ip)) == 0
         isp = 0;    % N is divisable by p(ip)
         break;
      end
   end
   
   if isp, return; end   % N is a prime
   N = N + 2;
end

This is faster than creating the temporary vector in ~all(rem(N, partialsieve)) and in some cases inside isprime again as ~all(rem(N, primes(sqrt(N))). The loop outperforms the vectorized method, because it stops soon usually.

The problem of this method is, that p is determined before the loop, but N is growing and requires more primes to be tested exhaustively. To handle this, p must be larger: primes(sqrt(N + M)) with a sufficient M. And if floor(sqrt(N)) > p(end), p must be updated accordingly. My tests with primes(sqrt(N * 1.1)) are promissing.

[Jan]

A speed comparison:

tic; for k = 1:1e5; nextprimeJohn(k); end; toc

Elapsed time is 4.307417 seconds.

tic; for k = 1:1e5; nextprimeJan(k); end; toc

Elapsed time is 1.090072 seconds.

tic; for k = 1e12:1e12 + 1000; nextprimeJohn(k); end; toc

Elapsed time is 13.324377 seconds.

tic; for k = 1e12:1e12 + 1000; nextprimeJan(k); end; toc

Elapsed time is 9.318602 seconds.

function N = nextprimeJan(N)
% Get next prime greater than input N
if N < 2, N = 2; return; end

if ~isfloat(N) || N > flintmax(class(N)) || ~isfinite(N) || ~isscalar(N)
   error('Jan:nextprime:BadN', 'Input must by a finite scalar.');
end

N = floor(N) + 1;    % Start search at next natural number
N = N + ~rem(N, 2);  % Ckeck odd numbers only

limN = -1;           % Any number < 0
while 1
   if N > limN       % Expand list of primes for checking
      % Create a list of primes until SQRT(N + Delta). Let Delta be small enough
      % to avoid producing too many primes, and large enough to avoid to expand
      % it too often.
      wp   = ceil(sqrt(N * 1.1));
      p    = primes(wp);
      np   = numel(p);
      limN = wp * wp;
   end
   
   % p(np) is only slightly larger than SQRT(N), such that limiting the search
   % for dividers in all p is cheaper than determine the exact p for each N:
   isp = 1;          % Assume N is a prime at first
   for ip = 2:np     % Omit the leading p(1)==2
      if rem(N, p(ip)) == 0
         isp = 0;    % N is divisable by p(ip)
         break;
      end
   end
   
   if isp            % Return if N is a prime
      return;
   end
   
   N = N + 2;        % Next odd N
end
end

[John D'Errico]

Intersting code by @Jan. However, it is only efficient for reasonably small values of N. The reason is, it actually generates the entire list of primes up to sqrt(1.1*N). How many primes is that? If N starts to become large enough, that set of primes will be huge. However, this code will work for doubles, since sqrt(flintmax) is not that huge.

sqrt(1.1*flintmax)

ans = 9.9539e+07

So approximately 1e8.

numel(primes(1e8))

ans = 5761455

which is not too large.

Answer 2

[Giancarlo milon]

You just want the next prime so you can just do

function k = next_prime(n)
k = n + 1
% now k = input + 1
while isprime(k) == 0
% if the k is not prime add 1 till its prime
    k = k+1;
    % when its prime thats the answer end the loop
end    
end

[Rik]

Direct complete solutions to homework questions are discouraged, as they tend to encourage cheating. Why did you post this?

[Giancarlo milon]

i'm posting other type of solutions, if you want to cheat that's you. I'm not forcing anyone to cheat, that's your desicion. Everyone is responsible for their own actions.

[Rik]

I would question how other this is from other solutions. This is still a while loop with a 1 step increment.

I don't question your intentions, but I do question what the most likely use would be. At least your solution is documented, which is more than can be said from most homework-solutions.

Answer 3

[Shubham Sangle]

This is link of answer to similar question https://stackoverflow.com/questions/2468412/finding-a-prime-number-after-a-given-number

Answer mentioned in above link mentions:

Some other methods have been suggested and I think that they are good, but it really depends on how much you want to have to store or compute on the spot. For instance if you are looking for the next prime after a very large number, then using the Sieve of Eratosthenes might not be so great because of the number of bits you would need to store.

Alternatively, you could check all odd integers between (and including) 3 and sqrt(N) on every number odd number N greater than the input number until you find the correct number. Of course you can stop checking when you find it is composite.

If you want a different method, then I would suggest using the Miller-Rabin primality test on all odd numbers above the input number (assuming the input is > 1) until a prime is found. If you follow the list, located at the bottom of the page, of numbers a to check for the given ranges, you can significantly cut down on the number of as you need to check. Of course, you might want to check at least a few of the smaller primes (3,5,7,11 for instance) before checking with Miller-Rabin.

[John D'Errico]

For a little expansion on the ideas mentioned here...

A rough number is a number that is divisible by no small primes. A k-rough number is a number that is divisible by no prime less than k. That allows us to define the extent that a number is known to be rough. (Some definiitions would have that inequality be less than, others would say less than or equal to. The difference is not really important here.)

Rough-ness is useful because it can help you to exclude numbers from testing them using isprime, because isprime is somewhat computationally intensive. Excluding all even numbers from a test essentially allows you to exclude half the numbers from testing to see if they are prime. But then you might also like to exclude all numbers divisible by 3, then 5, 7, 11, etc.

A simple test to identify 42-rough numbers is easy to implement using gcd.

prod(primes(41))
ans =
           304250263527210
log2(prod(primes(41)))
ans =
          48.1122518409569

So we see the product of all primes no larger than 41 is small enough that it fits perfectly into the dynamic range of a double precision integer. And if we could exclude all numbers that have at least one prime factor no larger than 41, we would manage to exclude roughly 85% of all integers.

1 - prod(1 - 1./primes(41))
ans =
         0.854906329597798

For example, is the number 1e9+1 a 41-rough number? GCD will tell us.

gcd(prod(primes(41)),1e9+1)
ans =
     19019

No. It has at least one small prime factor, here we see that to be 3. It may have multiple factors of 3, but we really don't care, since we know now that 1e9+1 cannot possibly be prime.

In fact, if we consider the next 10 integers that exceed 1e9, how many are 41 rough?

gcd(prod(primes(41)),1e9 + (1:10))
ans =
   19019   6   23   82   15    2    1   42   1  170

The idea is to look only for numbers that have a gcd of exactly 1 with our product of small prime, since that implies the number is not divisible by any of them.

So we immediately know that of the integers following 1e9, only 1e9+7 and 1e9+9 can possibly be prime, because only they are 41 rough.

isprime(1e9+[7 9])
ans =
  1×2 logical array
   1   1

In fact, they are both prime numbers. If we wanted to find more primes exceeding 1e9, we could just do:

kr41 = find(gcd(prod(primes(41)),1e9 + (1:50)) == 1)
kr41 =
     7     9    13    19    21    31    33    37
     
isprime(1e9 + kr41)
ans =
  1×8 logical array
   1   1   0   0   1   0   1   0

So only 8 explicit calls to isprime were needed, and half of those tests were successful in identifying a prime number.

Answer 4

[BHARAT JAWLE]

function k=next_prime(n)
k=n+1;

while  factor(k)~= k
    k=k+1;
  
end 

[Rik]

You can solve it like this, but do you actually understand why this works? Do you check the output of the factor function for non-primes? Did you check what while does for non-scalar inputs?

[Rik]

@BHARAT JAWLE You did not edit anything I could notice about your answer. Can you answer the questions in my previous comment?

Answer 5

[RP]

I was working on that exercise too and I tried

function next_prime(n)

x=n
while isprime(x) = 0
    x = n+1
end

however, Matlab returns "The expression to the left of the equals sign is not a valid target for an assignment."

[Stephen23]

You need to read the documentation for every operation that you use, no matter how trivial you think it is. Then you would learn the difference between = (allocation) and == (logical equivalence). In any case, a more MATLAB-ish solution is to use ~ :

if ~isprime(x)

You also forgot to specify any output argument: read the documentation to know how to define functions with input/output arguments

https://www.mathworks.com/help/matlab/ref/function.html

[RP]

Thank you for the tip! Unfortunately the function is still not doing what it's supposed to do and I see that it doesn't look right, but I don't really know how to fix it.

[Stephen23]

You forgot to define an output argument.

Also you never increment n, so x is always just n or n+1, nothing more. Try this:

function n = next_prime(n)
n = n+1;
while ~isprime(n)
    n = n+1;
end
end

[Xenium Adil]

thanks a lot...... actualy i was a little bit wrong... like this

function k=next_prime(n)
k=n;
while ~isprime(n)
    k=n+1;
end

So... thnks

[Stephen23]

@Xenium Adil: your function does not work:

>> next_prime(5) % wrong
ans = 5
>> next_prime(7) % wrong
ans = 7
>> next_prime(8) % infinite loop until I press ctrl+c

>>

Your function either returns an incorrect value or goes into an infinite loop. You really need to test your code and read the comments in this thread.

Answer 6

[VIGNESH B S]

function answer = next_prime(n)
%the function recieves a scalar 'n' 
flag = 1;
% a variable flag will act as break statement (helps to break out of while loop when its value is changed to 0). 
n = n+1;
% We are adding 1 to the scalar n because the question asks about the next prime to 'n'. 
% ( so if n is a prime this function returns the same n because this function is to find primes... but we actually want next prime to n not n itself). 
while flag == 1            % Then comes the while loop to check whether it is a prime or not  
    if isprime(n) == 1
        flag = 0;          % if it is prime it breaks out of while loop 
                           %  ( as flag is made from 1 to 0 and while is only if flag == 1).
    else
        n = n + 1;         %  check for prime if not add 1 .... till a prime next to n is obtained .
    end
end
answer = n;
end

[Rik]

Why did you post this answer? What does it teach? You're more than welcome to start answering questions, but why post a solution to a homework question without any explanation?

If you don't respond this answer will be deleted.

[Rik]

You didn't understand my point. I personally understand what you're doing. But I'm not the intended audience. You posted a solution for a homework question. That means your audience consists of fellow students. They can do two things with your post: learn and cheat.

To be able to learn from your code they must understand it. Your code doesn't contain any comments, so they can only understand it if they already understand what they should do. If they already understand, why would they bother reading this?

So the only thing they can do here is cheat.

If you don't comment your code, how can a student who can't already solve this on their own learn from what you post? (this is not a rhetorical question)

Also, your code is also fairly complex when you compare it to other posts in this thread that do essentially the same.

[VIGNESH B S]

Hope the code is fine now. BTW im new to community and got carried away in solvign and answering questions... so i thought the commenting part was not necessary .. Now its clear .Thank You !!!

[Rik]

You're most welcome to start answering questions, but I would suggest you don't answer any homework questions. You can find guidelines for posting homework on this forum here, I'll leave it to you to extend them to the answering of homework.

If you just want to practice solving puzzles with Matlab, I can recommend Cody. It will not reward good coding skills, but you will have many different problems you can solve.

Answer 7

[Zia Ur Rehman]

My code is running very well.

if there is any mistake plz mention me because I'm very new in coding.

function k = next_prime(n)
k =n+1;
while ~isprime(k)
    k=k+1;
end
end

[Rik]

Why did you post this answer? What does it teach? You're more than welcome to start answering questions, but why post a solution to a homework question without any explanation?

Your code works, but it is far from optimal. To understand what you can do to improve it, you should read the accepted answer.

Answer 8

[Aziz ur Rehman]

Simple as that

function Final_Prime=next_prime(n)
prime=primes(n+n)    ;    i=1
temp=prime(1);
while prime(i)<=n
    i=i+1;
    temp=prime(i);
end
Final_Prime=temp;

[Aziz ur Rehman]

it is just another simple approach for finding the answer.

[Rik]

I can't tell why this would work, be efficient, or how this would be simple. Can you explain it? Where is primes coming from? Why are you using an array?

[Jan]

@Aziz ur Rehman: This apporach is complicated and inefficient.

prime = primes(n+n);

If n is large, e.g. 4e14, calulating all primes until 2*n is expensive.

There is no need to introduce the temporary variable "temp". Simply use the output directly.

Checking all primes, if they are smaller equal than n is a waste of time also.

A simplified (but still less efficient) version of your code:

function m = next_prime(n)
prime = primes(n+n);  % ???
i     = 1;
while prime(i) <= n
    i = i + 1;
end
m = prime(i);

[John D'Errico]

Does this work? Well, yes, but in a terribly inefficient way. (I'm sorry, but that is true.)

First, for a given value of n, this computes ALL prime numbers up to 2*n. If n is at all large, for example, 1e9.

numel(primes(2e9))

ans = 98222287

So in that case, it generates just under 100 million primes, requiring a little under 1 gigabyte of RAM to store them all.

Then it tests EVERY one of those primes, starting at the first prime (which is 2), and compares that number to n. Is it greater than n? If not, then it looks at the next prime in the list of all primes.

The code relies on something that not all people know (ok, Bertrand and Chebyshev knew it)

https://en.wikipedia.org/wiki/Bertrand%27s_postulate

i.e., that there always exists at least one prime in the interval [n,2*n]. When n is small, this is not so problematic. However, I think you would not want to find the next prime larger than say 1e15 with the code by @Aziz ur Rehman.

nextprime(sym(1e15))

ans = 

1000000000000037

The code in this answer is massive overkill. The mathematical equivalent to using a Mack truck to carry a pea to Boston. Rather than just finding the next prime, it first generates potentially millions, or even billions of primes, and then looks to see which one is greater than n. And it does that search very ineffiiciently.

Could this code have been written more efficiently? If a while loop is an absolute imperative, AND the value of n is not too large, AND you cannot use a tool like find, then a while loop does work. But as I said, one can use a Mack truck to carry a pea to Boston. But that does not mean I want to do so as this code was written. For example, had the author used a simple bisection search in that vector of all primes, I would have given the code some grudging respect. And you could write that using a while loop.

Answer 9

[RAVI RANJAN]

function k = next_prime(n)
    total = 0;
    k = 0;
    if isprime(n) == 1
        k = n + 1;
    else
while isprime(n) == 0
      n = n + 1;
      k = n;
end
    end
fprintf('next prime number = %d',k)

[John D'Errico]

You define variables that are never used. What is the purpose of total in there for? Except to dump CPU cycles on the floor?

And your code will not work. If you are going to publish code that claims to give a solution to a problem, surely you would actually test the code? I've copied it directly into this comment. I even left total in there. I did need to add an extra end to the code so it would run though.

Now try that exact code, when n is a prime number. What is the next prime after 7?

k = next_prime(7)

k = 8

Oh, yes. Now I remember. 8 is a prime number. Wow! How could I have forgotten that 8 is prime? Your code works really well. NOT.

function k = next_prime(n)
total = 0;
k = 0;
if isprime(n) == 1
    k = n + 1;
else
    while isprime(n) == 0
        n = n + 1;
        k = n;
    end
end
end

Answer 10

[Aramis]

THE BEST ANSWER

function k = next_prime(n)
    k=n;
    while isprime(k)==false || k==n
            k = k+1;
    end
end

[Stephen23]

SOME IMPROVEMENTS

The loop condition can be simplified by simply incrementing k before the loop (just as some comments on this thread already show). Note that such an increment occurs just once whereas your approach calls EQ and && on every loop iteration.

Using EQ on logical values is not required. Use NOT instead.