Project Euler 659 - Largest Prime

Official link: https://projecteuler.net/problem=659

Thought Process

This problem taught me a lot about sieve type algorithms, but first a mathematical proof.

So this means we are looking for the largest prime q that divides 4k^2 + 1, here is where I upped my level on sieve understanding.

I found the following article: http://www.devalco.de/quadr_Sieb_4x%5E2+1.php which explains how to prime sieve 4n^2 + 1. You can read the article for the mathematical proof, I will explain the algorithm.

  1. Initialise a list f where f[x] = 4x^2 + 1 from x = 0 to x = limit (limit is defined by you)

  2. Initialise a second list max_prime_factor = [0] * (limit + 1), in this list max_prime_factor[x] = the maximum prime factors of 4x^2 + 1

  3. Go through the list f from the first element to the last element

    1. Initialise divisor = f[x]

    2. If f[x] > 1, then we continue with the sieve, otherwise we continue going through the loop

      1. From the article we need to check if, f[x + k*divisor] is divisible by divisor where (x + k*divisor) < limit

        1. While f[x + k*divisor] is divisible by divisor, we divide it and we need to assign max_prime_factor[x + k*divisor] = the maximum between it's current value and divisor, in python it looks like this: max_prime_factor[x + k*divisor] = max(max_prime_factor[x + k*divisor], fx)

      2. From the article we also need to check if, f[-x + k*divisor] is divisible by divisor where (-x + k*divisor) < limit

        1. While f[-x + k*divisor] is divisible by divisor, we divide it and we need to assign max_prime_factor[-x + k*divisor] = the maximum between it's current value and divisor, in python it looks like this: max_prime_factor[-x + k*divisor] = max(max_prime_factor[-x + k*divisor], fx)

  4. At the end the list f should be reduce to a list of only 1's, and max_prime_factor is a list where max_prime_factor[x] = the maximum prime factors of 4x^2 + 1, so all we need to do is sum the list max_prime_factor and mod 10^18 to find the last 18 digits

Interactive Code

Enter an integer (yourinput)

Code will output the last 18 digits of the sum from k = 1 to yourinput of P(k)

I have also added my code below, because this sieve is quite complicated