/*
Bezier Curve
Elegant addition-only differencing method
.
WARNING: untested code.
2002-04-13:DAV: Translated
  from lovely idiomatic PostScript
  to lowest-common-denominator C language.
1998-02-14: Don Lancaster wrote Version 3.2 in PostScript
  and posted online at
  http://www.tinaja.com/text/bezgen3.html
  .

How the method works:
Given a point on the spline curve, the next
point on the curve can be found
by simple addition of running terms.

The method is device and language independent,
letting you generate
a cubic spline **without** using the PostScript -curveto- operator!

TODO: Is it more efficient to leave like this,
breaking up a smooth Bezier curve into 30 or so straight line segments,
or would it be cool to
make this even more like the Bresenham circle algorithm,
*only* doing putpixel() and directly sets the correct pixels
without reference to a lineto() routine ?

*/

extern void moveto( int x, int y);

extern void lineto( int x, int y);

void
bezier(
	const int x0, const int y0, // start point
	const int x1, const int y1, // first control point
	const int x2, const int y2, // last control point
	const int x3, const int y3  // end point
	){


	const int numincs = 100;
	int i = numincs;
// % Plot the spline as a standard PostScript curveto...
// x0 y0 moveto x1 y1 x2 y2 x3 y3 curveto line1 stroke

// Find t-power coefficients a-d given control points. See
//   HACK62.PDF http://www.tinaja.com/glib/hack62.pdf
// for more details. This moves you from graph space to math space..

const float ax = x3 - x2*3 + x1*3 - x0;
const float ay = y3 - y2*3 + y1*3 - y0;

const float bx = x2*3 - x1*6 + x0*3;
const float by = y2*3 - y1*6 + y0*3;

const float cx = x1*3 - x0*3;
const float cy = y1*3 - y0*3;

// /dx x0 def
// /dy y0 def


// Elegant addition-only differencing method


	// /numincs 100 def % number of steps to generate the spline

	// % next, precalculate delta-t increment, square, and cube

	const float dt1 = 1/(float)numincs; // delta t
	const float dt2 = dt1*dt1; // delta t squared
	const float dt3 = dt1*dt2; // delta t cubed

	// % initialize some stuff...
	float xx = (float)x0;
	float yy = (float)y0;
	float ux = 0;
	float uy = 0;
	float vx = 0;
	float vy = 0;

	float mx1 = ax * dt3;
	float my1 = ay * dt3;

	float lx1 = bx * dt2;
	float ly1 = by * dt2;

	float kx = mx1 + lx1 + cx * dt1;
	float ky = my1 + ly1 + cy * dt1;

	float mx = mx1 * 6;
	float my = my1 * 6;

	float lx = mx + 2*lx1;
	float ly = my + 2*ly1;

	// % now generate the fake curveto...

	moveto( xx, yy ); // xx yy moveto

	for( i = numincs; i; i-- ){
		xx += ux + kx;
		yy += uy + ky;
		ux += vx + lx;
		uy += vy + ly;
		vx += mx;
		vy += my;

		lineto( xx, yy ); // xx yy lineto

	};

	/*
	Note that successive points along the curve require only *ten* repeated
	additions. There are no multiplies or other time-intensive operators.
	Note that this is **not** an approximation, but an exact solution!

	Try changing -numincs- to 3 to prove this to yourself.
*/

};