The circle algorithm seems to come from the trigonometric identities with small value approximations for b -> 0. ie. cos(b)~=1 and sin(b)~=b.
sin(a+b)=sin(a)cos(b)+cos(a)sin(b)
cos(a+b)=cos(a)cos(b)-sin(a)sin(b)
sin(a+b)=sin(a)cos(b)+cos(a)sin(b)
cos(a+b)=cos(a)cos(b)-sin(a)sin(b)
But the mystery is that experimentally it's more correct than that. Why does it seem to join up with itself exactly after 360°, and not spiral slightly inward or outward due to accumulated error? Why does it stop working if you use a temporary variable so that x and y are updated simultaneously instead of alternately?