Official link: https://projecteuler.net/problem=752
Note: My code takes ~230 seconds, 8 seconds with pypy. I have given details on why everything is true, I understand the more detailed explanations in the thread so I am happy with my understanding.
For a thorough understanding of this problem please read the thread
I found this problem really hard. It took me a long time to figure out basic things about the problem, and I found most of them by pure brute force and was only able to prove them after taking the time to understand other peoples solutions.
Through complete brute force for G(1000) and using Binomial expansions (which is very slow) I noticed some very crucial compenents of the problem.
I could not prove these and even then my program was way too slow. After some inspiration from https://mathlesstraveled.com/2017/07/06/the-curious-powers-of-1-sqrt-2-recurrences/ I realised I had wasted a lot of time and not realised the obvious way to approach this problem.
With these expression I could calculate (1 + √(7))^n extremely quickly using binary exponentiation. With this, the problem is solvable
Go through all primes, p
Go through all divisors, d, of p^2 - 1
Calculate (1 + √(7))^d and see if a(n), b(n) = 1, 0
If yes then we have found g(x)
If no then try next divisor
Go through all numbers x = 6k +- 1
Let g(x) = 1
Go through the prime factorization of n = p
For each pair (p, e) representing p^e dividing n, g(x) = lcm(g(x), g(p)*p^(e - 1))
I struggled to find a smart way to deal with the constant need to find the divisors of p^2 - 1 in the end I just used sympy divisors function as it is much faster than mine and I could solve the problem.
Below I give some justification for the points 1-4 listed above:
For actual proof we need some group/field theory which I can personally understand but my explanations will not be as good as those in the thread, hence I recommend you read the thread for a proper understanding of what is going on.
Enter an integer (yourinput)
Code will output G(yourinput)