Thanks to comments, it should be a plane but why does it look a bit like a fractal? Does my code overlook something or some err in plotting tool? I used Python and GNUplot.

Apparently an animated illusion.
And generated the points with Python and generated graph with GnuPlot like:
$ python copyPasteTheCodeToFile.py > .data
$ gnuplot -e "set terminal png; set grid; splot '.data'"
Python code
import itertools
def main():
density = 100.0
mySet= [x/density for x in range(int(density))]
points=""
for (x,y,z) in [(x,y,z) for x,y,z in itertools.product(mySet, repeat=3) if x+y+z==1]:
points = points + "%s\t%s\t%s\n"%(str(x),str(y),str(z))
print points
main()
can you spot an err? I cannot hence I thought A) mathematical err, B) algoritmic err or C) outside-party/program-err, cannot really say which one.
[Update]
It is a floating point error! The error is the conditional x+y+z==1. Since computers don't evaluate floats like 1.0 + 1.0 to 2.0 but something slightly different, the conditional fails with some right points. The fix is:
$abs(x+y+z-1) < \epsilon$
where the epsilon $\epsilon$ is some very small number like 1e-10. It is still open why there are just the holes in the grid. Why not other points? And why is it vertically symmetric? Is it so with every computer? Does the pattern vary between computers? Anyway I am still investigating why there this specific pattern with this conditional. We know now why there are the patters but we don't know why these patterns. I am running the commands with bulk i86 comp.
I suspect that the error here comes down to the use of floating-point arithmetic (which has a rather fascinating mathematical structure all its own). Your code is iterating over all of the points $(x, y, z)$ and doing a precise compare of the sum to 1.0; but, for instance, it may well be the case that the triad $(30.0/100.0, 30.0/100.0, 40.0/100.0)$ won't be in your set simply because the floating point representations of .3 and .4 aren't exact and so the sum $.3 + .3 + .4$ comes out to, say, $1.00000001$. While I can't speak specifically to the fractal structure you're seeing, I suspect it's loosely related to the representation of the Sierpinski Gasket in terms of the set of all $(x, y)$ with $x$ AND $y == 0$; see http://ecademy.agnesscott.edu/~lriddle/ifs/siertri/binary.htm for more details on this. The easiest way to fix this would be to just iterate over all points $(x, y)$ and define $z$ in terms of $x$ and $y$ by $z=1.0-(x+y)$; this also has the advantage of speeding up your generation a hundredfold (by eliminating the loop over $z$). More generally (for instance, in circumstances where there isn't an explicit equation for $z$ in terms of $x$ and $y$), floating point compares should almost never be exact, and a test for $(1.0-\epsilon)\lt x+y+z \lt (1.0+\epsilon)$ (for some small $\epsilon$) is a better bet.
EDIT to answer the new parts of the question that came in while I was answering: in general, floating-point representations are likely to be the same across platforms; the IEEE has a standard specification for floating-point numbers and the arithmetic using them, which most platforms support. There are corner cases where specific platforms might differ, but they're unlikely to be exercised by simple additions like this.