diff --git a/c++/vector/bendarea.cpp b/c++/vector/bendarea.cpp
index bdf160b..3bb2879 100644
--- a/c++/vector/bendarea.cpp
+++ b/c++/vector/bendarea.cpp
@@ -77,6 +77,21 @@ namespace {
 				jump = false;
 			}
 		}
+		
+		void smooth() {
+			if (size() < 2) return;
+			iterator i0, i2 = begin(), i1 = i2++;
+			i2->t0 = i2->p1 - i1->p1;
+			for(i0 = i1, i1 = i2++; i2 != end(); i0 = i1, i1 = i2++) {
+				if (equal(0, dot(i1->t1, i2->t0.perp())) && less(0, dot(i1->t1, i2->t0))) {
+					i1->t1 = i2->t0 = (i1->p1 - i0->p1)*0.5;
+				} else {
+					i1->t1 = i1->p1 - i0->p1;
+					i2->t0 = i2->p1 - i1->p1;
+				}
+			}
+			i1->t1 = i1->p1 - i0->p1;
+		}
 	};
 	
 	class BendContour: public std::vector<BendPoint> {
@@ -186,10 +201,12 @@ namespace {
 			}
 		}
 
-		void bend(Contour &dst, const Contour &src, const Matrix &matrix = Matrix()) const {
+		void bend(Contour &dst, const Contour &src, const Matrix &matrix, int segments) const {
 			if (empty() || src.empty())
 				dst.insert(dst.end(), src.begin(), src.end());
 
+			const int last_segment = segments - 1;
+			const Real step = Real(1)/segments;
 			Contour::const_iterator i = src.begin();
 			Vector p0 = matrix.transform(i->p1);
 			size_t dst_size = dst.size();
@@ -199,25 +216,25 @@ namespace {
 				Hermite h(p0, matrix.transform(i->p1), matrix.transform(i->t0, false), matrix.transform(i->t1, false));
 				p0 = h.p1;
 				
-				Real bends[2];
+				Real bends[3];
 				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;
+				bends_count += h.inflection_x(bends + bends_count);
 				
 				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;
+					for(Real *i = roots, *end = i + count; i != end; ++i)
 						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 )) );
+				for(Real *i = bends, *end = i + bends_count; i != end; ++i)
+					intersections.push_back( Intersection(*i, interpolate( h.p(*i).x )) );
+				Real bsl = step;
+				for(int i = 1; i < last_segment; ++i, bsl += step)
+					intersections.push_back( Intersection(bsl, interpolate( h.p(bsl).x )) );
 				sort(intersections.begin(), intersections.end());
 				
 				Intersection prev(0, interpolate(h.p0.x));
@@ -228,32 +245,35 @@ namespace {
 					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);
+						Real dl = next.point.l - prev.point.l;
+						Vector t0 = flip ? prev.point.t0 : prev.point.t1;
+						Vector t1 = flip ? next.point.t1 : next.point.t0;
 						
-						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;
-						}
+						Hermite bend_h( prev.point.p, next.point.p, t0*dl, t1*dl );
+						Bezier src_b = h.sub(prev.l, next.l).get_bezier();
 						
-						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);
+						Real a = std::min(src_b.p0.x, src_b.p1.x) - 0.1;
+						Real b = std::max(src_b.p0.x, src_b.p1.x) + 0.1;
+						assert( !less(src_b.pp0.x, a) && !less(src_b.pp1.x, a)
+						     && !less(b, src_b.pp0.x) && !less(b, src_b.pp1.x) );
+						
+						Real k = next.point.length - prev.point.length;
+						k = equal(k, 0) ? 0 : 1/k;
+						Real ll0 = (src_b.pp0.x - prev.point.length)*k;
+						Real ll1 = (src_b.pp1.x - next.point.length)*k + 1;
+						Real kk = flip ? -1 : 1;
+						Bezier dst_b(
+							bend_h.p0     + bend_h.d(  0).perp()*src_b.p0.y*kk,
+							bend_h.p1     + bend_h.d(  1).perp()*src_b.p1.y*kk,
+							bend_h.p(ll0) + bend_h.d(ll0).perp()*src_b.pp0.y*kk,
+							bend_h.p(ll1) + bend_h.d(ll1).perp()*src_b.pp1.y*kk );
+
+						half_corner(dst, prev.point, src_b.p0.y, flip, true);
+						dst.hermite( dst_b.get_hermite() );
+						half_corner(dst, next.point, src_b.p1.y, flip, false);
 					}
 					
 					prev = next;
@@ -301,7 +321,7 @@ namespace {
 	}
 
 	void build_bend(BendContour &dst, const ActiveCurve &curve, bool round, int segments) {
-		const Real dl = Real(1)/segments;
+		const Real step = Real(1)/segments;
 		size_t dst_size = dst.size();
 		BendPoint bp;
 		BendMode mode = round ? BendRound : BendCorner;
@@ -321,9 +341,9 @@ namespace {
 			dst.back().t1 = h.t0;
 			
 			Real l0 = index;
-			Real l = dl;
+			Real l = step;
 			bp.mode = BendNone;
-			for(int k = 2; k < segments; ++k, l += dl) {
+			for(int k = 2; k < segments; ++k, l += step) {
 				Vector p = h.p(l);
 				bp.l = l0 + l;
 				bp.length += (p - bp.p).len();
@@ -348,8 +368,8 @@ namespace {
 			last.t1 = first.t1;
 			first.e0 = last.e1 = false;
 		} else {
-			first.t0 = first.t1;
-			last.t1 = last.t0;
+			first.t0 = first.t1.norm();
+			last.t1 = last.t0.norm();
 			last.mode = BendNone;
 		}
 	}
@@ -390,7 +410,7 @@ BendArea::BendArea():
 	b1.x = p1->position.x;
 	b1.y = 780;
 	Vector bt(100, 100);
-	dst_curve.loop = true;
+	//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);
@@ -437,8 +457,9 @@ void BendArea::on_draw_content(const Context &context) {
 				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);
+						  -kx*p0->position.x, -ky*p0->position.y, 16 );
+				bend.bend(dst_contour, src_contour, m, 1);
+				//dst_contour.smooth();
 				
 				dst_contour.put(context, true);
 				context->close_path();
diff --git a/c++/vector/curve.cpp b/c++/vector/curve.cpp
index 2e9e0c3..0908c6b 100644
--- a/c++/vector/curve.cpp
+++ b/c++/vector/curve.cpp
@@ -2,12 +2,27 @@
 #include "curve.h"
 
 
+int Hermite::inflection(Real *l, const Real &p0, const Real &p1, const Real &t0, const Real &t1) {
+	Real root;
+	if (solve_equation(
+		&root,
+		-3*p0 + 3*p1 - 2*t0 -   t1,
+		 6*p0 - 6*p1 + 3*t0 + 3*t1 ))
+	{
+		if (less(0, root) && less(root, 1)) {
+			if (l) *l = root;
+			return 1;
+		}
+	}
+	return 0;
+}
+
 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,
+		                 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)
@@ -46,12 +61,27 @@ int Hermite::intersections(Real *l, const Real &p, const Real &p0, const Real &p
 }
 
 
+int Bezier::inflection(Real *l, const Real &p0, const Real &p1, const Real &pp0, const Real &pp1) {
+	Real root;
+	if (solve_equation(
+		&root,
+		 p0 - 2*pp0 +   pp1     ,
+		-p0 + 3*pp0 - 3*pp1 + p1 ))
+	{
+		if (less(0, root) && less(root, 1)) {
+			if (l) *l = root;
+			return 1;
+		}
+	}
+	return 0;
+}
+
 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     ,
+		2*p0 - 4*pp0 + 2*pp1     ,
 		 -p0 + 3*pp0 - 3*pp1 + p1 );
 	int valid_count = 0;	
 	for(Real *i = roots, *end = i + count; i != end; ++i)
diff --git a/c++/vector/curve.h b/c++/vector/curve.h
index 9b55cd6..9e8d239 100644
--- a/c++/vector/curve.h
+++ b/c++/vector/curve.h
@@ -38,10 +38,12 @@ public:
 
 	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 d(Real l) const {
+		Vector t = this->t(l);
+		Real lensqr = t.lensqr();
+		if (precision*precision < lensqr) return t/sqrt(lensqr);
+		return (p(l + differential) - p(l - differential)).norm();
+	}
 	
 	Vector pp0() const { return p0 + t0/3; }
 	Vector pp1() const { return p1 - t1/3; }
@@ -90,6 +92,11 @@ public:
 	Rect bounds_accurate() const
 		{ return Rect(bounds_accurate_x(), bounds_accurate_y()); }
 
+	int inflection_x(Real *l) const
+		{ return inflection(l, p0.x, p1.x, t0.x, t1.x); }
+	int inflection_y(Real *l) const
+		{ return inflection(l, p0.y, p1.y, t0.y, t1.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
@@ -121,6 +128,7 @@ public:
 	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 inflection(Real *l, const Real &p0, const Real &p1, const Real &t0, const Real &t1);
 	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);
@@ -157,10 +165,12 @@ public:
 
 	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 d(Real l) const {
+		Vector t = this->t(l);
+		Real lensqr = t.lensqr();
+		if (precision*precision < lensqr) return t/sqrt(lensqr);
+		return (p(l + differential) - p(l - differential)).norm();
+	}
 	
 	Vector t0() const { return 3*(pp0 - p0); }
 	Vector t1() const { return 3*(p1 - pp1); }
@@ -178,7 +188,7 @@ public:
 		Vector b1 =  a1*ll +  a2*l;
 		Vector c  =  b0*ll +  b1*l;
 		if (a) {
-			a->p0  = b->p0;
+			a->p0  = p0;
 			a->pp0 = a0;
 			a->pp1 = b0;
 			a->p1  = c;
@@ -187,7 +197,7 @@ public:
 			b->p0  = c;
 			b->pp0 = b1;
 			b->pp1 = a2;
-			b->p1  = b->p1;
+			b->p1  = p1;
 		}
 	}
 	
@@ -223,6 +233,11 @@ public:
 	Rect bounds_accurate() const
 		{ return Rect(bounds_accurate_x(), bounds_accurate_y()); }
 
+	int inflection_x(Real *l) const
+		{ return inflection(l, p0.x, p1.x, pp0.x, pp1.x); }
+	int inflection_y(Real *l) const
+		{ return inflection(l, p0.y, p1.y, pp0.y, pp1.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
@@ -260,6 +275,7 @@ public:
 	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 inflection(Real *l, const Real &p0, const Real &p1, const Real &pp0, const Real &pp1);
 	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);
diff --git a/c++/vector/curvearea.cpp b/c++/vector/curvearea.cpp
new file mode 100644
index 0000000..030bd8c
--- /dev/null
+++ b/c++/vector/curvearea.cpp
@@ -0,0 +1,177 @@
+
+
+#include "random.h"
+
+#include "curvearea.h"
+
+
+CurveArea::CurveArea():
+	hp0(new ActivePoint(Vector(100, 120), Color::white())),
+	hp1(new ActivePoint(Vector(400, 420), Color::white())),
+	ht0(new ActivePoint(Vector(400, 120), Color::yellow())),
+	ht1(new ActivePoint(Vector(800, 420), Color::yellow())),
+	bp0(new ActivePoint(Vector(600, 80), Color::white())),
+	bp1(new ActivePoint(Vector(900, 380), Color::white())),
+	bpp0(new ActivePoint(Vector(700, 80), Color::yellow())),
+	bpp1(new ActivePoint(Vector(800, 380), Color::yellow())),
+	cross(new ActivePoint(Vector(500, 250), Color::gray()))
+{
+	set_size_request(1000, 500);
+	add_point(hp0);
+	add_point(hp1);
+	add_point(ht0);
+	add_point(ht1);
+	add_point(bp0);
+	add_point(bp1);
+	add_point(bpp0);
+	add_point(bpp1);
+	add_point(cross);
+}
+
+void CurveArea::on_point_move(const ActivePoint::Handle &point, const Vector &oldposition, const Vector &newposition) {
+	if (point == hp0)
+		ht0->position += newposition - oldposition;
+	if (point == hp1)
+		ht1->position += newposition - oldposition;
+}
+
+void CurveArea::put_rect(const Context &context, const Rect &rect) {
+	context->move_to(rect.x0, rect.y0);
+	context->line_to(rect.x1, rect.y0);
+	context->line_to(rect.x1, rect.y1);
+	context->line_to(rect.x0, rect.y1);
+	context->close_path();
+}
+
+void CurveArea::put_hmark(const Context &context, const Vector &point) {
+	context->move_to(point.x, point.y - 10);
+	context->line_to(point.x, point.y + 10);
+}
+
+void CurveArea::put_vmark(const Context &context, const Vector &point) {
+	context->move_to(point.x - 10, point.y);
+	context->line_to(point.x + 10, point.y);
+}
+
+void CurveArea::put_hline(const Context &context, const Vector &point) {
+	Real width = get_allocated_width();
+	context->move_to(0, point.y);
+	context->line_to(width, point.y);
+	put_vmark(context, point);
+}
+
+void CurveArea::put_vline(const Context &context, const Vector &point) {
+	Real height = get_allocated_height();
+	context->move_to(point.x, 0);
+	context->line_to(point.x, height);
+	put_hmark(context, point);
+}
+
+void CurveArea::put_cross(const Context &context, const Vector &point) {
+	put_hmark(context, point);
+	put_vmark(context, point);
+}
+
+void CurveArea::on_draw_content(const Context &context) {
+	Real roots[3];
+	
+	context->save();
+
+	context.set_line_width_px(1);
+	context.set_source_rgba(Color::gray().mulalpha(0.5));
+	put_vline(context, cross->position);
+	put_hline(context, cross->position);
+	context->stroke();
+	
+	{
+		context.set_line_width_px(1);
+		Hermite h(
+			hp0->position,
+			hp1->position,
+			ht0->position - hp0->position,
+			ht1->position - hp1->position );
+		
+		context.set_source_rgba(Color::gray());
+		context.move_to(hp0->position);
+		context.line_to(ht0->position);
+		context.move_to(hp1->position);
+		context.line_to(ht1->position);
+		context->stroke();
+
+		context.set_source_rgba(Color::gray().mulalpha(0.25));
+		put_rect(context, h.bounds_accurate());
+		context->stroke();
+		
+		context.set_source_rgba(Color::gray().mulalpha(0.5));
+		for(int i = h.bends_x(roots) - 1; i >= 0; --i)
+			put_cross(context, h(roots[i]));
+		for(int i = h.bends_y(roots) - 1; i >= 0; --i)
+			put_cross(context, h(roots[i]));
+		for(int i = h.intersections_x(roots, cross->position.x) - 1; i >= 0; --i)
+			put_vmark(context, h(roots[i]));
+		for(int i = h.intersections_y(roots, cross->position.y) - 1; i >= 0; --i)
+			put_hmark(context, h(roots[i]));
+		for(int i = h.inflection_x(roots) - 1; i >= 0; --i)
+			put_cross(context, h(roots[i]));
+		for(int i = h.inflection_y(roots) - 1; i >= 0; --i)
+			put_cross(context, h(roots[i]));
+		context->stroke();
+		
+		context.set_source_rgba(Color::blue().withalpha(0.5));
+		context.hermite(h);
+		context->stroke();
+
+		context.set_line_width_px(3);
+		Bezier bb = h.sub(rnd::real(), rnd::real()).get_bezier();
+		context.set_source_rgba(Color::blue());
+		context.bezier(bb);
+		context->stroke();
+	}
+	
+	{
+		context.set_line_width_px(1);
+		Bezier b(
+			bp0->position,
+			bp1->position,
+			bpp0->position,
+			bpp1->position );
+		
+		context.set_source_rgba(Color::gray());
+		context.move_to(bp0->position);
+		context.line_to(bpp0->position);
+		context.line_to(bpp1->position);
+		context.line_to(bp1->position);
+		context->stroke();
+		
+		context.set_source_rgba(Color::gray().mulalpha(0.25));
+		put_rect(context, b.bounds_accurate());
+		context->stroke();
+		
+		context.set_source_rgba(Color::gray().mulalpha(0.5));
+		for(int i = b.bends_x(roots) - 1; i >= 0; --i)
+			put_cross(context, b(roots[i]));
+		for(int i = b.bends_y(roots) - 1; i >= 0; --i)
+			put_cross(context, b(roots[i]));
+		for(int i = b.intersections_x(roots, cross->position.x) - 1; i >= 0; --i)
+			put_vmark(context, b(roots[i]));
+		for(int i = b.intersections_y(roots, cross->position.y) - 1; i >= 0; --i)
+			put_hmark(context, b(roots[i]));
+		for(int i = b.inflection_x(roots) - 1; i >= 0; --i)
+			put_cross(context, b(roots[i]));
+		for(int i = b.inflection_y(roots) - 1; i >= 0; --i)
+			put_cross(context, b(roots[i]));
+		context->stroke();
+		
+		context.set_source_rgba(Color::blue().mulalpha(0.5));
+		context.bezier(b);
+		context->stroke();
+
+		context.set_line_width_px(3);
+		Hermite hh = b.sub(rnd::real(), rnd::real()).get_hermite();
+		context.set_source_rgba(Color::blue());
+		context.hermite(hh);
+		context->stroke();
+	}
+	
+	context->restore();
+}
diff --git a/c++/vector/curvearea.h b/c++/vector/curvearea.h
new file mode 100644
index 0000000..2cb0a4c
--- /dev/null
+++ b/c++/vector/curvearea.h
@@ -0,0 +1,30 @@
+#ifndef CURVEAREA_H
+#define CURVEAREA_H
+
+
+#include "activearea.h"
+#include "activecurve.h"
+
+
+class CurveArea: public ActiveArea {
+protected:
+	ActivePoint::Handle hp0, hp1, ht0, ht1;
+	ActivePoint::Handle bp0, bp1, bpp0, bpp1;
+	ActivePoint::Handle cross;
+	
+	void put_rect(const Context &context, const Rect &rect);
+	void put_hmark(const Context &context, const Vector &point);
+	void put_vmark(const Context &context, const Vector &point);
+	void put_hline(const Context &context, const Vector &point);
+	void put_vline(const Context &context, const Vector &point);
+	void put_cross(const Context &context, const Vector &point);
+
+	void on_point_move(const ActivePoint::Handle &point, const Vector &oldposition, const Vector &newposition) override;
+	void on_draw_content(const Context &context) override;
+	
+public:
+	CurveArea();
+};
+
+
+#endif
diff --git a/c++/vector/mainwindow.cpp b/c++/vector/mainwindow.cpp
index 2d9d1d2..47bce2e 100644
--- a/c++/vector/mainwindow.cpp
+++ b/c++/vector/mainwindow.cpp
@@ -2,6 +2,7 @@
 
 #include "outlinearea.h"
 #include "bendarea.h"
+#include "curvearea.h"
 
 #include "mainwindow.h"
 
@@ -9,6 +10,7 @@
 MainWindow::MainWindow() {
 	//OutlineArea *content = manage(new OutlineArea());
 	BendArea *content = manage(new BendArea());
+	//CurveArea *content = manage(new CurveArea());
 
 	content->show();
 	add(*content);
diff --git a/c++/vector/random.h b/c++/vector/random.h
index 749d501..1e00e86 100644
--- a/c++/vector/random.h
+++ b/c++/vector/random.h
@@ -1,3 +1,5 @@
+#ifndef RANDOM_H
+#define RANDOM_H
 
 
 #include <cstdlib>
@@ -7,21 +9,24 @@
 
 
 namespace rnd {
-	Real real()
+	inline Real real()
 		{ return Real(rand())/Real(RAND_MAX); }
-	Real real(const Real &max)
+	inline Real real(const Real &max)
 		{ return real()*max; }
-	Real real(const Real &min, const Real &max)
+	inline Real real(const Real &min, const Real &max)
 		{ return min + real(max - min); }
 
-	Vector vector()
+	inline Vector vector()
 		{ return Vector(real(), real()); }
-	Vector vector(const Real &maxx, const Real &maxy)
+	inline 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)
+	inline 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)
+	inline Vector vector(const Vector &max)
 		{ return vector(max.x, max.y); }
-	Vector vector(const Vector &min, const Vector &max)
+	inline Vector vector(const Vector &min, const Vector &max)
 		{ return vector(min.x, min.y, max.x, max.y); }
 }
+
+
+#endif
diff --git a/c++/vector/real.cpp b/c++/vector/real.cpp
index 66f3805..ead514a 100644
--- a/c++/vector/real.cpp
+++ b/c++/vector/real.cpp
@@ -1,5 +1,6 @@
 
 #include <cmath>
+#include <iostream>
 
 #include "real.h"
 
@@ -75,7 +76,7 @@ int solve_equation(Real* roots, const Real& k0, const Real& k1, const Real& k2, 
 	} else {
 		if (roots) {
 			Real ph = asinh( fabs(R)/sqrt(-Q3) )/3;
-			roots[0] = -2*sign(R)*sqrt(-Q)*cosh(ph) - a/3;
+			roots[0] = -2*sign(R)*sqrt(-Q)*sinh(ph) - a/3;
 		}
 		return 1;
 	}
diff --git a/c++/vector/real.h b/c++/vector/real.h
index e7a6101..c6e9c40 100644
--- a/c++/vector/real.h
+++ b/c++/vector/real.h
@@ -9,12 +9,14 @@ typedef double Real;
 
 
 const Real precision = 1e-10;
-const Real differential = 1e-5;
+const Real differential = 0.05;
 const Real pi = 3.1415926535897932384626433;
 
 
-inline bool equal(const Real &a, const Real &b)
+inline bool equal(const Real &a, const Real &b, const Real &precision)
 	{ return (a < b ? b - a : a - b) < precision; }
+inline bool equal(const Real &a, const Real &b)
+	{ return equal(a, b, precision); }
 inline bool less(const Real &a, const Real &b)
 	{ return !equal(a, b) && a < b; }
 inline bool greater(const Real &a, const Real &b)
@@ -46,7 +48,8 @@ T merge(const Real &l, T &a, T &b)
 
 inline Real sign(const Real &x)
 	{ return equal(x, 0) ? 0 : x < 0 ? -1 : 1; }
-
+inline Real clamp(const Real &x, const Real &min, const Real &max)
+	{ return x < min ? min : (max < x ? max : x); }
 
 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);
diff --git a/c++/vector/vector.h b/c++/vector/vector.h
index 7eec9f0..f0d8393 100644
--- a/c++/vector/vector.h
+++ b/c++/vector/vector.h
@@ -46,7 +46,7 @@ public:
 	Real lensqr() const { return dot(*this, *this); }
 	Real len() const { return sqrt(lensqr()); }
 
-	Vector norm() const { return iszero() ? Vector() : *this/len(); }
+	Vector norm() const { Real l = len(); return equal(l, 0) ? Vector() : *this/l; }
 	Vector perp() const { return Vector(-y, x); }
 	Vector scale(const Vector &a) const { return Vector(x*a.x, y*a.y); }
 	Vector rotate(const Real &angle) const {