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].
I should've summarized this as A001110 = A001333^2 * A000129^2 = A000129[2n]^2/4 = binom(A001108,2). Mike Hirschhorn apparently found most or all of this in 1996, and refers me to ". . .Tom Beldon and Tony Gardiner, ``Triangular numbers and perfect squares'', The Mathematical Gazette, 2002, pp423--431, esp pp424--426. {This is a British publication for senior high--school students and their teachers.]" I still can't prove the polylog puzzles--I'm just detecting them with Rich's old lattice reducer. Ramanujan apparently found a similar relation between Li_2(1/3) and Li_2(-1/3). --rwg