Bresenham Line and Circle Drawing
W. Randolph Franklin (WRF)
Bresenham Algorithm  Optimized Line Drawing Code
We want to draw a line from (0,0) to (x1,y1), where 0<=y1<=x1, by setting one pixel per column. For example, if x=10, y=7, we get this:
Here are several versions, ending with Bresenham's algorithm. The point of this is to use the simplest possible operations.

Simple program:
m=y1/x1; pixel(0,0); for(x=1;x<=x1;x++) { y=round(m*x); pixel(x,y); }

Track the deviation of the current y from the perfect line. When the deviation gets too large, increment y, and decrement the deviation.
m=y1/x1; d=0; pixel(0,0); y=0; for(x=1;x<=x1;x++) { d+= m; if (d>= 1/2) { d = 1; y++; } pixel(x,y); }
 Scale up m and d by 2x to convert fractions to
integers. This is ok since they are used only in the if;
and it still branches at the same times as before.
m=2*y1; d=0; pixel(0,0); y=0; for(x=1;x<=x1;x++) { d+= m; if (d>=x1) { d= 2*x1; y++;} pixel(x,y); }

Shift d by x1 to remove a subtraction.
m=2*y1; d= x1; pixel(0,0); y=0; for(x=1;x<=x1;x++) { d+= m; if (d>=0) { d= 2*x1; y++;} pixel(x,y); }

Assume that the pixels are stored in a framebuffer, and that pixel_{x y} is stored at location fb+x*maxy+y where fb is the address of the start of the framebuffer, and maxy is the max value of y in the framebuffer, i.e., 0<= y<= y1<= maxy. Then we get:
m=2*y1; d= x1; *fb=1; y=0; for(x=1;x<=x1;x++) { d+= m; if (d>=0) { d= 2*x1; y++;} *(fb+x*maxy+y)=1; }

Letting p be the current address, we can optimize.
m=2*y1; d= x1; p=fb; *p=1; for(x=1;x<=x1;x++) { d+= m; p+=maxy; if (d>=0) { d= 2*x1; p++;} *p=1; }

The next problem is that modern computers use a pipelined architecture, where different stages of several instructions are processed in parallel. However, when there is a conditional, the computer can't tell what the next instructions are, and so is forced to flush the pipeline. Therefore we want to eliminate the conditional in the inner loop. The following assumes that the right shift is arithmetic, which it is on the Sun.
Mask is all ones when d>=0, all zeroes otherwise.
We also pulled 2*x1 out of the loop.
m=2*y1; d= x1; p=fb; *p=1; doublex1= x1<1; for(x=1;x<=x1;x++) { d+= m; p+=maxy; mask = (~(d>>31)); d = (doublex1&mask); p += (1&mask); *p=1; }

We can also offset x by x1 so that the termination condition is a test against 0.
Bresenham Circle Algorithm
Here are several ways to draw a circle, ending with the Bresenham circle algorithm.

Taylor series for sqrt(1x^{2}), in the 2nd octant: 0<= x<= y, where the curve doesn't get too vertical. The following is expanded about the origin. The error at x=.71 is 0.0003.
y(x) = 1  (1/2) x^{2}  (1/8) x^{4}  (1/16) x^{6}  (5/128) x^{8}  (7/256) x^{10}  (21/1024) x ^{12}+ ...
This is popular but usually the wrong way.

If you must use Taylor, at least expand about the center of the interval:
z = x 0.3536
y(z) = ( 0.9353967287 0.3780214203 z 0.6109173568 z^{2}  0.2468897312 z^{3}  0.2992736736 z^{4}  0.2821915805 z^{5}  0.3420828869 z^{6}  0.4015387039 z^{7}  0.5080487872 z^{8}  0.6481416197 z^{9}  0.8517407025 z^{10}  1.133286086 z^{11}  1.531581071 z^{12}  2.091410948 z^{13}+...
The error at x=.71 is 0.000 004. This is better, but there are more nonzero coefficients, which means that it is slower.

A Chebyshev series, which tries to minimize the maximum error over an interval beats a Taylor series:
y(x) = + 1.000000093  0.00002511816061 x  0.4988901592 x^{2}  0.01882886918 x^{3} + 0.03498610076 x^{4}  0.7655182976 x^{5} + 2.101697099 x^{6}  3.603814683 x^{7} + 3.259950635 x^{8}  1.315059748 x^{9}
This has a max error of 1.5*10^{7}, which is much better than above, expecially as there are only 8 coefficients. We could get as good as either Taylor with many fewer terms than they.

Rotate the point (R,0) by Dt repeatedly (t is the angle). What should Dt be?

Approx the rotation matrix:
cos t sin t  sin t cos t 1 t t 1 with everything multiplied by R. But this circle may grow since the determinant >1.

Improve the approx matrix so the determinant =1:
1 t t 1t^{2} 
Binomial theorem on sqrt(1x^{2}). Does this do the same as Taylor expanded about 0?

Here are two parametric ways to draw circles.

x = R cos t
y = R sin tWhat should Dt be?

x = R (2t)/(t^{2}+1}
y = R (t^{2}1)/(t^{2}+1) ; infinity<t<infinityUnfortunately, the speed of the curve wrt t varies with t, so the stepsize must also.

There is no parametric polynomial, of any degree, that exactly represents a circle. However you can get good approximations.

x = R cos t

???? (You suggest)

Now let's optimize C code. Do the second octant only.
pixel(0,r); for(x=1;x<r/sqrt(2);x++) { y=sqrt(r^{2}x^{2}); pixel(x,y); }

Assume that the last point was (x1,y). Find whether (x,y½) is above or below the circle, and move SE to (x,y1) or E to (x,y), respectively.
y=r; pixel(0,r); for(x=1;x<r/sqrt(2);x++) { d=x^{2}+(y1/2)^{2}r^{2}; if (d>=0) y; pixel(x,y); }

Optimize: d=x^{2}+(y1/2)^{2}r^{2} becomes d=x^{2}+y^{2}y+1/4r^{2}. How does the d for a given x differ from d for x1? If y is unchanged, then d_{x} = d_{x1}+2x1. If y decreases, add 2y+2, where that is the old y. Also the only time that 1/4 affects the test d>0 is when d=1/4, so ignore it, but make the test d>=0.
y=r; d= r; pixel(0,r); for(x=1;x<r/sqrt(2);x++) { d+= 2x1; if (d>=0) { y; d = 2y; /* Must do this AFTER y */ } pixel(x,y); }
This is Bresenham's circle algorithm. Note that you can store 2x1 and 2y and update them by adding or subtracting 2.

Lolevel optimizations can also be done thus:
y=r; d= r; x2m1= 1; pixel(0,r); for(x=1;x<r/sqrt(2);x++) { x2m1 += 2; d+= x2m1; if (d>=0) { y; d = (y<<1); /* Must do this AFTER y */ } pixel(x,y); }

You can also do the same pipelining techniques used for the line algorithm.

Finally, on an Intel processor, there are ways to use the available instructions efficiently.