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

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

       1. If yes then we have found g(x)

       2. If no then try next divisor

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