struct point{
    double x, y;
    point () {}
    point (double a, double b){
        x = a; y = b;
    }
};

struct line{
   point p1;
   point p2;
};

double cross(point a, point b, point c){
   return (b.x-a.x)*(c.y-a.y)-(c.x-a.x)*(b.y-a.y);
}

bool cmp(point a, point b){
   if (a.x != b.x) return (a.x<b.x);
   else return (a.y < b.y);
}

void convexhull(int size){
   int k = 0;
   sort(pt,pt+size,cmp);
   s.resize(2*size);
   for (int i = 0; i < size; i++){
      while (k >= 2 && cross(s[k-2],s[k-1], pt[i])<=0){
         k--;
      }
      s[k++] = pt[i];
   }
   for (int i = size-2, t = k+1; i >= 0; i--){
      while (k >= t && cross(s[k-2],s[k-1],pt[i]) <= 0){
         k--;
      }
      s[k++] = pt[i];
   }
   s.resize(k);
   //s[s.size()-1] = s[0];
}

int dblcmp(double d){
   if (d > 0.000001)
      return 1;
   if (d < -0.000001)
      return -1;
   return 0;
}

int btwcmp(point a, point b, point c){
   return dblcmp(dot(a,b,c));
}

double dot(point a, point b, point c){
   return (b.x-a.x)*(c.x-a.x)+(b.y-a.y)*(c.y-a.y);
}

int intersect(line l1, line l2){
   double s1,s2,s3,s4;
   int d1,d2,d3,d4;
   s1 = cross(l1.p1,l1.p2,l2.p1);
   s2 = cross(l1.p1,l1.p2,l2.p2);
   s3 = cross(l2.p1,l2.p2,l1.p1);
   s4 = cross(l2.p1,l2.p2,l1.p2);
   d1 = dblcmp(s1); d2 = dblcmp(s2); d3 = dblcmp(s3); d4 = dblcmp(s4);
   if (d1*d2<0 && d3*d4<0){
      meet.x = (l2.p1.x*s2 - l2.p2.x*s1)/(s2-s1);
      meet.y = (l2.p1.y*s2 - l2.p2.y*s1)/(s2-s1);
      return 1;
   }
   if ((d1 == 0 && btwcmp(l2.p1,l1.p1,l1.p2)<=0) ||
      (d2 == 0 && btwcmp(l2.p2,l1.p1,l1.p2)<=0) ||
      (d3 == 0 && btwcmp(l1.p1,l2.p1,l2.p2)<=0) ||
      (d4 == 0 && btwcmp(l1.p2,l2.p1,l2.p2)<=0))
      return 1;
   return 0;
}

double angle(point a, point b, point c){
   double ux,uy,vx,vy;
   ux = b.x-a.x;
   uy = b.y-a.y;
   vx = c.x-a.x;
   vy = c.y-a.y;
   return acos((ux*vx+uy*vy)/sqrt((ux*ux+uy*uy)*(vx*vx+vy*vy)));
}

bool ptinpolygon(point p, vector <point> &pt){
   double sum = 0;
   for (int i = 0; i < pt.size()-1; i++){
      if (cross(p,pt[i],pt[i+1]) < 0) sum -= angle(p,pt[i],pt[i+1]);
      else sum+=angle(p,pt[i],pt[i+1]);
   }
   if (fabs(sum-2*pi)<1e-8 || fabs(sum+2*pi)<1e-8) return true;
   return false;
}

bool ptinpolygon(point p){
    int left = 1, right = s.size()-2, mid, ans = s.size();
    if (fabs(p.x-s[0].x) < 1e-9 && fabs(p.y-s[0].y) < 1e-9) return true;
    if (fabs(p.x-s[s.size()-2].x) < 1e-9 && fabs(p.y-s[s.size()-2].y) < 1e-9) return true;
    while (left <= right){
        mid = (left+right)/2;
        if (cross(s[0], p, s[mid]) >= 0){
            ans = min(ans,mid);
            right = mid-1;
        }else{
            left = mid+1;
        }
    }
    if (ans == 1){ 
        if (fabs(cross(s[0], s[1], p)) < 1e-9 && dist(s[0],p) <= dist(s[0],s[1])) return true;
        return false;
    }
    if (ans == s.size()){ 
        if (fabs(cross(s[0], s[s.size()-2], p)) < 1e-9 && dist(s[0],p) <= dist(s[0],s[s.size()-2])) return true;
        return false;
    }
    double sum = fabs(cross(s[0],s[ans-1],p)) + fabs(cross(s[0],p,s[ans])) + fabs(cross(p,s[ans-1],s[ans]));
    if (fabs(sum-cross(s[0], s[ans-1], s[ans])) < 1e-9) return true;
    return false;
}

//Line segment p-q intersect with line A-B.
point lineIntersectSeg(point p, point q, point A, point B){
    double a = B.y-A.y;
    double b = A.x-B.x;
    double c = B.x * A.y - A.x * B.y;
    double u = fabs(a*p.x + b*p.y + c);
    double v = fabs(a*q.x + b*q.y + c);
    return point((p.x*v + q.x*u)/(u+v), (p.y*v + q.y*u)/(u+v));
}

//Cut up the polygon q with the along the line formed by a->b.
//a & b: two points that form a line segment
//q: given vector of point (the last point must be the same as first point)
void cutPolygon(point a, point b, vector <point> &q){
    vector <point> P;
    for (int i = 0; i < q.size(); i++){
        double left1 = cross(a,b,q[i]);
        double left2 = 0.0;
        if (i != q.size()-1)
            left2 = cross(a,b,q[i+1]);
        if (left1 > -1e-9){
            P.push_back(q[i]);
        }
        if (left1 * left2 < -1e-9){
            P.push_back(lineIntersectSeg(q[i],q[i+1],a,b));
        }
    }
    q.clear();
    q = vector <point> (P);
    if (q.empty()) return;
    if (fabs(q.back().x - q.front().x) > 1e-9 || fabs(q.back().y - q.front().y) > 1e-9) q.push_back(P.front());
}

//input: point p[3]
//output: r, x, y
void circle3pts(){
   double A,X1,X2,Y1,Y2,L1,L2;
   X1 = p[1].x - p[0].x;
   X2 = p[2].x - p[0].x;
   Y1 = p[1].y - p[0].y;
   Y2 = p[2].y - p[0].y;
   L1 = X1*X1 + Y1*Y1;
   L2 = X2*X2 + Y2*Y2;
   A = (X1*Y2 - X2*Y1)/2;
   x = p[0].x + (Y2*L1-Y1*L2)/(4*A);
   y = p[0].y + (X1*L2-X2*L1)/(4*A);
   double dx = p[0].x - x;
   double dy = p[0].y - y;
   r = sqrt(dx*dx+dy*dy);
}

/*
Minimum bounding rectangle for a set of points
s = convex hull with s[N] = s[0]
O(n^2) solution
There exist O(n) solution using rotating calipers
*/
void boundingbox(){
   point u0,u1,D;
   double s0, t0, t1,s1,test;
   minarea = 1e10;
   for (int i = 1; i < s.size(); i++){
      u0.x = s[i].x - s[i-1].x;
      u0.y = s[i].y - s[i-1].y;
      double u0len = sqrt(u0.x*u0.x + u0.y*u0.y);
      u0.x /= u0len;
      u0.y /= u0len;
      u1.x = -u0.y;
      u1.y = u0.x;
      s0 = t0 = s1 = t1 = 0;
      for (int j = 1; j < s.size(); j++){
         D.x = s[j].x-s[0].x;
         D.y = s[j].y-s[0].y;
         test = dot(u0, D);
         if (test < s0) s0 = test; else if (test > s1) s1 = test;
         test = dot(u1, D);
         if (test < t0) t0 = test; else if (test > t1) t1 = test;
      }
      double area = (s1-s0)*(t1-t0);
      if (area < minarea)
         minarea = area;
   }
}

/*
Code for closestpair
O(n log n)
*/
double closestpair(int start, int end, vector <point> *y){
   double ans;
   if (end - start <= 3){
      ans = 1e10;
      for (int i = start; i <= end-1; i++)
         for (int j = i+1; j <= end; j++){
            ans = min(ans,dis(X[i],X[j]));
         }
   }else{
      int mid = (start+end)/2;
      vector <point> YR, YL;
      YR.clear(); YL.clear();
      for (int i = 0; i < y->size(); i++){
         if (y->at(i).x <= X[mid].x)
            YL.push_back(y->at(i));
         else
            YR.push_back(y->at(i));
      }
      double a = closestpair(start,mid,&YL);
      double b = closestpair(mid+1,end,&YR);
      ans = min(a,b);
      sol.clear();
      double midpt = X[mid].x+X[mid+1].x;
      midpt *= 0.5;
      for (int i = 0; i < y->size(); i++)
         if (midpt <= y->at(i).x+ans && y->at(i).x-ans <= midpt)
            sol.push_back(y->at(i));
      for (int i = 1; i < sol.size(); i++)
         ans = min(ans,dis(sol[i],sol[i-1]));
   }
   return ans;
}

bool cmp1(point a, point b){
   return a.x < b.x;
}

bool cmp2(point a, point b){
   return a.y < b.y;
}

//how to call the code.
...
//copy the points into both vector <point> X and Y
sort(X.begin(),X.end(),cmp1);
sort(Y.begin(),Y.end(),cmp2);
double ans = closestpair(0,n-1,&Y);
...


-----------------------------------------------
Number theory
-----------------------------------------------

int bigmod(int b, int p, int m){
    int a;
    if (p == 0)
        return 1;
    if (p%2 == 0){
        a = bigmod(b,p/2,m);
        return (a*a) % m;
    }else{
        return ((b % m) * bigmod(b,p-1,m)) % m;
    }
}
int GCD(int a,int b) {
    while (b > 0) {
        a = a % b;
        a ^= b;
        b ^= a;
        a ^= b;
    }
    return a;
}

typedef struct{
    int d;
    int x;
    int y;
}ee;

/*d = gcd(a,b) = ax + by*/
void extendedeuclid(int a, int b, ee *c){
    ee data1;
    if (b == 0){
        c->d = a; c->x = 1; c->y = 0;
        return;
    }
    extendedeuclid(b, a % b, &data1);
    c->d = data1.d; c->x = data1.y; c->y = data1.x - (a/b)*data1.y;
    return;
}

//Prime sieve
#define MAX 20000000
char p[MAX+5];

void sieve(){
    p[2] = 0;
    for (int i = 2; i <= MAX; i++) if (p[i] == 0){
        plist.push_back(i);
        for (int j = 2*i; j <= MAX; j += i) p[j] = 1;
    }
    return;
}

//Euler totient function
//requires prime list in plist and could be used with sieve
int totient(int a){
   int ans = a;
   for (int i = 0; plist[i]*plist[i] <= a; i++)
      if (a % plist[i] == 0){
         ans /= plist[i];
         ans *= plist[i]-1;
         while (a % plist[i] == 0)
            a /= plist[i];
      }
   if (a != 1){
      ans /= a;
      ans *= a-1;
   }
   return ans;
}

vector <pii> list;
void factor(int a){
   list.clear();
   for (int i=0; plist[i]*plist[i]<=a; i++){
      int cnt = 0;
      while (a % plist[i] == 0){
         a /= plist[i];
         cnt++;
      }
      if (cnt > 0){
         list.push_back(pii(plist[i],cnt));
      }
   }
   if (a != 1){
      list.push_back(pii(a, 1));
   }
}

//returns the multiplicative order of a modulo n
//ie the smallest h such that a^h = 1 (mod n)
int ord(int a, int n){
   if (a <= 1 || n <= 1 || GCD(a,n) != 1)
      return -1;
   int phi = totient(n);
   factor(phi);
   int ans = phi;
   for (int i = 0; i < list.size(); i++){
      for (int j = 0; j < list[i].second; j++)
         ans /= list[i].first;
      int g1 = bigmod(a,ans,n);
      while (g1 != 1){
         g1 = bigmod(g1,list[i].first,n);
         ans = ans * list[i].first;
      }
   }
   return ans;
}

-------------------------------------------------------------
Numerical Methods
-------------------------------------------------------------

#define EPS 1e-9
#define INF 1e9

double A[MAXN][MAXM], X[MAXN];
int basis[MAXN], out[MAXM];
int x[MAXM], y[MAXN];

void pivot (int n, int m, int a, int b) {
    int i, j;
    for (i = 0; i <= n; i++)
        if (i != a)
            for (j = 0; j <= m; j++)
                if (j != b)
                    A[i][j] -= A[a][j] * A[i][b] / A[a][b];
    for (j = 0; j <= m; j++)
        if (j != b) A[a][j] /= A[a][b];
    for (i = 0; i <= n; i++)
        if (i != a) A[i][b] = -A[i][b] / A[a][b];
    A[a][b] = 1 / A[a][b];
    i = basis[a]; basis[a] = out[b]; out[b] = i;
}

/*--------------------------------------------------------------
Input: n = last row in A, m = last col in A

A[0] = Objective function to minimise
The remaining entries are <= with the last col being b ie
Ax <= b

Output: INF if data is unbounded, -INF if data is inconsistent
----------------------------------------------------------------*/

double simplex (int n, int m) {
    int i, j, ii, jj;
    memset(X, 0, sizeof X);
    memset(basis, 0, sizeof basis);
    for (i = 0; i <= n; i++)
        basis[i] = -i;
    for (j = 0; j <= m; j++)
        out[j] = j;
    for (;;) {
        for (i = ii = 1; i <= n; i++)
            if (A[i][m] < A[ii][m] || (A[i][m] == A[ii][m] && basis[i] < basis[ii]))
                ii= i;
        if (A[ii][m] >= -EPS) break;
        for (j = jj = 0; j < m; j++)
            if (A[ii][j] < A[ii][jj] - EPS || (A[ii][j] < A[ii][jj] - EPS && out[i] < out[j]))
                jj = j;
        if (A[ii][jj] >= -EPS) return -INF;
        pivot (n, m, ii, jj);
    }
    for (;;) {
        for (j = jj = 0; j < m; j++)
            if (A[0][j] < A[0][jj] || (A[0][j] == A[0][jj] && out[j] < out[jj]))
                jj = j;
        if (A[0][jj] > -EPS) break;
        for (i = 1, ii = 0; i <= n; i++)
            if ((A[i][jj] > EPS) &&
            (!ii || (A[i][m] / A[i][jj] < A[ii][m] / A[ii][jj] - EPS) ||
            ((A[i][m] / A[i][jj] < A[ii][m] / A[ii][jj] + EPS) && (basis[i] < basis[ii]))))
                ii = i;
        if (A[ii][jj] <= EPS) return INF;
        pivot (n, m, ii, jj);
    }
    
    for (j = 0; j < m; j++)
        X[j] = 0;
    for (i = 1; i <= n; i++)
        if (basis[i] >= 0)
            X[basis[i]] = A[i][m];
    return A[0][m]; 
}

-----------------------------------------------------------------
polynomial class
-----------------------------------------------------------------

struct poly{
    vector <int> a;
    poly(){ a.clear(); }
    poly operator + (const poly &A){
        poly p = poly();
        int m = max(A.a.size(), a.size());
        for (int i = 0; i < m; i++){
            int h = (i < a.size() ? a[i] : 0);
            int k = (i < A.a.size() ? A.a[i] : 0);
            p.a.push_back(h+k);
        }
        return p;
    }
    poly operator - (const poly &A){
        poly p = poly();
        int m = max(A.a.size(), a.size());
        for (int i = 0; i < m; i++){
            int h = (i < a.size() ? a[i] : 0);
            int k = (i < A.a.size() ? A.a[i] : 0);
            p.a.push_back(h-k);
        }
        return p;
    }
    poly operator * (const poly & A){
        poly p = poly();
        int n = a.size(), m = A.a.size();
        p.a.resize(n*m);
        for (int i = 0; i < m; i++){
            for (int j = 0; j < n; j++){
                p.a[i+j] += A.a[i] * a[j];
            }
        }
        while (p.a.size() != 0 && p.a.back() == 0) p.a.pop_back();
        return p;
    }
    poly operator % (const poly &A){
        poly p = poly();
        for (int i = 0; i < a.size(); i++) p.a.push_back(a[i]);
        int n = A.a.size();
        for (int i = p.a.size()-1; i >= 0 && i+1 - n >= 0; i--){
            int rem = p.a[i] / A.a.back();
            for (int j = 0; j < n; j++){
                p.a[i+j-A.a.size()+1] -= rem*A.a[j];
            }
        }
        while (p.a.size() != 0 && p.a.back() == 0) p.a.pop_back();
        return p;
    }
    void print(){
        while (a.size() != 0 && a.back() == 0) a.pop_back();
        if (a.size() == 0){
            printf("0\n"); return;
        }
        for (int i = a.size()-1; i >= 0; i--){
            if (a[i] == 0) continue;
            if (i != a.size()-1 && a[i] > 0) printf("+");
            if (i == 0){
                printf("%d",a[i]);
            }else if (i == 1){
                if (a[i] == -1) printf("-");
                if (abs(a[i]) != 1) printf("%d",a[i]);
                printf("x");
            }else{
                if (a[i] == -1) printf("-");
                if (abs(a[i]) != 1) printf("%d",a[i]);
                printf("x^%d",i);
            }
        }
        printf("\n");
    }
};


-------------------------------------------------------------
Dijkstra
-------------------------------------------------------------

//PQ
int dis[90010];
int n;

int dijkstra(int start){
   set < pair<int, int> > pq;
   dis[start] = 0;
   pq.insert(pii(dis[start], start));
   while (!pq.empty()){
      int node = pq.begin()->second;
      if (node == end) return dis[node];
      pq.erase(pq.begin());
      for (int i = 0; i < adj[node].size(); i++){
         int u = adj[node][i].first;
         int weight = adj[node][i].second;
         if (weight + dis[node] < dis[u]){
            pq.erase(pii(dis[u],u));
            dis[u] = weight + dis[node];
            pq.insert(pii(dis[u],u));
         }
      }
   }
    return -1;
}

//..
adj[u].push_back(pii(v,cost));


-------------------------------------------------------------
Karp Minimum Mean Cycle Weight
-------------------------------------------------------------

      //add s to the graph and for all vertex x in the graph, add edge (s,x) with weight 0
      memset(dis,0x7f,sizeof(dis));
      dis[s][0] = 0;
      for (k = 1; k <= n; k++){
         for (int u = 0; u < n; u++){
            for (vector <pii>::iterator it = indeg[u].begin(); it != indeg[u].end(); it++){
               if (dis[it->first][k-1] != 0x7f7f7f7f && dis[it->first][k-1] + it->second < dis[u][k])
                  dis[u][k] = dis[it->first][k-1] + it->second;
            }
         }
      }
      bool flag = true;
      for (i = 0; i < n; i++) if (dis[i][n] != 0x7f7f7f7f){
         flag = false; break;
      }
      if (flag == true){
         //graph is acyclic
      }
      int mincycleweight = 0x7f7f7f7f;
      for (u = 0; u < n; u++){
         int cycleweight = 0;
         for (k = 0; k < n; k++) cycleweight = max(cycleweight, (dis[u][n]-dis[u][k])/(n-k));
         mincycleweight = min(mincycleweight, cycleweight);
      }


-------------------------------------------------------------
Articulation points
-------------------------------------------------------------

typedef struct{
    int cnt;
    int vertex[105];
}adj;

adj list[105];
int visit[105],d[105],low[105],pred[105],pts[105];

    int min(int a, int b){
        if (a < b)
            return a;
        return b;
    }

void ArtPt(int u){
    int ii,v,cnt=0;
    visit[u] = 1;
    low[u] = d[u] = ++time;
    for (ii=0;ii<list[u].cnt;ii++){
        v = list[u].vertex[ii];
        if (visit[v] == 0){
            pred[v] = u;
            ArtPt(v);
            low[u] = min(low[u],low[v]);
            cnt++;
            if (pred[u] == NULL){ /*Watch this part*/
                if (cnt == 2)
                    pts[u] = 1;
            }else if (low[v] >= d[u]){
                pts[u] = 1;
            }
        }else if (v != pred[u]){
            low[u] = min(low[u],d[v]);
        }
    }
    return;
}

-------------------------------------------------------------
Finding bridges
-------------------------------------------------------------

void Bridgefinding(int u){
    int ii,v;
    visit[u] = 1;
    low[u] = d[u] = ++time;
    for (ii=0;ii<list[u].cnt;ii++){
        v = list[u].vertex[ii];
        if (visit[v] == 0){
            pred[v] = u;
            Bridgefinding(v);
            low[u] = min(low[u],low[v]);
            if (low[v] > d[u]){
                bridge[cnt].small = min(u,v);
                bridge[cnt].big = max(u,v);
                cnt++;
            }
        }else if (v != pred[u]){
            low[u] = min(low[u],d[v]);
        }
    }
    return;
}

--------------------------------------------------------------
Network Flow
--------------------------------------------------------------

#define MAXN 105
#define MAXE 100000

int nsize, source, sink, nedge;
int head[MAXN], point[MAXE], cap[MAXE], flow[MAXE], next[MAXE];
int dis[MAXN], Q[MAXN], work[MAXN];

void init(int _size, int s, int t){
   nsize = _size; source = s; sink = t;
   memset(head,-1,sizeof(head));
   nedge = 0;
}

void add(int u, int v, int c1, int c2){
   point[nedge] = v, cap[nedge] = c1, flow[nedge] = 0, next[nedge] = head[u], head[u] = (nedge++);
   point[nedge] = u, cap[nedge] = c2, flow[nedge] = 0, next[nedge] = head[v], head[v] = (nedge++);
}

bool dinic_bfs(){
   memset(dis,-1,sizeof(dis));
   dis[source] = 0;
   queue <int> q;
   q.push(source);
   while (!q.empty()){
      int u = q.front(); q.pop();
      for (int i = head[u]; i >= 0; i = next[i])
         if (flow[i] < cap[i] && dis[point[i]] < 0){
            dis[point[i]] = dis[u]+1;
            q.push(point[i]);
         }
   }
   return (dis[sink]>=0);
}

int dinic_dfs(int node, int limit){
   if (node == sink) return limit;
   for (int &i = work[node]; i >= 0; i = next[i]){
      int v = point[i], tmp;
      if (flow[i] < cap[i] && dis[v] == dis[node]+1 && (tmp = dinic_dfs(v, min(limit,cap[i]-flow[i])))>0){
         flow[i] += tmp;
         flow[i^1] -= tmp;
         return tmp;
      }
   }
   return 0;
}

int dinic_flow(){
   int ans = 0;
   while (dinic_bfs()){
      for (int i = 0; i < nsize; i++) work[i] = head[i];
      while (1){
         int delta = dinic_dfs(source, INF);
         if (delta == 0) break;
         ans += delta;
      }
   }
   return ans;
}

///////////////////
//BFS VARIANT
///////////////////   
int totalflow;
int cap[205][205],pre[205],flow[205];
vector <int> adj[205];

int find(int source, int sink){
   memset(pre,-1,sizeof(pre));
   memset(flow,0x0,sizeof(flow));
   queue <int> q;
   q.push(source);
   flow[source] = MAXINT;
   while (!q.empty()){
      int u = q.front();
      q.pop();
      if (u == sink)
         break;
      for (int i = 0; i < adj[u].size(); i++){
         int next = adj[u][i];
         if (flow[next] || cap[u][next] <= 0) 
            continue;
         flow[next] = min(flow[u], cap[u][next]);
         pre[next] = u;
         q.push(next);
      }
   }
   int fl = flow[sink];
   while (pre[sink]!=-1){
      cap[sink][pre[sink]] += fl;
      cap[pre[sink]][sink] -= fl;
      sink = pre[sink];
   }
   return fl;
}

void networkflow(int source, int sink){
   int pathcap = 0;
   totalflow = 0;
   do {
      pathcap = find(source,sink);
      totalflow += pathcap;
   } while (pathcap);
   //totalflow = maxflow_val
}

//suppose a is adjacent to b with edge weight N on a directed graph
//then cap[a][b] = N;
//adj[a].push_back(b);
//adj[b].push_back(a);
//for residue flow

--------------------------------------------------------------
All Pairs Max Flow - Stoer Wagner
--------------------------------------------------------------
int adj[155][155];
int w[155], V[155];
bool A[155];

int mincut(){
    int cut = 1000000000;
    int N = n;
    FOR(i,1,n) V[i] = i;
    while (N > 1){
        //a = v[1]
        memset(A, 0, sizeof A);
        A[1] = true;
        FOR(i,2,N) w[V[i]] = adj[V[1]][V[i]];
        int prev = V[1];
        FOR(i,2,N){
            //find most tightly connected
            int loc = -1;
            FOR(j,1,N){
                if (A[j] == false && (loc < 0 || w[V[j]] > w[V[loc]]))
                    loc = j;
            }
            A[loc] = true;
            if (i == N){
                if (cut > w[V[loc]]) cut = w[V[loc]];
                //merge prev & loc
                FOR(j,1,N){
                    adj[V[j]][prev] += adj[V[loc]][V[j]];
                    adj[prev][V[j]] += adj[V[loc]][V[j]];
                }
                V[loc] = V[N];
                break;
            }
            prev = V[loc];
            FOR(j,1,N){
                if (A[j] == true) continue;
                w[V[j]] += adj[V[loc]][V[j]];
            }
        }
        N--;
    }
    return cut;
}


--------------------------------------------------------------
Minimum cost maximum flow (mcmf)
--------------------------------------------------------------

struct Edge{
int x,y;
long long cap, cost;
};
vector <Edge> E;
vector <int> adj[MAXN];
int N, prev[MAXN];
long long dis[MAXN], phi[MAXN],flowVal, flowCost;

void init(){
   E.clear();
   for (int i = 0; i < N; i++) adj[i].clear();
}

void add(int x, int y, long long cap, long long cost){
   Edge e1 = {x,y,cap,cost}, e2 = {y,x,0,-cost};
   adj[e1.x].push_back(E.size()); E.push_back(e1);
   adj[e2.x].push_back(E.size()); E.push_back(e2);
}

void mcmf(int s, int t){
   int x;
   flowVal = flowCost = 0;
   memset(phi,0,sizeof(phi));
   while (true){
      for (x = 0; x < N; x++) prev[x] = -1;
      for (x = 0; x < N; x++) dis[x] = INF;
      dis[s] = prev[s] = 0;
      set < pair<long long, int> > Q;
      Q.insert(make_pair(dis[s], s));
      while (!Q.empty()){
         x = Q.begin()->second; Q.erase(Q.begin());
         for (vector <int>::iterator it=adj[x].begin(); it!=adj[x].end(); it++){
            const Edge &e = E[*it];
            if (e.cap <= 0) continue;
            long long cc = e.cost + phi[x] - phi[e.y];
            if (dis[x] + cc < dis[e.y]){
               Q.erase(make_pair(dis[e.y],e.y));
               dis[e.y] = dis[x] + cc;
               prev[e.y] = *it;
               Q.insert(make_pair(dis[e.y],e.y));
            }
         }
      }
      if (prev[t] == -1) break;
      long long z = INF;
      for (x = t; x != s; x = E[prev[x]].x) z = min(z,E[prev[x]].cap);
      for (x = t; x != s; x = E[prev[x]].x){
         E[prev[x]].cap -= z;
         E[prev[x]^1].cap += z;
      }
      flowVal += z;
      flowCost += z * (dis[t]-phi[s]+phi[t]);
      for (x = 0; x < N; x++) if (prev[x] != -1) phi[x] += dis[x];
   }
}

--------------------------------------------------------------
Bipartite Matching O(mn^2)
--------------------------------------------------------------

int adj[M][N];
int seen[N];
int matchL[M], matchR[N];
int n, m;

bool bpm(int u){
   for (int v = 0; v < n; v++){
      if (adj[u][v]){
         if (seen[v])
            continue;
         seen[v] = true;
         if (matchR[v] < 0 || bpm(matchR[v])){
            matchL[u] = v;
            matchR[v] = u;
            return true;
         }
      }
   }
   return false;
}

main(){
   memset(matchL, -1, sizeof matchL);
   memset(matchR, -1, sizeof matchR);
   for (i = 0; i < m; i++){
      memset(seen,0,sizeof seen);
      if (bpm(i)) 
         matching++;
   }
   //matching = number of matching
   //matchL[m]: for each pigeon, a hole or -1
   //matchR[n]: for each hole, a pigeon or -1
}

--------------------------------------------------------------
Blossom shrinking - Maximum matching on general graph
--------------------------------------------------------------
#define FOR(i,a,b) for (int i = a; i <= (int) b; i++)
int n;
int match[105], visited[105];

bool dfs(int node, vector <vector <int> > &adj, vector <int> &blossom){
    visited[node] = 0;
    FOR(i, 0, n-1) if (adj[node][i]){
        if (visited[i] == -1){
            visited[i] = 1;
            if (match[i] == -1 || dfs(match[i], adj, blossom)){
                match[node] = i;
                match[i] = node;
                return true;
            }
        }
        if (visited[i] == 0 || !blossom.empty()){
            blossom.push_back(i);
            blossom.push_back(node);
            if (node == blossom[0]){
                match[node] = -1; return true;
            }
            return false;
        }
    }
    return false;
}

bool augment(vector <vector <int> > &adj){
    FOR(m, 0, n-1) if (match[m] == -1){
        vector <int> blossom;
        memset(visited, -1, sizeof visited);
        if (!dfs(m, adj, blossom)) continue;
        if (blossom.empty()) return true;
        int base = blossom[0];
        vector <vector <int> > newadj = adj;
        FOR(i, 1, blossom.size()-2) FOR(j, 0, n-1) newadj[base][j] = newadj[j][base] |= adj[blossom[i]][j];
        FOR(i, 1, blossom.size()-2) FOR(j, 0, n-1) newadj[blossom[i]][j] = newadj[j][blossom[i]] = 0;
        newadj[base][base] = 0;
        if (!augment(newadj)) return false;
        int k = match[base];
        if (k != -1) FOR(i, 0, blossom.size()-1) if (adj[blossom[i]][k]){
            match[blossom[i]] = k;
            match[k] = blossom[i];
            if (i & 1){
                for (int j = i+1; j < blossom.size(); j += 2){
                    match[blossom[j]] = blossom[j+1];
                    match[blossom[j+1]] = blossom[j];
                }
            }else{
                for (int j = 0; j < i; j += 2){
                    match[blossom[j]] = blossom[j+1];
                    match[blossom[j+1]] = blossom[j];
                }
            }
            break;
        }
        return true;
    }
    return false;
}

int edmonds(vector <vector <int> > &adj){
    int ans = 0;
    memset(match, -1, sizeof match);
    while (augment(adj)) ans++;
    return ans;
}

--------------------------------------------------------------
Strongly Connected Components (Tarjan's Algorithm)
--------------------------------------------------------------

vector <int> list[100100];
vector <int> indegree[100100];
int index2[100100],lowlink[100100],curindex;
vector <int> s;

    int min(int a, int b){
        if (a < b)
            return a;
        return b;
    }

int isnodeinstack(int node){
    for (int i = 0; i < s.size(); i++)
        if (s[i] == node)
            return 1;
    return 0;
}

void DFStarjan(int node){
    int i, u;
    index2[node] = curindex;
    lowlink[node] = curindex;
    curindex++;
    s.push_back(node);
    for (i = 0; i < list[node].size(); i++){
        u = list[node][i];
        if (index2[u]==0){
            DFStarjan(u);
            lowlink[node] = min(lowlink[node],lowlink[u]);
        }else if (isnodeinstack(u)){
            lowlink[node] = min(lowlink[node], lowlink[u]);
        }
    }
    if (lowlink[node] == index2[node]){
        printf("SCC: ");
        do {
            u = s.back();
            s.pop_back();
            printf("%d ",u);
        }while (u != node);
        printf("\n");
    }
}

//how to call it
for (i = 1; i <= V; i++){
    if (index2[i] == 0){
        s.clear();
        curindex = 1;
        DFStarjan(i);
    }
}

-------------------------------------------------------------
Prim O(V^2)
-------------------------------------------------------------

    int graph[105][105];
    int dis[105];
    char intree[105];

    int prim(int nvert){
        int ii,jj,mindis,treesize = 1, start = 0, total = 0;
        for (ii=0; ii < nvert; ii++){
            dis[ii] = 1000000000;
            intree[ii] = 0;
        }
        intree[0] = 1;
        for (ii=1; ii < nvert; ii++)
            dis[ii] = graph[0][ii];
        while (treesize < nvert){
            mindis = 1000000000;
            for (jj=0; jj<nvert; jj++){
                if (dis[jj] < mindis && intree[jj]==0){
                    mindis = dis[jj];
                    ii = jj;
                }
            }
            treesize++;
            total += mindis;
            intree[ii] = 1;
            for (jj = 0; jj < nvert; jj++)
                if (dis[jj] > graph[ii][jj] && intree[jj] == 0)
                    dis[jj] = graph[ii][jj];
        }
        return total;
    }

-------------------------------------------------------------
Prim O((V+E)lg V)
-------------------------------------------------------------

    //PQ
    int dis[90010];
    int numnode;

    struct my_less{
        bool operator () (int u, int v){
            if (dis[u] == dis[v])
                return (u < v);
            else
                return (dis[u] < dis[v]);
        }
    };

set <int, my_less> pq;

int prim(){
    int minloc, mindis,u,weight,total;
    total = 0;
    for (int ii = 0; ii < numnode; ii++){
        if (ii == start)
            continue;
        dis[ii] = MAXINT;
        pq.insert(ii);
    }
    dis[start] = 0;
    pq.insert(start);
    while (!pq.empty()){
        minloc = *pq.begin();
        pq.erase(minloc);
        total += dis[minloc];
        for (int ii = 0; ii < adj[minloc].size();ii++){
            u = adj[minloc][ii].first;
            weight = adj[minloc][ii].second;
            if (weight < dis[u]){
                pq.erase(u);
                dis[u] = weight;
                pq.insert(u);
            }
        }
    }
    return total;
}

-------------------------------------------------------------
Kruskal's algorithm
-------------------------------------------------------------

int parent[50100];
int rank[50100];

int findroot(int a){
   int root = a;
   while (root != parent[root])
      root = parent[root];
   int j = parent[a];
   while (j!=root){
      parent[a] = root; a = j; j = parent[a];
   }
   return root;
}

void unionset(int a, int b){
   a = findroot(a);
   b = findroot(b);
   if (a == b)
      return;
   if (rank[a] < rank[b])
      parent[a] = b;
   else if (rank[a] > rank[b])
      parent[b] = a;
   else{
      parent[a] = b;
      rank[b] ++;
   }
}

int kruskal(){
   int i,total = 0;
   sort(e.begin(), e.end());
   for (i = 0; i<nedges; i++){
      rank[i] = 0;
      parent[i] = i;
   }
   for (i = 0;i < nedges; i++){
      if (findroot(e[i].v1) != findroot(e[i].v2)){
         total += e[i].weight;
         unionset(e[i].v1, e[i].v2);
      }
   }
   return total
}

-------------------------------------------------------------
Bellman-Ford
-------------------------------------------------------------

/*Bellman-Ford*/
dis[source] = 0;
    for (i=0;i<n;i++)
        if (i!=source)
            dis[i] = 1000000000;
    for (i=0;i<n;i++)
        for (j=0;j<nedge;j++){
            u = edges[j].src;
            v = edges[j].dest;
            tmp = edges[j].weight + dis[u];
            if (dis[v] > tmp)
                dis[v] = tmp;
        }
/*Detect negative cycle*/
flag = 0;
for (i=0;i<nedge;i++){
    u = edges[i].src;
    v = edges[i].dest;
    if (dis[v] > dis[u] + edges[i].weight){
        flag = 1;
        break;
    }
}

-------------------------------------------------------------
Floyd Warshall
-------------------------------------------------------------

/*Initialisation*/
for (i=0;i<n;i++)
    for (j=0;j<n;j++){
        if (i==j)
            w[i][j] = 0;
        else
            w[i][j] = MAXINT;
    }

/*DP part*/

for (k=0;k<n;k++)
    for (i=0;i<n;i++)
        for (j=0;j<n;j++){
            if (w[i][j] > w[i][k] + w[k][j])
                w[i][j] = w[i][k] + w[k][j];
        }

----------------------------------------------------
Fleury's Algorithm
----------------------------------------------------

int conn[i][j] = number of edges between i and j
int plen;

void find_path(int loc){
    int lv;
    for (lv = 0; lv < nconn; lv++)
        if (conn[loc][lv]){
            /* delete edge */
            conn[loc][lv]--;
            conn[lv][loc]--;
            deg[lv]--;
            deg[loc]--;

            /* find path from new location */
            find_path(lv);
        }

    /* add this node to the `end' of the path */
    path[plen++] = loc;
}

--------------------------------
Topological Sorting
--------------------------------

list <int> L;

void TopSort(){
    memset(visit,0,sizeof visit)
    for (i=0;i<n;i++)
        if (visit[i] == 0)
            TopVisit(i);
}

void TopVisit(int node){
   visit[node] = 1;
   for (int j = 0; j < adj[node].size(); j++){
      int u = adj[node][j];
      if (visit[u] == 0) TopVisit(u);
   }
   L.push_front(node);
}

----------------------------------------------------
Stable Marriage Problem
----------------------------------------------------

/*
rankM[man][x] = xth most preferred woman for man
rankW[woman][man] = the ranking for man according to woman
*/

void stableMarriage(){
   memset(matchM,-1,sizeof(matchM));
   memset(matchW,-1,sizeof(matchW));
   memset(p,0,sizeof(p));
   for (int i = 0; i < m; i++){ //m = number of man
      int man = i;
      while (man >= 0){
         int woman;
         while (true){
            woman = rankM[man][p[man]++];
            if (matchW[woman] < 0 || rankW[woman][man] > rankW[woman][matchW[woman]]) break;
         }
         int hubby = matchW[woman];
         matchM[man] = woman;
         matchW[woman] = man;
         man = hubby;
      }
   }
}

--------------------------------
KMP string search
--------------------------------

int phi[10000];
char *KMP(char *text, char *pat){
   int i,k,n=strlen(pat);
   phi[0] = phi[1] = k = 0;
   for (i = 1; i < n ; i++){
      while (k > 0 && pat[k] != pat[i])
         k = phi[k];
      phi[i+1] = (pat[k]==pat[i] ? ++k : 0);
   }
   for (i = k = 0; text[i] && k < n; i++){
      while (k > 0 && pat[k] != text[i])
         k = phi[k];
      if (pat[k] == text[i]) k++;
   }
   return (k == n ? text+i-n: NULL);
}

--------------------------------
Suffix Array
--------------------------------

long long key[200000];
int sa[200000], rank[200000], *lcp = (int *)key;
struct Cmp { bool operator ()(int i, int j) const {return key[i]<key[j]; } };
void build(char *text, int n){
   for (int i = 0; i < n; i++){ sa[i] = i; key[i] = text[i]; }
   sort(sa, sa+n, Cmp());
   for (int k = 1; ; k *= 2){
      for (int i = 0; i < n; i++)
         rank[sa[i]] = ( i>0 && key[sa[i-1]] == key[sa[i]] ? rank[sa[i-1]] : i);
      if (k >= n)
         break;
      for (int i = 0; i < n; i++)
         key[i] = rank[i] * (n+1LL) + (i+k < n ? rank[i+k]+1 : 0);
      sort(sa,sa+n,Cmp());
   }
   for (int i = 0, k = 0; i < n; i++){
      if (k > 0) k--;
      if (rank[i] == n-1){ lcp[n-1] = -1; k = 0; continue;}
      int j = sa[rank[i]+1];
      while (text[i+k] == text[j+k]) k++;
      lcp[rank[i]] = k;
   }
}

//Binary search O(M + log N)
void query(char *Q, int Qlen){
   int L,R,M,l,r,m,mlr;
   L = 0; R = len*2-1;
   l = 0; r = 0;
   while (true){
      mlr = min(l,r);
      M = (L+R)/2;
      for (m = mlr;; m++)
         if (Q[m] != str[sa[M]+m])
            break;
      if (m == Qlen)
         break;
      if (M == L || M == R){
         //Q not found
         return;
      }
      if (str[sa[M]+m] > Q[m]){
         R = M;
         r = m;
      }else{
         L = M;
         l = m;
      }
   }
   printf("substring found at %d\n",M);
}

-------------------------------------------------------------
Trie
-------------------------------------------------------------
struct nodeT{
    int end;
    struct nodeT *next[26];
};

struct nodeT root;
struct nodeT nodes[10000];

void addWord(struct nodeT *p, char *word){
    int ii, id,len = strlen(word);
    struct nodeT *v;
    for (ii=0;ii<len;ii++){
        id = word[ii]-'a';
        if (p->next[id]==NULL){
            memset(&nodes[np],0x0,sizeof (struct nodeT));
            p->next[id] = &nodes[np];
            np++;
        }
        p = p->next[id];
        if (ii == len - 1)
            p->end = 1;
    }
    return;
}

void printtrie(struct nodeT *p, int depth){
    int ii,iszero;
    path[depth] = 0;
    if (p->end == 1){
        printf("%s\n",path);
    }
    for (ii=0,iszero=1;ii<26;ii++){
        if (p->next[ii] != 0){
            path[depth] = ii+'a';
            printtrie(p->next[ii],depth+1);
            iszero = 0;
        }
    }
    if (iszero != 0 && p->end == 0)
        printf("%s\n",path);
    return;
}

...
memset(&root,0x0,sizeof(struct nodeT));
...
AddWord(&root,word);
...
printtrie(&root,0);

-------------------------------------------------------------
Fenwick Tree
-------------------------------------------------------------
#define MAXN 10000005
int tree[MAXN];

void create(){
   memset(tree,0,sizeof(tree));
}

int query(int a, int b){
   if (a == 0){
      int sum = 0;
      for (; b >= 0; b = (b & (b+1))-1)
         sum += tree[b];
      return sum;
   }else
      return query(0,b)-query(0,a-1);
}

void increase(int k, int inc){
   for (; k < MAXN; k |= k+1) tree[k] += inc;
}

int kthelement(int k){ //ie kth smallest element (starting from 1)
   int left = -1;
   int right = MAXN - 1;
   while (right - left > 1){
      int mid = (left+right)/2;
      if (query(0,mid) < k)
         left = mid;
      else
         right = mid;
   }
   return right;
}

-------------------------------------------------------
RMQ - Range MAX Query
For Range MIN Query, change all the >= to <=
-------------------------------------------------------
int RMQ[100005][18]; //RMQ[MAXN][LOGMAXN+1]

int log_2(int a){
   return (int)floor(log(a)/log(2));
}

int queryRMQ(int a, int b){
   int k = log_2(b-a+1);
   if (A[RMQ[a][k]] >= A[RMQ[b-(1<<k)+1][k]])
      return RMQ[a][k];
   else
      return RMQ[b-(1<<k)+1][k];
}

void processRMQ(){
   for (int i = 0; i < n; i++) RMQ[i][0] = i;
   for (int j = 1; (1<<j) <= n; j++)
      for (int i = 0; i + (1<<j) - 1 < n; i++)
         if (A[RMQ[i][j-1]] >= A[RMQ[i+(1<<(j-1))][j-1]])
            RMQ[i][j] = RMQ[i][j-1];
         else
            RMQ[i][j] = RMQ[i+(1<< (j-1))][j-1];
}

-------------------------------------------------------
kd tree - k=2 here
-------------------------------------------------------

#define pii pair<double, int>

struct point{
double coord[2];
int id;
};

struct node{
int id;
double coord[2];
struct node *left;
struct node *right;
};

int axis;
vector <point> pt;
vector <pii> ans;

bool cmp(point a, point b){
   return (a.coord[axis] < b.coord[axis]);
}

node* build_kdtree(vector <point> const& pt, int depth){
   vector <point> a = pt;
   if (a.size() == 0) return NULL;
   axis = depth % 2;
   sort(a.begin(),a.end(),cmp);
   int median = a.size()/2;
   node *tmp = &nodeT[n];
   n++;
   tmp->id = a[median].id;
   tmp->coord[0] = a[median].coord[0];
   tmp->coord[1] = a[median].coord[1];
   vector <point> b;
   vector <point> c;
   for (int i = 0; i < median; i++) b.push_back(a[i]);
   for (int i = median+1; i < a.size(); i++) c.push_back(a[i]);
   tmp->left = build_kdtree(b,depth+1);
   tmp->right = build_kdtree(c,depth+1);
   return tmp;
}

double distance(double x1, double y1, double x2, double y2){
   double dx = x1-x2;
   double dy = y1-y2;
   return dx*dx+dy*dy;
}

void check_nearest(int k, node *curnode, double x, double y){
   double d = distance(curnode->coord[0], curnode->coord[1], x,y);
   if (ans.size() < k || d < ans[ans.size()-1].first){
      if (ans.size() >= k) ans.pop_back();
      ans.push_back(pii(d,curnode->id));
      sort(ans.begin(),ans.end());
   }
}

void knearest(int k,node *curnode, double x, double y, int depth){
   node *nearer, *further;
   double target[2];
   target[0] = x; target[1] = y;
   int axis = depth % 2;
   if (curnode->left == NULL && curnode->right == NULL){
      check_nearest(k,curnode,x,y);
      return;
   }
   if (curnode->right == NULL || (curnode->left != NULL &&target[axis] <=curnode->coord[axis])){
      nearer = curnode->left;
      further = curnode->right;
   }else{
      nearer = curnode->right;
      further = curnode->left;
   }
   knearest(k,nearer, x,y, depth+1);
   if (further != NULL){
      double tmp = target[axis]-curnode->coord[axis];
      if (ans.size() < k || tmp*tmp < ans[ans.size()-1].first)
         knearest(k,further, x,y,depth+1);
   }
   check_nearest(k,curnode, x,y);
}


///MISC FUNCTION
int getnode(string s){
   map<string,int>::iterator iter = mymap.find(s);
   if (iter != mymap.end())
      return iter->second;
   mymap.insert(make_pair(s,cnt));
   cnt++;
   return cnt-1;
}