How to find next Pythagorean prime

After reading a number from the keyboard, I need to find the next Pythagorean prime that can be expressed as 4n-1 and as the sum of two squares. I'm trying to find the next prime number that can be expressed as 4n-1 but I'm not having much luck.

prompt = 'Enter number';
x = input(prompt);
next = nextprime(sym(x));
while ((next-1)/4 ~= 0)
  next = nextprime(sym(next));
end
disp(next);

This is what I wrote but when I run the code, it never displays a number. Is there anyway to fix this problem?

[Josh Meyer]

How robust does the program need to be? Are you going to be feeding in very large numbers, for which the program might need to run for quite some time to find the next 4n-1 prime?

If the search is limited to some range, you could precompute the 4n-1 primes in that range and just look up the next 4n-1 prime by comparing them to the input number.

[John D'Errico]

Actually, you won't need to go too far. Roughly 50% of the odd primes are of the form 4*n+1.

p = primes(1000000);
sum(mod(p,4) == 1)/numel(p)
ans =
       0.499057300822951

Pretty close to 50%, even if we don't worry about Chebyshev's bias.

So a geometric series can predict the average number of tests necessary using this scheme.

1/(1-0.5)
ans =
     2

See the link I give. It states that a prime is pythagorean if it is of the form 4*n+1, not 4*n-1.

Thus 73 is a prime of the form 4*n+1.

isprime(73)
ans =
  logical
   1

73 == 3^2 + 8^2
ans =
  logical
     1

mod(73,4)
ans =
     1

And is clearly Pythagorean.

Answer 1

[John D'Errico]

Your problem is easy. The test you did was wrong. I see what you were thinking, but you had a bug in your thoughts. :)

You also have a bug in what you wrote here in your question. So reading this link:

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

To quote the first sentence:

"A Pythagorean prime is a prime number of the form 4n + 1. Pythagorean primes are exactly the odd prime numbers that are the sum of two squares..."

So you need to find primes that are of the form 4*n+1. Not 4*n-1 as you stated.

Anyway, the test that you did was wrong.

while ((next-1)/4 ~= 0)

That effectively tests only if (next-1) is zero. Dividing by 4 does not test for divisibility by 4. What you need to know is if (next-1) is a MULTIPLE of 4. You need to use mod in there. So you could have done this:

while mod(next,4) ~= 1

This continues the loop as long as next is NOT 1 more than a multiple of 4. So the loop continues until next is of the desired form.