#include <math.h>

typedef struct { double x, y; } PT;
typedef double F;

extern void comment(const char *message, ...);
extern void moveto(F x, F y);
extern void lineto(F x, F y);

#define epsilon (.01)

static inline F SQ(F a) { return a*a; }

static F dist_p2l(F x0, F y0, F x1, F y1, F xa, F ya) {
    x1 -= x0;
    y1 -= y0;
    xa -= x0;
    ya -= y0;

    return fabs(x1 * ya + xa * y1) / hypot(x1,y1);
}

static void bezier_inner(PT p1, PT p2, PT p3, F t0, F t1) {
    F ta = (t1+t0)/2;
    F
        x0 = SQ(1-t0) * p1.x + 2*t0*(1-t0) * p2.x + SQ(t0) * p3.x,
        y0 = SQ(1-t0) * p1.y + 2*t0*(1-t0) * p2.y + SQ(t0) * p3.y,
        x1 = SQ(1-t1) * p1.x + 2*t1*(1-t1) * p2.x + SQ(t1) * p3.x,
        y1 = SQ(1-t1) * p1.y + 2*t1*(1-t1) * p2.y + SQ(t1) * p3.y,
        xa = SQ(1-ta) * p1.x + 2*ta*(1-ta) * p2.x + SQ(ta) * p3.x,
        ya = SQ(1-ta) * p1.y + 2*ta*(1-ta) * p2.y + SQ(ta) * p3.y;
    F d;
    if(t1 != t0 && (d = dist_p2l(x0,y0,x1,y1,xa,ya)) > epsilon) {
        comment("subdivide d=%.2f", d);
        bezier_inner(p1, p2, p3, t0, ta);
        lineto(xa, ya);
        bezier_inner(p1, p2, p3, ta, t1);
    } else {
        if(t0 != 0) {
            lineto(x0, y0);
        }
    }
}

void bezier(PT p1, PT p2, PT p3, long int offset) {
    p1.x += offset;
    p2.x += offset;
    p3.x += offset;
    moveto(p1.x, p1.y);
    bezier_inner(p1, p2, p3, 0., 1.);
    lineto(p3.x, p3.y);
}