For comparison with other CAS choose from: Axiom Derive Macsyma Maple Reduce
    Factor[ a ^2 - b^2 ]
(a - b) (a + b)
Factor[ a
^2 + b^2 ]
     2    2
    a  + b
 Factor[a
^2 + b^2, GaussianIntegers -> True ]
(-I a + b) (I a + b)
fa := ( x
^2 z + y^4 z^2 + 5) * ( x y^3 + z^2) * (-x^3 y + z^2 + 3) * ( x^3 y^4 + z^2 )fa = Expand[fa]
        4  7      7  8      6  7      9  8           3  2
    15 x  y  - 5 x  y  + 3 x  y  z - x  y  z + 15 x y  z  +
         3  4  2      4  4  2      6  5  2      4  7  2
     15 x  y  z  - 5 x  y  z  - 5 x  y  z  + 5 x  y  z  +
        4  11  2    7  12  2      3  3  3      5  4  3
     3 x  y   z  - x  y   z  + 3 x  y  z  + 3 x  y  z  -
      6  4  3    8  5  3    6  7  3       4      3    4
     x  y  z  - x  y  z  + x  y  z  + 15 z  - 5 x  y z  +
          3  4      3  4  4        7  4      3  8  4
     5 x y  z  + 5 x  y  z  + 3 x y  z  + 3 x  y  z  -
      4  8  4    6  9  4    4  11  4      2  5    5    5
     x  y  z  - x  y  z  + x  y   z  + 3 x  z  - x  y z  +
      3  3  5    5  4  5      6      4  6    3  5  6
     x  y  z  + x  y  z  + 5 z  + 3 y  z  - x  y  z  +
        7  6    3  8  6    2  7    4  8
     x y  z  + x  y  z  + x  z  + y  z
 Factor[ fa ]
          3      2      3    2    3  4    2
    (3 - x  y + z ) (x y  + z ) (x  y  + z )
           2      4  2
     (5 + x  z + y  z )