Comment by Iv

I think that the key lies in the fact that we are creating a numerical solution to a differential equation: y = y - 1/16x actually means that we add an amount to y that is dependent on the value in x, i mathematical terms, it means that the derivative (approximated here) of y is equal to deltat(1/16)x

Assuming deltat=16, we have:

y' = -x

x' = y

Which is solved with y=cos and x=sin or with x=-cos and y=sin.

The fact that one line adds and the other substracts may be explain why rounding errors will compensate each other over time, but we need a more detailed working of the PDP arithmetics to be really sure about that.

Actually, many CS students may recognize a prey-and-predator model, that creates very simply sinusoids that are dephased of pi/2. That is, if one is a sine, the other is a cosine.