Project Euler 123: Determining the remainder when (pn − 1)^n + (pn + 1)^n is divided by pn^2

Problem 123 of Project Euler reads

Let pn be the nth prime: 2, 3, 5, 7, 11, …, and let r be the remainder when (pn-1)n + (pn+1)n is divided by pn2.

For example, when n = 3, p3 = 5, and 43 + 63 = 280 ≡ 5 mod 25.

The least value of n for which the remainder first exceeds 109 is 7037.

Find the least value of n for which the remainder first exceeds 1010.

A nice and short problem description which will also have a nice and short solution post since we solved Problem 120 which is very similar to this problem.

In problem 120, we established that the remainder r for this equation is given as

for odd values of n, and 2 for even values. So all we need to do now is to make a long enough list of primes and search through it until we reach an r which exceed the value we are looking for. So with our trusty old Sieve the code looks like

long limit = 10000000000;
BigInteger r = 0;
int[] primes = Sieve(2, 500000);
int n = 7037;

while (r < limit) {
    BigInteger p = primes[n-1];
    r = 2 * p * n;

There are two things to note here. The upper limit to the sieve is a guess that I made. It might very well be too much or too little. The second thing is that I pull a prime into the variable p inside the loop in line 8. I do that to convert the calculation of r into a biginteger, otherwise I will end up with overflow. I could have done that inline, but I think this is more readable.

I chose to start n at 7037 since we know from the problem that this is the first value where we exceed 109 so any value smaller than that will not give us the result.

The problem runs in less than 5ms where the majority of the time is spent generating the primes. The code gives us the following result

The least value of n = 21035

You can download the source code here. I know this was the short explanation, but please go check Problem 120 for a detailed version of how to derive the basics of this problem

Posted by Kristian


Jean-Marie Hachey

As shown in Fig.1, the parameters of Kristian’s code were increased in order to reach higher values of n and r.

The example illustrated here shows that the solution took nearly 15000 ms.


(One example).
1) long limit: 10^14
2) int[] primes = Sieve(2, 5×10^8)
3) Starting limit: int n = 579351

Results :

The least value of n = 1762729
Solution took 14881.8795 ms

Project Euler 123.
Variation of remainder r as a function of least value n (0 =< n <= 10^14).


Calculations of remainder r done with :
Microsoft Visual C# 2010
(Implementation : System Numerics)


Note :

Variation of the density of primes in the sequence of natural numbers.

Prime numbers distribution

Citation :
« The prime number theorem gives a general description of how the primes are distributed amongst the positive integers.
Informally speaking, the prime number theorem states that if a random integer is selected in the range of zero to some large integer N, the probability that the selected integer is prime is about 1 / ln(N), where ln(N) is the natural logarithm of N. »

So, at 10^10000, the density of prime numbers is about 1 PPM.

Babak Sairafi

Hello Dear Kristian
I think in your code must write r=(2*p*n)%(p*p)
of course that is right but % operator must be write
for example if n=3 then 280 is not correct answer
do you agree?

I am not sure I understand your question. If n= 3 I get that r = 2*5*3 = 30.

Babak Sairafi

yes,but it is not correct, 5 is correct answer .
you must calculate 30%25=5

Yes I see the problem now. You have found the only solution where it makes a difference. And let me explain why.

We are dealing with the fact that r = 2np mod p*p

for n = 5 we have

r = 2*5*11 = 110 mod 121

So now r is less than p^2. And we also have that p increases faster than 2n, since we for certain know that only odd numbers can be prime, so that makes it at least as fast as 2n. And we also know that there are many odd numbers that are not primes.

So what we end up with is a 2n that increases slower than p, and thus we will always have an r which is less than p^2, so any number n> 3 and thus we can skip the modulo part.

Babak Sairafi

Yes,that’s true
thank you

Carl J Appellof

Small observation – in C#, a “long” variable is 64 bits big. Much bigger than your limit of 10^10, so there’s no need to convert to the expensive BigInteger.

Leave a Reply