diff --git a/c++/vector/activearea.h b/c++/vector/activearea.h
index 9a32217..f1d8af4 100644
--- a/c++/vector/activearea.h
+++ b/c++/vector/activearea.h
@@ -20,7 +20,9 @@ private:
 	ActivePoint::Handle drag_point;
 	Vector drag_offset;
 
-protected:
+public:
+	ActiveArea();
+
 	const PointList& get_points() const { return points; }
 	bool is_point_added(const ActivePoint::Handle &point) const;
 	void add_point(const ActivePoint::Handle &point);
@@ -29,6 +31,7 @@ protected:
 
 	ActivePoint::Handle get_point_at(const Vector &position) const;
 
+protected:
 	virtual void on_point_move(const ActivePoint::Handle &point, const Vector &oldposition, const Vector &newposition);
 	virtual void on_draw(const Context &context);
 	virtual void on_draw_content(const Context &context);
@@ -37,9 +40,6 @@ protected:
 	bool on_button_press_event(GdkEventButton* event) override;
 	bool on_button_release_event(GdkEventButton* event) override;
 	bool on_motion_notify_event(GdkEventMotion* event) override;
-
-public:
-	ActiveArea();
 };
 
 
diff --git a/c++/vector/activecurve.h b/c++/vector/activecurve.h
new file mode 100644
index 0000000..3e2160e
--- /dev/null
+++ b/c++/vector/activecurve.h
@@ -0,0 +1,161 @@
+#ifndef ACTIVECURVE_H
+#define ACTIVECURVE_H
+
+
+#include <vector>
+
+#include "curve.h"
+#include "activearea.h"
+
+class ActiveCurvePoint {
+public:
+	ActivePoint::Handle p, t0, t1;
+	
+	ActiveCurvePoint(const Vector &p, const Vector &t0, const Vector &t1):
+		p(new ActivePoint(p)),
+		t0(new ActivePoint(p + t0, Color::yellow())),
+		t1(new ActivePoint(p + t1, Color::yellow()))
+		{ }
+	explicit ActiveCurvePoint(const Vector &p = Vector(), const Vector &t = Vector()):
+		ActiveCurvePoint(p, t, t) { }
+
+	const Vector& get_p() const { return p->position; }
+	Vector get_t0() const { return t0->position - get_p(); }
+	Vector get_t1() const { return t1->position - get_p(); }
+	Vector get_pp0() const { return get_t0()/3 + get_p(); }
+	Vector get_pp1() const { return get_t1()/3 + get_p(); }
+	
+	void add_to(ActiveArea &area) const {
+		area.add_point(p);
+		area.add_point(t0);
+		area.add_point(t1);
+	}
+	void remove_from(ActiveArea &area) const {
+		area.remove_point(p);
+		area.remove_point(t0);
+		area.remove_point(t1);
+	}
+};
+
+class ActiveCurve: public std::vector<ActiveCurvePoint> {
+public:
+	template<typename T>
+	class SegIterBase {
+	public:
+		typedef T Type;
+	private:
+		Type begin, end, last, i0, i1;
+		bool loop, valid;
+	protected:
+		SegIterBase():
+			begin(), end(), i0(), i1(), loop(), valid() { }
+		SegIterBase(const Type &begin, const Type &end, bool loop):
+			begin(begin), end(end), last(end), i0(begin), i1(begin), loop(loop), valid(true)
+		{
+			--last; ++i1;
+			if (i0 == end) i1 = end;
+			if (loop && i1 == end) i1 = begin;
+		}
+	public:
+		operator bool () const
+			{ return valid && i0 != end && i1 != end; }
+		bool is_first() const
+			{ assert(*this); return i0 == begin; }
+		bool is_last() const
+			{ assert(*this); return i1 == (loop ? begin : last); }
+	
+		bool next() {
+			assert(*this);
+			++i0; ++i1;
+			if (loop && i1 == end) i1 = begin;
+			return *this;
+		}
+		
+		const ActiveCurvePoint& point0() const { assert(*this); return *i0; }
+		const ActiveCurvePoint& point1() const { assert(*this); return *i1; }
+		const Vector& p0() const { return point0().get_p(); }
+		const Vector& p1() const { return point1().get_p(); }
+	};
+
+	class SegIter: public SegIterBase<const_iterator> {
+	public:
+		SegIter() { }
+		explicit SegIter(const ActiveCurve &curve):
+			SegIterBase(curve.begin(), curve.end(), curve.loop) { }
+		Vector t0()  const { return point0().get_t1(); }
+		Vector t1()  const { return point1().get_t0(); }
+		Vector pp0() const { return point0().get_pp1(); }
+		Vector pp1() const { return point1().get_pp0(); }
+		Hermite hermite() const { return Hermite(p0(), p1(), t0(), -t1()); }
+		Bezier  bezier()  const { return Bezier (p0(), p1(), pp0(), pp1()); }
+	};
+
+	class SegRIter: public SegIterBase<const_reverse_iterator> {
+	public:
+		SegRIter() { }
+		explicit SegRIter(const ActiveCurve &curve):
+			SegIterBase(curve.rbegin(), curve.rend(), curve.loop) { }
+		Vector t0()  const { return point0().get_t0(); }
+		Vector t1()  const { return point1().get_t1(); }
+		Vector pp0() const { return point0().get_pp0(); }
+		Vector pp1() const { return point1().get_pp1(); }
+		Hermite hermite() const { return Hermite(p0(), p1(), t0(), -t1()); }
+		Bezier  bezier()  const { return Bezier (p0(), p1(), pp0(), pp1()); }
+	};
+
+	bool loop;
+	Color color;
+	Color tangent_color;
+	
+	ActiveCurve(): loop(), color(Color::black()), tangent_color(Color::gray()) { }
+	
+	void set_tail_tangents_visible(bool x)
+		{ if (!empty()) front().t0->visible = back().t1->visible = x; }
+	void set_tail_tangents_visible()
+		{ set_tail_tangents_visible(loop); }
+	
+	void put(const Context &context, bool jump = false) const {
+		for(SegIter i(*this); i; i.next())
+			context.hermite( i.hermite(), jump && i.is_first() );
+		if (loop && jump) context->close_path();
+	}
+
+	void put_tangents(const Context &context) const {
+		for(SegIter i(*this); i; i.next()) {
+			context.move_to(i.p0());
+			context.line_to(i.p0() + i.t0());
+			context.move_to(i.p1());
+			context.line_to(i.p1() + i.t1());
+		}
+	}
+	
+	void on_point_move(const ActivePoint::Handle &point, const Vector &oldposition, const Vector &/*newposition*/) {
+		for(const_iterator i = begin(); i != end(); ++i) {
+			if (i->p == point) {
+				Vector d = i->p->position - oldposition;
+				i->t0->position += d;
+				i->t1->position += d;
+				break;
+			}
+		}
+	}
+
+	void draw(const Context &context) const {
+		context->save();
+		
+		context.set_line_width_px(1);
+		
+		context.set_source_rgba(color);
+		put_tangents(context);
+		context->stroke();
+
+		context.set_source_rgba(tangent_color);
+		put(context, true);
+		context->stroke();
+		
+		context->restore();
+	}
+};
+
+
+#endif
diff --git a/c++/vector/bendarea.cpp b/c++/vector/bendarea.cpp
new file mode 100644
index 0000000..bdf160b
--- /dev/null
+++ b/c++/vector/bendarea.cpp
@@ -0,0 +1,457 @@
+
+#include <vector>
+
+#include "random.h"
+#include "matrix.h"
+
+#include "bendarea.h"
+
+
+namespace {
+	enum BendMode {
+		BendNone,
+		BendRound,
+		BendCorner
+	};
+	
+	class ContourChunk {
+	public:
+		Vector p1, t0, t1;
+		explicit ContourChunk(const Vector &p1 = Vector(), const Vector &t0 = Vector(), const Vector &t1 = Vector()):
+			p1(p1), t0(t0), t1(t1) { }
+	};
+
+	class BendPoint {
+	public:
+		Vector p;
+		Vector t0, t1;
+		bool e0, e1;
+		Real l;
+		Real length;
+		BendMode mode;
+		BendPoint(): e0(), e1(), l(), length(), mode(BendNone) { }
+	};
+	
+	class Intersection {
+	public:
+		Real l;
+		BendPoint point;
+		explicit Intersection(const Real &l = Real(), const BendPoint &point = BendPoint()):
+			l(l), point(point) { }
+		bool operator<(const Intersection &other) const { return l < other.l; }
+	};
+	
+	typedef std::vector<Intersection> IntersectionList;
+	
+	class Contour: public std::vector<ContourChunk> {
+	public:
+		void hermite_to(const Vector &p1, const Vector &t0, const Vector &t1)
+			{ push_back( ContourChunk(p1, t0, t1) ); }
+		void line_to(const Vector &p1) {
+			if (empty()) {
+				hermite_to(p1, Vector(), Vector());
+			} else
+			if (p1 != back().p1) {
+				Vector t = p1 - back().p1;
+				hermite_to(p1, t, t);
+			}
+		}
+
+		void line(const Vector &p0, const Vector &p1)
+			{ line_to(p0); line_to(p1); }
+		void hermite(const Vector &p0, const Vector &p1, const Vector &t0, const Vector &t1)
+			{ line_to(p0); hermite_to(p1, t0, t1); }
+		void hermite(const Hermite &h)
+			{ hermite(h.p0, h.p1, h.t0, h.t1); }
+		
+		void put(const Context &context, bool jump = false) const {
+			if (empty()) return;
+			const_iterator i = begin();
+			Hermite h(i->p1, i->p1);
+			for(++i; i != end(); ++i) {
+				h.p0 = h.p1;
+				h.p1 = i->p1;
+				h.t0 = i->t0;
+				h.t1 = i->t1;
+				context.hermite(h, jump);
+				jump = false;
+			}
+		}
+	};
+	
+	class BendContour: public std::vector<BendPoint> {
+	public:
+		const_iterator find(const Real &length) const {
+			if (empty()) return end();
+			const_iterator a = begin(), b = end() - 1, c;
+			if (!less(a->length, length)) return a;
+			if (!less(length, b->length)) return b;
+			while( (c = a + (b - a)/2) != a ) {
+				if (equal(length, c->length)) return c;
+				(length < c->length ? b : a) = c;
+			}
+			return c;
+		}
+		
+		const BendPoint* find_exact(const Real &length) const {
+			const_iterator i = find(length);
+			return i == end() || !equal(length, i->length) ? nullptr : &*i;
+		}
+
+		BendPoint interpolate(const Real &length) const {
+			const_iterator i = find(length);
+			if (i == end()) return BendPoint();
+			
+			if (equal(length, i->length)) return *i;
+			
+			if (length < i->length) {
+				// extarpolation before
+				Real dlen = length - i->length;
+				BendPoint s;
+				s.p = i->p + i->t0.norm()*dlen;
+				s.t0 = s.t1 = i->t0.norm();
+				s.l = i->l + dlen;
+				s.length = length;
+				s.e0 = s.e1 = i->e0;
+				return s;
+			}
+			
+			const_iterator j = i + 1;
+			if (j == end()) {
+				// extarpolation after
+				Real dlen = length - i->length;
+				BendPoint s;
+				s.p = i->p + i->t1.norm()*dlen;
+				s.t0 = s.t1 = i->t1.norm();
+				s.l = i->l + dlen;
+				s.length = length;
+				s.e0 = s.e1 = i->e1;
+				return s;
+			}
+			
+			// interpolation
+			Real l = (length - i->length)/(j->length - i->length);
+			Real dl = j->l - i->l;
+			Hermite h(i->p, j->p, i->t1*dl, j->t0*dl);
+			
+			BendPoint s;
+			s.p = h.p(l);
+			s.t0 = s.t1 = h.t(l);
+			s.l = i->l + l*(j->l - i->l);
+			s.length = length;
+			s.e0 = s.e1 = (i->e1 && j->e0);
+			
+			return s;
+		}
+		
+		void half_corner(Contour &dst, const BendPoint &point, const Real &radius, bool flip, bool out) const {
+			if (point.mode == BendNone || equal(radius, 0))
+				return;
+
+			const Vector &t0 = flip ? point.t1 : point.t0;
+			const Vector &t1 = flip ? point.t0 : point.t1;
+			const Vector p0 = t0.perp().norm() * radius;
+			const Vector p1 = t1.perp().norm() * radius;
+			const Vector &p = out ? p1 : p0;
+			const Vector &center = point.p;
+			const Real rmod = fabs(radius);
+			
+			Real d = dot( t0, t1.perp() );
+			if (!equal(0, d) && (d*radius < 0) == flip) {
+				if (point.mode == BendRound) {
+					const Vector pp = (p0 + p1).norm()*rmod;
+					const Vector ppp = (p + pp).norm()*rmod;
+					Real k = (ppp - p).len()/rmod;
+					if (dot(p, ppp.perp()) < 0) k = -k;
+					if (out) {
+						dst.line_to(center + pp);
+						dst.hermite_to(center + ppp, pp.perp()*k, ppp.perp()*k);
+						dst.hermite_to(center + p, ppp.perp()*k, p.perp()*k);
+					} else {
+						k = -k;
+						dst.line_to(center + p);
+						dst.hermite_to(center + ppp, p.perp()*k, ppp.perp()*k);
+						dst.hermite_to(center + pp, ppp.perp()*k, pp.perp()*k);
+					}
+				} else
+				if (point.mode == BendCorner) {
+					Vector pp = p0 + t0*(dot(p1 - p0, t1.perp())/d);
+					if (out) dst.line(center + pp, center + p);
+						else dst.line(center + p, center + pp);
+				}
+			} else {
+				Vector pp = (p0 + p1)*0.5;
+				if (out) dst.line(center + pp, center + p);
+				    else dst.line(center + p, center + pp);
+			}
+		}
+
+		void bend(Contour &dst, const Contour &src, const Matrix &matrix = Matrix()) const {
+			if (empty() || src.empty())
+				dst.insert(dst.end(), src.begin(), src.end());
+
+			Contour::const_iterator i = src.begin();
+			Vector p0 = matrix.transform(i->p1);
+			size_t dst_size = dst.size();
+			
+			IntersectionList intersections;
+			while(++i != src.end()) {
+				Hermite h(p0, matrix.transform(i->p1), matrix.transform(i->t0, false), matrix.transform(i->t1, false));
+				p0 = h.p1;
+				
+				Real bends[2];
+				int bends_count = h.bends_x(bends);
+				Range r(h.p0.x);
+				r.expand(h.p1.x);
+				for(Real *i = bends, *end = i + bends_count; i != end; ++i)
+					r.expand( Hermite::p(*i, h.p0.x, h.p1.x, h.t0.x, h.t1.x) );
+				bool b0 = bends_count > 0, b1 = bends_count > 1;
+				
+				for(const_iterator bi = find(r.a0); bi != end() && !less(r.a1, bi->length); ++bi) {
+					Real roots[3];
+					int count = h.intersections_x(roots, bi->length);
+					for(Real *i = roots, *end = i + count; i != end; ++i) {
+						if (b0 && equal(*i, bends[0])) b0 = false;
+						if (b1 && equal(*i, bends[1])) b1 = false;
+						intersections.push_back( Intersection(*i, *bi) );
+					}
+				}
+				if (b0) intersections.push_back( Intersection(bends[0], interpolate( h.p(bends[0]).x )) );
+				if (b1) intersections.push_back( Intersection(bends[1], interpolate( h.p(bends[1]).x )) );
+				sort(intersections.begin(), intersections.end());
+				
+				Intersection prev(0, interpolate(h.p0.x));
+				intersections.push_back( Intersection(1, interpolate(h.p1.x)) );
+				
+				for(IntersectionList::const_iterator bi = intersections.begin(); bi != intersections.end(); ++bi) {
+					const Intersection &next = *bi;
+					if (equal(prev.l, next.l)) continue;
+
+					bool flip = next.point.l < prev.point.l;
+					Real dl = next.point.l - prev.point.l;
+					bool e0 = flip ? prev.point.e0 : prev.point.e1;
+					bool e1 = flip ? next.point.e1 : next.point.e0;
+					if (e0 && e1) {
+						Vector bt0 = flip ? prev.point.t0 : prev.point.t1;
+						Vector bt1 = flip ? next.point.t1 : next.point.t0;
+						Hermite bh(prev.point.p, next.point.p, bt0*dl, bt1*dl);
+						
+						Hermite src_h(h.sub(prev.l, next.l));
+						Vector x0 = bt0.norm(), y0 = x0.perp();
+						Vector x1 = bt1.norm(), y1 = x1.perp();
+						Hermite dst_h(
+							prev.point.p + y0*src_h.p0.y,
+							next.point.p + y1*src_h.p1.y,
+							x0*src_h.t0.x + y0*src_h.t0.y,
+							x1*src_h.t1.x + y1*src_h.t1.y );
+						Real k = (src_h.p1 - src_h.p0).lensqr();
+						if (!equal(k, 0)) {
+							k = sqrt( (dst_h.p1 - dst_h.p0).lensqr()/k );
+							dst_h.t0 *= k;
+							dst_h.t1 *= k;
+						}
+						
+						half_corner(dst, prev.point, src_h.p0.y, flip, true);
+						dst.hermite(dst_h);
+						half_corner(dst, next.point, src_h.p1.y, flip, false);
+					}
+					
+					prev = next;
+				}
+				
+				intersections.clear();
+			}
+			
+			if (dst_size < dst.size())
+				dst.line_to(dst[dst_size].p1);
+		}
+		
+		void draw(const Context &context) const {
+			if (empty()) return;
+			
+			context->save();
+			
+			context.set_line_width_px(10);
+			context.set_source_rgba(Color::blue().mulalpha(0.1));
+
+			for(const_iterator i = begin(); i != end(); ++i) {
+				context.move_to(i->p + i->t0);
+				context.line_to(i->p);
+				context.line_to(i->p + i->t1);
+			}
+			context->stroke();
+
+			context.set_source_rgba(Color::blue().mulalpha(0.25));
+			for(const_iterator j = begin(), i = j++; j != end(); i = j++) {
+				Real dl = j->l - i->l;
+				context.hermite(Hermite(i->p, j->p, i->t1*dl, j->t0*dl), true);
+			}
+			context->stroke();
+
+			context->restore();
+		}
+	};
+	
+	void build_contour(Contour &dst, const ActiveCurve &curve) {
+		size_t dst_size = dst.size();
+		for(ActiveCurve::SegIter i(curve); i; i.next())
+			dst.hermite( i.hermite() );
+		if (dst_size < dst.size())
+			dst.line_to( dst[dst_size].p1 );
+	}
+
+	void build_bend(BendContour &dst, const ActiveCurve &curve, bool round, int segments) {
+		const Real dl = Real(1)/segments;
+		size_t dst_size = dst.size();
+		BendPoint bp;
+		BendMode mode = round ? BendRound : BendCorner;
+		
+		bp.p = curve.front().get_p();
+		bp.e0 = bp.e1 = true;
+		bp.l = 0;
+		bp.length = 0;
+		bp.mode = BendNone;
+		dst.push_back( bp );
+		
+		int index = 0;
+		for(ActiveCurve::SegIter i(curve); i; i.next()) {
+			Hermite h = i.hermite();
+			
+			dst.back().p = h.p0;
+			dst.back().t1 = h.t0;
+			
+			Real l0 = index;
+			Real l = dl;
+			bp.mode = BendNone;
+			for(int k = 2; k < segments; ++k, l += dl) {
+				Vector p = h.p(l);
+				bp.l = l0 + l;
+				bp.length += (p - bp.p).len();
+				bp.p = p;
+				bp.t0 = bp.t1 = h.t(l);
+				dst.push_back( bp );
+			}
+			
+			bp.l = ++index;
+			bp.length += (h.p1 - bp.p).len();
+			bp.p = h.p1;
+			bp.t0 = h.t1;
+			bp.mode = mode;
+			dst.push_back( bp );
+		}
+		
+		BendPoint &first = dst[dst_size];
+		BendPoint &last = dst.back();
+		if (curve.loop) {
+			first.mode = mode;
+			first.t0 = last.t0;
+			last.t1 = first.t1;
+			first.e0 = last.e1 = false;
+		} else {
+			first.t0 = first.t1;
+			last.t1 = last.t0;
+			last.mode = BendNone;
+		}
+	}
+}
+
+
+
+
+BendArea::BendArea():
+	p0(new ActivePoint(Vector(100, 100), Color::blue())),
+	p1(new ActivePoint(Vector(900, 100), Color::blue()))
+{
+	set_size_request(1000, 900);
+
+	add_point(p0);
+	add_point(p1);
+	
+	// generate src curve
+	Vector b0 = p0->position - Vector(20, 80);
+	Vector b1 = p1->position + Vector(20, 0);
+	Vector pp(b0.x, b1.y);
+	while(true) {
+		Vector pt0 = pp + Vector(-50, 0); //rnd::vector(pp.x - 20, b0.y, pp.x + 20, b1.y);
+		Vector pt1 = pp + Vector( 50, 0); //rnd::vector(pp.x - 20, b0.y, pp.x + 20, b1.y);
+
+		src_curve.push_back(ActiveCurvePoint(pp, pt0 - pp, pt1 - pp));
+		src_curve.back().add_to(*this);
+
+		if (!less(pp.x, b1.x)) break;
+		pp = Vector(pp.x + 60, 0.5*(b0.y + b1.y)); //rnd::vector(pp.x + 20, b0.y, pp.x + 100, b1.y);
+		if (!less(pp.x, b1.x)) pp = b1;
+	}
+	src_curve.set_tail_tangents_visible();
+
+	// generate dst curve
+	b0.x = p0->position.x;
+	b0.y = p0->position.y + 200;
+	b1.x = p1->position.x;
+	b1.y = 780;
+	Vector bt(100, 100);
+	dst_curve.loop = true;
+	for(int i = 0; i < 5; ++i) {
+		dst_curve.push_back(ActiveCurvePoint(rnd::vector(b0, b1), rnd::vector(-bt, bt), rnd::vector(-bt, bt)));
+		dst_curve.back().add_to(*this);
+	}
+	dst_curve.set_tail_tangents_visible();
+}
+
+void BendArea::on_point_move(const ActivePoint::Handle &point, const Vector&oldposition, const Vector &newposition) {
+	src_curve.on_point_move(point, oldposition, newposition);
+	dst_curve.on_point_move(point, oldposition, newposition);
+	if (point == p0 || point == p1)
+		p0->position.y = p1->position.y = newposition.y;
+}
+
+void BendArea::on_draw_content(const Context &context) {
+	context->save();
+	
+	context.set_line_width_px(1);
+	context.set_source_rgba(Color::blue());
+	context.move_to(p0->position);
+	context.line_to(p1->position);
+	context->stroke();
+	
+	if (!src_curve.empty()) {
+		context.set_source_rgba(Color::gray().mulalpha(0.5));
+		src_curve.put(context);
+		context->close_path();
+		context->fill();
+
+		if (!dst_curve.empty()) {
+			Contour src_contour;
+			BendContour bend;
+			Contour dst_contour;
+			
+			build_contour(src_contour, src_curve);
+			build_bend(bend, dst_curve, false, 16);
+			
+			bend.draw(context);
+			
+			Real dst_length = bend.back().length;
+			Real src_length = p1->position.x - p0->position.x;
+			if (!equal(src_length, 0) && !equal(dst_length, 0)) {
+				Real kx = dst_length/src_length;
+				Real ky = 1;
+				Matrix m( kx, 0, 0,
+						  0, ky, 0,
+						  -kx*p0->position.x, -ky*p0->position.y, 1 );
+				bend.bend(dst_contour, src_contour, m);
+				
+				dst_contour.put(context, true);
+				context->close_path();
+				context->fill_preserve();
+				context.set_source_rgba(Color::blue());
+				context->stroke();
+			}
+		}
+	}
+	
+	context->restore();
+	
+	src_curve.draw(context);
+	dst_curve.draw(context);
+}
+
diff --git a/c++/vector/bendarea.h b/c++/vector/bendarea.h
new file mode 100644
index 0000000..40cc9fc
--- /dev/null
+++ b/c++/vector/bendarea.h
@@ -0,0 +1,24 @@
+#ifndef BENDAREA_H
+#define BENDAREA_H
+
+
+#include "activearea.h"
+#include "activecurve.h"
+
+
+class BendArea: public ActiveArea {
+protected:
+	ActivePoint::Handle p0;
+	ActivePoint::Handle p1;
+	ActiveCurve src_curve;
+	ActiveCurve dst_curve;
+
+	void on_point_move(const ActivePoint::Handle &point, const Vector &oldposition, const Vector &newposition) override;
+	void on_draw_content(const Context &context) override;
+	
+public:
+	BendArea();
+};
+
+
+#endif
diff --git a/c++/vector/context.cpp b/c++/vector/context.cpp
index b2d3242..9d990bb 100644
--- a/c++/vector/context.cpp
+++ b/c++/vector/context.cpp
@@ -38,9 +38,7 @@ void Context::point(const ActivePoint &point) const {
 	ctx()->restore();
 }
 
-void Context::segment(const Segment &segment) const {
-	line_to(segment.p0);
-	Vector b0 = segment.p0 + segment.t0/3;
-	Vector b1 = segment.p1 - segment.t1/3;
-	ctx()->curve_to(b0.x, b0.y, b1.x, b1.y, segment.p1.x, segment.p1.y);
+void Context::bezier(const Bezier &bezier, bool jump) const {
+	if (jump) move_to(bezier.p0); else line_to(bezier.p0);
+	ctx()->curve_to(bezier.pp0.x, bezier.pp0.y, bezier.pp1.x, bezier.pp1.y, bezier.p1.x, bezier.p1.y);
 }
diff --git a/c++/vector/context.h b/c++/vector/context.h
index 1fd9f0b..6b6cff1 100644
--- a/c++/vector/context.h
+++ b/c++/vector/context.h
@@ -6,8 +6,9 @@
 
 #include <cairomm/cairomm.h>
 
-#include "vector.h"
 #include "color.h"
+#include "vector.h"
+#include "curve.h"
 
 class ActivePoint;
 
@@ -39,7 +40,9 @@ public:
 		{ ctx()->line_to(p.x, p.y); }
 	
 	void point(const ActivePoint &point) const;
-	void segment(const Segment &segment) const;
+	void bezier(const Bezier &bezier, bool jump = false) const;
+	void hermite(const Hermite &hermite, bool jump = false) const
+		{ bezier(hermite.get_bezier(), jump); }
 };
 
 
diff --git a/c++/vector/curve.cpp b/c++/vector/curve.cpp
new file mode 100644
index 0000000..2e9e0c3
--- /dev/null
+++ b/c++/vector/curve.cpp
@@ -0,0 +1,90 @@
+
+#include "curve.h"
+
+
+int Hermite::bends(Real *l, const Real &p0, const Real &p1, const Real &t0, const Real &t1) {
+	Real roots[2];
+	int count = solve_equation(
+		roots,
+		   p0 +          t0       ,
+		-6*p0 + 6*p1 - 4*t0 + 2*t1,
+		 6*p0 - 6*p1 + 3*t0 + 3*t1 );
+	int valid_count = 0;	
+	for(Real *i = roots, *end = i + count; i != end; ++i)
+		if (less(0, *i) && less(*i, 1)) {
+			if (l) l[valid_count] = *i;
+			++valid_count;
+		}
+	return valid_count;
+}
+
+Range Hermite::bounds_accurate(const Real &p0, const Real &p1, const Real &t0, const Real &t1) {
+	Range range(p0);
+	range.expand(p1);
+	Real roots[2];
+	int count = bends(roots, p0, p1, t0, t1);
+	for(Real *i = roots, *end = i + count; i != end; ++i)
+		range.expand( p(*i, p0, p1, t0, t1) );
+	return range;
+}
+
+int Hermite::intersections(Real *l, const Real &p, const Real &p0, const Real &p1, const Real &t0, const Real &t1) {
+	Real roots[3];
+	int count = solve_equation(
+		roots,
+		   p0 - p               ,
+		                 t0     ,
+		-3*p0 + 3*p1 - 2*t0 - t1,
+		 2*p0 - 2*p1 +   t0 + t1 );
+	int valid_count = 0;
+	for(Real *i = roots, *end = i + count; i != end; ++i)
+		if (less(0, *i) && less(*i, 1)) {
+			if (l) l[valid_count] = *i;
+			++valid_count;
+		}
+	return valid_count;
+}
+
+
+int Bezier::bends(Real *l, const Real &p0, const Real &p1, const Real &pp0, const Real &pp1) {
+	Real roots[2];
+	int count = solve_equation(
+		roots,
+		 -p0 +   pp0             ,
+		2*p0 - 6*pp0 + 2*pp1     ,
+		 -p0 + 3*pp0 - 3*pp1 + p1 );
+	int valid_count = 0;	
+	for(Real *i = roots, *end = i + count; i != end; ++i)
+		if (less(0, *i) && less(*i, 1)) {
+			if (l) l[valid_count] = *i;
+			++valid_count;
+		}
+	return valid_count;
+}
+
+Range Bezier::bounds_accurate(const Real &p0, const Real &p1, const Real &pp0, const Real &pp1) {
+	Range range(p0);
+	range.expand(p1);
+	Real roots[2];
+	int count = bends(roots, p0, p1, pp0, pp1);
+	for(Real *i = roots, *end = i + count; i != end; ++i)
+		range.expand( p(*i, p0, p1, pp0, pp1) );
+	return range;
+}
+
+int Bezier::intersections(Real *l, const Real &p, const Real &p0, const Real &p1, const Real &pp0, const Real &pp1) {
+	Real roots[3];
+	int count = solve_equation(
+		roots,
+		   p0 - p                 ,
+		-3*p0 + 3*pp0             ,
+		 3*p0 - 6*pp0 + 3*pp1     ,
+		  -p0 + 3*pp0 - 3*pp1 + p1 );
+	int valid_count = 0;
+	for(Real *i = roots, *end = i + count; i != end; ++i)
+		if (less(0, *i) && less(*i, 1)) {
+			if (l) l[valid_count] = *i;
+			++valid_count;
+		}
+	return valid_count;
+}
diff --git a/c++/vector/curve.h b/c++/vector/curve.h
new file mode 100644
index 0000000..9b55cd6
--- /dev/null
+++ b/c++/vector/curve.h
@@ -0,0 +1,272 @@
+#ifndef CURVE_H
+#define CURVE_H
+
+
+#include <cassert>
+
+#include "vector.h"
+
+
+class Bezier;
+
+class Hermite {
+public:
+	Vector p0, p1, t0, t1;
+
+	explicit Hermite(
+		const Vector &p0 = Vector(),
+		const Vector &p1 = Vector(),
+		const Vector &t0 = Vector(),
+		const Vector &t1 = Vector()
+	):
+		p0(p0), p1(p1), t0(t0), t1(t1) { }
+
+	bool operator< (const Hermite &a) const {
+		return p0 < a.p0 ? true : a.p0 < p0 ? false
+			 : p1 < a.p1 ? true : a.p1 < p1 ? false
+			 : t0 < a.t0 ? true : a.t0 < t0 ? false
+			 : t1 < a.t1;
+	}
+	bool operator>(const Hermite &a) const { return a < *this; }
+	bool operator==(const Hermite &a) const { return p0 == a.p0 && p1 == a.p1 && t0 == a.t0 && t1 == a.t1; }
+	bool operator!=(const Hermite &a) const { return !(*this == a); }
+
+	Vector p(const Real &l) const
+		{ return p(l, p0, p1, t0, t1); }
+	Vector t(const Real &l) const
+		{ return t(l, p0, p1, t0, t1); }
+
+	Vector operator() (const Real &l) const  { return p(l); }
+	
+	Vector d0(const Real &dl = differential)
+		{ return (p(dl) - p0).norm(); }
+	Vector d1(const Real &dl = differential)
+		{ return (p1 - p(1 - dl)).norm(); }
+	
+	Vector pp0() const { return p0 + t0/3; }
+	Vector pp1() const { return p1 - t1/3; }
+	
+	Hermite flip() const
+		{ return Hermite(p1, p0, -t1, -t0); }
+	
+	Hermite sub(const Real &a, const Real &b) const {
+		Real k = b - a;
+		return Hermite(p(a), p(b), t(a)*k, t(b)*k);
+	}
+	
+	void split(const Real &l, Hermite *a, Hermite *b) const {
+		if (!a && !b) return;
+		assert(a != b);
+		Vector point = p(l);
+		Vector tangent = t(l);
+		Hermite h = *this;
+		if (a) {
+			a->p0 = h.p0;
+			a->t0 = h.t0*l;
+			a->p1 = point;
+			a->t1 = tangent*l;
+		}
+		if (b && b != a) {
+			Real ll = 1 - l;
+			b->p0 = point;
+			b->t0 = tangent*ll;
+			b->p1 = h.p1;
+			b->t1 = h.t1*ll;
+		}
+	}
+
+	Rect bounds() const
+		{ return Rect(p0).expand(p1).expand(pp0()).expand(pp1()); }
+
+	int bends_x(Real *l) const
+		{ return bends(l, p0.x, p1.x, t0.x, t1.x); }
+	int bends_y(Real *l) const
+		{ return bends(l, p0.y, p1.y, t0.y, t1.y); }
+
+	Range bounds_accurate_x() const
+		{ return bounds_accurate(p0.x, p1.x, t0.x, t1.x); }
+	Range bounds_accurate_y() const
+		{ return bounds_accurate(p0.y, p1.y, t0.y, t1.y); }
+	Rect bounds_accurate() const
+		{ return Rect(bounds_accurate_x(), bounds_accurate_y()); }
+
+	int intersections_x(Real *l, const Real &x) const
+		{ return intersections(l, x, p0.x, p1.x, t0.x, t1.x); }
+	int intersections_y(Real *l, const Real &y) const
+		{ return intersections(l, y, p0.y, p1.y, t0.y, t1.y); }
+	
+	inline Bezier get_bezier() const;
+	
+	static Hermite linear(const Vector &p, const Vector &t)
+		{ return Hermite(p, p + t, t, t);  }
+	static Hermite point(const Vector &p)
+		{ return Hermite(p, p);  }
+
+	static Real p(const Real &l, const Real &p0, const Real &p1, const Real &t0, const Real &t1) {
+		return p0*((( 2*l - 3)*l    )*l + 1)
+			 + p1*(((-2*l + 3)*l    )*l    )
+			 + t0*(((   l - 2)*l + 1)*l    )
+			 + t1*(((   l - 1)*l    )*l    );
+	}
+
+	static Real t(const Real &l, const Real &p0, const Real &p1, const Real &t0, const Real &t1) {
+		return p0*(( 6*l - 6)*l    )
+			 + p1*((-6*l + 6)*l    )
+			 + t0*(( 3*l - 4)*l + 1)
+			 + t1*(( 3*l - 2)*l    );
+	}
+
+	static Vector p(const Real &l, const Vector &p0, const Vector &p1, const Vector &t0, const Vector &t1)
+		{ return Vector( p(l, p0.x, p1.x, t0.x, t1.x), p(l, p0.y, p1.y, t0.y, t1.y) ); }
+	static Vector t(const Real &l, const Vector &p0, const Vector &p1, const Vector &t0, const Vector &t1)
+		{ return Vector( t(l, p0.x, p1.x, t0.x, t1.x), t(l, p0.y, p1.y, t0.y, t1.y) ); }
+	
+	static int bends(Real *l, const Real &p0, const Real &p1, const Real &t0, const Real &t1);
+	static int intersections(Real *l, const Real &p, const Real &p0, const Real &p1, const Real &t0, const Real &t1);
+	static Range bounds_accurate(const Real &p0, const Real &p1, const Real &t0, const Real &t1);
+};
+
+
+class Bezier {
+public:
+	Vector p0, pp0, pp1, p1;
+
+	Bezier(const Vector &p0, const Vector &p1, const Vector &pp0, const Vector &pp1):
+		p0(p0), pp0(pp0), pp1(pp1), p1(p1) { }
+	Bezier(const Vector &p0, const Vector &p1, const Vector &pp):
+		Bezier(p0, p1, pp, pp) { }
+	Bezier(const Vector &p0, const Vector &p1):
+		Bezier(p0, p1, (p1 - p0)/3 + p0, (p0 - p1)/3 + p1) { }
+	explicit Bezier(const Vector &p = Vector()):
+		Bezier(p, p, p, p) { }
+
+	bool operator< (const Bezier &a) const {
+		return p0  < a.p0  ? true : a.p0  < p0  ? false
+			 : p1  < a.p1  ? true : a.p1  < p1  ? false
+			 : pp0 < a.pp0 ? true : a.pp0 < pp0 ? false
+			 : pp1 < a.pp1;
+	}
+	bool operator>(const Bezier &a) const { return a < *this; }
+	bool operator==(const Bezier &a) const { return p0 == a.p0 && p1 == a.p1 && pp0 == a.pp0 && pp1 == a.pp1; }
+	bool operator!=(const Bezier &a) const { return !(*this == a); }
+
+	Vector p(const Real &l) const
+		{ return p(l, p0, p1, pp0, pp1); }
+	Vector t(const Real &l) const
+		{ return t(l, p0, p1, pp0, pp1); }
+
+	Vector operator() (const Real &l) const  { return p(l); }
+	
+	Vector d0(const Real &dl = differential)
+		{ return (p(dl) - p0).norm(); }
+	Vector d1(const Real &dl = differential)
+		{ return (p1 - p(1 - dl)).norm(); }
+	
+	Vector t0() const { return 3*(pp0 - p0); }
+	Vector t1() const { return 3*(p1 - pp1); }
+	
+	void split(const Real &l, Bezier *a, Bezier *b) const {
+		if (!a && !b) return;
+		assert(a != b);
+		Real ll = 1 - l;
+		Vector p0 = this->p0;
+		Vector p1 = this->p1;
+		Vector a0 =  p0*ll + pp0*l;
+		Vector a1 = pp0*ll + pp1*l;
+		Vector a2 = pp1*ll +  p1*l;
+		Vector b0 =  a0*ll +  a1*l;
+		Vector b1 =  a1*ll +  a2*l;
+		Vector c  =  b0*ll +  b1*l;
+		if (a) {
+			a->p0  = b->p0;
+			a->pp0 = a0;
+			a->pp1 = b0;
+			a->p1  = c;
+		}
+		if (b && b != a) {
+			b->p0  = c;
+			b->pp0 = b1;
+			b->pp1 = a2;
+			b->p1  = b->p1;
+		}
+	}
+	
+	Bezier flip() const
+		{ return Bezier(p1, p0, pp1, pp0); }
+
+	Bezier sub(const Real &a, const Real &b) const {
+		if (equal(a, b)) return Bezier( p(0.5*(a + b)) );
+
+		Bezier s;
+		if (equal(a, 1)) {
+			split(b, nullptr, &s);
+			return s.flip();
+		}
+		
+		split(a, nullptr, &s);
+		s.split((b - a)/(1 - a), &s, nullptr);
+		return s;
+	}
+	
+	Rect bounds() const
+		{ return Rect(p0).expand(p1).expand(pp0).expand(pp1); }
+
+	int bends_x(Real *l) const
+		{ return bends(l, p0.x, p1.x, pp0.x, pp1.x); }
+	int bends_y(Real *l) const
+		{ return bends(l, p0.y, p1.y, pp0.y, pp1.y); }
+
+	Range bounds_accurate_x() const
+		{ return bounds_accurate(p0.x, p1.x, pp0.x, pp1.x); }
+	Range bounds_accurate_y() const
+		{ return bounds_accurate(p0.y, p1.y, pp0.y, pp1.y); }
+	Rect bounds_accurate() const
+		{ return Rect(bounds_accurate_x(), bounds_accurate_y()); }
+
+	int intersections_x(Real *l, const Real &x) const
+		{ return intersections(l, x, p0.x, p1.x, pp0.x, pp1.x); }
+	int intersections_y(Real *l, const Real &y) const
+		{ return intersections(l, y, p0.y, p1.y, pp0.y, pp1.y); }
+
+	inline Hermite get_hermite() const;
+
+	static Bezier linear(const Vector &p0, const Vector &p1)
+		{ return Bezier(p0, p1, p0 + (p1 - p0)/3, p1 + (p0 - p1)/3);  }
+	static Bezier point(const Vector &p)
+		{ return Bezier(p, p);  }
+
+	static Real p(const Real &l, const Real &p0, const Real &p1, const Real &pp0, const Real &pp1) {
+		Real ll = 1 - l;
+		const Real a0 =  p0*ll + pp0*l;
+		const Real a1 = pp0*ll + pp1*l;
+		const Real a2 = pp1*ll +  p1*l;
+		const Real b0 =  a0*ll +  a1*l;
+		const Real b1 =  a1*ll +  a2*l;
+		return b0*ll +  b1*l;
+	}
+
+	static Real t(const Real &l, const Real &p0, const Real &p1, const Real &pp0, const Real &pp1) {
+		Real ll = 1 - l;
+		const Real a0 =  p0*ll + pp0*l;
+		const Real a1 = pp0*ll + pp1*l;
+		const Real a2 = pp1*ll +  p1*l;
+		const Real b0 =  a0*ll +  a1*l;
+		const Real b1 =  a1*ll +  a2*l;
+		return 3*(b1 - b0)*l;
+	}
+
+	static Vector p(const Real &l, const Vector &p0, const Vector &p1, const Vector &pp0, const Vector &pp1)
+		{ return Vector( p(l, p0.x, p1.x, pp0.x, pp1.x), p(l, p0.y, p1.y, pp0.y, pp1.y) ); }
+	static Vector t(const Real &l, const Vector &p0, const Vector &p1, const Vector &pp0, const Vector &pp1)
+		{ return Vector( t(l, p0.x, p1.x, pp0.x, pp1.x), t(l, p0.y, p1.y, pp0.y, pp1.y) ); }
+	
+	static int bends(Real *l, const Real &p0, const Real &p1, const Real &pp0, const Real &pp1);
+	static int intersections(Real *l, const Real &p, const Real &p0, const Real &p1, const Real &pp0, const Real &pp1);
+	static Range bounds_accurate(const Real &p0, const Real &p1, const Real &pp0, const Real &pp1);
+};
+
+
+inline Bezier Hermite::get_bezier() const { return Bezier(p0, p1, pp0(), pp1()); }
+inline Hermite Bezier::get_hermite() const { return Hermite(p0, p1, t0(), t1()); }
+
+#endif
diff --git a/c++/vector/curvesarea.cpp b/c++/vector/curvesarea.cpp
deleted file mode 100644
index 10d8b56..0000000
--- a/c++/vector/curvesarea.cpp
+++ /dev/null
@@ -1,128 +0,0 @@
-
-
-#include "curvesarea.h"
-
-
-namespace {
-	struct DrawSegment {
-	private:
-		const Context &context;
-		const Segment &segment;
-		const Real width;
-		const Real precision_sqr;
-		const int max_level = 16;
-		int level;
-		Vector p0;
-		Real l0;
-
-		Vector func(const Real &l)
-			{ return segment(l) + segment.t(l).perp().norm()*width; }
-
-		void step(const Real &l1) {
-			Vector p1 = segment(l1) + segment.t(l1).perp().norm()*width;
-			if (++level < max_level)
-				while (greater((p1 - p0).lensqr(), precision_sqr))
-					step(l0 + 0.5*(l1 - l0));
-			--level;
-			context.line_to(p1);
-			p0 = p1; l0 = l1;
-		}
-
-	public:
-		DrawSegment(const Context &context, const Segment &segment, const Real &width):
-			context(context), segment(segment), width(width), precision_sqr(context.pixelsize_sqr()), level(), l0()
-		{
-			context.line_to(p0 = func(0));
-			step(1);
-		}
-	};
-}
-
-
-CurvesArea::CurvesArea():
-	w(new ActivePoint(Vector(), Color::red()))
-{
-	set_size_request(500, 500);
-	
-	Vector center(250, 250);
-	const Real radius = 100;
-	const Real tangent_len = 50;
-	const int count = 5;
-	for(int i = 0; i < count; ++i) {
-		Real a = 2*pi*i/count;
-		Vector v = Vector::from_angle(a);
-		curve.push_back(CurvePoint(center + v*radius, v*tangent_len));
-		curve.back().add_to(*this);
-	}
-	curve.loop = true;
-	add_point(w);
-	set_width(50);
-}
-
-Real CurvesArea::get_width() const
-	{ return fabs(dot(curve.front().get_t().perp().norm(), w->position - curve.front().get_p())); }
-
-void CurvesArea::set_width(const Real &width)
-	{ w->position = curve.front().get_p() + curve.front().get_t().perp().norm()*width; }
-
-void CurvesArea::on_point_move(const ActivePoint::Handle &point, const Vector &oldposition, const Vector &newposition) {
-	Real width = -1;
-	for(Curve::const_iterator i = curve.begin(); i != curve.end(); ++i) {
-		if (i == curve.begin() && (i->p == point || i->t == point)) {
-			point->position = oldposition;
-			width = get_width();
-			point->position = newposition;
-		}
-		if (i->p == point)
-			{ i->t->position += i->p->position - oldposition; break; }
-	}
-	if (point == w)
-		width = get_width();
-	if (width >= 0)
-		set_width(width);
-}
-
-void CurvesArea::put_segment(const Context &context, const Segment &segment, const Real width) {
-	if (equal(width, 0)) context.segment(segment);
-		            else DrawSegment(context, segment, width);
-}
-
-void CurvesArea::put_curve(const Context &context, const Curve &curve, const Real width) {
-	for(Curve::const_iterator curr = curve.begin(); curr != curve.end(); ++curr) {
-		Curve::const_iterator next = curve.next(curr);
-		put_segment(context, Segment(curr->get_p(), next->get_p(), curr->get_t(), next->get_t()), width);
-	}
-}
-
-void CurvesArea::draw_curve(const Context &context, const Curve &curve, const Real width) {
-	context->save();
-	
-	context.set_line_width_px(1);
-
-	// tangents
-	context.set_source_rgba(Color::gray());
-	for(Curve::const_iterator i = curve.begin(); i != curve.end(); ++i) {
-		context.move_to(i->p->position);
-		context.line_to(i->t->position);
-		context->stroke();
-	}
-
-	// middle line
-	context.set_source_rgba(Color::black());
-	put_curve(context, curve);
-	context->stroke();
-	
-	// borders
-	context.set_source_rgba(Color::blue());
-	put_curve(context, curve, width);
-	context->stroke();
-	put_curve(context, curve, -width);
-	context->stroke();
-	
-	context->restore();
-}
-
-
-void CurvesArea::on_draw_content(const Context &context) {
-	draw_curve(context, curve, get_width());
-}
diff --git a/c++/vector/curvesarea.h b/c++/vector/curvesarea.h
deleted file mode 100644
index c233511..0000000
--- a/c++/vector/curvesarea.h
+++ /dev/null
@@ -1,57 +0,0 @@
-#ifndef CURVEAREA_H
-#define CURVEAREA_H
-
-
-#include <vector>
-
-#include "activearea.h"
-
-
-class CurvesArea: public ActiveArea {
-public:
-	class CurvePoint {
-	public:
-		ActivePoint::Handle p, t;
-		
-		explicit CurvePoint(const Vector &p = Vector(), const Vector &t = Vector()):
-			p(new ActivePoint(p)),
-			t(new ActivePoint(p + t, Color::yellow()))
-			{ }
-		
-		const Vector& get_p() const { return p->position; }
-		Vector get_t() const { return t->position - p->position; }
-		
-		void add_to(CurvesArea &area) const { area.add_point(p); area.add_point(t); }
-		void remove_from(CurvesArea &area) const { area.remove_point(p); area.remove_point(t); }
-	};
-
-	class Curve: public std::vector<CurvePoint> {
-	public:
-		bool loop;
-		Curve(): loop() { }
-		Curve::const_iterator next(Curve::const_iterator i) const
-			{ return ++i != end() ? i : loop ? begin() : --i; }
-		Curve::const_iterator prev(Curve::const_iterator i) const
-			{ return --i != begin() ? i : loop ? end() : ++i; }
-	};
-
-protected:
-	Curve curve;
-	ActivePoint::Handle w;
-
-	Real get_width() const;
-	void set_width(const Real &width);
-
-	void put_segment(const Context &context, const Segment &segment, const Real width = 0);
-	void put_curve(const Context &context, const Curve &curve, const Real width = 0);
-	void draw_curve(const Context &context, const Curve &curve, const Real width);
-
-	void on_point_move(const ActivePoint::Handle &point, const Vector &oldposition, const Vector &newposition) override;
-	void on_draw_content(const Context &context) override;
-	
-public:
-	CurvesArea();
-};
-
-
-#endif
diff --git a/c++/vector/mainwindow.cpp b/c++/vector/mainwindow.cpp
index 8fd8055..2d9d1d2 100644
--- a/c++/vector/mainwindow.cpp
+++ b/c++/vector/mainwindow.cpp
@@ -1,12 +1,16 @@
 
 
-#include "curvesarea.h"
+#include "outlinearea.h"
+#include "bendarea.h"
 
 #include "mainwindow.h"
 
 
 MainWindow::MainWindow() {
-	CurvesArea *curves = manage(new CurvesArea());
-	curves->show();
-	add(*curves);
+	//OutlineArea *content = manage(new OutlineArea());
+	BendArea *content = manage(new BendArea());
+
+	content->show();
+	add(*content);
+
 }
diff --git a/c++/vector/matrix.cpp b/c++/vector/matrix.cpp
new file mode 100644
index 0000000..53bbec4
--- /dev/null
+++ b/c++/vector/matrix.cpp
@@ -0,0 +1,51 @@
+
+#include "matrix.h"
+
+
+Matrix
+Matrix::get_inverted() const {
+	Real d = det();
+	if (equal(d, 0)) {
+		// try to make back transform for X-axis
+		if (equal(m02, 0) && equal(m10, 0) && equal(m11, 0) && equal(m12, 0) && equal(m22, 1)) {
+			d = axis_x().lensqr();
+			if (d > precision) {
+				d = 1.0/d;
+				return Matrix(
+				    d*m00,                 0, 0,
+				    d*m01,                 0, 0,
+				   -d*(m20*m00 + m21*m01), 0, 1 );
+			}
+		}
+
+		// try to make back transform for Y-axis
+		if (equal(m00, 0) && equal(m01, 0) && equal(m02, 0) && equal(m12, 0) && equal(m22, 1)) {
+			d = axis_y().lensqr();
+			if (d > precision) {
+				d = 1/d;
+				return Matrix(
+				    0,  d*m10,                 0,
+					0,  d*m11,                 0,
+					0, -d*(m20*m10 + m21*m11), 1 );
+			}
+		}
+
+		// give up
+		return Matrix( 0, 0, 0,
+					   0, 0, 0,
+					   0, 0, 0 );
+	}
+
+	// proper inversion
+	Real p = 1/d, m = -p;
+	return Matrix(
+		p*(m11*m22 - m12*m21), // row0
+		m*(m01*m22 - m02*m21),
+		p*(m01*m12 - m02*m11),
+		m*(m10*m22 - m12*m20), // row1
+		p*(m00*m22 - m02*m20),
+		m*(m00*m12 - m02*m10),
+		p*(m10*m21 - m11*m20), // row2
+		m*(m00*m21 - m01*m20),
+		p*(m00*m11 - m01*m10) );
+}
diff --git a/c++/vector/matrix.h b/c++/vector/matrix.h
new file mode 100644
index 0000000..5512ab0
--- /dev/null
+++ b/c++/vector/matrix.h
@@ -0,0 +1,132 @@
+#ifndef MATRIX_H
+#define MATRIX_H
+
+
+#include <cmath>
+
+#include "vector.h"
+
+
+class Matrix {
+public:
+	union {
+		Real m[3][3];
+		struct {
+			Real m00, m01, m02;
+			Real m10, m11, m12;
+			Real m20, m21, m22;
+		};
+	};
+
+	Matrix(
+		Real m00, Real m01, Real m02,
+		Real m10, Real m11, Real m12,
+		Real m20, Real m21, Real m22
+	):
+		m00(m00), m01(m01), m02(m02),
+		m10(m10), m11(m11), m12(m12),
+		m20(m20), m21(m21), m22(m22)
+	{ }
+
+	Matrix(): Matrix(
+		1, 0, 0,
+		0, 1, 0,
+		0, 0, 1 ) { }
+
+	Matrix(
+		const Vector &axis_x,
+		const Vector &axis_y,
+		const Vector &offset = Vector()
+	): Matrix(
+		axis_x.x, axis_x.y, 0,
+		axis_y.x, axis_y.y, 0,
+		offset.x, offset.y, 1 ) { }
+
+	const Vector& axis(int index) const { return *(const Vector*)(m[index]); }
+	const Vector& axis_x() const { return axis(0); }
+	const Vector& axis_y() const { return axis(1); }
+	const Vector& offset() const { return axis(2); }
+
+	Vector& axis(int index) { return *(Vector*)(m[index]); }
+	Vector& axis_x() { return axis(0); }
+	Vector& axis_y() { return axis(1); }
+	Vector& offset() { return axis(2); }
+
+	Matrix& set_identity()
+		{ return *this = Matrix(); }
+	bool is_identity() const
+		{ return *this == Matrix(); }
+
+	Matrix& set_scale(const Real &x, const Real &y) {
+		return *this = Matrix( x, 0, 0,
+							   0, y, 0,
+							   0, 0, 1 );
+	}
+	Matrix& set_scale(const Real &s)
+		{ return set_scale(s, s); }
+	Matrix& set_scale(const Vector &s)
+		{ return set_scale(s.x, s.y); }
+
+	Matrix& set_rotate(const Real &angle) {
+		Real s = sin(angle), c = cos(angle);
+		return *this = Matrix( c, s, 0,
+							  -s, c, 0,
+							   0, 0, 1 );
+	}
+
+	Matrix& set_translate(const Real &x, const Real &y) {
+		return *this = Matrix( 1, 0, 0,
+							   0, 1, 0,
+							   x, y, 1 );
+	}
+	Matrix& set_translate(const Vector &t)
+		{ return set_translate(t.x, t.y); }
+
+	void transform(
+		Real &out_x, Real &out_y, Real &out_z,
+		const Real &x, const Real &y, const Real &z ) const
+	{
+		out_x = x*m00 + y*m10 + z*m20;
+		out_y = x*m01 + y*m11 + z*m21;
+		out_z = x*m02 + y*m12 + z*m22;
+	}
+	void transform(Real &out_x, Real &out_y, const Real &x, const Real &y, bool translate = true) const
+		{ Real z; transform(out_x, out_y, z, x, y, translate ? 1.0 : 0.0); }
+	Vector transform(const Vector &v, bool translate = true)const
+		{ Vector vv; transform(vv.x, vv.y, v.x, v.y, translate); return vv; }
+
+	bool operator==(const Matrix &other) const {
+		return equal(m00, other.m00) && equal(m01, other.m01) && equal(m02, other.m02)
+			&& equal(m10, other.m10) && equal(m11, other.m11) && equal(m12, other.m12)
+			&& equal(m20, other.m20) && equal(m21, other.m21) && equal(m22, other.m22);
+	}
+	bool operator!=(const Matrix &other) const
+		{ return !(*this == other); }
+
+	Matrix& operator*=(const Matrix &other)
+		{ return *this = *this * other; }
+
+	Matrix operator*(const Matrix &other) const {
+		Matrix m;
+		transform(m.m00, m.m01, m.m02, other.m00, other.m01, other.m02);
+		transform(m.m10, m.m11, m.m12, other.m10, other.m11, other.m12);
+		transform(m.m20, m.m21, m.m22, other.m20, other.m21, other.m22);
+		return m;
+	}
+
+	Real det() const {
+		return m00*(m11*m22 - m12*m21)
+			 - m01*(m10*m22 - m12*m20)
+			 + m02*(m10*m21 - m11*m20);
+	}
+	
+	bool is_invertible() const
+		{ return !equal(det(), 0); }
+	
+	Matrix get_inverted() const;
+
+	Matrix& invert()
+		{ return *this = get_inverted(); }
+};
+
+#endif
diff --git a/c++/vector/outlinearea.cpp b/c++/vector/outlinearea.cpp
new file mode 100644
index 0000000..47bfbdd
--- /dev/null
+++ b/c++/vector/outlinearea.cpp
@@ -0,0 +1,115 @@
+
+
+#include "outlinearea.h"
+
+
+namespace {
+	struct DrawBezier {
+	private:
+		const Context &context;
+		const Bezier &bezier;
+		const Real width;
+		const Real precision_sqr;
+		const int max_level = 16;
+		int level;
+		Vector p0;
+		Real l0;
+
+		Vector func(const Real &l)
+			{ return bezier(l) + bezier.t(l).perp().norm()*width; }
+
+		void step(const Real &l1) {
+			Vector p1 = bezier(l1) + bezier.t(l1).perp().norm()*width;
+			if (++level < max_level)
+				while (greater((p1 - p0).lensqr(), precision_sqr))
+					step(l0 + 0.5*(l1 - l0));
+			--level;
+			context.line_to(p1);
+			p0 = p1; l0 = l1;
+		}
+
+	public:
+		DrawBezier(const Context &context, const Bezier &bezier, const Real &width = 0, bool jump = false):
+			context(context), bezier(bezier), width(width), precision_sqr(context.pixelsize_sqr()), level(), l0()
+		{
+			p0 = func(0);
+			if (jump) context.move_to(p0); else context.line_to(p0);
+			step(1);
+		}
+	};
+}
+
+
+OutlineArea::OutlineArea():
+	w(new ActivePoint(Vector(), Color::red()))
+{
+	set_size_request(500, 500);
+	
+	Vector center(250, 250);
+	const Real radius = 100;
+	const Real tangent_len = 50;
+	const int count = 5;
+	for(int i = 0; i < count; ++i) {
+		Real a = 2*pi*i/count;
+		Vector v = Vector::from_angle(a);
+		curve.push_back(ActiveCurvePoint(center + v*radius, v*tangent_len));
+		curve.back().add_to(*this);
+	}
+	//curve.loop = true;
+	curve.set_tail_tangents_visible();
+	add_point(w);
+	set_width(50);
+}
+
+Real OutlineArea::get_width() const
+	{ return fabs(dot(curve.front().get_t1().perp().norm(), w->position - curve.front().get_p())); }
+
+void OutlineArea::set_width(const Real &width)
+	{ w->position = curve.front().get_p() + curve.front().get_t1().perp().norm()*width; }
+
+void OutlineArea::on_point_move(const ActivePoint::Handle &point, const Vector &oldposition, const Vector &newposition) {
+	Real width = -1;
+	if (point == w)
+		width = get_width();
+	if (!curve.empty())
+		if (curve.front().p == point || curve.front().t1 == point)
+			width = get_width();
+
+	curve.on_point_move(point, oldposition, newposition);
+
+	if (width >= 0)
+		set_width(width);
+}
+
+void OutlineArea::put_bezier(const Context &context, const Bezier &bezier, const Real &width, bool jump) {
+	if (equal(width, 0))
+		context.bezier(bezier, jump);
+	else
+		DrawBezier(context, bezier, width, jump);
+}
+
+void OutlineArea::put_curve(const Context &context, const ActiveCurve &curve, const Real &width) {
+	for(ActiveCurve::SegIter i(curve); i; i.next())
+		put_bezier(context, i.bezier(), width);
+	if (curve.loop) context->close_path();
+	for(ActiveCurve::SegRIter i(curve); i; i.next())
+		put_bezier(context, i.bezier(), width);
+	context->close_path();
+}
+
+void OutlineArea::draw_curve(const Context &context, const ActiveCurve &curve, const Real &width) {
+	context->save();
+	
+	curve.draw(context);
+	
+	context.set_line_width_px(1);
+	context.set_source_rgba(Color::blue());
+	put_curve(context, curve, width);
+	context->stroke();
+	
+	context->restore();
+}
+
+void OutlineArea::on_draw_content(const Context &context) {
+	draw_curve(context, curve, get_width());
+}
diff --git a/c++/vector/outlinearea.h b/c++/vector/outlinearea.h
new file mode 100644
index 0000000..69e60fb
--- /dev/null
+++ b/c++/vector/outlinearea.h
@@ -0,0 +1,31 @@
+#ifndef OUTLINEAREA_H
+#define OUTLINEAREA_H
+
+
+#include "activearea.h"
+#include "activecurve.h"
+
+
+class OutlineArea: public ActiveArea {
+protected:
+	ActiveCurve curve;
+	ActivePoint::Handle w;
+
+	Real get_width() const;
+	void set_width(const Real &width);
+
+	void put_bezier(const Context &context, const Bezier &bezier, const Real &width = 0, bool jump = false);
+	void put_hermite(const Context &context, const Hermite &hermite, const Real &width = 0, bool jump = false)
+		{ put_bezier(context, hermite.get_bezier(), width, jump); }
+	void put_curve(const Context &context, const ActiveCurve &curve, const Real &width);
+	void draw_curve(const Context &context, const ActiveCurve &curve, const Real &width);
+
+	void on_point_move(const ActivePoint::Handle &point, const Vector &oldposition, const Vector &newposition) override;
+	void on_draw_content(const Context &context) override;
+	
+public:
+	OutlineArea();
+};
+
+
+#endif
diff --git a/c++/vector/random.h b/c++/vector/random.h
new file mode 100644
index 0000000..749d501
--- /dev/null
+++ b/c++/vector/random.h
@@ -0,0 +1,27 @@
+
+
+#include <cstdlib>
+
+#include "real.h"
+#include "vector.h"
+
+
+namespace rnd {
+	Real real()
+		{ return Real(rand())/Real(RAND_MAX); }
+	Real real(const Real &max)
+		{ return real()*max; }
+	Real real(const Real &min, const Real &max)
+		{ return min + real(max - min); }
+
+	Vector vector()
+		{ return Vector(real(), real()); }
+	Vector vector(const Real &maxx, const Real &maxy)
+		{ return Vector(real(maxx), real(maxy)); }
+	Vector vector(const Real &minx, const Real &miny, const Real &maxx, const Real &maxy)
+		{ return Vector(real(minx, maxx), real(miny, maxy)); }
+	Vector vector(const Vector &max)
+		{ return vector(max.x, max.y); }
+	Vector vector(const Vector &min, const Vector &max)
+		{ return vector(min.x, min.y, max.x, max.y); }
+}
diff --git a/c++/vector/real.cpp b/c++/vector/real.cpp
new file mode 100644
index 0000000..66f3805
--- /dev/null
+++ b/c++/vector/real.cpp
@@ -0,0 +1,83 @@
+
+#include <cmath>
+
+#include "real.h"
+
+
+int solve_equation(Real* roots, const Real& k0, const Real& k1) {
+	if (equal(k1, 0)) return 0;
+	if (roots) roots[0] = -k0/k1;
+	return 1;
+}
+
+int solve_equation(Real* roots, const Real& k0, const Real& k1, const Real& k2) {
+	if (equal(k2, 0)) return solve_equation(roots,k0, k1);
+	Real D = k1*k1 - 4*k2*k0;
+	if (equal(D, 0)) {
+		if (roots) roots[0] = -0.5*k1/k2;
+		return 1;
+	} else
+	if (D > 0) {
+		if (roots) {
+			Real a = sqrt(D);
+			Real b = -0.5/k2;
+			roots[0] = (k1 - a)*b;
+			roots[1] = (k1 + a)*b;
+		}
+		return 2;
+	}
+	return 0;
+}
+
+int solve_equation(Real* roots, const Real& k0, const Real& k1, const Real& k2, const Real& k3) {
+	if (equal(k3, 0)) return solve_equation(roots, k0, k1, k2);
+
+	Real k = 1/k3;
+	Real a = k2*k;
+	Real b = k1*k;
+	Real c = k0*k;
+	
+	Real Q = (a*a - 3*b)/9;
+	Real Q3 = Q*Q*Q;
+	Real R = (2*a*a*a - 9*a*b + 27*c)/54;
+	Real S = Q3 - R*R;
+	
+	if (equal(S, 0)) {
+		if (roots) {
+			Real rr = cbrt(R);
+			Real aa = -a/3;
+			roots[0] = aa - 2*rr;
+			roots[1] = aa + rr;
+		}
+		return 2;
+	} else
+	if (S > 0) {
+		if (roots) {
+			Real ph = acos(R/sqrt(Q3))/3;
+			Real qq = -2*sqrt(Q);
+			Real aa = -a/3;
+			roots[0] = qq*cos(ph) + aa;
+			roots[1] = qq*cos(ph + 2*pi/3) + aa;
+			roots[2] = qq*cos(ph - 2*pi/3) + aa;
+		}
+		return 3;
+	} else
+	if (equal(Q, 0)) {
+		if (roots) roots[0] = -cbrt(c - a*a*a/27) - a/3;
+		return 1;
+	} else
+	if (Q > 0) {
+		if (roots) {
+			Real ph = acosh( fabs(R)/sqrt(Q3) )/3;
+			roots[0] = -2*sign(R)*sqrt(Q)*cosh(ph) - a/3;
+		}
+		return 1;
+	} else {
+		if (roots) {
+			Real ph = asinh( fabs(R)/sqrt(-Q3) )/3;
+			roots[0] = -2*sign(R)*sqrt(-Q)*cosh(ph) - a/3;
+		}
+		return 1;
+	}
+	return 0;
+}
diff --git a/c++/vector/real.h b/c++/vector/real.h
index 6919d30..e7a6101 100644
--- a/c++/vector/real.h
+++ b/c++/vector/real.h
@@ -2,10 +2,14 @@
 #define REAL_H
 
 
+#include <utility>
+
+
 typedef double Real;
 
 
-const Real precision = 1e-5;
+const Real precision = 1e-10;
+const Real differential = 1e-5;
 const Real pi = 3.1415926535897932384626433;
 
 
@@ -16,9 +20,36 @@ inline bool less(const Real &a, const Real &b)
 inline bool greater(const Real &a, const Real &b)
 	{ return less(b, a); }
 inline bool lesseq(const Real &a, const Real &b)
-	{ return equal(a, b) || a < b; }
+	{ return !less(b, a); }
 inline bool greatereq(const Real &a, const Real &b)
 	{ return lesseq(b, a); }
 
 
+template<typename T>
+bool sort(T &a, T &b) {
+	if (!(a < b)) { std::swap(a, b); return true; }
+	return false;
+}
+
+template<typename T>
+bool sort(T &a, T &b, T &c) {
+	bool x = false;
+	if (sort(a, b)) x = true;
+	if (sort(b, c)) x = true;
+	if (x) sort(a, b);
+	return x;
+}
+
+template<typename T>
+T merge(const Real &l, T &a, T &b)
+	{ return a*(1 - l) + b*l; }
+
+inline Real sign(const Real &x)
+	{ return equal(x, 0) ? 0 : x < 0 ? -1 : 1; }
+
+
+int solve_equation(Real *roots, const Real &k0, const Real &k1);
+int solve_equation(Real *roots, const Real &k0, const Real &k1, const Real &k2);
+int solve_equation(Real *roots, const Real &k0, const Real &k1, const Real &k2, const Real &k3);
+
 #endif
diff --git a/c++/vector/vector.h b/c++/vector/vector.h
index ade2d2d..7eec9f0 100644
--- a/c++/vector/vector.h
+++ b/c++/vector/vector.h
@@ -3,6 +3,9 @@
 
 
 #include <cmath>
+
+#include <algorithm>
+
 #include "real.h"
 
 
@@ -18,19 +21,21 @@ public:
 	const Real& operator[] (int i) const { return coords[i]; }
 	Real& operator[] (int i) { return coords[i]; }
 
-	bool operator< (const Vector &a) const { return less(x, a.x) || (!less(a.x, x) && less(y, a.y)); }
-	bool operator== (const Vector &a) const { return equal(x, a.x) && equal(y, a.y); }
-	bool operator!= (const Vector &a) const { return !(*this == a); }
+	bool operator<(const Vector &a) const { return less(x, a.x) || (!less(a.x, x) && less(y, a.y)); }
+	bool operator==(const Vector &a) const { return equal(x, a.x) && equal(y, a.y); }
+	bool operator!=(const Vector &a) const { return !(*this == a); }
 
-	Vector operator+ (const Vector &a) const { return Vector(x + a.x, y + a.y); }
-	Vector operator- (const Vector &a) const { return Vector(x - a.x, y - a.y); }
-	Vector operator* (const Real &a) const { return Vector(x*a, y*a); }
-	Vector operator/ (const Real &a) const { return *this * (Real(1)/a); }
+	Vector operator-() const { return Vector(-x, -y); }
+	
+	Vector operator+(const Vector &a) const { return Vector(x + a.x, y + a.y); }
+	Vector operator-(const Vector &a) const { return Vector(x - a.x, y - a.y); }
+	Vector operator*(const Real &a) const { return Vector(x*a, y*a); }
+	Vector operator/(const Real &a) const { return *this * (Real(1)/a); }
 	
-	Vector& operator+= (const Vector &a) { x += a.x; y += a.y; return *this; }
-	Vector& operator-= (const Vector &a) { x -= a.x; y -= a.y; return *this; }
-	Vector& operator*= (const Real &a) { x *= a; y *= a; return *this; }
-	Vector& operator/= (const Real &a) { return *this *= (Real(1)/a); }
+	Vector& operator+=(const Vector &a) { x += a.x; y += a.y; return *this; }
+	Vector& operator-=(const Vector &a) { x -= a.x; y -= a.y; return *this; }
+	Vector& operator*=(const Real &a) { x *= a; y *= a; return *this; }
+	Vector& operator/=(const Real &a) { return *this *= (Real(1)/a); }
 
 	friend inline Vector operator*(const Real &a, const Vector &b) { return b*a; }
 	friend inline Real dot(const Vector &a, const Vector &b) { return a.x*b.x + a.y*b.y; }
@@ -54,44 +59,87 @@ public:
 };
 
 
-class Segment {
+class RealPair {
 public:
-	Vector p0, p1, t0, t1;
-
-	explicit Segment(
-		const Vector &p0 = Vector(),
-		const Vector &p1 = Vector(),
-		const Vector &t0 = Vector(),
-		const Vector &t1 = Vector()
-	):
-		p0(p0), p1(p1), t0(t0), t1(t1) { }
-
-	bool operator< (const Segment &a) const {
-		return p0 < a.p0 ? true : a.p0 < p0 ? false
-			 : p1 < a.p1 ? true : a.p1 < p1 ? false
-			 : t0 < a.t0 ? true : a.t0 < t0 ? false
-			 : t1 < a.t1;
-	}
-	bool operator> (const Segment &a) const { return a < *this; }
-	bool operator== (const Segment &a) const { return p0 == a.p0 && p1 == a.p1 && t0 == a.t0 && t1 == a.t1; }
-	bool operator!= (const Segment &a) const { return !(*this == a); }
-
-	Vector p(const Real &l) const {
-		return p0*((( 2.0*l - 3.0)*l      )*l + 1.0)
-			 + p1*(((-2.0*l + 3.0)*l      )*l      )
-			 + t0*(((     l - 2.0)*l + 1.0)*l      )
-			 + t1*(((     l - 1.0)*l      )*l      );
+	Real a0, a1;
+	RealPair(const Real &a0, const Real &a1):
+		a0(a0), a1(a1) { }
+	explicit RealPair(const Real &a = Real()):
+		RealPair(a, a) { }
+	Real dist() const { return a1 - a0; }
+	void reset() { a0 = a1 = Real(); }
+};
+
+
+class Range: public RealPair {
+public:
+	using RealPair::RealPair;
+	Real size() const { return dist(); }
+	bool isempty() const { return !less(a0, a1); }
+	Range& expand(const Real &a) {
+		if (a < a0) a0 = a;
+		if (a1 < a) a1 = a;
+		return *this;
 	}
+	Range get_expanded(const Real &a) const
+		{ return Range(std::min(a, a0), std::max(a1, a)); }
+};
 
-	Vector t(const Real &l) const {
-		return p0*(( 6.0*l - 6.0)*l      )
-			 + p1*((-6.0*l + 6.0)*l      )
-			 + t0*(( 3.0*l - 4.0)*l + 1.0)
-			 + t1*(( 3.0*l - 2.0)*l      );
+
+class VectorPair {
+public:
+	Real x0, y0, x1, y1;
+	
+	VectorPair(const Real &x0, const Real &y0, const Real &x1, const Real &y1):
+		x0(x0), y0(y0), x1(x1), y1(y1) { }
+	explicit VectorPair(const Real &x = Real(), const Real &y = Real()):
+		VectorPair(x, y, x, y) { }
+
+	VectorPair(const Vector &p0, const Vector &p1):
+		VectorPair(p0.x, p0.y, p1.x, p1.y) { }
+	explicit VectorPair(const Vector &p):
+		VectorPair(p, p) { }
+
+	VectorPair(const RealPair &x, const RealPair &y):
+		VectorPair(x.a0, y.a0, x.a1, y.a1) { }
+
+	Vector& p0() { return *(Vector*)&x0; };
+	Vector& p1() { return *(Vector*)&x1; };
+	const Vector& p0() const { return *(const Vector*)&x0; };
+	const Vector& p1() const { return *(const Vector*)&x1; };
+	
+	Vector dist() const { return p1() - p0(); }
+	Real dist_x() const { return x1 - x0; }
+	Real dist_y() const { return y1 - y0; }
+	
+	void reset() { p0().reset(); p1().reset(); }
+};
+
+
+class Rect: public VectorPair {
+public:
+	using VectorPair::VectorPair;
+	
+	Vector size() const { return dist(); }
+	Real width() const { return dist_x(); }
+	Real height() const { return dist_y(); }
+	
+	bool isempty() const { return !less(x0, x1) || !less(y0, y1); }
+	
+	Rect& expand(const Real &x, const Real &y) {
+		if (x0 < x) x0 = x;
+		if (y0 < y) y0 = y;
+		if (x1 < x) x1 = x;
+		if (y1 < y) y1 = y;
+		return *this;
 	}
+	Rect& expand(const Vector &p)
+		{ return expand(p.x, p.y); }
 
-	Vector operator() (const Real &l) const  { return p(l); }
-	Segment sub(const Real &a, const Real &b) const { return Segment(p(a), p(b), t(a), t(b)); }
+	Rect get_expanded(const Real &x, const Real &y) const
+		{ return Rect(std::min(x0, x), std::min(y0, y), std::max(x1, x), std::max(y1, y)); }
+	Rect get_expanded(const Vector &p) const
+		{ return get_expanded(p.x, p.y); }
 };