1 #include <iostream>
   2 #include <vector>
   3 #include <map>
   4 #include <algorithm>
   5 #include <queue>
   6 #include <cstring>
   7 #include <complex>
   8 
   9 using namespace std;
  10 
  11 #define rep(i,n) for(int i=0;i<n;++i)
  12 #define pb(a) push_back(a)
  13 #define mp(a,b) make_pair(a,b)
  14 #define F first
  15 #define S second
  16 #define SZ(a) (int)((a).size())
  17 

   1 // primitives
   2 typedef long double Real;
   3 typedef complex<Real> Pt;
   4 typedef pair<Pt, Pt> Line;
   5 typedef pair<Pt, Real> Circle;
   6 typedef vector<Pt> Poly;
   7 
   8 Real eps = 1e-9;
   9 
  10 inline istream& operator>>(istream& s, Pt& p) {return s >> p.real() >> p.imag();}
  11 inline bool near(const Pt& p, const Pt& q){return abs(p - q) < eps;}

   1 inline Real dot(const Pt& p, const Pt& q){return p.real() * q.real() + p.imag() * q.imag();}
   2 inline Real cross(const Pt& p, const Pt& q){return p.real() * q.imag() - p.imag() * q.real();}

   1 // ccw_b : old version of ccw.
   2 inline Real ccw_b(const Pt& p, const Pt& q, const Pt& r) {return cross(q - p, r - p);}
   3 
   4 /* ccw :
   5 CD  : counter direction
   6 CW  : clock wise
   7 OS  : on segment
   8 CCW : counter clock wise
   9 D   : direction
  10  */
  11 enum LPposit { P_CD = -2, P_CW = -1, P_OS = 0, P_CCW = 1, P_D = 2};
  12 LPposit ccw(const Pt& p, const Pt& q, const Pt& r) {
  13   Real c = ccw_b(p, q, r);
  14   if (c < -eps) return P_CW;
  15   if (c >  eps) return P_CCW;
  16   if (dot(q - p, r - p) < -eps) return P_CD;
  17   if (dot(p - q, r - q) < -eps) return P_D;
  18   return P_OS;
  19 }

   1 inline Real Sabs(const Line& l) {return abs(l.first - l.second); }
   2 inline Real LPdist(const Line& l, const Pt& p) {return abs(ccw_b(l.first, l.second, p)) / Sabs(l); }
   3 inline Real SPdist(const Line& l, const Pt& p) {
   4   if ((Labs(l) > abs(l.first - p) && Labs(l) > abs(l.second - p)) ||
   5       cos(arg(p - l.first) - arg(p - l.second)) > 0){
   6 
   7     return LPdist(l, p);
   8   }
   9 
  10   return min(abs(p - l.first), abs(p - l.second));
  11 }

   1 bool crossS(const Line& p, const Line& q){
   2   return
   3     ccw(p.first, p.second, q.first) * ccw(p.first, p.second, q.second) <= 0 &&
   4     ccw(q.first, q.second, p.first) * ccw(q.first, q.second, p.second) <= 0;
   5 }

   1 Pt intersect(const Line& p, const Line& q) {
   2   Pt vp = p.second - p.first;
   3   Pt vq = q.second - q.first;
   4   Pt c(cross(vp, p.first), cross(vq, q.first));
   5   return Pt(cross(c, Pt(vp.real(), vq.real())), cross(c, Pt(vp.imag(), vq.imag()))) / cross(vp, vq);
   6 }

   1 Real area(const Poly& p) {
   2   Real ret = 0.0;
   3   for(Poly::iterator it = p.begin(); it != p.end(); it++) ret += cross(*it, *(it + 1));
   4   return ret / 2;
   5 }

   1 vector<pt> circle_line_intersect(line l,circle c){
   2   vector<pt> ret;
   3   LD di = dist(l,c.first);
   4   LD r=c.second;
   5   if(di+EPS > r) return ret;
   6   pt v=(l.second-l.first);
   7   v/=abs(v);  
   8   pt rv=v*pt(0,1);
   9   rv*=di;  
  10   if(dist(l,c.first+rv) > di+EPS) rv = -rv;
  11   v*=sqrt(r*r-di*di);
  12   ret.push_back(c.first+rv-v);
  13   ret.push_back(c.first+rv+v);
  14   return ret;
  15 }

   1 vector<line> contact(circle p, circle q){
   2   vector<line> ret;
   3   if(p.second < q.second) swap(p, q);
   4   long double d = abs(p.first - q.first);
   5   pt n = q.first - p.first;
   6   n /= abs(n);
   7 
   8   if(d + eps < abs(p.second - q.second)){
   9     ret.clear();
  10   } else if(eq(d, abs(p.second - q.second))){
  11     pt t, u;
  12     t = p.first + p.second * n;
  13     u = t + n * pt(0, 1);
  14     ret.push_back(make_pair(t, u));
  15   } else {
  16     if(!eq(p.second, q.second)){
  17       pt t = p.first + (p.second * d / (p.second - q.second)) * n;
  18       long double theta = asin((p.second - q.second) / d);
  19       pt u = n * pt(cos(theta), sin(theta));
  20       pt v = n * pt(cos(-theta), sin(-theta));
  21       u += t;
  22       v += t;
  23       ret.push_back(make_pair(t, u));
  24       ret.push_back(make_pair(t, v));
  25     } else {
  26       pt t = p.first + n * pt(0, 1) * p.second;
  27       pt u = p.first - n * pt(0, 1) * p.second;
  28       ret.push_back(make_pair(t, t+n));
  29       ret.push_back(make_pair(u, u+n));
  30     }
  31 
  32     if(eq(d, p.second + q.second)){
  33       pt t, u;
  34       t = p.first + p.second * n;
  35       u = t + n * pt(0, 1);
  36       ret.push_back(make_pair(t, u));
  37     } else if(d > p.second + q.second){
  38       pt t = p.first + (p.second * d / (p.second + q.second)) * n;
  39       long double theta = asin((p.second + q.second) / d);
  40       pt u = n * pt(cos(theta), sin(theta));
  41       pt v = n * pt(cos(-theta), sin(-theta));
  42       u += t;
  43       v += t;
  44       ret.push_back(make_pair(t, u));
  45       ret.push_back(make_pair(t, v));
  46     }
  47   }
  48   return ret;
  49 }

   1 // convex-hull                                                                  
   2 // return minimum convex polygon that contains every point in vs.               
   3 // depend on ccw                                                                
   4 namespace std{
   5   bool operator<(const Pt& a, const Pt& b) {
   6     return a.real() != b.real() ? a.real() < b.real() : a.imag() < b.imag();
   7   }
   8 }
   9 
  10 Poly convexHull(vector<Pt> ps) {
  11   int n = ps.size();
  12   sort(ps.begin(), ps.end());
  13 
  14   Poly ret(2 * n);
  15   int m = 0;
  16   for (int i = 0; i < n; i++) {
  17     while (m >= 2 && ccw(ret[m-2], ret[m-1], ps[i]) < 0) m--;
  18     ret[m++] = ps[i];
  19   }
  20 
  21   int t = m;
  22   for (int i = n-2; i >= 0; i--) {
  23     while (m >= t && ccw(ret[m-2], ret[m-1], ps[i]) < 0) m--;
  24     ret[m++] = ps[i];
  25   }
  26 
  27   ret.resize(m - 1);
  28   return ret;
  29 }

   1 long long powmod(long long x, long long p, long long m){
   2   if(p == 0) return 1;
   3   long long rt = powmod(x, p/2, m);
   4   if(p % 2 == 0){
   5     return rt * rt % m;
   6   } else {
   7     return (rt * rt % m) * x % m;
   8   }
   9 }
  10 
  11 int miller_test[9] = {2,3,5,7,11,13,17,19,23};
  12 bool isprime(long long n){
  13   if(n <= ISPR_MAX) return !isnpr[n];
  14 
  15   long long d = n-1, s = 0;
  16   while(d%2 == 0){
  17     d /= 2;
  18     s++;
  19   }
  20 
  21   for(int i=0; i<9; i++){
  22     bool iscomp = true;
  23     li x = powmod(miller_test[i], d, n);
  24     iscomp = iscomp && (x % n != 1);
  25 
  26     for(int r=0; r<s; r++){
  27       iscomp = iscomp && (x % n != n-1);
  28       x = x * x % n;
  29     }
  30     if(iscomp) return false;
  31   }
  32 
  33   return true;
  34 }

   1 long long myrand(long long c, long long n, long long x){
   2   return (x * x + c) % n;
   3 }
   4  
   5 long long pollard(long long n){
   6   long long x = 2, y = 2, d = 1, c = rand();
   7   while(d == 1){
   8     x = myrand(c, n, x);
   9     y = myrand(c, n, myrand(c, n, y));
  10     d = __gcd(abs(x-y), n);
  11   }
  12   return d;
  13 }

   1 typedef vector<vector<li> > matrix;
   2 const li mod = 1000000007;
   3 
   4 matrix ident(const int& n) {
   5     matrix ret(n, vi(n, 0));
   6     rep(i, n) ret[i][i] = 1;
   7     return ret;
   8 }
   9 
  10 matrix matadd(const matrix& p, const matrix& q) {
  11     int n = p.size(), m = p[0].size();
  12     matrix ret = p;
  13     rep(i, n) rep(j, m) ret[i][j] = (p[i][j] + q[i][j]) % mod;
  14     return ret;
  15 }
  16 
  17 matrix matmul(const matrix& p, const matrix& q) {
  18     int n = p.size(), m = q[0].size(), l = p[0].size();
  19     matrix ret = matrix(n, vi(m, 0));
  20     rep(i, n) rep(j, m) rep(k, l) ret[i][j] = (ret[i][j] + p[i][k] * q[k][j] % mod) % mod;
  21     return ret;
  22 }
  23 
  24 matrix matpow(const matrix& p, const li& x) {
  25     if (x == 0) return ident(p.size());
  26     if (x == 1) return p;
  27     matrix ret = matpow(p, x / 2);
  28     ret = matmul(ret, ret);
  29     if (x % 2) ret = matmul(ret, p);
  30     return ret;
  31 }

   1 typedef double Real;
   2 typedef complex<Real> Complex;
   3 
   4 vector<Complex> fastFourierTransform(const vector<Complex>& x) {
   5     if (__builtin_popcount(x.size()) > 1) {
   6         cerr << "fastFourierTransform: x.size() should be power of 2." << endl;
   7         return x;
   8     }
   9 
  10     const int n = x.size();
  11     if (n == 1) {
  12         return x;
  13     }
  14 
  15     vector<Complex> even(n / 2), odd(n / 2), ret(n);
  16     for (int i = 0; i < n; ++i) {
  17         if (i % 2 == 0) {
  18             even[i / 2] = x[i];
  19         } else {
  20             odd[i / 2] = x[i];
  21         }
  22     }
  23 
  24     even = fastFourierTransform(even);
  25     odd = fastFourierTransform(odd);
  26 
  27     for (int i = 0; i < n; ++i) {
  28         int j = i % (n / 2);
  29         ret[i] = even[j] + odd[j] * exp(Complex(0, -2 * M_PI * i / n));
  30     }
  31 
  32     return ret;
  33 }
  34 
  35 vector<Complex> inverseFastFourierTransform(vector<Complex> x) {
  36     const int n = x.size();
  37     for (int i = 0; i < n; ++i) {
  38         x[i] = conj(x[i]);
  39     }
  40     x = fastFourierTransform(x);
  41     for (int i = 0; i < n; ++i) {
  42         x[i] = conj(x[i]) / Complex(n);
  43     }
  44     return x;
  45 }
  46 
  47 int conv_minpow2(int x) {
  48     int ret = 1;
  49     while (ret < x) ret <<= 1;
  50     return ret;
  51 }
  52 
  53 vector<Complex> convolution(vector<Complex> x,
  54                             vector<Complex> y) {
  55     int l = 2 * max(conv_minpow2(x.size()), conv_minpow2(y.size()));
  56     x.resize(l);
  57     y.resize(l);
  58 
  59     x = fastFourierTransform(x);
  60     y = fastFourierTransform(y);
  61 
  62     vector<Complex> ret(l);
  63     for (int i = 0; i < l; i++) {
  64         ret[i] = x[i] * y[i];
  65     }
  66 
  67     ret = inverseFastFourierTransform(ret);
  68     while (not ret.empty() && ret.back() == Complex(0)) ret.pop_back();
  69     return ret;
  70 }

   1 // linear programming by simplex method
   2 // solve maximize(cx) s.t. Ax <= a
   3 // if not bounded, return -1.
   4 double simplex(vector<double>& c, 
   5                vector<vector<double> >& A,
   6                vector<double>& a) {
   7   // introduce slack variable
   8   int n = c.size(), m = a.size(), dim = n + m;
   9   vector<int> N(n), B(m);
  10   for (int i = 0; i < n; ++i) N[i] = i;
  11   for (int i = 0; i < m; ++i) {
  12     B[i] = i + n;
  13     c.push_back(0);
  14     for (int j = 0; j < n; j++) A[i][j] *= -1;
  15     for (int j = 0; j < m; j++) A[i].push_back(0);
  16   }
  17 
  18   double ans = 0;
  19 
  20   while (true) {
  21     // check optimized or not
  22     int s = -1;
  23     for (int i = 0; i < n; i++) if (c[N[i]] > 0) s = N[i];
  24     if (s < 0) break;
  25 
  26     // check bounded or not
  27     double bound = 1e300;
  28     int r = -1;
  29     for (int i = 0; i < m; i++) {
  30       if (A[i][N[s]] < 0) {
  31         double nbound = -a[i] / A[i][N[s]];
  32         if (nbound < bound) {
  33           bound = nbound;
  34           r = i;
  35         }
  36       }
  37     }
  38 
  39     if (r < 0) return -1;
  40 
  41     // pivotting
  42     for (int i = 0; i < dim; i++) if (i != N[s]) A[r][i] /= -A[r][N[s]];
  43     a[r] /= -A[r][N[s]];
  44     A[r][B[r]] = 1.0 / A[r][N[s]];
  45     A[r][N[s]] = 0;
  46 
  47     for (int i = 0; i < m; i++) {
  48       if (i == r) continue;
  49 
  50       for (int j = 0; j < dim; j++) if (j != N[s]) A[i][j] += A[i][N[s]] * A[r][j];
  51       a[i] += A[i][N[s]] * a[r];
  52       A[i][N[s]] = 0;
  53     }
  54 
  55     ans += c[N[s]] * a[r];
  56     for (int i = 0; i < dim; i++) if (i != N[s]) c[i] += c[N[s]] * A[r][i];
  57     c[N[s]] = 0;
  58 
  59     swap(N[s], B[r]);
  60   }
  61 
  62   return ans;
  63 }