Donde dice:

Hay que resolver

  N^2 - 3p^2 -3p - 1 = 0

Multiplicando por 4 esto queda

  (2N)^2 - 3(2p+1)^2 = 1

que es la ecuación de Pell, con la condición de que x sea par e y impar.

La solución general de esta ecuación de Pell es

  x(0) = 2

  y(0) = 1

  x(n) = 2x(n-1) + 3y(n-1)

  y(n) = x(n-1) + 2y(n-1)

De esta solución vemos que las paridades de x e y van cambiando en cada
iteración, por lo que necesitamos las soluciones de índice par, que cumplen

  x(2n) = 7x(2n-2) + 12y(2n-2)

  y(2n) = 4x(2n-2) + 7y(2n-2)

Haciendo x(2n) = 2N(n), y(2n) = 2p(n)+1, resulta

   N(n) = 7N(n-1) + 12p(n-1)+6

   p(n) = 4N(n-1) + 7p(n-1)+3

Los primeros valores para N son

  1 = 0^2 + 1^2

  13 = 2^2 + 3^2

  181 = 9^2 + 10^2

  2521 = 35^2 + 36^2

Resolviendo la ecuación obtenemos la expresión general

  N(n) = ((2+rq(3))^(2n+1) + (2-rq(3))^(2n+1))/4

que todavía no es lo que buscamos... Luego sigo

--

Tenemos entonces

(m + 1)^3 - m^3 =3m^2 +3m + 1 = n^2

Para completar cuadrados, multipliquemos por 4,

3(2m)^2 + 3*2*(2m) + 3 + 1 = (2n)^2

(2n)^2 - 3(2m + 1)^2 = 1

Hagamos 2n = p, 2m + 1 = q y nos queda la ecuación de Pell

p^2 - 3q^2 = 1

de la que nos interesan las soluciones con p par (y por tanto q impar). La
mínima solución no trivial de esta ecuación es (p, q) = (2, 1), por lo que
todas son de la forma:

p(k) + rq(3)q(k) = (2 + rq(3))^k

Las relaciones que verifican p(k) y q(k) se obtienen de

p(k+1)  + rq(3)q(k+1) = (2 + rq(3))(p(k) + rq(3)q(k)) ==>

p(k+1) = 2p(k) + 3q(k)
q(k+1) =  p(k) + 2q(k)

Como p(0) = 1 y q(0) = 0, de estas relaciones vemos que para k impar es p
par y q impar. Podemos considerar entonces las relaciones

p(2k+1)  + rq(3)q(2k+1) = (2 + rq(3))^2(p(2k-1) + rq(3)q(2k-1))

                                       = (7 + 4rq(3))(p(2k-1) +
rq(3)q(2k-1))

p(2k+1) = 7p(2k-1) + 12q(2k-1)
q(2k+1) = 4p(2k-1) + 7q(2k-1)

Consoderando ahora r(k) = p(2k+1), s(k) = q(2k+1), tenemos que r(0) = 2,
q(0) = 1, y

r(k+1) = 7r(k) + 12s(k)
s(k+1) = 4r(k) + 7s(k)

Para obtener relaciones independientes, escribimos

r(k+2) = 7r(k+1) + 12s(k+1) = 7r(k+1) + 48r(k) + 84s(k) = 7r(k+1) + 48r(k) +
7(r(k+1) - 7r(k))

r(k+2) = 14r(k+1) - r(k),  r(0) = 2, r(1) = p(3) = 26

similarmente, se obtiene

s(k+2) = 14s(k+1) - s(k),  s(0) = 1, s(1) = q(3) = 15

La ecuación característica es entonces x^2 - 14x - 1 = 0, con las raíces  7
+/- 4rq(3). Entonces,

r(k) = A(7 + 4rq(3))^k + B(7 - 4rq(3))^k

que con r(0) = 2, r(1) = 26, nos da

r(k) = ((2 + rq(3))/2)(7 + 4rq(3))^k + ((2 - rq(3))/2)(7 - 4rq(3))^k

        = (1/2)((2 + rq(3))^(2k+1) + ((2 - rq(3))^(2k+1)

Estos son los valores de 2n. Los de n(k) serán

n(k) = (1/4)((2 + rq(3))^(2k+1) + (2 - rq(3))^(2k+1))

Hemos de ver que n(k) = 1 (mod 4), para todo k. Es decir, que

4n(k) = (2 + rq(3))^(2k+1) + (2 - rq(3))^(2k+1)

es múltiplo de 4 pero no de 8.

Pero

(2 + rq(3))^(2k+1) + (2 - rq(3))^(2k+1) =  2(2k+1)2(rq(3))^(2k) = 4*3^k = 4
(mod 8)

(q.e.d)

--
Saludos,

--

   N(n) = 14N(n-1) - N(n-2)

Construimos ahora la recurrencia

   b(0) = 0

   b(1) = 2

   b(n) = 4b(n-1) - b(n-2)

(que no sale por arte de magia, proviene de la relación de primer orden
que escribí antes para x e y).

Se trata ahora de ver que

  N(n) = 2b(n)^2 + 2b(n) + 1 = b(n)^2 + (b(n)+1)^2

Tras la cena sigo.
--

Ahora ya en serio, esto solo es una parte. Qudaría por ver que que todos los
factores de n(k) iguales a 3 (mod 4) están elevados a potencias pares. Lo
anterior es solo una condición necesaria, pero no suficiente.

--
Saludos,

     b(n) = 4b(n-1) - b(n-2) + 1

   N(0) = b(0)^2 + (b(0)+1)^2

   N(1) = b(1)^2 + (b(1)+1)^2

Ahora aplicamos inducción. Suponemos que

   N(n) = b(n)^2 + (b(n)+1)^2  (#1)

y demostraremos que

   N(n) - 14N(n-1) + N(n-2) = 0

Sustituimos la relación (#1) y nos queda el residuo

   R(n-1) = N(n) - 14N(n-1) + N(n-2) =

   = 4(b(n-1)^2 - b(n-1) +b(n-2)^2 - b(n-2)-4b(n-1)b(n-2)-2)

Veamos que esta cantidad es nula. tenemos que

   R(1) = 4(0^2 - 0 + 4 - 2 - 0 - 2) = 0

y que restando

   R(n)-R(n-1) = 4(b(n)^2 - b(n) - b(n-2)^2 + b(n-2) -

    -4b(n)b(n-1) - 4 b(n-1)b(n-2)) =

   = 4(b(n)-b(n-2))(b(n) - 4b(n-1) + b(n-2)-1) = 0    c.q.d.

ya que el paréntesis es la relación de recurrencia que cumple b(n).

Por tanto, resumiendo,

Los valores de N cumplen:

   N(0) = 1

   N(1) = 13

   N(n) = 14N(n-1) - N(n-2)

o, explícitamente,

   N(n) = ((2+rq(3))^(2n+1) + (2-rq(3))^(2n+1))/4

Estos valores se pueden poner como suma de cuadrados de dos números
consecutivos, b(n) y b(n+1), con b(n) cumpliendo

   b(0) = 0

   b(1) = 2

   b(n) = 4b(n-1) - b(n-2) + 1

o, explícitamente

   b(n) = ((2+rq(3))^(n+1/2) - (2-rq(3))^(n-1/2) - 1)/2

--