refactor NetworkFlow, add a faster convolution

zimpha · zimpha · commit 2578f8dfdb01 · 2017-09-26T10:17:18.000+08:00
diff --git a/graph-utility/NetworkFlow.cc b/graph-utility/NetworkFlow.cc
@@ -2,79 +2,102 @@
 #include <algorithm>
 
 namespace NetFlow {
-  const int N=100000,MAXM=100000,inf=1e9;
-  struct Edge {
-    int v,c,f,nx;//c:capcity, f:flow
-    Edge() {}
-    Edge(int v,int c,int f,int nx):v(v),c(c),f(f),nx(nx) {}
-  } E[MAXM];
-  int G[N],cur[N],pre[N],dis[N],gap[N],n,sz;
+  using flow_t = int;
+  const int N = 1e5 + 10, M = 1e5 + 10;
+  const flow_t inf = 1e9;
+  struct edge_t {
+    int to, nx;
+    flow_t cap, flow;
+    edge_t() {}
+    edge_t(int to, int nx, flow_t cap): 
+      to(to), nx(nx), cap(cap), flow(0) {}
+  } edges[M];
+  int G[N], cur[N], pre[N], gap[N], n, sz;
+  flow_t dis[N];
   void init(int _n) {
-    n=_n,sz=0; memset(G,-1,sizeof(G[0])*n);
+    n = _n, sz = 0;
+    memset(G, -1, sizeof(*G) * n);
   }
-  void link(int u,int v,int c) {
-    E[sz]=Edge(v,c,0,G[u]); G[u]=sz++;
-    E[sz]=Edge(u,0,0,G[v]); G[v]=sz++;
+  void link(int u, int v, flow_t c) {
+    edges[sz] = edge_t(v, G[u], c); G[u] = sz++;
+    edges[sz] = edge_t(u, G[v], 0); G[v] = sz++;
   }
-  int ISAP(int S,int T) {//S -> T
-    int maxflow=0,aug=inf,flag=false,u,v;
-    for (int i=0;i<n;++i)cur[i]=G[i],gap[i]=dis[i]=0;
-    for (gap[S]=n,u=pre[S]=S;dis[S]<n;flag=false) {
-      for (int &it=cur[u];~it;it=E[it].nx) {
-        if (E[it].c>E[it].f&&dis[u]==dis[v=E[it].v]+1) {
-          if (aug>E[it].c-E[it].f) aug=E[it].c-E[it].f;
-          pre[v]=u,u=v; flag=true;
-          if (u==T) {
-            for (maxflow+=aug;u!=S;) {
-              E[cur[u=pre[u]]].f+=aug;
-              E[cur[u]^1].f-=aug;
+  flow_t ISAP(int S, int T) {//S -> T
+    flow_t maxflow = 0, aug = inf;
+    memcpy(cur, G, sizeof(*G) * n);
+    memset(gap, 0, sizeof(*gap) * n);
+    memset(dis, 0, sizeof(*dis) * n);
+    gap[S] = n, pre[S] = S;
+    for (int u = S, flag = 0; dis[S] < n; flag = 0) {
+      for (int &it = cur[u]; ~it; it = edges[it].nx) {
+        int v = edges[it].to;
+        if (edges[it].cap > edges[it].flow && dis[u] == dis[v] + 1) {
+          aug = std::min(aug, edges[it].cap - edges[it].flow);
+          pre[v] = u, u = v, flag = true;
+          if (u == T) {
+            for (maxflow += aug; u != S; ) {
+              u = pre[u];
+              edges[cur[u]].flow += aug;
+              edges[cur[u] ^ 1].flow -= aug;
             }
-            aug=inf;
+            aug = inf;
           }
           break;
         }
       }
       if (flag) continue;
-      int mx=n;
-      for (int it=G[u];~it;it=E[it].nx) {
-        if (E[it].c>E[it].f&&dis[E[it].v]<mx) {
-          mx=dis[E[it].v]; cur[u]=it;
+      int mx = n;
+      for (int it = G[u]; ~it; it = edges[it].nx) {
+        if (edges[it].cap > edges[it].flow && dis[edges[it].to] < mx) {
+          mx = dis[edges[it].to];
+          cur[u] = it;
         }
       }
-      if ((--gap[dis[u]])==0) break;
-      ++gap[dis[u]=mx+1]; u=pre[u];
+      if (--gap[dis[u]] == 0) break;
+      ++gap[dis[u] = mx + 1];
+      u = pre[u];
     }
     return maxflow;
   }
-  bool bfs(int S,int T) {
-    static int Q[N]; memset(dis,-1,sizeof(dis[0])*n);
-    dis[S]=0; Q[0]=S;
-    for (int h=0,t=1,u,v,it;h<t;++h) {
-      for (u=Q[h],it=G[u];~it;it=E[it].nx) {
-        if (dis[v=E[it].v]==-1&&E[it].c>E[it].f) {
-          dis[v]=dis[u]+1; Q[t++]=v;
+  bool bfs(int S, int T) {
+    static int Q[N];
+    memset(dis, -1, sizeof(*dis) * n);
+    dis[S] = 0, Q[0] = S;
+    for (int h = 0, t = 1; h < t; ++h) {
+      for (int u = Q[h], it = G[u]; ~it; it = edges[it].nx) {
+        int v = edges[it].to;
+        if (dis[v] == -1 && edges[it].cap > edges[it].flow) {
+          dis[v] = dis[u] + 1;
+          Q[t++] = v;
         }
       }
     }
-    return dis[T]!=-1;
+    return dis[T] != -1;
   }
-  int dfs(int u,int T,int low) {
-    if (u==T) return low;
-    int ret=0,tmp,v;
-    for (int &it=cur[u];~it&&ret<low;it=E[it].nx) {
-      if (dis[v=E[it].v]==dis[u]+1&&E[it].c>E[it].f) {
-        if (tmp=dfs(v,T,std::min(low-ret,E[it].c-E[it].f))) {
-          ret+=tmp; E[it].f+=tmp; E[it^1].f-=tmp;
+  flow_t dfs(int u, int T, flow_t low) {
+    if (u == T) return low;
+    flow_t ret = 0;
+    for (int &it = cur[u]; ~it && ret < low; it = edges[it].nx) {
+      int v = edges[it].to;
+      if (dis[v] == dis[u] + 1 && edges[it].cap > edges[it].flow) {
+        flow_t tmp = dfs(v, T, std::min(low - ret, edges[it].cap - edges[it].flow));
+        if (tmp > 0) {
+          ret += tmp;
+          edges[it].flow += tmp;
+          edges[it ^ 1].flow -= tmp;
         }
       }
     }
-    if (!ret) dis[u]=-1; return ret;
+    if (!ret) dis[u] = -1;
+    return ret;
   }
-  int dinic(int S,int T) {
-    int maxflow=0,tmp;
-    while (bfs(S,T)) {
-      memcpy(cur,G,sizeof(G[0])*n);
-      while (tmp=dfs(S,T,inf)) maxflow+=tmp;
+  flow_t dinic(int S, int T) {
+    flow_t maxflow = 0, tmp;
+    while (bfs(S, T)) {
+      memcpy(cur, G, sizeof(*G) * n);
+      while (tmp = dfs(S, T, inf)) {
+        maxflow += tmp;
+      }
     }
     return maxflow;
   }
diff --git a/mathematics/fast-fourier-transform.cc b/mathematics/fast-fourier-transform.cc
@@ -35,12 +35,12 @@ namespace fft {
   const double pi = acos(-1.0);
   Complex w[N];
   int rev[N];
-  void init(int n) {
-    int LN = __builtin_ctz(n);
+  void init(int L) {
+    int n = 1 << L;
     for (int i = 0 ; i < n ; ++ i) {
       double ang = 2 * pi * i / n;
       w[i] = Complex(cos(ang) , sin(ang));
-      rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (LN - 1));
+      rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (L - 1));
     }
   }
   void trans(Complex P[], int n, int oper) {
@@ -68,9 +68,9 @@ namespace fft {
   }
   Complex A[N] , B[N] , C1[N] , C2[N];
   std::vector<int64> conv(const std::vector<int> &a, const std::vector<int> &b) {
-    int n = a.size(), m = b.size(), s = 1;
-    while (s <= n + m - 2) s <<= 1;
-    init(s);
+    int n = a.size(), m = b.size(), L = 0, s = 1;
+    while (s <= n + m - 2) s <<= 1, ++L;
+    init(L);
     for (int i = 0; i < s; ++i) {
       A[i] = i < n ? Complex(a[i], 0) : Complex();
       B[i] = i < m ? Complex(b[i], 0) : Complex();
@@ -90,17 +90,41 @@ namespace fft {
     }
     return res;
   }
+  std::vector<int64> fast_conv(const std::vector<int> &a, const std::vector<int> &b) {
+    int n = a.size(), m = b.size(), L = 0, s = 1;
+    for (; s <= n + m - 2; s <<= 1, ++L);
+    s >>= 1, --L;
+    init(L);
+    for (int i = 0; i < s; ++i) {
+      A[i].x = (i << 1) < n ? a[i << 1] : 0;
+      A[i].y = (i << 1 | 1) < n ? a[i << 1 | 1] : 0;
+      B[i].x = (i << 1) <  m ? b[i << 1] : 0;
+      B[i].y = (i << 1 | 1) < m ? b[i << 1 | 1] : 0;
+    }
+    trans(A, s, 1); trans(B, s, 1);
+    for (int i = 0; i < s; ++i) {
+      int j = (s - i) & (s - 1);
+      A[i] = (Complex(4, 0) * (A[j] * B[j]).conj() - (A[j].conj() - A[i]) * (B[j].conj() - B[i]) * (w[i] + Complex(1, 0))) * Complex(0, 0.25);
+    }
+    std::reverse(w + 1, w + s);
+    trans(A, s, -1);
+    std::vector<int64> res(n + m - 1);
+    for (int i = 0; i <= (n + m - 1) / 2; ++i) {
+      res[i << 1] = int64(A[i].y + 0.5);
+      res[i << 1 | 1] = int64(A[i].x + 0.5);
+    }
+    return res;
+  }
   // arbitrary modulo convolution
   void conv(int a[], int b[], int n, int m, int mod, int res[]) {
-    int s = 1;
-    while (s <= n + m - 2) s <<= 1;
-    init(s);
+    int s = 1, L = 0;
+    while (s <= n + m - 2) s <<= 1, ++L;
+    init(L);
     for (int i = 0; i < s; ++i) {
       A[i] = i < n ? Complex(a[i] / M, a[i] % M) : Complex();
       B[i] = i < m ? Complex(b[i] / M, b[i] % M) : Complex();
     }
-    trans(A, s, 1);
-    trans(B, s, 1);
+    trans(A, s, 1); trans(B, s, 1);
     for (int i = 0; i < s; ++i) {
       int j = i ? s - i : i;
       Complex a1 = (A[i] + A[j].conj()) * Complex(0.5, 0);
@@ -112,8 +136,7 @@ namespace fft {
       C1[j] = c11 + c12 * Complex(0, 1);
       C2[j] = c21 + c22 * Complex(0, 1);
     }
-    trans(C1, s, -1);
-    trans(C2, s, -1);
+    trans(C1, s, -1); trans(C2, s, -1);
     for (int i = 0 ; i < n + m - 1; ++i) {
       int x = int64(C1[i].x + 0.5) % mod;
       int y1 = int64(C1[i].y + 0.5) % mod;
