Weird solve() problem

edit

As for me, on 8.2.beta6, i got the same NotImplementedError as @dan_fulea. Without the assume statements, i got:

sage: var('a,b,c')
(a, b, c)
sage: solve([c c^2],[a,b,c],solution_dict=True)
[{a: r1, c: sqrt((2*r2^2 + 1)*r1^2 + r2^2), b: r2},
 {a: r3, c: -sqrt((2*r4^2 + 1)*r3^2 + r4^2), b: r4}]

which seems perfectly valid. That said, if you do not care about the integer nature of the solution, there is no benefit in using of a computer algebra system to solve such an equation since the variable c is isolated, it is clear that for each a and each b, there are two solutions sqrt(a^2 + b^2 + 2*a^2*b^2) and sqrt(a^2 + b^2 + 2*a^2*b^2).

The question is when a^2 + b^2 + 2*a^2*b^2) is a square. You can try to guess some pattern for the integer solutions as follows (try with size=10,100,1000,10000):

sage: from itertools import product
sage: p = product(range(size),range(size))
sage: L = []
sage: for a,b in p:
....:     if (a^2 + b^2 + 2*a^2*b^2).is_square():
....:         L.append((a,b))
sage: points(L, aspect_ratio=1)

Given those pictures, no obvious paterns seems to appear, and i would bet that there is no simple formula (like a fnite union of discrete hyperbolas or so) to describe such solution set.