In[]:= AbsoluteTiming[
 sol = Solve[
   {x^2 + y^2 == z^2, x + y + z == 1000, x > 0, y > x, z > y},
   {x, y, z}, Integers];
 a = x /. sol[[1, 1]];
 b = y /. sol[[1, 2]];
 c = z /. sol[[1, 3]];
 ans = a* b* c]

Out[]= {0.0548628, 31875000}