I don't recall seeing the simple observation that this sequence is trivially constructed from the successive approximants to sqrt(2): Given a rational approximation a/b, the fraction 2 b/a errs on the opposite side, and their average would be the next Newton approximation. (Starting with a=b=1.) These form an exponentially sparse subsequence of the continued fraction (1+/2+/2+/2+/...) approximants resulting from taking the mediant, (a + 2 b)/(b + a), instead of the average. E.g.,
a 1 3 7 17 41 ... b 1 2 5 12 29 ...
Then the square triangular numbers are just a^2 b^2 = binom(a^2+(n mod 2),2):
= (b[2n]/2)^2, where the nth a and b are a[n] and b[n]. I.e., b[2n] = 2 a[n] b[n].
(c607) define(inversetriang(n),rhs(solve(binomial(x+1,2) = n,x)[2]))
sqrt(8 n + 1) - 1 (d607) inversetriang(n) := ----------------- 2
(c608) block([m : matrix([1,1],[1,0]),fancy_display : false], thru 9 do (print(m,(m[1,1]/m[2,1])^2,inversetriang((m[1,1]*m[2,1])^2)),m:m.matrix([2,1],[1,0])))$
[a a'] -1 [ ] a^2/b^2 triang (a^2 b^2) = a^2-(n+1 mod 2) [b b'] --------------------------------------------------------
[ 1 1 ] 1 [ ] - 1 [ 1 0 ] 1
[ 3 1 ] 9 [ ] - 8 [ 2 1 ] 4
[ 7 3 ] 49 [ ] -- 49 [ 5 2 ] 25
[ 17 7 ] 289 [ ] --- 288 [ 12 5 ] 144
[ 41 17 ] 1681 [ ] ---- 1681 [ 29 12 ] 841
[ 99 41 ] 9801 [ ] ---- 9800 [ 70 29 ] 4900
[ 239 99 ] 57121 [ ] ----- 57121 [ 169 70 ] 28561
[ 577 239 ] 332929 [ ] ------ 332928 [ 408 169 ] 166464
[ 1393 577 ] 1940449 [ ] ------- 1940449 [ 985 408 ] 970225
where
[a a'] 1 1 2 1 n [ ] = ( )( ) and a'/b' is the previous a/b. [b b'] 1 0 1 0
--rwg