Problem: Problem F: Magnum Tornado

We have a toy that consists of a small racing circuit and a tiny car. For simplicity you can regard the circuit as a 2-dimensional closed loop, made of line segments and circular arcs. The circuit has no branchings. All segments and arcs are connected smoothly, i.e. there are no sharp corners.

The car travels on this circuit with one distinct feature: it is capable of jumping, which enables short cuts. It can jump at any time for any distance. The constraints are that 1) the traveling direction will never change during each jump, 2) the jumping direction needs to match the traveling direction at the time of take-off, and 3) the car must land on the circuit in parallel to the tangent at the landing point. Note that, however, the traveling direction at the landing point may be the opposite of the forward direction (we define forward direction as the direction from the starting point to the ending point of each line segment in the circuit.) That is, the car can traverse part of the circuit in the reverse direction.

Your job is to write a program which, given the shape of the circuit, calculates the per-lap length of the shortest route. You can ignore the height of the jump, i.e. just project the route onto the plane on which the circuit resides. The car must start from the starting point of the first line segment, heading towards the ending point of that segment, and must come back to the same starting point heading in the same direction.

Figure 1 shows the solution for the first sample input.

<image>
Figure 1: The solution for the first sample input



Input

The input begins with a line that solely consists of an integer N (2 <= N <= 100), the number of line segments in the circuit. This is followed by N lines, where the i-th line corresponds to the i-th line segment (1 <= i <= N). Each of these N lines contains 4 integers x0, y0, x1 and y1 (-100 <= x0, y0, x1, y1 <= 100) in this order, separated by a space. Here (x0, y0) is the starting point and (x1, y1) is the ending point of the i-th line segment. For each i, the i-th and (i+1)-th line segments are connected smoothly by a circular arc, which you may assume exists uniquely (for simplicity, we consider the (N+1)-th line as the 1st line). You may also assume that, while two line segments or circular arcs may cross each other, they will never overlap nor be tangent to each other.

Output

For each test case, output one line that solely consists of a decimal value, representing the per-lap length. The output value should be in a decimal fraction and should not contain an error greater than 0.001.

Examples

Input

5
0 1 0 2
1 3 2 3
2 2 1 2
1 1 2 1
2 0 1 0


Output

9.712


Input

12
4 5 4 6
3 7 1 7
0 8 0 10
1 11 3 11
4 10 4 9
5 8 99 8
100 7 100 4
99 3 4 3
3 2 3 1
2 0 1 0
0 1 0 3
1 4 3 4


Output

27.406

Solution (in CPP):

Solution:
#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <climits>
#include <ctime>
#include <queue>
#include <stack>
#include <algorithm>
#include <list>
#include <vector>
#include <set>
#include <map>
#include <iostream>
#include <deque>
#include <complex>
#include <string>
#include <iomanip>
#include <sstream>
#include <bitset>
#include <valarray>
#include <unordered_map>
#include <iterator>
#include <assert.h>
using namespace std;
typedef long long int ll;
typedef unsigned int uint;
typedef unsigned char uchar;
typedef unsigned long long ull;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
typedef vector<int> vi;

#define REP(i,x) for(int i=0;i<(int)(x);i++)
#define REPS(i,x) for(int i=1;i<=(int)(x);i++)
#define RREP(i,x) for(int i=((int)(x)-1);i>=0;i--)
#define RREPS(i,x) for(int i=((int)(x));i>0;i--)
#define FOR(i,c) for(__typeof((c).begin())i=(c).begin();i!=(c).end();i++)
#define RFOR(i,c) for(__typeof((c).rbegin())i=(c).rbegin();i!=(c).rend();i++)
#define ALL(container) (container).begin(), (container).end()
#define RALL(container) (container).rbegin(), (container).rend()
#define SZ(container) ((int)container.size())
#define mp(a,b) make_pair(a, b)
#define UNIQUE(v) v.erase( unique(v.begin(), v.end()), v.end() );

template<class T> bool chmax(T &a, const T &b) { if (a<b) { a=b; return 1; } return 0; }
template<class T> bool chmin(T &a, const T &b) { if (a>b) { a=b; return 1; } return 0; }
template<class T> ostream& operator<<(ostream &os, const vector<T> &t) {
os<<"["; FOR(it,t) {if(it!=t.begin()) os<<","; os<<*it;} os<<"]"; return os;
}
template<class T> ostream& operator<<(ostream &os, const set<T> &t) {
os<<"{"; FOR(it,t) {if(it!=t.begin()) os<<","; os<<*it;} os<<"}"; return os;
}
template<class S, class T> ostream& operator<<(ostream &os, const pair<S,T> &t) { return os<<"("<<t.first<<","<<t.second<<")";}
template<class S, class T> pair<S,T> operator+(const pair<S,T> &s, const pair<S,T> &t){ return pair<S,T>(s.first+t.first, s.second+t.second);}
template<class S, class T> pair<S,T> operator-(const pair<S,T> &s, const pair<S,T> &t){ return pair<S,T>(s.first-t.first, s.second-t.second);}

namespace geom{
#define X real()
#define Y imag()
#define at(i) ((*this)[i])
#define SELF (*this)
	enum {TRUE = 1, FALSE = 0, BORDER = -1};
	typedef int BOOL;
	typedef double R;
	const R INF = 1e8;
	const R EPS = 1e-8;
	const R PI = 3.1415926535897932384626;
	inline int sig(const R &x) { return (abs(x) < EPS ? 0 : x > 0 ? 1 : -1); }
	inline BOOL less(const R &x, const R &y) {return sig(x-y) ? x < y : BORDER;}
	typedef complex<R> P;
	inline R norm(const P &p){return p.X*p.X+p.Y*p.Y;}
	inline R inp(const P& a, const P& b){return (conj(a)*b).X;}
	inline R outp(const P& a, const P& b){return (conj(a)*b).Y;}
	inline P unit(const P& p){return p/abs(p);}
	inline P proj(const P &s, const P &t){return t*inp(s, t)/norm(t);}
	inline int ccw(const P &s, const P &t, const P &p, int adv=0){
		int res = sig(outp(t-s, p-s));
		if(res || !adv) return res;
		if(sig(inp(t-s, p-s)) < 0) return -2;	// p-s-t
		if(sig(inp(s-t, p-t)) < 0) return 2;	// s-t-p
		return 0;								// s-p-t
	}
	
	
	struct L : public vector<P>{	// line
		L(const P &p1, const P &p2){this->push_back(p1);this->push_back(p2);}
		L(){}
		P dir()const {return at(1) - at(0);}
		BOOL online(const P &p)const {return !sig(outp(p-at(0), dir()));}
	};
	struct S : public L{	// segment
		S(const P &p1, const P &p2):L(p1, p2){}
		S(){}
		BOOL online(const P &p)const {
			if(!sig(norm(p - at(0))) || !sig(norm(p - at(1)))) return BORDER;
			return !sig(abs(at(0)-p) + abs(at(1) - p) - abs(at(0) - at(1)));
		}
	};
	struct C : public P{
		C(){}
		C(const P& p, const R r):P(p), r(r){}
		R r;
		BOOL inside(const P& p)const { return less(norm(p-SELF), r*r);}
	};
	struct F : public C{
		R s, t;
		int cc;
		F(const C &c, R ss, R tt, int cc=1):C(c), s(ss), t(tt),cc(cc){
			if(PI < s) s -= 2*PI;
			if(PI < t) t -= 2*PI;
		}
		R length(){
			return ((t < s) ? t - (s-2*PI) : t-s) * r;
		}
		BOOL inside(const P& p)const {
			P v = p - SELF;
			if(!sig(norm(v))) return BORDER;
			R a = arg(v);
			if(t < s){
				if((!less(s, a) && !less(a, t)) || !less(norm(v), r*r)) return FALSE;
				return less(s, a) | less(a, t) | less(norm(v), r*r);
			}else{
				if(!less(s, a) || !less(a, t) || !less(norm(v), r*r)) return FALSE;
				return less(s, a) | less(a, t) | less(norm(v), r*r);
			}
		}
	};
	P crosspoint(const L &l, const L &m);
	struct G : public vector<P>{
		G(size_type size=0):vector(size){}
		S edge(int i)const {return S(at(i), at(i+1 == size() ? 0 : i+1));}
		BOOL contains(const P &p)const {
			R sum = .0;
			REP(i, size()){
				if(S(at(i), at((i+1)%size())).online(p)) return BORDER;	// online
				sum += arg((at(i) - p) / (at((i+1)%size()) - p));
			}
			return !!sig(sum);
		}
		R area()const {
			R sum = 0;
			REP(i, size()) sum += outp(at(i), at((i+1)%size()));
			return abs(sum / 2.);
		}
		
		G convex_hull(bool online = false) {
			if(size() < 2) return *this;
			sort(ALL(*this));
			G r;
			r.resize((int)size()*2);
			int k=0;
			for(int i=0;i<size();r[k++]=at(i++))
				while(k>1 && ccw(r[k-2], r[k-1], at(i)) < 1-online) k--;
			int t = k;
			for(int i=(int)size()-1;i>=0;r[k++]=at(i--))
				while(k>t && ccw(r[k-2], r[k-1], at(i)) < 1-online) k--;
			r.resize(k-1);
			return r;
		}
		G cut(const L &l)const {
			G g;
			REP(i, size()){
				const S &s = edge(i);
				if(ccw(l[0], l[1], s[0], 0) >= 0) g.push_back(s[0]);
				if(ccw(l[0], l[1], s[0], 0) * ccw(l[0], l[1], s[1], 0) < 0)
					g.push_back(crosspoint(s, l));
			}
			return g;
		}
		G Voronoi(const vector<P> &p, const int t)const {
			G g = *this;
			REP(i, p.size())if(i!=t){
				const P m = (p[t]+p[i])*0.5;
				g = g.cut(L(m, m+(p[i]-p[t])*P(0, 1)));
			}
			return g;
		}
	};

	inline P proj(const P &s, const L &t){return t[0] + proj(s-t[0], t[1]-t[0]);}
	inline P reflect(const P &s, const L &t){return 2.*proj(s, t) - s;}
	inline S reflect(const S &s, const L &t){return S(reflect(s[0], t), reflect(s[1], t));}
	BOOL intersect(const S &s, const S &t){
		const int p = ccw(t[0], t[1], s[0], 1) * ccw(t[0], t[1], s[1], 1);
		const int q = ccw(s[0], s[1], t[0], 1) * ccw(s[0], s[1], t[1], 1);
		return (p>0||q>0) ? FALSE : (!p||!q) ? BORDER : TRUE;
	}
	BOOL intersect(const S &s, const L &l){
		if(l.online(s[0]) || l.online(s[1])) return BORDER;
		return (sig(outp(l.dir(), s[0]-l[0])) * sig(outp(l.dir(), s[1]-l[0])) <= 0);
	}
	R dist2(const L &l, const P &p){return norm(outp(l.dir(), p - l[0])) / norm(l.dir());}
	R dist2(const S &s, const P &p){
		if(inp(p-s[0],  s.dir()) < EPS) return norm(p - s[0]);
		if(inp(p-s[1], -s.dir()) < EPS) return norm(p - s[1]);
		return dist2((const L &)s, p);
	}
	R dist2(const S &s, const L &l){
		return intersect(s, l) ? .0 : min(dist2(l, s[0]), dist2(l, s[1]));
	}
	R dist2(const S &s, const S &t){
		return intersect(s, t) ? .0 : min(min(dist2(s, t[0]), dist2(t, s[0])), 
									   	  min(dist2(s, t[1]), dist2(t, s[1])));
	}
	template <class T> R dist2(const G &g, const T& t){ // todo: 内部に完全に含まれる場合
		R res = INF;
		REP(i, g.size()) res = min(res, dist2(g.edge(i), t));
		return res;
	}
	template<class S, class T> R dist(const S& s, const T& t){return sqrt(dist2(s, t));}
	inline BOOL intersect(const C &a, const C &b){
		return less((a.r-b.r)*(a.r-b.r), norm(a-b)) + less(norm(a-b), (a.r+b.r)*(a.r+b.r)) - 1;
	}
	inline BOOL intersect(const C &c, const L &l){
		return less(dist2(l, c), c.r*c.r);
	}
	inline BOOL intersect(const C &c, const S &s){
		int d = less(dist2(s, c), c.r*c.r);
		if(d != TRUE) return d;
		int p = c.inside(s[0]), q = c.inside(s[1]);
		return (p<0 || q<0) ? BORDER : p&q;
	}
	
	/*
	 * CCWでs[0]->s[1]にかけてが共通部分
	 */
	inline S crosspoint(const C &c1, const C &c2){
		if(!intersect(c1, c2)) return S();
		R d = abs(c1 - c2);
		R x = (c1.r*c1.r - c2.r*c2.r + d*d) / (2*d);
		R h = sqrt(c1.r*c1.r - x*x);
		P u = unit(c2-c1);
		return S(c1 + u*x + u*P(0,-1)*h, c1 + u*x + u*P(0,1)*h);
	}
	inline P crosspoint(const L &l, const L &m){
		R A = outp(l.dir(), m.dir()), B = outp(l.dir(), l[1] - m[0]);
		if(!sig(abs(A)) && !sig(abs(B))) return m[0]; // same line
		if(abs(A) < EPS) assert(false); // !!!PRECONDITION NOT SATISFIED!!!
		return m[0] + B / A * (m[1] - m[0]);
	}
	inline S crosspoint(const C &c, const L &l){
		R d2 = dist2(l, c);
		if(c.r*c.r+EPS < d2) return S();
		P m = proj(c, l);
		P u = unit(l[1]-l[0]);
		R d = sqrt(c.r*c.r - d2);
		return S(m+u*d, m-u*d);
	}
	inline vector<P> crosspoint(const C &c, const S &s){
		vector<P> res = crosspoint(c, (const L&)s);
		RREP(i, res.size()){
			if(inp(res[i]-s[0], s.dir())*inp(res[i]-s[1], -s.dir())<EPS)
				res.erase(res.begin() + i);
		}
		return res;
	}
	inline R commonarea(const C &a, const C &b){
		if(less(norm(a-b), (a.r-b.r)*(a.r-b.r)) == TRUE) return min(a.r*a.r, b.r*b.r)*PI;
		if(less((a.r+b.r)*(a.r+b.r), norm(a-b)) == TRUE) return .0;
		double d = abs(a-b);
		double rc = (d*d + a.r*a.r - b.r*b.r) / (2*d);
		double theta = acos(rc / a.r);
		double phi = acos((d - rc) / b.r);
		return a.r*a.r*theta + b.r*b.r*phi - d*a.r*sin(theta);
	}
	vector<L> CommonTangent(C c1, C c2){
		if(c1.r > c2.r) swap(c1, c2);
		double d = abs(c1-c2);
		vector<L> res;
		if(d < EPS) return res;
		if(d + EPS > c1.r + c2.r){
			// 内接線
			P crs = (c1*c2.r + c2*c1.r) / (c1.r + c2.r);
			double rad = asin(c1.r/abs(crs-c1));
			res.push_back(L(crs, crs + (c1-crs)*polar(1., rad)));
			res.push_back(L(crs, crs + (c1-crs)*polar(1., -rad)));
		}
		if(c1.r + d + EPS > c2.r){
			// 外接線
			double rad = 0.5*PI+asin((c2.r-c1.r) / d);
			P v = unit(c2-c1)*polar(1., rad);
			if(c1.r + d - EPS < c2.r){
				res.push_back(L(c1+v*c1.r, c1+v*c1.r+(c1-c2)*P(0, 1)));
			}else{
				res.push_back(L(c1+v*c1.r, c2+v*c2.r));
				v = 2.*proj(v, c2-c1) - v;
				res.push_back(L(c1+v*c1.r, c2+v*c2.r));
			}
		}
		return res;
	}
	
	
	struct Arrangement{
		struct AEdge{
			int u, v, t;
			R cost;
			AEdge(int u=0, int v=0, int t=0, R cost=0)
				:u(u), v(v), t(t), cost(cost){}
		};
		typedef vector<vector<AEdge>> AGraph;
		vector<P> p;
		AGraph g;
		Arrangement(){}
		Arrangement(vector<S> seg){
			int m = seg.size();
			REP(i, m){
				p.push_back(seg[i][0]);
				p.push_back(seg[i][1]);
				REP(j, i) if(sig(outp(seg[i].dir(), seg[j].dir())) && intersect(seg[i], seg[j]) == TRUE)
					p.push_back(crosspoint(seg[i], seg[j]));
			}
			sort(ALL(p)); UNIQUE(p);
			int n=p.size();
			g.resize(n);
			REP(i, m){
				S &s = seg[i];
				vector<pair<R, int>> ps;
				REP(j, n) if(s.online(p[j])) ps.emplace_back(norm(p[j] - s[0]), j);
				sort(ALL(ps));
				REP(j, (int)ps.size()-1){
					const int u=ps[j].second;
					const int v=ps[j+1].second;
					g[u].emplace_back(u, v, 0, abs(p[u] - p[v]));
					g[v].emplace_back(v, u, 0, abs(p[u] - p[v]));
				}
			}
		}
		
		int getIdx(P q){
			auto it = lower_bound(ALL(p), q);
			if(it == p.end() || *it != q) return -1;
			return it - p.begin();
		}
	};

	struct min_ball {
		P center;
		R radius2;
		min_ball(const vector<P>& p) {
			FOR(it, p) ps.push_back(*it);
		}
		min_ball& compile() {
			m = 0;
			center = P(0, 0);
			radius2 = -1;
			make_ball(ps.end());
			return *this;
		}
	private:
		list<P> ps;
		list<P>::iterator supp_end;
		int m;
		P v[3], c[3];
		R z[3], r[3];
		void pop() { --m; }
		void push(const P& p) {
			if (m == 0) {
				c[0] = p; r[0] = 0;
			} else {
				R e = norm(p-c[m-1]) - r[m-1];
				P delta = p - c[0];
				v[m] = p - c[0];
				for (int i = 1; i < m; ++i)
					v[m] -= v[i] * inp(v[i], delta) / z[i];
				z[m] = inp(v[m], v[m]);
				c[m] = c[m-1] + e*v[m]/z[m]*.5;
				r[m] = r[m-1] + e*e/z[m]*.25;
			}
			center	= c[m];
			radius2 = r[m]; ++m;
		}
		void make_ball(list<P>::iterator i) {
			supp_end = ps.begin();
			if (m == 3) return;
			for (list<P>::iterator k = ps.begin(); k != i; ) {
				list<P>::iterator j = k++;
				if (norm(*j-center) > radius2) {
					push(*j); make_ball(j); pop();
					move_to_front(j);
				}
			}
		}
		void move_to_front(list<P>::iterator j) {
			if (supp_end == j) ++supp_end;
			ps.splice (ps.begin(), ps, j);
		}
	};
	struct DualGraph{
		struct DEdge{
			int u, v, f, l;
			R a;
			DEdge(int u, int v, R a):u(u),v(v),f(0),l(0){
				while(PI < a) a -= 2*PI;
				while(a < -PI) a += 2*PI;
				this->a = a;
			}
			bool operator<(const DEdge &opp) const{
				return a>opp.a;
			}
			bool operator<(const R &opp) const{
				return a>opp;
			}
			friend ostream& operator<<(ostream &os, const DEdge &t) { return os<<"("<<t.u<<","<<t.v<<","<<t.a*180/PI<<")";}
		};
		
		int n;
		vector<P> p;
		vector<vector<DEdge>> g;
		DualGraph(const vector<P> &p):p(p),g(p.size()),n(p.size()){}
		void add_edge(int s, int t){
			R a = arg(p[t]-p[s]);
			g[s].emplace_back(s, t, a);
			g[t].emplace_back(t, s, a+PI);
		}
		vector<G> poly;
		
		void add_polygon(int s, int t, R a){
			auto e = lower_bound(ALL(g[s]), a-EPS);
			if(e == g[s].end()) e = g[s].begin();
			if(e->f) return;
			e->f = 1;
			e->l = t;
			poly[t].push_back(p[s]);
			add_polygon(e->v, t, e->a > 0 ? e->a-PI : e->a+PI);
		}
		
		vector<vi> dual(){
			REP(i, n) sort(ALL(g[i]));
			int s = min_element(ALL(p)) - p.begin();
			poly.emplace_back();
			add_polygon(s, poly.size()-1, -PI*.5);
			REP(i, n)REP(j, g[i].size())if(!g[i][j].f){
				poly.emplace_back();
				add_polygon(i, poly.size()-1, g[i][j].a+EPS);
			}
			int m = poly.size();
			vector<vi> dg(m);
			vector<unordered_map<int, int>> rev(n);
			REP(i, n)REP(j, g[i].size()){
				int u = i, v = g[i][j].v, l = g[i][j].l;
				if(u>v) swap(u, v);
				if(rev[u].find(v) != rev[u].end()){
					dg[l].push_back(rev[u][v]);
					dg[rev[u][v]].push_back(l);
				}
				else rev[u][v] = l;
			}
			return dg;
		}
	};

#undef SELF
#undef at
}

using namespace geom;


namespace std{
	bool operator<(const P &a, const P &b){return sig(a.X-b.X) ? a.X < b.X : a.Y+EPS < b.Y;}
	bool operator==(const P &a, const P &b){return !sig(norm(a-b));}
	istream& operator>>(istream &is, P &p){R x,y;is>>x>>y;p=P(x, y);return is;}
	istream& operator>>(istream &is, L &l){l.resize(2);return is >> l[0] >> l[1];}
	istream& operator>>(istream &is, C &c){return is >> (P &)c >> c.r;}
/*	const double B = 200;
	const double Z = 20;
	ostream& operator<<(ostream &os, const C &c){return os << "circle("<<B+Z*(c.X)<<", "<<1000-B-Z*(c.Y)<<", "<<Z*(c.r)<<")";}
	ostream& operator<<(ostream &os, const P &p){return os << C(p, 2.);}
	ostream& operator<<(ostream &os, const S &s){return os << "line("<<B+Z*(s[0].X)<<", "<<1000-B-Z*(s[0].Y)<<", "<<B+Z*(s[1].X)<<", "<<1000-B-Z*(s[1].Y)<<")";}
	ostream& operator<<(ostream &os, const G &g){REP(i, g.size()) cout << g.edge(i) << endl;return os;}
*/
}

int n;

R dijkstra(const vector< vector<pair<int, R>> > &g, int start){
	int n = g.size();
	vector<R> ret(n, INF);
	priority_queue<pair<R, int>, vector<pair<R, int>>, greater<pair<R, int>> > dp;
	dp.push(pair<R, int>(0., start));
	while(!dp.empty()){
		int s = dp.top().second;
		R c = dp.top().first;
		dp.pop();
		if(s == start && c != 0) return c;
		if(s != start && ret[s] < c - EPS) continue;
		FOR(it, g[s]){
			int t = it->first;
			if(ret[t] - EPS > c+it->second){
				ret[t] = c+it->second;
				dp.push(pair<R, int>(ret[t], t));
			}
		}
	}
	return INF;
}


int main(){
	ios::sync_with_stdio(false);
	cin >> n;
	vector<pair<pair<P, P>, double>> e;
	vector<P> p;
	P fst, prev;
	vector<S> s(n);
	REP(i, n)cin >> s[i];
	vector<F> f;
	REP(i, n){
		const S &cur = s[i];
		const S &nxt = s[(i+1) % n];
		C c;
		(P&)c = sig(outp(cur.dir(), nxt.dir())) ? 
			   crosspoint(L(cur[1], cur[1]+cur.dir()*P(0, 1)), L(nxt[0], nxt[0]+nxt.dir()*P(0, 1)))
			: (cur[1] + nxt[0]) * .5;
		c.r = abs(cur[1] - c);
		R ca = arg(cur[1] - c);
		R cb = arg(nxt[0] - c);
		int cc = ccw(cur[0], cur[1], c);
		if(cc == -1) swap(ca, cb);
		f.emplace_back(c, ca, cb, cc);
		p.push_back(cur[0]);
		p.push_back(cur[1]);
		p.push_back(cur[0]+P(1000, 0));
		p.push_back(cur[1]+P(1000, 0));
		e.emplace_back(mp(cur[0], cur[1]), abs(cur.dir()));
		e.emplace_back(mp(cur[1]+P(1000, 0), cur[0]+P(1000, 0)), abs(cur.dir()));
		e.emplace_back(mp(cur[1], nxt[0]), f.back().length());
//		cout << f.back() << endl;
//		cout << s[i] << endl;
	}
	REP(i, n)REP(j, n)if(i!=j){
		const F &u = f[i];
		const F &v = f[j];
		vector<L> ls = CommonTangent(u, v);
		FOR(s, ls){
			S root(proj(u, *s), proj(v, *s));
			if(!u.inside(root[0]) || !v.inside(root[1])) continue;
			int f1 = ccw(root[0], root[1], u) != u.cc;
			int f2 = ccw(root[0], root[1], v) != v.cc;
//			if(ccw(root[0], root[1], u) != u.cc || ccw(root[0], root[1], v) != v.cc) continue;
//			cerr << root[0].X << "," << root[0].Y << " -> " << root[1].X << "," << root[1].Y << endl;
			p.push_back(root[0]);
			p.push_back(root[1]);
			p.push_back(root[0]+P(1000, 0));
			p.push_back(root[1]+P(1000, 0));
			e.emplace_back(mp(root[0]+(f1?P(1000, 0):P(0,0)), 
							  root[1]+(f2?P(1000, 0):P(0,0))), abs(root.dir()));
			e.emplace_back(mp(root[1]+(!f2?P(1000, 0):P(0,0)), 
							  root[0]+(!f1?P(1000, 0):P(0,0))), abs(root.dir()));
//			if(i == 3 && j == 4) cout << "!" << root[0] << "," << root[1] << endl;
		}
	}
	sort(ALL(p));UNIQUE(p);
	
	REP(i, n){
		vector<pair<R, int>> q;
		P st = s[i][1] - f[i];
		REP(j, p.size()) if(f[i].inside(p[j]) && !sig(norm(f[i]-p[j]) - f[i].r*f[i].r)){
			R ar = sig(norm(p[j]-s[i][1])) ? arg(((p[j]-f[i])/st) * P(-1, 0)) : f[i].cc == 1 ? -PI : PI;
			q.emplace_back(ar, j);
		}
		if(f[i].cc == 1) sort(ALL(q));
		else sort(RALL(q));
//		cout << q << endl;
		REP(j, (int)q.size() - 1){
//			cout << "!" << p[q[j].second] << "->" << p[q[j+1].second] << ": " << abs(q[j].first - q[j+1].first)*180/PI << endl;
			e.emplace_back(mp(p[q[j].second], p[q[j+1].second]), f[i].r * abs(q[j].first - q[j+1].first));
			e.emplace_back(mp(p[q[j+1].second]+P(1000, 0), p[q[j].second]+P(1000, 0)), f[i].r * abs(q[j].first - q[j+1].first));
		}
	}
	
	n = p.size();
	vector<vector<pair<int, R>>> g(n);
	int start = lower_bound(ALL(p), s[0][0]) - p.begin();
	sort(ALL(e));UNIQUE(e);
	FOR(it, e){
//		if(it->second < 21)cout << it->first.first << " -> " << it->first.second << ": " << it->second << endl; 
		int u = lower_bound(ALL(p), it->first.first) - p.begin();
		int v = lower_bound(ALL(p), it->first.second) - p.begin();
		g[u].emplace_back(v, it->second);
	}
	
	printf("%.14f\n", dijkstra(g, start));
	return 0;
}