can all triangle numbers that are squares be expressed as sum of squares I'm not sure if this is just a subset of Which integers can be expressed as a sum of squares of two coprime integers? which in turn points to http://en.wikipedia.org/wiki/Brahmagupta%E2%80%93Fibonacci_identity, but if so, I'm not seeing it.  
Basically: looking at the numbers between 0 and 1000, if (but not iff) n is a square, then the nth triangle number (i.e. $$\frac{n(n+1)}{2}$$) can be expressed as the sum of two perfect squares .  Does this hold for all squares, and can someone point me in the direction of why this is?  (For what it's worth, I is an engineer but haven't really touched number theoryish stuff since college.)  
(Also saw Prove that there are infinitely many natural numbers $n$, such that $n(n+1)$ can be expressed as sum of two positive squares in two distinct ways. which proves there are infinitely many for the similar $n(n+1)$ case but not that all squares will work unless I missed some aspect)
(Brute force approach to checking numbers, because brute force always works)
import math
maxsquare = 1001
squares = [i*i for i in xrange(maxsquare)]
for j in xrange(int(math.sqrt(maxsquare))):
    i = j * j 
    n = i * (i+1) / 2
    for s in squares:
        f = n - s
        if f in squares:
            print "%d^2 + %d^2 = %d * %d+1 / 2" % (int(math.sqrt(f)),int(math.sqrt(s)),i,i)
            break
    else:
        print "Could not find anything for i = %d" % i

 A: The answer is yes, it is true for all $n$.
Let me reformulate your question. You want to know if, whenever $n$ is a perfect square, it is the case that $\frac12 n(n+1)$ is a sum of two squares.  If $n$ is a perfect square, then it has an integer square root, which we can call $m$, and write $n = m^2$.  Then your question becomes whether $\frac12 m^2(m^2+1)$ is a sum of two squares for all integers $m$.
There is a very useful theorem about sums of two squares which says that a number $N$  is a sum of two squares if and only if every prime of the form $4k+3$ (such as $3, 7, 11, 19,$ etc.) appears in the prime factorization of $N$ an even number of times. 
Now $m^2(m^2 + 1) = {\left(m^2\right)}^2 + m^2$ is obviously a sum of two squares.  And dividing $m^2(m^2+1)$ by 2 can't possibly affect the number of times any prime of the form $4k+3$ appears in its factorization, so $\frac12 m^2(m^2+1)$ is also a sum of two squares.
(Addendum: Erick Wong points out below that we do not even need to invoke the Fermat $4k+3$ theorem.)
A: If you run your program you should find
$$
0^2+1^2 = \frac{1\times 2}{2}\\
1^2+3^2 = \frac{4\times 5}{2}\\
3^2+6^2 = \frac{9\times 10}{2} \\
6^2+10^2 = \frac{16\times 17}{2} \\
10^2+15^2 = \frac{25\times 26}{2}
$$
from which you might recognize the numbers in the sums as the triangular numbers again, in which case you can guess that
$$
\frac{j^2(j^2+1)}{2} = \left(\frac{j(j-1)}{2}\right)^2+\left(\frac{j(j+1)}{2}\right)^2
$$
which is true.
