AOJ 2691: Cost Performance Flow

aoj2691.cpp

#include <bits/stdc++.h>

using namespace std;

using Int = __int128_t;
Int abs128(Int val){return val<0?-val:val;}

ostream &operator<<(ostream &os,Int val){
  if(ostream::sentry(os)){
    __uint128_t tmp=abs128(val);
    char buf[64];
    char *d=end(buf);
    do{
      --d;
      *d=char(tmp%10+'0');
      tmp/=10;
    }while(tmp);
    if(val<0) *--d='-';
    int len=end(buf)-d;
    if(os.rdbuf()->sputn(d,len)!=len){
      os.setstate(ios_base::badbit);
    }
  }
  return os;
}

istream &operator>>(istream &is,Int &val){
  string s;
  is>>s;
  val=0;
  for(int i=0;i<(int)s.size();i++)
    if(isdigit(s[i])) val=val*10+s[i]-'0';
  if(s[0]=='-') val*=-1;
  return is;
}

#define all(v) (v).begin(), (v).end()
#define resz(v, ...) (v).clear(), (v).resize(__VA_ARGS__)
#define reps(i, m, n) for(Int i = (Int)(m); i < (Int)(n); i++)
#define rep(i, n) reps(i, 0, n)

template<class T1, class T2> void chmin(T1 &a, T2 b){if(a>b)a=b;}
template<class T1, class T2> void chmax(T1 &a, T2 b){if(a<b)a=b;}

using Pi = pair<Int, Int>;
using Tapris = tuple<Int, Int, Int>;
using vint = vector<Int>;

const Int inf = 1LL << 55;
const Int mod = 1e9 + 7;

Int sq(Int x) {
  return x*x;
}

// Sccessive Shortest Path(Primal Dual): minimum cost maximum flow
struct PrimalDual {
  struct edge {
    Int to, capacity, cost, rev;
    edge(){}
    edge(Int to, Int capacity, Int cost, Int rev):to(to), capacity(capacity), cost(cost), rev(rev){}
  };

  vector< vector<edge> > graph;
  vector<Int> potential, mincost, prevv, preve;
  map<Int, Int> mp;
  PrimalDual(Int V):graph(V), potential(V), mincost(V), prevv(V), preve(V){}
  void add_edge(Int from, Int to, Int capacity, Int cost) {
    graph[from].push_back(edge(to, capacity, cost, (Int)graph[to].size()));
    graph[to].push_back(edge(from, 0, -cost, (Int)graph[from].size()-1));
  }
  Pi min_cost_flow(Int source, Int sink, Int f) {
    Int minc = 0;
    Int maxf = 0;
    mp[0] = 0;
    fill(potential.begin(), potential.end(), 0);
    fill(prevv.begin(), prevv.end(), -1);
    fill(preve.begin(), preve.end(), -1);

    typedef pair<Int, Int> Pi;
    while(f > 0) {
      priority_queue<Pi, vector<Pi>, greater<Pi> > que;
      fill(mincost.begin(), mincost.end(), inf);
      mincost[source] = 0;
      que.push(Pi(0, source));
      while(!que.empty()) {
	Pi p = que.top(); que.pop();
	Int v = p.second;
	if(mincost[v] < p.first) continue;
	for(Int i = 0; i < (Int)graph[v].size(); i++) {
	  edge& e = graph[v][i];
	  Int dual_cost = mincost[v] + e.cost + potential[v] - potential[e.to];
	  if(e.capacity > 0 && dual_cost < mincost[e.to]) {
	    mincost[e.to] = dual_cost;
	    prevv[e.to] = v; preve[e.to] = i;
	    que.push(Pi(mincost[e.to], e.to));
	  }
	}
      }

      if(mincost[sink] == inf) return Pi(minc, maxf);
      for(Int v = 0; v < (Int)graph.size(); v++) potential[v] += mincost[v];
      Int d = f;
      for(Int v = sink; v != source; v = prevv[v]) d = min(d, graph[prevv[v]][preve[v]].capacity);
      f -= d;
      minc += d * potential[sink];
      maxf += d;
      mp[maxf] = minc;
      for(Int v = sink; v != source; v = prevv[v]) {
	edge& e = graph[prevv[v]][preve[v]];
	e.capacity -= d;
	graph[v][e.rev].capacity += d;
      }
    }
    return Pi(minc, maxf);
  }
  Pi solve(Int source, Int sink) {
    Pi res = Pi(inf, 1);
    Pi mcf = min_cost_flow(source, sink, inf);
    Int M = mcf.second;
    //cout<<mcf.first<<" "<<mcf.second<<endl;
    for(Pi p : mp) {
      chmin(res.first, sq(p.second)+sq(M-p.first));
    }
    //cout<<res.first<<endl;
    auto last = mp.end(); --last;
    for(auto it = mp.begin(); it != last; ++it) {
      auto nit = next(it);
      //cout<<it->first<<" "<<it->second << " " << nit->first << " " << nit->second << endl;
      Int a = nit->first - it->first;
      Int b = nit->second - it->second;
      {
	Int GCD = __gcd(a, b);
	a /= GCD;
	b /= GCD;
      }
      Int c = it->second;
      Int alpha = it->first;
      //cout<<a<<" "<<b<<" "<<c<<endl;
      Int F_den = b*b + a*a;
      Int F_num = b*b*alpha-a*b*c+M*a*a;
      {
	Int GCD = __gcd(F_den, F_num);
	F_den /= GCD;
	F_num /= GCD;
      }
      //cout<<F_den<<" "<<F_num<<endl;
      //cout<<it->first<<"<"<<F_num<<"/"<<F_den<<"<"<<nit->first<<endl;
      if(it->first*F_den <= F_num && F_num <= nit->first*F_den) {
	Int B_den = sq(a*F_den);
	Int B_num = sq(b*(F_num-alpha*F_den)+c*a*F_den)+sq(M*F_den-F_num)*a*a;
	{
	  Int GCD = __gcd(B_den, B_num);
	  B_den /= GCD;
	  B_num /= GCD;
	}
	//cout<<B_den<<" "<<B_num<<endl;
	if(B_num*res.second < res.first*B_den) res = Pi(B_num, B_den);
      }
    }
    return res;
  }
};


signed main() {
  cin.tie(0);
  ios_base::sync_with_stdio(0);
  cout << fixed << setprecision(12);

  Int N, M, s, t;
  cin >> N >> M >> s >> t;
  --s, --t;
  PrimalDual graph(N);
  rep(i, M) {
    Int a, b, c, d;
    cin >> a >> b >> c >> d;
    --a, --b;
    graph.add_edge(a, b, c, d);
  }
  Pi res = graph.solve(s, t);
  cout << res.first << "/" << res.second << endl;

  return 0;
}