Is there a way to determine a formula giving all *integer* values of $x$ for which the value of a polynomial $P(x)$ with *integer coefficients* is a square?


That is, is there a closed formula for:


$X = \{ x \in \mathbb{N} : \exists \ n \in \mathbb{N} : P(x) = n^2 \}$ ?


I'm interested in particular in $P'(x) = 8x^2-8x+1$, but am wondering about the general case as well.


For $P'(x)$ a sample of $X$ is $\{ 1, 3, 15, 85, 493, 2871, 16731, 97513, \ldots \}$.



Integer polynomials taking square values
There's a fairly detailed explanation of the solution to a similar equation [here](https://web.archive.org/web/20150908090650/http://mathforum.org/library/drmath/view/73118.html). See also [this page](https://www.alpertron.com.ar/QUAD.HTM), which can give you an automated step-by-step solution to such quadratic diophantine equations.


I'll also add that the command Reduce[8 x^2 - 8 x + 1 - y^2 == 0 && Element[x | y, Integers], {x, y}] will produce the answer to your particular problem in Mathematica fairly quickly. I'm making this an answer because the output is too huge to fit into the comments.



```
(C[1] \[Element] Integers && C[1] >= 0 && 
   x == 1/32 (16 + 
       4 (-2 (17 - 12 Sqrt[2])^C[1] + 
          Sqrt[2] (17 - 12 Sqrt[2])^C[1] - 2 (17 + 12 Sqrt[2])^C[1] - 
          Sqrt[2] (17 + 12 Sqrt[2])^C[1])) && 
   y == 1/2 ((17 - 12 Sqrt[2])^C[1] - 
       Sqrt[2] (17 - 12 Sqrt[2])^C[1] + (17 + 12 Sqrt[2])^C[1] + 
       Sqrt[2] (17 + 12 Sqrt[2])^C[1])) || (C[1] \[Element] Integers &&
    C[1] >= 0 && 
   x == 1/32 (16 + 
       4 (-2 (17 - 12 Sqrt[2])^C[1] + 
          Sqrt[2] (17 - 12 Sqrt[2])^C[1] - 2 (17 + 12 Sqrt[2])^C[1] - 
          Sqrt[2] (17 + 12 Sqrt[2])^C[1])) && 
   y == 1/2 (-(17 - 12 Sqrt[2])^C[1] + 
       Sqrt[2] (17 - 12 Sqrt[2])^C[1] - (17 + 12 Sqrt[2])^C[1] - 
       Sqrt[2] (17 + 12 Sqrt[2])^C[1])) || (C[1] \[Element] Integers &&
    C[1] >= 0 && 
   x == 1/32 (16 - 
       4 (-2 (17 - 12 Sqrt[2])^C[1] + 
          Sqrt[2] (17 - 12 Sqrt[2])^C[1] - 2 (17 + 12 Sqrt[2])^C[1] - 
          Sqrt[2] (17 + 12 Sqrt[2])^C[1])) && 
   y == 1/2 ((17 - 12 Sqrt[2])^C[1] - 
       Sqrt[2] (17 - 12 Sqrt[2])^C[1] + (17 + 12 Sqrt[2])^C[1] + 
       Sqrt[2] (17 + 12 Sqrt[2])^C[1])) || (C[1] \[Element] Integers &&
    C[1] >= 0 && 
   x == 1/32 (16 - 
       4 (-2 (17 - 12 Sqrt[2])^C[1] + 
          Sqrt[2] (17 - 12 Sqrt[2])^C[1] - 2 (17 + 12 Sqrt[2])^C[1] - 
          Sqrt[2] (17 + 12 Sqrt[2])^C[1])) && 
   y == 1/2 (-(17 - 12 Sqrt[2])^C[1] + 
       Sqrt[2] (17 - 12 Sqrt[2])^C[1] - (17 + 12 Sqrt[2])^C[1] - 
       Sqrt[2] (17 + 12 Sqrt[2])^C[1])) || (C[1] \[Element] Integers &&
    C[1] >= 0 && 
   x == 1/32 (16 + 
       4 (2 (17 - 12 Sqrt[2])^C[1] + Sqrt[2] (17 - 12 Sqrt[2])^C[1] + 
          2 (17 + 12 Sqrt[2])^C[1] - 
          Sqrt[2] (17 + 12 Sqrt[2])^C[1])) && 
   y == 1/2 (-(17 - 12 Sqrt[2])^C[1] - 
       Sqrt[2] (17 - 12 Sqrt[2])^C[1] - (17 + 12 Sqrt[2])^C[1] + 
       Sqrt[2] (17 + 12 Sqrt[2])^C[1])) || (C[1] \[Element] Integers &&
    C[1] >= 0 && 
   x == 1/32 (16 + 
       4 (2 (17 - 12 Sqrt[2])^C[1] + Sqrt[2] (17 - 12 Sqrt[2])^C[1] + 
          2 (17 + 12 Sqrt[2])^C[1] - 
          Sqrt[2] (17 + 12 Sqrt[2])^C[1])) && 
   y == 1/2 ((17 - 12 Sqrt[2])^C[1] + 
       Sqrt[2] (17 - 12 Sqrt[2])^C[1] + (17 + 12 Sqrt[2])^C[1] - 
       Sqrt[2] (17 + 12 Sqrt[2])^C[1])) || (C[1] \[Element] Integers &&
    C[1] >= 0 && 
   x == 1/32 (16 - 
       4 (2 (17 - 12 Sqrt[2])^C[1] + Sqrt[2] (17 - 12 Sqrt[2])^C[1] + 
          2 (17 + 12 Sqrt[2])^C[1] - 
          Sqrt[2] (17 + 12 Sqrt[2])^C[1])) && 
   y == 1/2 (-(17 - 12 Sqrt[2])^C[1] - 
       Sqrt[2] (17 - 12 Sqrt[2])^C[1] - (17 + 12 Sqrt[2])^C[1] + 
       Sqrt[2] (17 + 12 Sqrt[2])^C[1])) || (C[1] \[Element] Integers &&
    C[1] >= 0 && 
   x == 1/32 (16 - 
       4 (2 (17 - 12 Sqrt[2])^C[1] + Sqrt[2] (17 - 12 Sqrt[2])^C[1] + 
          2 (17 + 12 Sqrt[2])^C[1] - 
          Sqrt[2] (17 + 12 Sqrt[2])^C[1])) && 
   y == 1/2 ((17 - 12 Sqrt[2])^C[1] + 
       Sqrt[2] (17 - 12 Sqrt[2])^C[1] + (17 + 12 Sqrt[2])^C[1] - 
       Sqrt[2] (17 + 12 Sqrt[2])^C[1]))
```

