#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"
#include<cassert>
#include<limits>
#include<random>
#include<utility>
#include<vector>#include"../../../graph/flow/flow.hpp"
#include"../../../utilities/fast_io.hpp"voidtest_max_flow(){m1une::flow::MaxFlow<longlong>mf(4);inte0=mf.add_edge(0,1,2);inte1=mf.add_edge(0,2,1);inte2=mf.add_edge(1,2,1);inte3=mf.add_edge(1,3,1);inte4=mf.add_edge(2,3,2);(void)e1;(void)e2;(void)e3;(void)e4;assert(mf.size()==4);assert(mf.edge_count()==5);assert(mf.max_flow(0,3)==3);autoedges=mf.edges();longlongoutgoing=0;for(constauto&e:edges){if(e.from==0)outgoing+=e.flow;assert(0<=e.flow&&e.flow<=e.cap);}assert(outgoing==3);assert(mf.get_edge(e0).cap==2);autocut=mf.min_cut(0);assert(cut[0]);assert(!cut[3]);mf.change_edge(e0,3,1);autochanged=mf.get_edge(e0);assert(changed.cap==3);assert(changed.flow==1);m1une::flow::MaxFlow<longlong>undirected(2);undirected.reserve_edges(1,std::vector<int>{1,1});intundirected_id=undirected.add_undirected_edge(0,1,7);undirected.change_edge(undirected_id,7,-3);autoinitial_undirected=undirected.get_edge(undirected_id);assert(initial_undirected.cap==7);assert(initial_undirected.flow==-3);assert(undirected.max_flow(0,1)==10);assert(undirected.get_edge(undirected_id).flow==7);m1une::flow::MaxFlow<longlong>limited(2);intlimited_id=limited.add_edge(0,1,10);assert(limited.max_flow(0,1,4)==4);assert(limited.get_edge(limited_id).flow==4);assert(limited.max_flow(0,1)==6);assert(limited.get_edge(limited_id).flow==10);// Different path lengths exercise repeated highest-label changes.m1une::flow::MaxFlow<longlong>layered(17);intnext_vertex=1;for(intlength=2;length<=6;length++){intfrom=0;for(intedge=0;edge<length;edge++){intto=edge+1==length?16:next_vertex++;layered.add_edge(from,to,3);from=to;}}assert(next_vertex==16);assert(layered.max_flow(0,16)==15);assert(layered.max_flow(0,16)==0);structInputEdge{intfrom;intto;longlongcap;boolundirected;};std::mt19937random(19260817);for(intiteration=0;iteration<500;iteration++){intn=2+int(random()%6);intm=int(random()%13);std::vector<InputEdge>input_edges;m1une::flow::MaxFlow<longlong>flow(n);m1une::flow::MaxFlow<longlong>push_relabel_flow(n);m1une::flow::MaxFlow<longlong>dinic_flow(n);flow.reserve_edges(m);push_relabel_flow.reserve_edges(m);dinic_flow.reserve_edges(m);for(intedge=0;edge<m;edge++){InputEdgeinput{int(random()%n),int(random()%n),1+static_cast<longlong>(random()%10),bool(random()&1)};input_edges.push_back(input);if(input.undirected){flow.add_undirected_edge(input.from,input.to,input.cap);push_relabel_flow.add_undirected_edge(input.from,input.to,input.cap);dinic_flow.add_undirected_edge(input.from,input.to,input.cap);}else{flow.add_edge(input.from,input.to,input.cap);push_relabel_flow.add_edge(input.from,input.to,input.cap);dinic_flow.add_edge(input.from,input.to,input.cap);}}longlongexpected=std::numeric_limits<longlong>::max();for(intmask=0;mask<(1<<n);mask++){if((mask&1)==0||(mask>>(n-1)&1)!=0)continue;longlongcapacity=0;for(constauto&edge:input_edges){boolfrom_side=mask>>edge.from&1;boolto_side=mask>>edge.to&1;if(from_side&&!to_side)capacity+=edge.cap;if(edge.undirected&&!from_side&&to_side){capacity+=edge.cap;}}expected=std::min(expected,capacity);}longlongresult=flow.max_flow(0,n-1);assert(result==expected);longlongpush_relabel_result=push_relabel_flow.max_flow_push_relabel(0,n-1);assert(push_relabel_result==expected);longlongdinic_result=dinic_flow.max_flow_dinic(0,n-1);assert(dinic_result==expected);autovalidate=[&](constauto&solved_flow,longlongsolved_value){std::vector<longlong>net_flow(n,0);for(intedge=0;edge<m;edge++){autoresult_edge=solved_flow.get_edge(edge);assert(result_edge.cap==input_edges[edge].cap);if(input_edges[edge].undirected){assert(-result_edge.cap<=result_edge.flow);}else{assert(0<=result_edge.flow);}assert(result_edge.flow<=result_edge.cap);net_flow[result_edge.from]+=result_edge.flow;net_flow[result_edge.to]-=result_edge.flow;}assert(net_flow[0]==solved_value);assert(net_flow[n-1]==-solved_value);for(intvertex=1;vertex+1<n;vertex++){assert(net_flow[vertex]==0);}};validate(flow,result);validate(push_relabel_flow,push_relabel_result);validate(dinic_flow,dinic_result);assert(flow.max_flow_push_relabel(0,n-1)==0);assert(push_relabel_flow.max_flow_push_relabel(0,n-1)==0);assert(push_relabel_flow.max_flow(0,n-1)==0);assert(dinic_flow.max_flow_dinic(0,n-1)==0);}}voidtest_gomory_hu(){m1une::flow::GomoryHu<longlong>gh(4);gh.add_edge(0,1,3);gh.add_edge(1,2,2);gh.add_edge(0,2,1);gh.add_edge(2,3,4);gh.build();assert(gh.size()==4);assert(gh.edge_count()==4);assert(gh.tree_edges().size()==3);assert(gh.min_cut(0,1)==4);assert(gh.min_cut(0,2)==3);assert(gh.min_cut(0,3)==3);assert(gh.min_cut(2,3)==4);m1une::flow::GomoryHu<longlong>disconnected(3);disconnected.add_edge(0,1,5);disconnected.build();assert(disconnected.min_cut(0,1)==5);assert(disconnected.min_cut(0,2)==0);m1une::flow::GomoryHu<longlong>rebuilt(2);rebuilt.build();assert(rebuilt.min_cut(0,1)==0);rebuilt.add_edge(0,1,7);rebuilt.add_edge(0,0,100);rebuilt.build();assert(rebuilt.min_cut(0,1)==7);m1une::flow::GomoryHu<longlong>singleton(1);singleton.build();assert(singleton.tree_edges().empty());std::mt19937random(123456789);for(intiteration=0;iteration<200;iteration++){intn=2+int(random()%8);structInputEdge{intu;intv;longlongcap;};std::vector<InputEdge>edges;m1une::flow::GomoryHu<longlong>tree(n);intm=random()%(2*n*n+1);for(inti=0;i<m;i++){intu=random()%n;intv=random()%n;longlongcap=random()%1000001;edges.push_back(InputEdge{u,v,cap});tree.add_edge(u,v,cap);}tree.build();for(ints=0;s<n;s++){for(intt=s+1;t<n;t++){m1une::flow::MaxFlow<longlong>mf(n);for(constauto&edge:edges){mf.add_edge(edge.u,edge.v,edge.cap);mf.add_edge(edge.v,edge.u,edge.cap);}assert(tree.min_cut(s,t)==mf.max_flow(s,t));}}}}voidtest_bounded_flow(){m1une::flow::BoundedFlow<longlong>st(4);inta=st.add_edge(0,1,1,3);intb=st.add_edge(0,2,0,2);intc=st.add_edge(1,3,1,2);intd=st.add_edge(2,3,0,2);inte=st.add_edge(1,2,0,1);(void)a;(void)b;(void)c;(void)d;(void)e;autoexact=st.feasible_st_flow(0,3,3);assert(exact.has_value());std::vector<longlong>balance(4,0);for(constauto&edge:exact->edges){assert(edge.lower<=edge.flow&&edge.flow<=edge.upper);balance[edge.from]+=edge.flow;balance[edge.to]-=edge.flow;}assert((balance==std::vector<longlong>{3,0,0,-3}));autotoo_much=st.feasible_st_flow(0,3,6);assert(!too_much.has_value());m1une::flow::BoundedFlow<longlong>bf(3);intf01=bf.add_edge(0,1,1,3);intf02=bf.add_edge(0,2,0,4);bf.add_edge(1,2,0,2);bf.add_supply(0,4);bf.add_demand(1,1);bf.add_demand(2,3);assert(bf.balance(0)==4);autobflow=bf.feasible_flow();assert(bflow.has_value());assert(bflow->get_edge(f01).flow>=1);assert(bflow->get_edge(f02).flow>=0);std::vector<longlong>b_balance(3,0);for(constauto&edge:bflow->edges){assert(edge.lower<=edge.flow&&edge.flow<=edge.upper);b_balance[edge.from]+=edge.flow;b_balance[edge.to]-=edge.flow;}assert((b_balance==std::vector<longlong>{4,-1,-3}));m1une::flow::BoundedFlow<longlong>negative(2);intneg=negative.add_edge(0,1,-5,5);negative.add_demand(0,3);negative.add_supply(1,3);autonegative_flow=negative.feasible_flow();assert(negative_flow.has_value());assert(negative_flow->flow(neg)==-3);m1une::flow::BoundedFlow<longlong>impossible(2);impossible.add_edge(0,1,0,1);impossible.add_supply(0,2);impossible.add_demand(1,2);assert(!impossible.feasible_flow().has_value());m1une::flow::BFlow<longlong>alias(1);assert(alias.size()==1);}voidtest_bounded_min_cost_flow(){m1une::flow::BoundedMinCostFlow<longlong,longlong>st(3);st.reserve_edges(3);inte01=st.add_edge(0,1,1,3,2);inte12=st.add_edge(1,2,1,3,1);inte02=st.add_edge(0,2,0,3,10);autoexact=st.min_cost_st_flow(0,2,3);assert(exact.has_value());assert(exact->cost==9);assert(exact->flow(e01)==3);assert(exact->flow(e12)==3);assert(exact->flow(e02)==0);autoexact_polynomial=st.min_cost_st_flow_polynomial(0,2,3);assert(exact_polynomial.has_value());assert(exact_polynomial->cost==9);std::vector<longlong>balance(3,0);for(constauto&edge:exact->edges){assert(edge.lower<=edge.flow&&edge.flow<=edge.upper);balance[edge.from]+=edge.flow;balance[edge.to]-=edge.flow;}assert((balance==std::vector<longlong>{3,0,-3}));m1une::flow::BoundedMinCostFlow<longlong,longlong>bf(3);intp01=bf.add_edge(0,1,0,2,1);intp12=bf.add_edge(1,2,0,2,1);intp02=bf.add_edge(0,2,0,2,5);bf.add_supply(0,2);bf.add_demand(2,2);autobflow=bf.min_cost_flow();assert(bflow.has_value());assert(bflow->cost==4);assert(bflow->flow(p01)==2);assert(bflow->flow(p12)==2);assert(bflow->flow(p02)==0);m1une::flow::BoundedMinCostFlow<longlong,longlong>negative(2);intneg=negative.add_edge(0,1,-5,5,2);negative.add_demand(0,3);negative.add_supply(1,3);autonegative_flow=negative.min_cost_flow();assert(negative_flow.has_value());assert(negative_flow->flow(neg)==-3);assert(negative_flow->cost==-6);m1une::flow::BoundedMinCostFlow<longlong,longlong>cycle(2);intc01=cycle.add_edge(0,1,0,1,-5);intc10=cycle.add_edge(1,0,0,1,3);autocirculation=cycle.min_cost_flow();assert(circulation.has_value());assert(circulation->flow(c01)==1);assert(circulation->flow(c10)==1);assert(circulation->cost==-2);autopolynomial_circulation=cycle.min_cost_flow_polynomial();assert(polynomial_circulation.has_value());assert(polynomial_circulation->cost==-2);usingImmediateFallback=m1une::flow::BoundedMinCostFlow<longlong,longlong,longlong,0>;ImmediateFallbackfallback_cycle(2);fallback_cycle.add_edge(0,1,0,1,-5);fallback_cycle.add_edge(1,0,0,1,3);autofallback_circulation=fallback_cycle.min_cost_flow();assert(fallback_circulation.has_value());assert(fallback_circulation->cost==-2);m1une::flow::BoundedMinCostFlow<longlong,longlong>impossible(2);impossible.add_edge(0,1,0,1,0);impossible.add_supply(0,2);impossible.add_demand(1,2);assert(!impossible.min_cost_flow().has_value());usingWideCostFlow=m1une::flow::BoundedMinCostFlow<longlong,longlong,__int128_t>;WideCostFlowwide_cost(1);wide_cost.reserve_edges(1);constexprlonglongtrillion=1000000000000LL;wide_cost.add_edge(0,0,trillion,trillion,trillion);autowide_result=wide_cost.min_cost_flow();assert(wide_result.has_value());assert(wide_result->cost==__int128_t(trillion)*trillion);std::mt19937random(987654321);for(intiteration=0;iteration<500;iteration++){intn=1+int(random()%4);intm=int(random()%7);structSmallEdge{intfrom;intto;longlonglower;longlongupper;longlongcost;};std::vector<SmallEdge>edges;m1une::flow::BoundedMinCostFlow<longlong,longlong>solver(n);for(inti=0;i<m;i++){intfrom=int(random()%n);intto=int(random()%n);longlonglower=static_cast<longlong>(random()%5)-2;longlongupper=lower+static_cast<longlong>(random()%4);longlongcost=static_cast<longlong>(random()%9)-4;edges.push_back(SmallEdge{from,to,lower,upper,cost});solver.add_edge(from,to,lower,upper,cost);}std::vector<longlong>required_balance(n,0);longlongbalance_sum=0;for(intvertex=0;vertex+1<n;vertex++){required_balance[vertex]=static_cast<longlong>(random()%7)-3;balance_sum+=required_balance[vertex];}required_balance.back()=-balance_sum;boolfeasible=false;longlongbest_cost=0;std::vector<longlong>flow(m);autoenumerate=[&](auto&&self,intedge_id)->void{if(edge_id!=m){for(flow[edge_id]=edges[edge_id].lower;flow[edge_id]<=edges[edge_id].upper;flow[edge_id]++){self(self,edge_id+1);}return;}std::vector<longlong>actual_balance(n,0);longlongcost=0;for(inti=0;i<m;i++){actual_balance[edges[i].from]+=flow[i];actual_balance[edges[i].to]-=flow[i];cost+=flow[i]*edges[i].cost;}if(actual_balance!=required_balance)return;if(!feasible||cost<best_cost)best_cost=cost;feasible=true;};enumerate(enumerate,0);autovalidate_result=[&](constauto&result){assert(result.has_value()==feasible);if(!result.has_value())return;assert(result->cost==best_cost);std::vector<longlong>actual_balance(n,0);for(constauto&edge:result->edges){assert(edge.lower<=edge.flow&&edge.flow<=edge.upper);actual_balance[edge.from]+=edge.flow;actual_balance[edge.to]-=edge.flow;longlongreduced_cost=edge.cost+result->potential[edge.from]-result->potential[edge.to];if(edge.flow<edge.upper)assert(0<=reduced_cost);if(edge.lower<edge.flow)assert(reduced_cost<=0);}assert(actual_balance==required_balance);};validate_result(solver.min_cost_flow(required_balance));validate_result(solver.min_cost_flow_polynomial(required_balance));}m1une::flow::BMinCostFlow<longlong,longlong>alias(1);assert(alias.size()==1);m1une::flow::BMinCostFlow<longlong,longlong,__int128_t>wide_alias(1);assert(wide_alias.size()==1);m1une::flow::BMinCostFlow<longlong,longlong,longlong,0>fallback_alias(1);assert(fallback_alias.size()==1);}voidtest_min_cost_flow(){m1une::flow::MinCostFlow<longlong,longlong>mcf(4);mcf.add_edge(0,1,2,1);mcf.add_edge(0,2,1,2);mcf.add_edge(1,2,1,0);mcf.add_edge(1,3,1,3);mcf.add_edge(2,3,2,1);autoresult=mcf.flow(0,3,2);assert(result.first==2);assert(result.second==5);autoedges=mcf.edges();longlongtotal_source_flow=0;for(constauto&e:edges){if(e.from==0)total_source_flow+=e.flow;assert(0<=e.flow&&e.flow<=e.cap);}assert(total_source_flow==2);m1une::flow::MinCostFlow<longlong,longlong>negative(3);negative.add_edge(0,1,1,-5);negative.add_edge(1,2,1,2);negative.add_edge(0,2,1,10);autoslope=negative.slope(0,2,2);std::vector<std::pair<longlong,longlong>>expected_slope={std::pair<longlong,longlong>{0,0},std::pair<longlong,longlong>{1,-3},std::pair<longlong,longlong>{2,7},};assert(slope==expected_slope);std::mt19937random(31415926);for(intiteration=0;iteration<200;iteration++){intn=2+int(random()%7);intm=int(random()%30);longlongflow_limit=random()%16;m1une::flow::MinCostFlow<longlong,longlong>tested(n);m1une::flow::BoundedMinCostFlow<longlong,longlong>expected(n);m1une::flow::MaxFlow<longlong>maximum(n);tested.reserve_edges(m);expected.reserve_edges(m);maximum.reserve_edges(m);for(intedge=0;edge<m;edge++){intfrom=int(random()%(n-1));intto=from+1+int(random()%(n-from-1));longlongcap=random()%6;longlongcost=int(random()%21)-10;tested.add_edge(from,to,cap,cost);expected.add_edge(from,to,0,cap,cost);maximum.add_edge(from,to,cap);}longlongsent=maximum.max_flow(0,n-1,flow_limit);autoexpected_result=expected.min_cost_st_flow(0,n-1,sent);assert(expected_result.has_value());autotested_result=tested.flow(0,n-1,flow_limit);assert(tested_result.first==sent);assert(tested_result.second==expected_result->cost);}// At least eight terminal arcs are required, exercising the guarded// one-shot network-simplex path rather than successive shortest paths.m1une::flow::MinCostFlow<longlong,longlong>fast(10);m1une::flow::MinCostFlow<longlong,longlong>reference(10);fast.reserve_edges(64);reference.reserve_edges(64);for(intv=1;v<=8;v++){fast.add_edge(0,v,1,v);fast.add_edge(v,9,1,10-v);reference.add_edge(0,v,1,v);reference.add_edge(v,9,1,10-v);}while(fast.edge_count()<64){intfrom=1+fast.edge_count()%8;intto=1+(fast.edge_count()*3+1)%8;if(from==to)to=to==8?1:to+1;fast.add_edge(from,to,1,100);reference.add_edge(from,to,1,100);}autofast_result=fast.flow(0,9,10);autoreference_result=reference.slope(0,9,10).back();assert(fast_result==reference_result);std::pair<longlong,longlong>expected_fast_result(8,80);assert(fast_result==expected_fast_result);assert(fast.flow(0,9).first==0);// After one shortest-path augmentation, a limit below both remaining// terminal capacities exercises the exact residual s-t solve without// terminal contraction, including negative-cost reverse arcs.m1une::flow::MinCostFlow<longlong,longlong>partial(12);m1une::flow::MinCostFlow<longlong,longlong>partial_reference(12);for(intv=1;v<=10;v++){partial.add_edge(0,v,1,v);partial.add_edge(v,11,1,0);partial_reference.add_edge(0,v,1,v);partial_reference.add_edge(v,11,1,0);}while(partial.edge_count()<64){intfrom=1+partial.edge_count()%10;intto=1+(partial.edge_count()*7+1)%10;if(from==to)to=to==10?1:to+1;partial.add_edge(from,to,1,100);partial_reference.add_edge(from,to,1,100);}assert(partial.flow(0,11,1)==partial_reference.slope(0,11,1).back());assert(partial.flow(0,11,8)==partial_reference.slope(0,11,8).back());assert(partial.flow(0,11)==partial_reference.slope(0,11).back());// The terminal cut permits eight units but an internal bottleneck permits// only four, exercising the max-flow fallback before the exact-cost solve.m1une::flow::MinCostFlow<longlong,longlong>bottleneck(20);for(intv=1;v<=8;v++){bottleneck.add_edge(0,v,1,0);bottleneck.add_edge(v,9,1,0);bottleneck.add_edge(10,v+10,1,0);bottleneck.add_edge(v+10,19,1,0);}bottleneck.add_edge(9,10,4,1);while(bottleneck.edge_count()<64){intfrom=1+bottleneck.edge_count()%8;intto=1+(bottleneck.edge_count()*5+1)%8;if(from==to)to=to==8?1:to+1;bottleneck.add_edge(from,to,1,100);}std::pair<longlong,longlong>expected_bottleneck(4,4);assert(bottleneck.flow(0,19,8)==expected_bottleneck);// Randomized comparisons cover the terminal-contraction optimization on// feasible non-negative-cost instances.for(intiteration=0;iteration<100;iteration++){constexprintn=18;m1une::flow::MinCostFlow<longlong,longlong>adaptive(n);m1une::flow::MinCostFlow<longlong,longlong>shortest_paths(n);adaptive.reserve_edges(80);shortest_paths.reserve_edges(80);for(intv=1;v<=8;v++){longlongfirst_cost=random()%11;longlongmiddle_cost=random()%11;longlonglast_cost=random()%11;adaptive.add_edge(0,v,1,first_cost);adaptive.add_edge(v,v+8,1,middle_cost);adaptive.add_edge(v+8,17,1,last_cost);shortest_paths.add_edge(0,v,1,first_cost);shortest_paths.add_edge(v,v+8,1,middle_cost);shortest_paths.add_edge(v+8,17,1,last_cost);}while(adaptive.edge_count()<80){intfrom=1+int(random()%16);intto=1+int(random()%16);if(from==to)to=to==16?1:to+1;longlongcap=1+random()%3;longlongcost=random()%21;adaptive.add_edge(from,to,cap,cost);shortest_paths.add_edge(from,to,cap,cost);}assert(adaptive.flow(0,17,8)==shortest_paths.slope(0,17,8).back());}}intmain(){m1une::utilities::FastInputfast_input;m1une::utilities::FastOutputfast_output;test_max_flow();test_gomory_hu();test_bounded_flow();test_bounded_min_cost_flow();test_min_cost_flow();longlonga,b;fast_input>>a>>b;fast_output<<a+b<<'\n';}
#line 1 "verify/graph/flow/flow_algorithms.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/aplusb"
#include<cassert>
#include<limits>
#include<random>
#include<utility>
#include<vector>#line 1 "graph/flow/flow.hpp"
#line 1 "graph/flow/bounded_flow.hpp"
#line 5 "graph/flow/bounded_flow.hpp"
#include<optional>
#line 7 "graph/flow/bounded_flow.hpp"
#line 1 "graph/flow/max_flow.hpp"
#include<algorithm>
#line 6 "graph/flow/max_flow.hpp"
#include<cstddef>
#line 9 "graph/flow/max_flow.hpp"
namespacem1une{namespaceflow{template<classCap>structMaxFlow{structEdge{intfrom;intto;Capcap;Capflow;};private:structInternalEdge{intto;intrev;Capcap;};structPosition{intfrom;intedge;};int_n;std::vector<Position>_pos;std::vector<std::vector<InternalEdge>>_g;Caphighest_label_preflow_push(ints,intt){constintdead=2*_n;constintunreachable=_n+1;std::vector<Cap>excess(_n,Cap(0));std::vector<int>state(8*std::size_t(_n)+2);int*height=state.data();int*height_count=height+_n;int*current=height_count+dead+1;int*queue=current+_n;int*next=queue+_n;int*bucket_head=next+_n;std::vector<char>active(_n,false);inthighest=-1;longlongwork=0;constlonglongarc_count=2LL*static_cast<longlong>(_pos.size());constlonglongwork_limit=std::max(1LL,4*arc_count+_n);autoactivate=[&](intv){if(v==s||v==t||active[v]||excess[v]==Cap(0)||height[v]>=dead){return;}active[v]=true;next[v]=bucket_head[height[v]];bucket_head[height[v]]=v;highest=std::max(highest,height[v]);};autorebuild_buckets=[&](){std::fill(bucket_head,bucket_head+dead+1,-1);std::fill(active.begin(),active.end(),false);highest=-1;for(intv=0;v<_n;v++)activate(v);};autoglobal_relabel=[&](){std::fill(height,height+_n,unreachable);std::fill(height_count,height_count+dead+1,0);std::fill(current,current+_n,0);inthead=0;inttail=0;height[t]=0;height[s]=_n;queue[tail++]=t;while(head!=tail){intv=queue[head++];for(constauto&e:_g[v]){if(e.to==s||height[e.to]!=unreachable)continue;constauto&reverse=_g[e.to][e.rev];if(reverse.cap==Cap(0))continue;height[e.to]=height[v]+1;queue[tail++]=e.to;}}for(intv=0;v<_n;v++)height_count[height[v]]++;rebuild_buckets();work=0;};autogap=[&](intempty_height){for(intv=0;v<_n;v++){if(v==s||v==t||height[v]<=empty_height||height[v]>=_n){continue;}height_count[height[v]]--;height[v]=unreachable;height_count[height[v]]++;current[v]=0;}rebuild_buckets();};autorelabel=[&](intv)->bool{intold_height=height[v];intnew_height=dead;work+=int(_g[v].size());for(constauto&e:_g[v]){if(e.cap!=Cap(0)){new_height=std::min(new_height,height[e.to]+1);}}height_count[old_height]--;height[v]=std::min(new_height,dead);height_count[height[v]]++;current[v]=0;if(old_height<_n&&height_count[old_height]==0){gap(old_height);returntrue;}returnfalse;};autopush=[&](intv,InternalEdge&e){Capsent=std::min(excess[v],e.cap);boolwas_zero=excess[e.to]==Cap(0);e.cap-=sent;_g[e.to][e.rev].cap+=sent;excess[v]-=sent;excess[e.to]+=sent;if(was_zero)activate(e.to);};autodischarge=[&](intv){while(excess[v]!=Cap(0)&&height[v]<dead){if(current[v]==int(_g[v].size())){if(relabel(v))return;continue;}auto&e=_g[v][current[v]];work++;if(e.cap!=Cap(0)&&height[v]==height[e.to]+1){push(v,e);}else{current[v]++;}}activate(v);};for(auto&e:_g[s]){if(e.to==s||e.cap==Cap(0))continue;Capsent=e.cap;e.cap=Cap(0);_g[e.to][e.rev].cap+=sent;excess[e.to]+=sent;}global_relabel();while(highest>=0){if(bucket_head[highest]==-1){highest--;continue;}intv=bucket_head[highest];bucket_head[highest]=next[v];if(!active[v]||height[v]!=highest)continue;active[v]=false;discharge(v);if(work>=work_limit)global_relabel();}returnexcess[t];}public:MaxFlow():MaxFlow(0){}explicitMaxFlow(intn):_n(n),_g(n){assert(0<=n);}intsize()const{return_n;}intedge_count()const{returnint(_pos.size());}voidreserve_edges(intedge_count){assert(0<=edge_count);_pos.reserve(edge_count);if(_n==0||edge_count==0||2*std::size_t(edge_count)<std::size_t(_n)){return;}conststd::size_taverage_degree=(3*std::size_t(edge_count)+std::size_t(_n)-1)/std::size_t(_n);for(auto&edges:_g)edges.reserve(average_degree);}voidreserve_edges(intedge_count,conststd::vector<int>°rees){assert(0<=edge_count);assert(int(degrees.size())==_n);_pos.reserve(edge_count);for(intv=0;v<_n;v++){assert(0<=degrees[v]);_g[v].reserve(degrees[v]);}}intadd_edge(intfrom,intto,Capcap){assert(0<=from&&from<_n);assert(0<=to&&to<_n);assert(Cap(0)<=cap);intid=int(_pos.size());intfrom_id=int(_g[from].size());intto_id=int(_g[to].size());if(from==to)to_id++;_pos.push_back(Position{from,from_id});_g[from].push_back(InternalEdge{to,to_id,cap});_g[to].push_back(InternalEdge{from,from_id,Cap(0)});returnid;}intadd_undirected_edge(intfirst,intsecond,Capcap){static_assert(std::numeric_limits<Cap>::is_signed);assert(0<=first&&first<_n);assert(0<=second&&second<_n);assert(Cap(0)<=cap);assert(cap<=std::numeric_limits<Cap>::max()/Cap(2));intid=int(_pos.size());intfirst_id=int(_g[first].size());intsecond_id=int(_g[second].size());if(first==second)second_id++;_pos.push_back(Position{first,~first_id});_g[first].push_back(InternalEdge{second,second_id,cap});_g[second].push_back(InternalEdge{first,first_id,cap});returnid;}Edgeget_edge(inti)const{assert(0<=i&&i<int(_pos.size()));constauto&position=_pos[i];intfrom=position.from;boolundirected=position.edge<0;intidx=undirected?~position.edge:position.edge;constauto&e=_g[from][idx];constauto&re=_g[e.to][e.rev];if(undirected){returnEdge{from,e.to,(e.cap+re.cap)/Cap(2),(re.cap-e.cap)/Cap(2)};}returnEdge{from,e.to,e.cap+re.cap,re.cap};}std::vector<Edge>edges()const{std::vector<Edge>result;result.reserve(_pos.size());for(inti=0;i<int(_pos.size());i++)result.push_back(get_edge(i));returnresult;}voidchange_edge(inti,Capnew_cap,Capnew_flow){assert(0<=i&&i<int(_pos.size()));assert(Cap(0)<=new_cap);auto&position=_pos[i];intfrom=position.from;boolundirected=position.edge<0;intidx=undirected?~position.edge:position.edge;auto&e=_g[from][idx];auto&re=_g[e.to][e.rev];if(undirected){assert(new_cap<=std::numeric_limits<Cap>::max()/Cap(2));assert(-new_cap<=new_flow&&new_flow<=new_cap);e.cap=new_cap-new_flow;re.cap=new_cap+new_flow;}else{assert(Cap(0)<=new_flow&&new_flow<=new_cap);e.cap=new_cap-new_flow;re.cap=new_flow;}}Capmax_flow(ints,intt){assert(0<=s&&s<_n);assert(0<=t&&t<_n);assert(s!=t);returnhighest_label_preflow_push(s,t);}Capmax_flow_push_relabel(ints,intt){assert(0<=s&&s<_n);assert(0<=t&&t<_n);assert(s!=t);returnhighest_label_preflow_push(s,t);}Capmax_flow_dinic(ints,intt){returnmax_flow(s,t,std::numeric_limits<Cap>::max());}Capmax_flow(ints,intt,Capflow_limit){assert(0<=s&&s<_n);assert(0<=t&&t<_n);assert(s!=t);std::vector<int>work(3*std::size_t(_n));int*level=work.data();int*iter=level+_n;int*queue=iter+_n;autobfs=[&]()->bool{std::fill(level,level+_n,-1);inthead=0;inttail=0;level[s]=0;queue[tail++]=s;while(head!=tail){intv=queue[head++];for(constauto&e:_g[v]){if(level[e.to]!=-1||e.cap==Cap(0))continue;level[e.to]=level[v]+1;if(e.to==t)returntrue;queue[tail++]=e.to;}}returnlevel[t]!=-1;};autodfs=[&](auto&&self,intv,Capup)->Cap{if(v==s)returnup;Capresult=Cap(0);constintcurrent_level=level[v];auto&edges=_g[v];constintedge_count=int(edges.size());for(int&i=iter[v];i<edge_count;i++){auto&e=edges[i];if(level[e.to]+1!=current_level)continue;auto&reverse=_g[e.to][e.rev];if(reverse.cap==Cap(0))continue;Capd=self(self,e.to,std::min(up-result,reverse.cap));if(d==Cap(0))continue;e.cap+=d;reverse.cap-=d;result+=d;if(result==up)returnresult;}level[v]=_n;returnresult;};Capflow=0;while(flow<flow_limit&&bfs()){std::fill(iter,iter+_n,0);flow+=dfs(dfs,t,flow_limit-flow);}returnflow;}std::vector<bool>min_cut(ints)const{assert(0<=s&&s<_n);std::vector<bool>visited(_n,false);std::vector<int>queue(_n);inthead=0;inttail=0;visited[s]=true;queue[tail++]=s;while(head!=tail){intv=queue[head++];for(constauto&e:_g[v]){if(e.cap==Cap(0)||visited[e.to])continue;visited[e.to]=true;queue[tail++]=e.to;}}returnvisited;}};}// namespace flow}// namespace m1une#line 9 "graph/flow/bounded_flow.hpp"
namespacem1une{namespaceflow{template<classCap>structBoundedFlow{structEdge{intfrom;intto;Caplower;Capupper;};structResultEdge{intfrom;intto;Caplower;Capupper;Capflow;};structResult{std::vector<ResultEdge>edges;std::vector<Cap>balance;ResultEdgeget_edge(inti)const{assert(0<=i&&i<int(edges.size()));returnedges[i];}Capflow(inti)const{assert(0<=i&&i<int(edges.size()));returnedges[i].flow;}};private:int_n;std::vector<Edge>_edges;std::vector<Cap>_balance;public:BoundedFlow():BoundedFlow(0){}explicitBoundedFlow(intn):_n(n),_balance(n,Cap(0)){assert(0<=n);}intsize()const{return_n;}intedge_count()const{returnint(_edges.size());}intadd_edge(intfrom,intto,Caplower,Capupper){assert(0<=from&&from<_n);assert(0<=to&&to<_n);assert(lower<=upper);intid=int(_edges.size());_edges.push_back(Edge{from,to,lower,upper});returnid;}Edgeget_edge(inti)const{assert(0<=i&&i<int(_edges.size()));return_edges[i];}std::vector<Edge>edges()const{return_edges;}voidset_balance(intv,Capb){assert(0<=v&&v<_n);_balance[v]=b;}voidadd_balance(intv,Capb){assert(0<=v&&v<_n);_balance[v]+=b;}voidadd_supply(intv,Capsupply){assert(Cap(0)<=supply);add_balance(v,supply);}voidadd_demand(intv,Capdemand){assert(Cap(0)<=demand);add_balance(v,-demand);}Capbalance(intv)const{assert(0<=v&&v<_n);return_balance[v];}conststd::vector<Cap>&balances()const{return_balance;}std::optional<Result>feasible_flow()const{returnfeasible_flow(_balance);}std::optional<Result>feasible_flow(conststd::vector<Cap>&balance)const{assert(int(balance.size())==_n);intss=_n,tt=_n+1;MaxFlow<Cap>mf(_n+2);std::vector<int>edge_ids;edge_ids.reserve(_edges.size());std::vector<Cap>need=balance;for(constauto&e:_edges){edge_ids.push_back(mf.add_edge(e.from,e.to,e.upper-e.lower));need[e.from]-=e.lower;need[e.to]+=e.lower;}Cappositive_sum=Cap(0),negative_sum=Cap(0);for(intv=0;v<_n;v++){if(need[v]>Cap(0)){positive_sum+=need[v];mf.add_edge(ss,v,need[v]);}elseif(need[v]<Cap(0)){negative_sum+=-need[v];mf.add_edge(v,tt,-need[v]);}}if(positive_sum!=negative_sum)returnstd::nullopt;if(mf.max_flow(ss,tt)!=positive_sum)returnstd::nullopt;Resultresult;result.balance=balance;result.edges.reserve(_edges.size());for(inti=0;i<int(_edges.size());i++){autoused=mf.get_edge(edge_ids[i]).flow;constauto&e=_edges[i];result.edges.push_back(ResultEdge{e.from,e.to,e.lower,e.upper,e.lower+used});}returnresult;}std::optional<Result>feasible_st_flow(ints,intt,Capflow_value)const{assert(0<=s&&s<_n);assert(0<=t&&t<_n);assert(s!=t);std::vector<Cap>balance=_balance;balance[s]+=flow_value;balance[t]-=flow_value;returnfeasible_flow(balance);}};template<classCap>usingBFlow=BoundedFlow<Cap>;}// namespace flow}// namespace m1une#line 1 "graph/flow/bounded_min_cost_flow.hpp"
#line 7 "graph/flow/bounded_min_cost_flow.hpp"
#include<cmath>
#include<functional>
#line 11 "graph/flow/bounded_min_cost_flow.hpp"
#include<queue>
#line 14 "graph/flow/bounded_min_cost_flow.hpp"
namespacem1une{namespaceflow{template<classCap,classCost,classTotalCost=Cost,std::size_tPivotLimitFactor=8>structBoundedMinCostFlow{static_assert(std::numeric_limits<Cap>::is_integer);static_assert(std::numeric_limits<Cap>::is_signed);static_assert(std::numeric_limits<Cost>::is_specialized);static_assert(std::numeric_limits<Cost>::is_signed);structEdge{intfrom;intto;Caplower;Capupper;Costcost;};structResultEdge{intfrom;intto;Caplower;Capupper;Capflow;Costcost;};structResult{std::vector<ResultEdge>edges;std::vector<Cap>balance;std::vector<Cost>potential;TotalCostcost;ResultEdgeget_edge(inti)const{assert(0<=i&&i<int(edges.size()));returnedges[i];}Capflow(inti)const{assert(0<=i&&i<int(edges.size()));returnedges[i].flow;}};private:structNetworkEdge{intto;Capcap;Costcost;};structNetworkSimplexSolver{enumclassStatus{optimal,infeasible,pivot_limit_reached,};structParent{intvertex;intedge;Capup;Capdown;};intn;std::vector<NetworkEdge>edges;std::vector<Cap>excess;std::vector<Cost>potential;std::size_tpivot_count=0;NetworkSimplexSolver(intvertex_count,conststd::vector<Cap>&balance):n(vertex_count),excess(balance){}voidreserve_edges(intedge_count){edges.reserve(2*(edge_count+n));}intadd_edge(intfrom,intto,Caplower,Capupper,Costcost){intid=int(edges.size())/2;edges.push_back(NetworkEdge{to,upper-lower,cost});edges.push_back(NetworkEdge{from,Cap(0),-cost});excess[from]-=lower;excess[to]+=lower;returnid;}Statussolve(std::size_tpivot_limit){pivot_count=0;constintoriginal_edge_count=int(edges.size());potential.assign(n+1,Cost(0));Costartificial_cost=Cost(1);for(intedge=0;edge<original_edge_count;edge+=2){artificial_cost+=edges[edge].cost<Cost(0)?-edges[edge].cost:edges[edge].cost;}std::vector<Parent>parent(n);edges.reserve(original_edge_count+2*n);for(intvertex=0;vertex<n;vertex++){if(excess[vertex]>=Cap(0)){edges.push_back(NetworkEdge{n,Cap(0),artificial_cost});edges.push_back(NetworkEdge{vertex,excess[vertex],-artificial_cost});potential[vertex]=-artificial_cost;}else{edges.push_back(NetworkEdge{n,-excess[vertex],-artificial_cost});edges.push_back(NetworkEdge{vertex,Cap(0),artificial_cost});potential[vertex]=artificial_cost;}intedge=int(edges.size())-2;parent[vertex]=Parent{n,edge,edges[edge].cap,edges[edge^1].cap};}std::vector<int>depth(n+1,1);depth[n]=0;std::vector<int>next(2*(n+1));std::vector<int>previous(2*(n+1));autoconnect=[&](intfirst,intsecond){next[first]=second;previous[second]=first;};for(intvertex=0;vertex<=n;vertex++){connect(2*vertex,2*vertex+1);}for(intvertex=0;vertex<n;vertex++){connect(2*vertex+1,next[2*n]);connect(2*n,2*vertex);}autopush_flow=[&](intentering_edge){constintfirst=edges[entering_edge^1].to;constintsecond=edges[entering_edge].to;constCostcycle_cost=edges[entering_edge].cost+potential[first]-potential[second];Capamount=edges[entering_edge].cap;boolleave_first_side=true;intleaving_vertex=second;intfirst_ancestor=first;intsecond_ancestor=second;automove_first_up=[&]{if(parent[first_ancestor].down<amount){amount=parent[first_ancestor].down;leaving_vertex=first_ancestor;leave_first_side=true;}first_ancestor=parent[first_ancestor].vertex;};automove_second_up=[&]{if(parent[second_ancestor].up<=amount){amount=parent[second_ancestor].up;leaving_vertex=second_ancestor;leave_first_side=false;}second_ancestor=parent[second_ancestor].vertex;};if(depth[first_ancestor]>=depth[second_ancestor]){intdifference=depth[first_ancestor]-depth[second_ancestor];for(inti=0;i<difference;i++)move_first_up();}else{intdifference=depth[second_ancestor]-depth[first_ancestor];for(inti=0;i<difference;i++)move_second_up();}while(first_ancestor!=second_ancestor){move_first_up();move_second_up();}constintancestor=first_ancestor;if(amount!=Cap(0)){intvertex=first;while(vertex!=ancestor){parent[vertex].up+=amount;parent[vertex].down-=amount;vertex=parent[vertex].vertex;}vertex=second;while(vertex!=ancestor){parent[vertex].up-=amount;parent[vertex].down+=amount;vertex=parent[vertex].vertex;}}intvertex=first;intnew_parent=second;std::pair<Cap,Cap>parent_capacities{edges[entering_edge].cap-amount,edges[entering_edge^1].cap+amount};Costpotential_difference=-cycle_cost;if(!leave_first_side){std::swap(vertex,new_parent);std::swap(parent_capacities.first,parent_capacities.second);potential_difference=-potential_difference;}intparent_edge=entering_edge^(leave_first_side?0:1);while(new_parent!=leaving_vertex){intnew_depth=depth[new_parent];inttour_index=2*vertex;while(tour_index!=2*vertex+1){if((tour_index&1)==0){new_depth++;potential[tour_index/2]+=potential_difference;depth[tour_index/2]=new_depth;}else{new_depth--;}tour_index=next[tour_index];}connect(previous[2*vertex],next[2*vertex+1]);connect(2*vertex+1,next[2*new_parent]);connect(2*new_parent,2*vertex);std::swap(parent[vertex].edge,parent_edge);parent_edge^=1;std::swap(parent[vertex].up,parent_capacities.first);std::swap(parent[vertex].down,parent_capacities.second);std::swap(parent_capacities.first,parent_capacities.second);intold_parent=parent[vertex].vertex;parent[vertex].vertex=new_parent;new_parent=vertex;vertex=old_parent;}edges[parent_edge].cap=parent_capacities.first;edges[parent_edge^1].cap=parent_capacities.second;};boolpivot_limit_reached=false;autopivot=[&](intentering_edge){if(pivot_count==pivot_limit){pivot_limit_reached=true;returnfalse;}push_flow(entering_edge);pivot_count++;returntrue;};constintcandidate_limit=std::max(int(0.2*std::sqrt(double(original_edge_count))),10);constintminor_limit=std::max(candidate_limit/10,3);std::vector<int>candidates;candidates.reserve(candidate_limit);autominor_pivot=[&]{Costbest_cost=Cost(0);intbest_edge=-1;intindex=0;while(index<int(candidates.size())){intedge=candidates[index];if(edges[edge].cap==Cap(0)){candidates[index]=candidates.back();candidates.pop_back();continue;}Costreduced_cost=edges[edge].cost+potential[edges[edge^1].to]-potential[edges[edge].to];if(reduced_cost>=Cost(0)){candidates[index]=candidates.back();candidates.pop_back();continue;}if(reduced_cost<best_cost){best_cost=reduced_cost;best_edge=edge;}index++;}if(best_edge==-1)returnfalse;returnpivot(best_edge);};intedge=0;while(true){for(intiteration=0;iteration<minor_limit;iteration++){if(!minor_pivot())break;}if(pivot_limit_reached)returnStatus::pivot_limit_reached;Costbest_cost=Cost(0);intbest_edge=-1;candidates.clear();for(intscanned=0;scanned<int(edges.size());scanned++){if(edges[edge].cap!=Cap(0)){Costreduced_cost=edges[edge].cost+potential[edges[edge^1].to]-potential[edges[edge].to];if(reduced_cost<Cost(0)){if(reduced_cost<best_cost){best_cost=reduced_cost;best_edge=edge;}candidates.push_back(edge);if(int(candidates.size())==candidate_limit)break;}}edge++;if(edge==int(edges.size()))edge=0;}if(candidates.empty())break;if(!pivot(best_edge))returnStatus::pivot_limit_reached;}for(intvertex=0;vertex<n;vertex++){edges[parent[vertex].edge].cap=parent[vertex].up;edges[parent[vertex].edge^1].cap=parent[vertex].down;}boolfeasible=true;for(intvertex=0;vertex<n;vertex++){intartificial_edge=original_edge_count+2*vertex;if((excess[vertex]>=Cap(0)&&edges[artificial_edge^1].cap!=Cap(0))||(excess[vertex]<Cap(0)&&edges[artificial_edge].cap!=Cap(0))){feasible=false;break;}}potential.pop_back();returnfeasible?Status::optimal:Status::infeasible;}Capedge_flow(intedge_id,Caplower)const{returnlower+edges[2*edge_id+1].cap;}};structScalingEdge{intto;intreverse;Capcap;Capflow;Costcost;};structScalingSolver{intn;std::vector<std::vector<ScalingEdge>>graph;std::vector<std::pair<int,int>>positions;std::vector<Cap>excess;std::vector<Cost>potential;std::vector<Cost>distance;std::vector<int>parent_vertex;std::vector<int>parent_edge;std::vector<int>excess_vertices;std::vector<int>deficit_vertices;Costfarthest=Cost(0);ScalingSolver(intvertex_count,conststd::vector<Cap>&balance):n(vertex_count),graph(vertex_count),excess(balance),potential(vertex_count,Cost(0)){}voidreserve_edges(intedge_count){positions.reserve(edge_count);}intadd_edge(intfrom,intto,Caplower,Capupper,Costcost){intid=int(positions.size());intfrom_edge=int(graph[from].size());intto_edge=int(graph[to].size());if(from==to)to_edge++;positions.emplace_back(from,from_edge);graph[from].push_back(ScalingEdge{to,to_edge,upper,Cap(0),cost});graph[to].push_back(ScalingEdge{from,from_edge,-lower,Cap(0),-cost});returnid;}Capresidual_capacity(intfrom,intedge_id)const{constauto&edge=graph[from][edge_id];returnedge.cap-edge.flow;}Costresidual_cost(intfrom,constScalingEdge&edge)const{returnedge.cost+potential[from]-potential[edge.to];}voidpush(intfrom,intedge_id,Capamount){auto&edge=graph[from][edge_id];edge.flow+=amount;graph[edge.to][edge.reverse].flow-=amount;}voidsaturate_negative(Capdelta){excess_vertices.clear();deficit_vertices.clear();for(intfrom=0;from<n;from++){for(intedge_id=0;edge_id<int(graph[from].size());edge_id++){constauto&edge=graph[from][edge_id];Capresidual=edge.cap-edge.flow;residual-=residual%delta;if(residual_cost(from,edge)<Cost(0)||residual<Cap(0)){intto=edge.to;push(from,edge_id,residual);excess[from]-=residual;excess[to]+=residual;}}}for(intvertex=0;vertex<n;vertex++){if(excess[vertex]>Cap(0)){excess_vertices.push_back(vertex);}elseif(excess[vertex]<Cap(0)){deficit_vertices.push_back(vertex);}}}booldual(Capdelta){excess_vertices.erase(std::remove_if(excess_vertices.begin(),excess_vertices.end(),[&](intvertex){returnexcess[vertex]<delta;}),excess_vertices.end());deficit_vertices.erase(std::remove_if(deficit_vertices.begin(),deficit_vertices.end(),[&](intvertex){returnexcess[vertex]>-delta;}),deficit_vertices.end());constCostunreachable=std::numeric_limits<Cost>::max();distance.assign(n,unreachable);parent_vertex.assign(n,-1);parent_edge.assign(n,-1);usingQueueEntry=std::pair<Cost,int>;std::priority_queue<QueueEntry,std::vector<QueueEntry>,std::greater<QueueEntry>>queue;for(intvertex:excess_vertices){distance[vertex]=Cost(0);queue.emplace(Cost(0),vertex);}farthest=Cost(0);intreached_deficits=0;while(!queue.empty()){auto[current_distance,from]=queue.top();queue.pop();if(distance[from]!=current_distance)continue;farthest=current_distance;if(excess[from]<=-delta)reached_deficits++;if(reached_deficits>=int(deficit_vertices.size()))break;for(intedge_id=0;edge_id<int(graph[from].size());edge_id++){constauto&edge=graph[from][edge_id];if(edge.cap-edge.flow<delta)continue;Costnext_distance=current_distance+residual_cost(from,edge);if(next_distance>=distance[edge.to])continue;distance[edge.to]=next_distance;parent_vertex[edge.to]=from;parent_edge[edge.to]=edge_id;queue.emplace(next_distance,edge.to);}}for(intvertex=0;vertex<n;vertex++){potential[vertex]+=std::min(distance[vertex],farthest);}returnreached_deficits>0;}voidprimal(Capdelta){for(intsink:deficit_vertices){if(distance[sink]>farthest)continue;Capamount=-excess[sink];introot=sink;while(parent_edge[root]!=-1){intfrom=parent_vertex[root];amount=std::min(amount,residual_capacity(from,parent_edge[root]));root=from;}amount=std::min(amount,excess[root]);amount-=amount%delta;if(amount<=Cap(0))continue;intvertex=sink;while(parent_edge[vertex]!=-1){intfrom=parent_vertex[vertex];intedge_id=parent_edge[vertex];push(from,edge_id,amount);if(residual_capacity(from,edge_id)==Cap(0)){parent_edge[vertex]=-1;}vertex=from;}excess[sink]+=amount;excess[root]-=amount;}}boolsolve(){Capscale_bound=Cap(1);for(Capvalue:excess){scale_bound=std::max(scale_bound,value);scale_bound=std::max(scale_bound,-value);}for(constauto&edges:graph){for(constauto&edge:edges){Capresidual=edge.cap-edge.flow;scale_bound=std::max(scale_bound,residual);scale_bound=std::max(scale_bound,-residual);}}Capdelta=Cap(1);while(delta<=scale_bound/Cap(2))delta*=Cap(2);while(true){saturate_negative(delta);while(dual(delta))primal(delta);if(delta==Cap(1))break;delta/=Cap(2);}returnexcess_vertices.empty()&&deficit_vertices.empty();}Capedge_flow(intedge_id,Cap)const{auto[from,index]=positions[edge_id];returngraph[from][index].flow;}};int_n;std::vector<Edge>_edges;std::vector<Cap>_balance;template<classSolver>Resultmake_result(conststd::vector<Cap>&balance,constSolver&solver,std::vector<Cost>potential)const{Resultresult;result.balance=balance;result.cost=TotalCost(0);result.edges.reserve(_edges.size());for(inti=0;i<int(_edges.size());i++){constauto&edge=_edges[i];Capflow=solver.edge_flow(i,edge.lower);result.cost+=TotalCost(flow)*TotalCost(edge.cost);result.edges.push_back(ResultEdge{edge.from,edge.to,edge.lower,edge.upper,flow,edge.cost});}result.potential=std::move(potential);returnresult;}std::vector<Cost>residual_potential(conststd::vector<ResultEdge>&edges)const{std::vector<Cost>potential(_n,Cost(0));boolupdated=false;for(intiteration=0;iteration<_n;iteration++){updated=false;for(constResultEdge&edge:edges){if(edge.flow<edge.upper&&potential[edge.to]>potential[edge.from]+edge.cost){potential[edge.to]=potential[edge.from]+edge.cost;updated=true;}if(edge.lower<edge.flow&&potential[edge.from]>potential[edge.to]-edge.cost){potential[edge.from]=potential[edge.to]-edge.cost;updated=true;}}if(!updated)break;}assert(!updated);returnpotential;}std::optional<Result>polynomial_min_cost_flow_impl(conststd::vector<Cap>&balance)const{ScalingSolversolver(_n,balance);solver.reserve_edges(int(_edges.size()));for(constauto&edge:_edges){solver.add_edge(edge.from,edge.to,edge.lower,edge.upper,edge.cost);}if(!solver.solve())returnstd::nullopt;Resultresult=make_result(balance,solver,{});result.potential=residual_potential(result.edges);returnresult;}public:BoundedMinCostFlow():BoundedMinCostFlow(0){}explicitBoundedMinCostFlow(intn):_n(n),_balance(n,Cap(0)){assert(0<=n);}intsize()const{return_n;}intedge_count()const{returnint(_edges.size());}voidreserve_edges(intedge_count){assert(0<=edge_count);_edges.reserve(edge_count);}intadd_edge(intfrom,intto,Caplower,Capupper,Costcost){assert(0<=from&&from<_n);assert(0<=to&&to<_n);assert(lower<=upper);intid=int(_edges.size());_edges.push_back(Edge{from,to,lower,upper,cost});returnid;}Edgeget_edge(inti)const{assert(0<=i&&i<int(_edges.size()));return_edges[i];}std::vector<Edge>edges()const{return_edges;}voidset_balance(intv,Capb){assert(0<=v&&v<_n);_balance[v]=b;}voidadd_balance(intv,Capb){assert(0<=v&&v<_n);_balance[v]+=b;}voidadd_supply(intv,Capsupply){assert(Cap(0)<=supply);add_balance(v,supply);}voidadd_demand(intv,Capdemand){assert(Cap(0)<=demand);add_balance(v,-demand);}Capbalance(intv)const{assert(0<=v&&v<_n);return_balance[v];}conststd::vector<Cap>&balances()const{return_balance;}std::optional<Result>min_cost_flow()const{returnmin_cost_flow(_balance);}std::optional<Result>min_cost_flow(conststd::vector<Cap>&balance)const{assert(int(balance.size())==_n);Capbalance_sum=Cap(0);for(Capvalue:balance)balance_sum+=value;if(balance_sum!=Cap(0))returnstd::nullopt;NetworkSimplexSolversolver(_n,balance);solver.reserve_edges(int(_edges.size()));for(constauto&edge:_edges){solver.add_edge(edge.from,edge.to,edge.lower,edge.upper,edge.cost);}conststd::size_tgraph_size=std::size_t(_n)+_edges.size()+1;std::size_tpivot_limit=0;ifconstexpr(PivotLimitFactor!=0){conststd::size_tmaximum=std::numeric_limits<std::size_t>::max();pivot_limit=graph_size>maximum/PivotLimitFactor?maximum:PivotLimitFactor*graph_size;}autostatus=solver.solve(pivot_limit);if(status==NetworkSimplexSolver::Status::infeasible){returnstd::nullopt;}if(status==NetworkSimplexSolver::Status::pivot_limit_reached){returnpolynomial_min_cost_flow_impl(balance);}returnmake_result(balance,solver,std::move(solver.potential));}std::optional<Result>min_cost_flow_polynomial()const{returnmin_cost_flow_polynomial(_balance);}std::optional<Result>min_cost_flow_polynomial(conststd::vector<Cap>&balance)const{assert(int(balance.size())==_n);Capbalance_sum=Cap(0);for(Capvalue:balance)balance_sum+=value;if(balance_sum!=Cap(0))returnstd::nullopt;returnpolynomial_min_cost_flow_impl(balance);}std::optional<Result>min_cost_st_flow(ints,intt,Capflow_value)const{assert(0<=s&&s<_n);assert(0<=t&&t<_n);assert(s!=t);std::vector<Cap>balance=_balance;balance[s]+=flow_value;balance[t]-=flow_value;returnmin_cost_flow(balance);}std::optional<Result>min_cost_st_flow_polynomial(ints,intt,Capflow_value)const{assert(0<=s&&s<_n);assert(0<=t&&t<_n);assert(s!=t);std::vector<Cap>balance=_balance;balance[s]+=flow_value;balance[t]-=flow_value;returnmin_cost_flow_polynomial(balance);}};template<classCap,classCost,classTotalCost=Cost,std::size_tPivotLimitFactor=8>usingBMinCostFlow=BoundedMinCostFlow<Cap,Cost,TotalCost,PivotLimitFactor>;}// namespace flow}// namespace m1une#line 1 "graph/flow/gomory_hu.hpp"
#line 9 "graph/flow/gomory_hu.hpp"
namespacem1une{namespaceflow{template<classCap>structGomoryHu{structEdge{intu;intv;Capcap;};private:structFlowEdge{intto;intrev;Capcap;Capinitial_cap;};int_n;bool_built=false;std::vector<Edge>_edges;std::vector<Edge>_tree_edges;std::vector<int>_parent;std::vector<Cap>_cut_value;std::vector<std::vector<std::pair<int,Cap>>>_tree;std::vector<std::vector<int>>_up;std::vector<std::vector<Cap>>_minimum;std::vector<int>_depth;std::vector<std::vector<FlowEdge>>_graph;std::vector<Cap>_excess;std::vector<int>_height;std::vector<int>_height_count;std::vector<int>_current;std::vector<bool>_active;std::vector<std::vector<int>>_buckets;std::vector<int>_queue;int_highest;longlong_work;longlong_work_limit;voidadd_flow_edge(intu,intv,Capcap){if(u==v||cap==Cap(0))return;intui=int(_graph[u].size());intvi=int(_graph[v].size());_graph[u].push_back(FlowEdge{v,vi,cap,cap});_graph[v].push_back(FlowEdge{u,ui,cap,cap});}voidreset_flow(){for(auto&edges:_graph){for(auto&edge:edges)edge.cap=edge.initial_cap;}}voidactivate(intv,ints,intt){intdead=2*_n;if(v==s||v==t||_active[v]||_excess[v]==Cap(0)||_height[v]>=dead)return;_active[v]=true;_buckets[_height[v]].push_back(v);_highest=std::max(_highest,_height[v]);}voidrebuild_buckets(ints,intt){for(auto&bucket:_buckets)bucket.clear();std::fill(_active.begin(),_active.end(),false);_highest=-1;for(intv=0;v<_n;v++)activate(v,s,t);}voidglobal_relabel(ints,intt){intdead=2*_n;intunreachable=_n+1;std::fill(_height.begin(),_height.end(),unreachable);std::fill(_height_count.begin(),_height_count.end(),0);std::fill(_current.begin(),_current.end(),0);inthead=0;inttail=0;_height[t]=0;_height[s]=_n;_queue[tail++]=t;while(head<tail){intv=_queue[head++];for(constauto&edge:_graph[v]){constFlowEdge&reverse=_graph[edge.to][edge.rev];if(reverse.cap==Cap(0)||_height[edge.to]!=unreachable)continue;_height[edge.to]=_height[v]+1;_queue[tail++]=edge.to;}}for(intv=0;v<_n;v++){_height[v]=std::min(_height[v],dead);_height_count[_height[v]]++;}rebuild_buckets(s,t);_work=0;}voidpush(intv,FlowEdge&edge,ints,intt){if(edge.cap==Cap(0)||_height[v]!=_height[edge.to]+1)return;Capsent=std::min(_excess[v],edge.cap);if(sent==Cap(0))return;boolwas_zero=_excess[edge.to]==Cap(0);edge.cap-=sent;_graph[edge.to][edge.rev].cap+=sent;_excess[v]-=sent;_excess[edge.to]+=sent;if(was_zero)activate(edge.to,s,t);}voidgap(intheight,ints,intt){intunreachable=_n+1;for(intv=0;v<_n;v++){if(v==s||v==t||_height[v]<=height||_height[v]>=_n)continue;_height_count[_height[v]]--;_height[v]=unreachable;_height_count[_height[v]]++;_current[v]=0;}rebuild_buckets(s,t);}boolrelabel(intv,ints,intt){intdead=2*_n;intold_height=_height[v];intnew_height=dead;_work+=int(_graph[v].size());for(constauto&edge:_graph[v]){if(edge.cap!=Cap(0))new_height=std::min(new_height,_height[edge.to]+1);}_height_count[old_height]--;_height[v]=std::min(new_height,dead);_height_count[_height[v]]++;_current[v]=0;if(old_height<_n&&_height_count[old_height]==0){gap(old_height,s,t);returntrue;}returnfalse;}voiddischarge(intv,ints,intt){while(_excess[v]!=Cap(0)&&_height[v]<2*_n){if(_current[v]==int(_graph[v].size())){if(relabel(v,s,t))return;continue;}FlowEdge&edge=_graph[v][_current[v]];_work++;if(edge.cap!=Cap(0)&&_height[v]==_height[edge.to]+1){push(v,edge,s,t);}else{_current[v]++;}}activate(v,s,t);}Capmax_flow(ints,intt){reset_flow();std::fill(_excess.begin(),_excess.end(),Cap(0));for(auto&edge:_graph[s]){Capsent=edge.cap;if(sent==Cap(0))continue;edge.cap=Cap(0);_graph[edge.to][edge.rev].cap+=sent;_excess[edge.to]+=sent;}global_relabel(s,t);while(_highest>=0){if(_buckets[_highest].empty()){_highest--;continue;}intv=_buckets[_highest].back();_buckets[_highest].pop_back();if(!_active[v]||_height[v]!=_highest)continue;_active[v]=false;discharge(v,s,t);if(_work>=_work_limit)global_relabel(s,t);}return_excess[t];}std::vector<bool>source_side(ints){std::vector<bool>visited(_n,false);inthead=0;inttail=0;visited[s]=true;_queue[tail++]=s;while(head<tail){intv=_queue[head++];for(constauto&edge:_graph[v]){if(edge.cap==Cap(0)||visited[edge.to])continue;visited[edge.to]=true;_queue[tail++]=edge.to;}}returnvisited;}voidbuild_query_table(){intlog=1;while((1LL<<log)<=std::max(1,_n))log++;constCapinfinity=std::numeric_limits<Cap>::max();_up.assign(log,std::vector<int>(_n,0));_minimum.assign(log,std::vector<Cap>(_n,infinity));_depth.assign(_n,0);if(_n==0)return;std::vector<int>order;order.reserve(_n);order.push_back(0);for(inti=0;i<int(order.size());i++){intv=order[i];for(auto[to,cap]:_tree[v]){if(to==_up[0][v]&&v!=0)continue;_up[0][to]=v;_minimum[0][to]=cap;_depth[to]=_depth[v]+1;order.push_back(to);}}for(intk=1;k<log;k++){for(intv=0;v<_n;v++){intmiddle=_up[k-1][v];_up[k][v]=_up[k-1][middle];_minimum[k][v]=std::min(_minimum[k-1][v],_minimum[k-1][middle]);}}}public:GomoryHu():GomoryHu(0){}explicitGomoryHu(intn):_n(n){assert(0<=n);}intsize()const{return_n;}intedge_count()const{returnint(_edges.size());}intadd_edge(intu,intv,Capcap){assert(0<=u&&u<_n);assert(0<=v&&v<_n);assert(Cap(0)<=cap);_built=false;intid=int(_edges.size());_edges.push_back(Edge{u,v,cap});returnid;}voidbuild(){std::vector<Edge>flow_edges;flow_edges.reserve(_edges.size());for(autoedge:_edges){if(edge.u==edge.v||edge.cap==Cap(0))continue;if(edge.u>edge.v)std::swap(edge.u,edge.v);flow_edges.push_back(edge);}std::sort(flow_edges.begin(),flow_edges.end(),[](constEdge&lhs,constEdge&rhs){returnstd::pair<int,int>(lhs.u,lhs.v)<std::pair<int,int>(rhs.u,rhs.v);});intunique_edges=0;for(constauto&edge:flow_edges){if(unique_edges>0&&flow_edges[unique_edges-1].u==edge.u&&flow_edges[unique_edges-1].v==edge.v){flow_edges[unique_edges-1].cap+=edge.cap;}else{flow_edges[unique_edges++]=edge;}}flow_edges.resize(unique_edges);_graph.assign(_n,{});std::vector<int>degree(_n,0);for(constauto&edge:flow_edges){degree[edge.u]++;degree[edge.v]++;}for(intv=0;v<_n;v++)_graph[v].reserve(degree[v]);for(constauto&edge:flow_edges)add_flow_edge(edge.u,edge.v,edge.cap);_excess.resize(_n);_height.resize(_n);_height_count.resize(2*_n+1);_current.resize(_n);_active.resize(_n);_buckets.resize(2*_n+1);_queue.resize(_n);longlongarc_count=0;for(constauto&edges:_graph)arc_count+=int(edges.size());_work_limit=std::max(1LL,4*arc_count+_n);_parent.assign(_n,0);_cut_value.assign(_n,std::numeric_limits<Cap>::max());for(ints=1;s<_n;s++){intt=_parent[s];Capflow=max_flow(s,t);std::vector<bool>cut=source_side(s);for(intv=s+1;v<_n;v++){if(_parent[v]==t&&cut[v])_parent[v]=s;}if(cut[_parent[t]]){_parent[s]=_parent[t];_parent[t]=s;_cut_value[s]=_cut_value[t];_cut_value[t]=flow;}else{_cut_value[s]=flow;}}_tree.assign(_n,{});_tree_edges.clear();if(_n>0)_tree_edges.reserve(_n-1);for(intv=1;v<_n;v++){intp=_parent[v];Capcap=_cut_value[v];_tree_edges.push_back(Edge{v,p,cap});_tree[v].emplace_back(p,cap);_tree[p].emplace_back(v,cap);}build_query_table();_built=true;}conststd::vector<Edge>&tree_edges()const{assert(_built);return_tree_edges;}conststd::vector<int>&parent()const{assert(_built);return_parent;}conststd::vector<Cap>&cut_values()const{assert(_built);return_cut_value;}Capmin_cut(intu,intv)const{assert(_built);assert(0<=u&&u<_n);assert(0<=v&&v<_n);assert(u!=v);Capresult=std::numeric_limits<Cap>::max();if(_depth[u]<_depth[v])std::swap(u,v);intdifference=_depth[u]-_depth[v];for(intk=0;difference>0;k++,difference>>=1){if((difference&1)==0)continue;result=std::min(result,_minimum[k][u]);u=_up[k][u];}if(u==v)returnresult;for(intk=int(_up.size())-1;k>=0;k--){if(_up[k][u]==_up[k][v])continue;result=std::min(result,_minimum[k][u]);result=std::min(result,_minimum[k][v]);u=_up[k][u];v=_up[k][v];}result=std::min(result,_minimum[0][u]);result=std::min(result,_minimum[0][v]);returnresult;}};}// namespace flow}// namespace m1une#line 1 "graph/flow/min_cost_flow.hpp"
#line 5 "graph/flow/min_cost_flow.hpp"
#include<array>
#include<bit>
#line 12 "graph/flow/min_cost_flow.hpp"
#include<type_traits>
#line 15 "graph/flow/min_cost_flow.hpp"
#line 18 "graph/flow/min_cost_flow.hpp"
namespacem1une{namespaceflow{template<classCap,classCost>structMinCostFlow{structEdge{intfrom;intto;Capcap;Capflow;Costcost;};private:structInternalEdge{intto;intrev;Capcap;Costcost;};int_n;std::vector<std::pair<int,int>>_pos;std::vector<std::vector<InternalEdge>>_g;bool_has_negative_cost;bool_has_flow;template<classKey>structRadixHeap{usingUnsigned=std::make_unsigned_t<Key>;staticconstexprintbits=std::numeric_limits<Unsigned>::digits;std::array<std::vector<std::pair<Unsigned,int>>,bits+1>bucket;Unsignedlast=0;std::size_tcount=0;staticintindex(Unsignedfirst,Unsignedsecond){returnint(std::bit_width(first^second));}voidclear(){for(auto&values:bucket)values.clear();last=0;count=0;}boolempty()const{returncount==0;}voidpush(Keykey,intvertex){Unsignedvalue=static_cast<Unsigned>(key);assert(last<=value);bucket[index(value,last)].emplace_back(value,vertex);count++;}std::pair<Key,int>pop(){if(bucket[0].empty()){inti=1;while(bucket[i].empty())i++;last=bucket[i][0].first;for(constauto&value:bucket[i]){last=std::min(last,value.first);}for(constauto&value:bucket[i]){bucket[index(value.first,last)].push_back(value);}bucket[i].clear();}auto[key,vertex]=bucket[0].back();bucket[0].pop_back();count--;return{static_cast<Key>(key),vertex};}};template<classKey>structBinaryHeap{usingValue=std::pair<Key,int>;std::vector<Value>heap;voidclear(){heap.clear();}boolempty()const{returnheap.empty();}voidpush(Keykey,intvertex){heap.emplace_back(key,vertex);std::push_heap(heap.begin(),heap.end(),std::greater<Value>());}Valuepop(){std::pop_heap(heap.begin(),heap.end(),std::greater<Value>());Valueresult=heap.back();heap.pop_back();returnresult;}};template<classKey,boolUseRadix=std::numeric_limits<Key>::is_integer&&sizeof(Key)<=8>structHeapSelector{usingType=BinaryHeap<Key>;};template<classKey>structHeapSelector<Key,true>{usingType=RadixHeap<Key>;};booluse_network_simplex(ints,intt,Capflow_limit)const{if(_has_negative_cost)returnfalse;if(_pos.size()<64)returnfalse;autoadd_saturated=[](Capfirst,Capsecond){constCapmaximum=std::numeric_limits<Cap>::max();returnmaximum-first<second?maximum:first+second;};structTerminalCapacity{Captotal=Cap(0);std::array<Cap,7>largest{};};autoadd_capacity=[&](TerminalCapacity&terminal,Capcap){terminal.total=add_saturated(terminal.total,cap);for(Cap¤t:terminal.largest){if(cap<=current)break;std::swap(cap,current);}};TerminalCapacitysource;for(constauto&e:_g[s]){if(e.to==s)continue;add_capacity(source,e.cap);}TerminalCapacitysink;for(constauto&e:_g[t]){if(e.to==t)continue;Capcap=_g[e.to][e.rev].cap;add_capacity(sink,cap);}Captarget=std::min(flow_limit,std::min(source.total,sink.total));if(target==Cap(0))returnfalse;autorequires_eight_arcs=[&](constTerminalCapacity&terminal){Capsum=Cap(0);for(Capcap:terminal.largest){sum=add_saturated(sum,cap);}returnsum<target;};returnrequires_eight_arcs(source)&&requires_eight_arcs(sink);}std::pair<Cap,Cost>network_simplex_flow(ints,intt,Capflow_limit){structResidualArc{intedge;boolreverse;};usingSolver=BoundedMinCostFlow<Cap,Cost,Cost>;std::vector<ResidualArc>arcs;arcs.reserve(2*_pos.size());for(inti=0;i<int(_pos.size());i++){auto[from,idx]=_pos[i];constauto&e=_g[from][idx];constauto&reverse=_g[e.to][e.rev];if(e.cap!=Cap(0)){arcs.push_back(ResidualArc{i,false});}if(reverse.cap!=Cap(0)){arcs.push_back(ResidualArc{i,true});}}autoadd_saturated=[](Capfirst,Capsecond,bool&exact){constCapmaximum=std::numeric_limits<Cap>::max();if(maximum-first<second){exact=false;returnmaximum;}returnfirst+second;};boolsource_capacity_exact=true;Capsource_capacity=Cap(0);for(constauto&e:_g[s]){if(e.to==s)continue;source_capacity=add_saturated(source_capacity,e.cap,source_capacity_exact);}boolsink_capacity_exact=true;Capsink_capacity=Cap(0);for(constauto&e:_g[t]){if(e.to==t)continue;sink_capacity=add_saturated(sink_capacity,_g[e.to][e.rev].cap,sink_capacity_exact);}Captarget=std::min(flow_limit,std::min(source_capacity,sink_capacity));if(target==Cap(0))return{Cap(0),Cost(0)};structArcData{intfrom;intto;Capcap;Costcost;};autoarc_data=[&](constResidualArc&arc){auto[from,idx]=_pos[arc.edge];constauto&e=_g[from][idx];constauto&reverse=_g[e.to][e.rev];returnarc.reverse?ArcData{e.to,from,reverse.cap,reverse.cost}:ArcData{from,e.to,e.cap,e.cost};};autoapply_flow=[&](constResidualArc&arc,Capamount){auto[from,idx]=_pos[arc.edge];auto&e=_g[from][idx];auto&reverse=_g[e.to][e.rev];if(arc.reverse){reverse.cap-=amount;e.cap+=amount;}else{e.cap-=amount;reverse.cap+=amount;}};booltarget_infeasible=false;if(source_capacity_exact&&sink_capacity_exact&&target==source_capacity&&target==sink_capacity){Solverterminal_solver(_n);terminal_solver.reserve_edges(int(arcs.size()));std::vector<Cap>balance(_n,Cap(0));std::vector<int>internal_arcs;std::vector<int>fixed_arcs;internal_arcs.reserve(arcs.size());fixed_arcs.reserve(_g[s].size()+_g[t].size());Costfixed_cost=Cost(0);for(inti=0;i<int(arcs.size());i++){ArcDatadata=arc_data(arcs[i]);if(data.from==s){if(data.to==s)continue;fixed_arcs.push_back(i);fixed_cost+=Cost(data.cap)*data.cost;if(data.to!=t)balance[data.to]+=data.cap;}elseif(data.to==t){if(data.from==t)continue;fixed_arcs.push_back(i);fixed_cost+=Cost(data.cap)*data.cost;balance[data.from]-=data.cap;}elseif(data.to!=s&&data.from!=t){terminal_solver.add_edge(data.from,data.to,Cap(0),data.cap,data.cost);internal_arcs.push_back(i);}}autoterminal_result=terminal_solver.min_cost_flow(balance);if(terminal_result){for(inti:fixed_arcs){apply_flow(arcs[i],arc_data(arcs[i]).cap);}for(inti=0;i<int(internal_arcs.size());i++){apply_flow(arcs[internal_arcs[i]],terminal_result->flow(i));}_has_flow=true;return{target,fixed_cost+terminal_result->cost};}target_infeasible=true;}Solversolver(_n);solver.reserve_edges(int(arcs.size()));for(constauto&arc:arcs){ArcDatadata=arc_data(arc);solver.add_edge(data.from,data.to,Cap(0),data.cap,data.cost);}Capsent=target;std::optional<typenameSolver::Result>result;if(!target_infeasible&&target!=std::numeric_limits<Cap>::max()){result=solver.min_cost_st_flow(s,t,target);}if(!result){MaxFlow<Cap>feasible(_n);feasible.reserve_edges(int(arcs.size()));for(constauto&arc:arcs){auto[from,idx]=_pos[arc.edge];constauto&e=_g[from][idx];constauto&reverse=_g[e.to][e.rev];if(arc.reverse){feasible.add_edge(e.to,from,reverse.cap);}else{feasible.add_edge(from,e.to,e.cap);}}sent=feasible.max_flow(s,t,target);if(sent==Cap(0))return{Cap(0),Cost(0)};result=solver.min_cost_st_flow(s,t,sent);}assert(result.has_value());for(inti=0;i<int(arcs.size());i++){auto[from,idx]=_pos[arcs[i].edge];auto&e=_g[from][idx];auto&reverse=_g[e.to][e.rev];Capamount=result->flow(i);if(arcs[i].reverse){reverse.cap-=amount;e.cap+=amount;}else{e.cap-=amount;reverse.cap+=amount;}}_has_flow=true;return{sent,result->cost};}voidinit_potential(ints,std::vector<Cost>&potential,Costcost_inf)const{if(!_has_negative_cost&&!_has_flow){potential.assign(_n,Cost(0));return;}potential.assign(_n,cost_inf);potential[s]=Cost(0);for(intiter=0;iter<_n-1;iter++){boolupdated=false;for(intv=0;v<_n;v++){if(potential[v]==cost_inf)continue;for(constauto&e:_g[v]){if(e.cap==Cap(0))continue;Costnd=potential[v]+e.cost;if(nd<potential[e.to]){potential[e.to]=nd;updated=true;}}}if(!updated)break;}for(intv=0;v<_n;v++){if(potential[v]==cost_inf)potential[v]=Cost(0);}}public:MinCostFlow():MinCostFlow(0){}explicitMinCostFlow(intn):_n(n),_g(n),_has_negative_cost(false),_has_flow(false){assert(0<=n);}intsize()const{return_n;}intedge_count()const{returnint(_pos.size());}voidreserve_edges(intedge_count){assert(0<=edge_count);_pos.reserve(edge_count);if(_n==0||edge_count==0||2*std::size_t(edge_count)<std::size_t(_n)){return;}conststd::size_taverage_degree=(3*std::size_t(edge_count)+std::size_t(_n)-1)/std::size_t(_n);for(auto&edges:_g)edges.reserve(average_degree);}voidreserve_edges(intedge_count,conststd::vector<int>°rees){assert(0<=edge_count);assert(int(degrees.size())==_n);_pos.reserve(edge_count);for(intv=0;v<_n;v++){assert(0<=degrees[v]);_g[v].reserve(degrees[v]);}}intadd_edge(intfrom,intto,Capcap,Costcost){assert(0<=from&&from<_n);assert(0<=to&&to<_n);assert(Cap(0)<=cap);_has_negative_cost=_has_negative_cost||cost<Cost(0);intid=int(_pos.size());intfrom_id=int(_g[from].size());intto_id=int(_g[to].size());if(from==to)to_id++;_pos.emplace_back(from,from_id);_g[from].push_back(InternalEdge{to,to_id,cap,cost});_g[to].push_back(InternalEdge{from,from_id,Cap(0),-cost});returnid;}Edgeget_edge(inti)const{assert(0<=i&&i<int(_pos.size()));auto[from,idx]=_pos[i];constauto&e=_g[from][idx];constauto&re=_g[e.to][e.rev];returnEdge{from,e.to,e.cap+re.cap,re.cap,e.cost};}std::vector<Edge>edges()const{std::vector<Edge>result;result.reserve(_pos.size());for(inti=0;i<int(_pos.size());i++)result.push_back(get_edge(i));returnresult;}std::pair<Cap,Cost>flow(ints,intt){returnflow(s,t,std::numeric_limits<Cap>::max());}std::pair<Cap,Cost>flow(ints,intt,Capflow_limit){assert(0<=s&&s<_n);assert(0<=t&&t<_n);assert(s!=t);assert(Cap(0)<=flow_limit);if(flow_limit==Cap(0))return{Cap(0),Cost(0)};ifconstexpr(std::numeric_limits<Cap>::is_integer&&std::numeric_limits<Cap>::is_signed&&std::numeric_limits<Cost>::is_signed){if(use_network_simplex(s,t,flow_limit)){returnnetwork_simplex_flow(s,t,flow_limit);}}autoresult=slope(s,t,flow_limit);returnresult.back();}std::vector<std::pair<Cap,Cost>>slope(ints,intt){returnslope(s,t,std::numeric_limits<Cap>::max());}std::vector<std::pair<Cap,Cost>>slope(ints,intt,Capflow_limit){assert(0<=s&&s<_n);assert(0<=t&&t<_n);assert(s!=t);assert(Cap(0)<=flow_limit);constCostcost_inf=std::numeric_limits<Cost>::max()/Cost(4);std::vector<Cost>potential,dist(_n);std::vector<int>prev_v(_n),prev_e(_n);std::vector<int>settled;settled.reserve(_n);typenameHeapSelector<Cost>::Typeque;init_potential(s,potential,cost_inf);std::vector<std::pair<Cap,Cost>>result;result.emplace_back(Cap(0),Cost(0));Capflow=0;Costcost=0;while(flow<flow_limit){std::fill(dist.begin(),dist.end(),cost_inf);dist[s]=Cost(0);settled.clear();que.clear();que.push(Cost(0),s);while(!que.empty()){auto[d,v]=que.pop();if(dist[v]!=d)continue;settled.push_back(v);if(v==t)break;for(inti=0;i<int(_g[v].size());i++){constauto&e=_g[v][i];if(e.cap==Cap(0))continue;Costnd=d+e.cost+potential[v]-potential[e.to];if(nd>=dist[e.to])continue;dist[e.to]=nd;prev_v[e.to]=v;prev_e[e.to]=i;que.push(nd,e.to);}}if(dist[t]==cost_inf)break;for(intv:settled){potential[v]+=dist[v]-dist[t];}Capadd=flow_limit-flow;for(intv=t;v!=s;v=prev_v[v]){add=std::min(add,_g[prev_v[v]][prev_e[v]].cap);}Costpath_cost=potential[t]-potential[s];for(intv=t;v!=s;v=prev_v[v]){auto&e=_g[prev_v[v]][prev_e[v]];e.cap-=add;_g[e.to][e.rev].cap+=add;}flow+=add;cost+=Cost(add)*path_cost;result.emplace_back(flow,cost);}_has_flow=_has_flow||flow!=Cap(0);returnresult;}};}// namespace flow}// namespace m1une#line 9 "graph/flow/flow.hpp"
#line 1 "utilities/fast_io.hpp"
#line 6 "utilities/fast_io.hpp"
#include<cerrno>
#include<charconv>
#line 9 "utilities/fast_io.hpp"
#include<cstdio>
#include<cstdlib>
#include<cstdint>
#include<cstring>
#include<iterator>
#include<string>
#include<sys/stat.h>
#line 18 "utilities/fast_io.hpp"
#include<unistd.h>
#line 20 "utilities/fast_io.hpp"
namespacem1une{namespaceutilities{structFastOutput;namespaceinternal{// Shared with the convenience helpers in template.hpp.inlineFastOutput*standard_output_instance=nullptr;// Detect std::begin(x), std::end(x).template<classT,class=void>structis_range:std::false_type{};template<classT>structis_range<T,std::void_t<decltype(std::begin(std::declval<T&>())),decltype(std::end(std::declval<T&>()))>>:std::true_type{};template<classT>inlineconstexprboolis_range_v=is_range<T>::value;template<classT>usingrange_reference_t=decltype(*std::begin(std::declval<T&>()));template<classT>usingrange_value_t=std::remove_cv_t<std::remove_reference_t<range_reference_t<T>>>;template<classT,class=void>structrange_stored_value{usingtype=range_value_t<T>;};template<classT>structrange_stored_value<T,std::void_t<typenamestd::remove_cv_t<std::remove_reference_t<T>>::value_type>>{usingtype=typenamestd::remove_cv_t<std::remove_reference_t<T>>::value_type;};template<classT>usingrange_stored_value_t=typenamerange_stored_value<T>::type;// Treat strings and C strings as scalar output objects, not as ranges.template<classT>structis_char_array:std::false_type{};template<classT,std::size_tN>structis_char_array<T[N]>:std::bool_constant<std::is_same_v<std::remove_cv_t<T>,char>>{};template<classT>structis_string_like:std::bool_constant<std::is_same_v<std::decay_t<T>,std::string>||std::is_same_v<std::decay_t<T>,constchar*>||std::is_same_v<std::decay_t<T>,char*>||is_char_array<std::remove_reference_t<T>>::value>{};template<classT>inlineconstexprboolis_string_like_v=is_string_like<T>::value;// ModInt-like type: x.val() is printable, and x can be assigned from long long.template<classT,class=void>structhas_val_method:std::false_type{};template<classT>structhas_val_method<T,std::void_t<decltype(std::declval<constT&>().val())>>:std::true_type{};template<classT>inlineconstexprboolhas_val_method_v=has_val_method<T>::value;template<classT,class=void>structhas_static_mod_raw:std::false_type{};template<classT>structhas_static_mod_raw<T,std::void_t<decltype(T::mod()),decltype(T::raw(std::declval<uint32_t>()))>>:std::true_type{};template<classT>inlineconstexprboolhas_static_mod_raw_v=has_static_mod_raw<T>::value;// libstdc++ before GCC 16 does not classify __int128 as an integral type in// strict ISO modes such as -std=c++23. Keep the fast-I/O interface independent// of that implementation detail.template<classT>inlineconstexprboolis_integral_v=std::is_integral_v<T>||std::is_same_v<std::remove_cv_t<T>,__int128_t>||std::is_same_v<std::remove_cv_t<T>,__uint128_t>;template<classT>inlineconstexprboolis_signed_v=std::is_signed_v<T>||std::is_same_v<std::remove_cv_t<T>,__int128_t>;template<classT>structmake_unsigned{usingtype=std::make_unsigned_t<T>;};template<>structmake_unsigned<__int128_t>{usingtype=__uint128_t;};template<>structmake_unsigned<__uint128_t>{usingtype=__uint128_t;};template<classT>usingmake_unsigned_t=typenamemake_unsigned<std::remove_cv_t<T>>::type;}// namespace internalstructFastInput{staticconstexprintbuffer_size=1<<20;private:std::FILE*_stream;char_buffer[buffer_size];int_position;int_length;int_file_descriptor;bool_streaming;boolrefill(){_position=0;if(_streaming){ssize_tlength;do{length=::read(_file_descriptor,_buffer,buffer_size);}while(length<0&&errno==EINTR);if(length<=0){_length=0;returnfalse;}_length=int(length);}else{_length=int(std::fread(_buffer,1,buffer_size,_stream));}return_length!=0;}template<classT>boolread_integer_from_stream(T&value){if(!skip_spaces())returnfalse;intc=read_char_raw();boolnegative=false;if(c=='-'){negative=true;c=read_char_raw();}ifconstexpr(internal::is_signed_v<T>){Tresult=0;while('0'<=c&&c<='9'){result=negative?result*10-(c-'0'):result*10+(c-'0');c=read_char_raw();}value=result;}else{Tresult=0;while('0'<=c&&c<='9'){result=result*10+T(c-'0');c=read_char_raw();}value=negative?T(0)-result:result;}returntrue;}boolprepare_number(){if(_length-_position>=64)returntrue;constintremaining=_length-_position;if(remaining>0)std::memmove(_buffer,_buffer+_position,remaining);constintadded=int(std::fread(_buffer+remaining,1,buffer_size-remaining,_stream));_position=0;_length=remaining+added;if(_length<buffer_size)_buffer[_length]='\0';return_length!=0;}public:explicitFastInput(std::FILE*stream=stdin):_stream(stream),_position(0),_length(0),_file_descriptor(::fileno(stream)),_streaming([&]{structstatstatus;return_file_descriptor>=0&&::fstat(_file_descriptor,&status)==0&&!S_ISREG(status.st_mode);}()){}FastInput(constFastInput&)=delete;FastInput&operator=(constFastInput&)=delete;intread_char_raw(){if(_position==_length&&!refill())returnEOF;return_buffer[_position++];}boolskip_spaces(){intc=read_char_raw();while(c!=EOF&&c<=' ')c=read_char_raw();if(c==EOF)returnfalse;--_position;returntrue;}boolread(char&value){if(!skip_spaces())returnfalse;value=char(read_char_raw());returntrue;}boolread(std::string&value){if(!skip_spaces())returnfalse;value.clear();while(true){constintbegin=_position;while(_position<_length&&static_cast<unsignedchar>(_buffer[_position])>' '){++_position;}value.append(_buffer+begin,_position-begin);if(_position<_length){++_position;returntrue;}if(!refill())returntrue;}}boolread(bool&value){intx;if(!read(x))returnfalse;value=x!=0;returntrue;}template<classT>std::enable_if_t<internal::is_integral_v<T>&&!std::is_same_v<std::remove_cv_t<T>,bool>&&!std::is_same_v<std::remove_cv_t<T>,char>,bool>read(T&value){if(_streaming)returnread_integer_from_stream(value);if(!prepare_number())returnfalse;intc=static_cast<unsignedchar>(_buffer[_position++]);while(c<=' ')c=static_cast<unsignedchar>(_buffer[_position++]);boolnegative=false;if(c=='-'){negative=true;c=static_cast<unsignedchar>(_buffer[_position++]);}ifconstexpr(internal::is_signed_v<T>){Tresult=0;while('0'<=c&&c<='9'){constintfirst=c-'0';constintsecond=static_cast<unsignedchar>(_buffer[_position])-'0';if(0<=second&&second<=9){result=negative?result*100-(first*10+second):result*100+(first*10+second);++_position;}else{result=negative?result*10-first:result*10+first;}c=static_cast<unsignedchar>(_buffer[_position++]);}value=result;}else{Tresult=0;while('0'<=c&&c<='9'){constunsignedfirst=unsigned(c-'0');constintsecond=static_cast<unsignedchar>(_buffer[_position])-'0';if(0<=second&&second<=9){result=result*100+T(first*10+unsigned(second));++_position;}else{result=result*10+T(first);}c=static_cast<unsignedchar>(_buffer[_position++]);}value=negative?T(0)-result:result;}if(_position>_length)_position=_length;returntrue;}template<classT>std::enable_if_t<std::is_floating_point_v<T>,bool>read(T&value){if(!skip_spaces())returnfalse;intc=read_char_raw();boolnegative=false;if(c=='-'||c=='+'){negative=c=='-';c=read_char_raw();}longdoubleresult=0;while('0'<=c&&c<='9'){result=result*10+(c-'0');c=read_char_raw();}if(c=='.'){longdoubleplace=0.1L;c=read_char_raw();while('0'<=c&&c<='9'){result+=(c-'0')*place;place*=0.1L;c=read_char_raw();}}if(c=='e'||c=='E'){c=read_char_raw();boolexponent_negative=false;if(c=='-'||c=='+'){exponent_negative=c=='-';c=read_char_raw();}intexponent=0;while('0'<=c&&c<='9'){exponent=exponent*10+(c-'0');c=read_char_raw();}longdoublescale=1;longdoublepower=10;while(exponent>0){if(exponent&1)scale*=power;power*=power;exponent>>=1;}result=exponent_negative?result/scale:result*scale;}value=static_cast<T>(negative?-result:result);returntrue;}template<classT>std::enable_if_t<internal::has_val_method_v<T>&&!internal::is_integral_v<T>&&!internal::is_range_v<T>,bool>read(T&value){longlongx;if(!read(x))returnfalse;ifconstexpr(internal::has_static_mod_raw_v<T>){if(x>=0&&uint64_t(x)<uint64_t(T::mod())){value=T::raw(uint32_t(x));}else{value=T(x);}}else{value=T(x);}returntrue;}template<classFirst,classSecond>boolread(std::pair<First,Second>&value){if(!read(value.first))returnfalse;returnread(value.second);}template<classRange>std::enable_if_t<internal::is_range_v<Range>&&!internal::is_string_like_v<Range>,bool>read(Range&range){usingStoredValue=internal::range_stored_value_t<Range>;constexprboolnested=internal::is_range_v<StoredValue>&&!internal::is_string_like_v<StoredValue>;for(auto&&value:range){ifconstexpr(std::is_same_v<StoredValue,bool>&&!nested){boolx;if(!read(x))returnfalse;value=x;}else{if(!read(value))returnfalse;}}returntrue;}template<classFirst,classSecond,class...Rest>boolread(First&first,Second&second,Rest&...rest){if(!read(first))returnfalse;returnread(second,rest...);}template<classT>FastInput&operator>>(T&value){if(!read(value))std::abort();return*this;}};structFastOutput{staticconstexprintbuffer_size=1<<20;private:inlinestaticconstautodigit_quads=[]{std::array<char,40000>result{};for(inti=0;i<10000;i++){intvalue=i;for(intj=3;j>=0;j--){result[4*i+j]=char('0'+value%10);value/=10;}}returnresult;}();std::FILE*_stream;char_buffer[buffer_size];int_position;int_precision;std::chars_format_float_format;char_range_separator;std::string*_capture=nullptr;template<classT>std::stringformat_cell(constT&value){std::stringresult;structCaptureGuard{std::string*⌖std::string*previous;~CaptureGuard(){target=previous;}}guard{_capture,_capture};_capture=&result;write(value);returnresult;}template<classMatrix>voidwrite_aligned_matrix(constMatrix&matrix){std::vector<std::vector<std::string>>rows;std::vector<std::size_t>widths;for(constauto&row:matrix){auto&cells=rows.emplace_back();std::size_tcolumn=0;for(constauto&value:row){cells.push_back(format_cell(value));if(column==widths.size())widths.push_back(0);widths[column]=std::max(widths[column],cells.back().size());++column;}}boolfirst=true;for(constauto&row:rows){if(!first)write_char('\n');first=false;for(std::size_tcolumn=0;column<row.size();++column){if(column!=0)write_char(_range_separator);for(std::size_tpadding=row[column].size();padding<widths[column];++padding){write_char(' ');}write(row[column]);}}}public:explicitFastOutput(std::FILE*stream=stdout):_stream(stream),_position(0),_precision(6),_float_format(std::chars_format::general),_range_separator(' '){if(_stream==stdout&&internal::standard_output_instance==nullptr){internal::standard_output_instance=this;}}FastOutput(constFastOutput&)=delete;FastOutput&operator=(constFastOutput&)=delete;~FastOutput(){flush();if(internal::standard_output_instance==this){internal::standard_output_instance=nullptr;}}voidflush(){if(_position!=0){std::fwrite(_buffer,1,_position,_stream);_position=0;}std::fflush(_stream);}voidwrite_char(charc){if(_capture!=nullptr){_capture->push_back(c);return;}if(_position==buffer_size)flush();_buffer[_position++]=c;}voidwrite(constchar*s){while(*s!='\0')write_char(*s++);}voidwrite(conststd::string&s){if(_capture!=nullptr){_capture->append(s);return;}std::size_tposition=0;while(position<s.size()){if(_position==buffer_size)flush();conststd::size_tcopied=std::min<std::size_t>(buffer_size-_position,s.size()-position);std::memcpy(_buffer+_position,s.data()+position,copied);_position+=int(copied);position+=copied;}}voidwrite(charc){write_char(c);}voidwrite(boolvalue){write_char(value?'1':'0');}template<classT>std::enable_if_t<std::is_floating_point_v<T>>write(Tvalue){chardigits[128];auto[end,error]=std::to_chars(digits,digits+sizeof(digits),value,_float_format,_precision);if(error!=std::errc())std::abort();for(constchar*pointer=digits;pointer!=end;pointer++){write_char(*pointer);}}template<classT>std::enable_if_t<internal::is_integral_v<T>&&!std::is_same_v<std::remove_cv_t<T>,bool>&&!std::is_same_v<std::remove_cv_t<T>,char>>write(Tvalue){usingRaw=std::remove_cv_t<T>;usingUnsigned=internal::make_unsigned_t<Raw>;Unsignedmagnitude;ifconstexpr(internal::is_signed_v<Raw>){if(value<0){write_char('-');magnitude=Unsigned(0)-Unsigned(value);}else{magnitude=Unsigned(value);}}else{magnitude=value;}if(magnitude==0){write_char('0');return;}unsignedchunks[16];intcount=0;while(magnitude>=10000){constUnsignedquotient=magnitude/10000;chunks[count++]=unsigned(magnitude-quotient*10000);magnitude=quotient;}if(_capture==nullptr&&_position>buffer_size-64)flush();charcaptured[64];char*constbegin=_capture!=nullptr?captured:_buffer+_position;char*destination=begin;constunsignedleading=unsigned(magnitude);constchar*first=digit_quads.data()+4*leading;intskip=leading<10?3:leading<100?2:leading<1000?1:0;for(;skip<4;skip++)*destination++=first[skip];while(count--){constchar*digits=digit_quads.data()+4*chunks[count];std::memcpy(destination,digits,4);destination+=4;}if(_capture!=nullptr){_capture->append(begin,destination-begin);}else{_position+=int(destination-begin);}}template<classT>std::enable_if_t<internal::has_val_method_v<T>&&!internal::is_integral_v<T>&&!internal::is_range_v<T>>write(constT&value){write(value.val());}template<classFirst,classSecond>voidwrite(conststd::pair<First,Second>&value){write(value.first);write_char(' ');write(value.second);}template<classRange>std::enable_if_t<internal::is_range_v<Range>&&!internal::is_string_like_v<Range>>write(constRange&range){usingStoredValue=internal::range_stored_value_t<constRange>;constexprboolnested=internal::is_range_v<StoredValue>&&!internal::is_string_like_v<StoredValue>;boolfirst=true;for(constauto&value:range){if(!first)write_char(nested?'\n':_range_separator);first=false;ifconstexpr(std::is_same_v<StoredValue,bool>&&!nested){write(static_cast<bool>(value));}else{write(value);}}}template<classFirst,class...Rest>voidprint(constFirst&first,constRest&...rest){write(first);((write_char(' '),write(rest)),...);}voidprintln(){write_char('\n');}voidset_precision(intprecision){_precision=precision;}voidset_fixed(intprecision=6){_float_format=std::chars_format::fixed;_precision=precision;}voidset_general(intprecision=6){_float_format=std::chars_format::general;_precision=precision;}voidset_range_separator(charseparator){_range_separator=separator;}template<classMatrix>voidwrite_aligned(constMatrix&matrix){usingRow=internal::range_stored_value_t<constMatrix>;usingCell=internal::range_stored_value_t<constRow>;static_assert(internal::is_range_v<Row>&&!internal::is_string_like_v<Row>,"write_aligned requires a two-dimensional range");static_assert(!internal::is_range_v<Cell>||internal::is_string_like_v<Cell>,"write_aligned requires scalar cells");write_aligned_matrix(matrix);}template<classMatrix>voidprintln_aligned(constMatrix&matrix){write_aligned(matrix);write_char('\n');}template<class...Args>voidprintln(constArgs&...args){print(args...);write_char('\n');}template<classT>FastOutput&operator<<(constT&value){write(value);return*this;}};}// namespace utilities}// namespace m1une#line 11 "verify/graph/flow/flow_algorithms.test.cpp"
voidtest_max_flow(){m1une::flow::MaxFlow<longlong>mf(4);inte0=mf.add_edge(0,1,2);inte1=mf.add_edge(0,2,1);inte2=mf.add_edge(1,2,1);inte3=mf.add_edge(1,3,1);inte4=mf.add_edge(2,3,2);(void)e1;(void)e2;(void)e3;(void)e4;assert(mf.size()==4);assert(mf.edge_count()==5);assert(mf.max_flow(0,3)==3);autoedges=mf.edges();longlongoutgoing=0;for(constauto&e:edges){if(e.from==0)outgoing+=e.flow;assert(0<=e.flow&&e.flow<=e.cap);}assert(outgoing==3);assert(mf.get_edge(e0).cap==2);autocut=mf.min_cut(0);assert(cut[0]);assert(!cut[3]);mf.change_edge(e0,3,1);autochanged=mf.get_edge(e0);assert(changed.cap==3);assert(changed.flow==1);m1une::flow::MaxFlow<longlong>undirected(2);undirected.reserve_edges(1,std::vector<int>{1,1});intundirected_id=undirected.add_undirected_edge(0,1,7);undirected.change_edge(undirected_id,7,-3);autoinitial_undirected=undirected.get_edge(undirected_id);assert(initial_undirected.cap==7);assert(initial_undirected.flow==-3);assert(undirected.max_flow(0,1)==10);assert(undirected.get_edge(undirected_id).flow==7);m1une::flow::MaxFlow<longlong>limited(2);intlimited_id=limited.add_edge(0,1,10);assert(limited.max_flow(0,1,4)==4);assert(limited.get_edge(limited_id).flow==4);assert(limited.max_flow(0,1)==6);assert(limited.get_edge(limited_id).flow==10);// Different path lengths exercise repeated highest-label changes.m1une::flow::MaxFlow<longlong>layered(17);intnext_vertex=1;for(intlength=2;length<=6;length++){intfrom=0;for(intedge=0;edge<length;edge++){intto=edge+1==length?16:next_vertex++;layered.add_edge(from,to,3);from=to;}}assert(next_vertex==16);assert(layered.max_flow(0,16)==15);assert(layered.max_flow(0,16)==0);structInputEdge{intfrom;intto;longlongcap;boolundirected;};std::mt19937random(19260817);for(intiteration=0;iteration<500;iteration++){intn=2+int(random()%6);intm=int(random()%13);std::vector<InputEdge>input_edges;m1une::flow::MaxFlow<longlong>flow(n);m1une::flow::MaxFlow<longlong>push_relabel_flow(n);m1une::flow::MaxFlow<longlong>dinic_flow(n);flow.reserve_edges(m);push_relabel_flow.reserve_edges(m);dinic_flow.reserve_edges(m);for(intedge=0;edge<m;edge++){InputEdgeinput{int(random()%n),int(random()%n),1+static_cast<longlong>(random()%10),bool(random()&1)};input_edges.push_back(input);if(input.undirected){flow.add_undirected_edge(input.from,input.to,input.cap);push_relabel_flow.add_undirected_edge(input.from,input.to,input.cap);dinic_flow.add_undirected_edge(input.from,input.to,input.cap);}else{flow.add_edge(input.from,input.to,input.cap);push_relabel_flow.add_edge(input.from,input.to,input.cap);dinic_flow.add_edge(input.from,input.to,input.cap);}}longlongexpected=std::numeric_limits<longlong>::max();for(intmask=0;mask<(1<<n);mask++){if((mask&1)==0||(mask>>(n-1)&1)!=0)continue;longlongcapacity=0;for(constauto&edge:input_edges){boolfrom_side=mask>>edge.from&1;boolto_side=mask>>edge.to&1;if(from_side&&!to_side)capacity+=edge.cap;if(edge.undirected&&!from_side&&to_side){capacity+=edge.cap;}}expected=std::min(expected,capacity);}longlongresult=flow.max_flow(0,n-1);assert(result==expected);longlongpush_relabel_result=push_relabel_flow.max_flow_push_relabel(0,n-1);assert(push_relabel_result==expected);longlongdinic_result=dinic_flow.max_flow_dinic(0,n-1);assert(dinic_result==expected);autovalidate=[&](constauto&solved_flow,longlongsolved_value){std::vector<longlong>net_flow(n,0);for(intedge=0;edge<m;edge++){autoresult_edge=solved_flow.get_edge(edge);assert(result_edge.cap==input_edges[edge].cap);if(input_edges[edge].undirected){assert(-result_edge.cap<=result_edge.flow);}else{assert(0<=result_edge.flow);}assert(result_edge.flow<=result_edge.cap);net_flow[result_edge.from]+=result_edge.flow;net_flow[result_edge.to]-=result_edge.flow;}assert(net_flow[0]==solved_value);assert(net_flow[n-1]==-solved_value);for(intvertex=1;vertex+1<n;vertex++){assert(net_flow[vertex]==0);}};validate(flow,result);validate(push_relabel_flow,push_relabel_result);validate(dinic_flow,dinic_result);assert(flow.max_flow_push_relabel(0,n-1)==0);assert(push_relabel_flow.max_flow_push_relabel(0,n-1)==0);assert(push_relabel_flow.max_flow(0,n-1)==0);assert(dinic_flow.max_flow_dinic(0,n-1)==0);}}voidtest_gomory_hu(){m1une::flow::GomoryHu<longlong>gh(4);gh.add_edge(0,1,3);gh.add_edge(1,2,2);gh.add_edge(0,2,1);gh.add_edge(2,3,4);gh.build();assert(gh.size()==4);assert(gh.edge_count()==4);assert(gh.tree_edges().size()==3);assert(gh.min_cut(0,1)==4);assert(gh.min_cut(0,2)==3);assert(gh.min_cut(0,3)==3);assert(gh.min_cut(2,3)==4);m1une::flow::GomoryHu<longlong>disconnected(3);disconnected.add_edge(0,1,5);disconnected.build();assert(disconnected.min_cut(0,1)==5);assert(disconnected.min_cut(0,2)==0);m1une::flow::GomoryHu<longlong>rebuilt(2);rebuilt.build();assert(rebuilt.min_cut(0,1)==0);rebuilt.add_edge(0,1,7);rebuilt.add_edge(0,0,100);rebuilt.build();assert(rebuilt.min_cut(0,1)==7);m1une::flow::GomoryHu<longlong>singleton(1);singleton.build();assert(singleton.tree_edges().empty());std::mt19937random(123456789);for(intiteration=0;iteration<200;iteration++){intn=2+int(random()%8);structInputEdge{intu;intv;longlongcap;};std::vector<InputEdge>edges;m1une::flow::GomoryHu<longlong>tree(n);intm=random()%(2*n*n+1);for(inti=0;i<m;i++){intu=random()%n;intv=random()%n;longlongcap=random()%1000001;edges.push_back(InputEdge{u,v,cap});tree.add_edge(u,v,cap);}tree.build();for(ints=0;s<n;s++){for(intt=s+1;t<n;t++){m1une::flow::MaxFlow<longlong>mf(n);for(constauto&edge:edges){mf.add_edge(edge.u,edge.v,edge.cap);mf.add_edge(edge.v,edge.u,edge.cap);}assert(tree.min_cut(s,t)==mf.max_flow(s,t));}}}}voidtest_bounded_flow(){m1une::flow::BoundedFlow<longlong>st(4);inta=st.add_edge(0,1,1,3);intb=st.add_edge(0,2,0,2);intc=st.add_edge(1,3,1,2);intd=st.add_edge(2,3,0,2);inte=st.add_edge(1,2,0,1);(void)a;(void)b;(void)c;(void)d;(void)e;autoexact=st.feasible_st_flow(0,3,3);assert(exact.has_value());std::vector<longlong>balance(4,0);for(constauto&edge:exact->edges){assert(edge.lower<=edge.flow&&edge.flow<=edge.upper);balance[edge.from]+=edge.flow;balance[edge.to]-=edge.flow;}assert((balance==std::vector<longlong>{3,0,0,-3}));autotoo_much=st.feasible_st_flow(0,3,6);assert(!too_much.has_value());m1une::flow::BoundedFlow<longlong>bf(3);intf01=bf.add_edge(0,1,1,3);intf02=bf.add_edge(0,2,0,4);bf.add_edge(1,2,0,2);bf.add_supply(0,4);bf.add_demand(1,1);bf.add_demand(2,3);assert(bf.balance(0)==4);autobflow=bf.feasible_flow();assert(bflow.has_value());assert(bflow->get_edge(f01).flow>=1);assert(bflow->get_edge(f02).flow>=0);std::vector<longlong>b_balance(3,0);for(constauto&edge:bflow->edges){assert(edge.lower<=edge.flow&&edge.flow<=edge.upper);b_balance[edge.from]+=edge.flow;b_balance[edge.to]-=edge.flow;}assert((b_balance==std::vector<longlong>{4,-1,-3}));m1une::flow::BoundedFlow<longlong>negative(2);intneg=negative.add_edge(0,1,-5,5);negative.add_demand(0,3);negative.add_supply(1,3);autonegative_flow=negative.feasible_flow();assert(negative_flow.has_value());assert(negative_flow->flow(neg)==-3);m1une::flow::BoundedFlow<longlong>impossible(2);impossible.add_edge(0,1,0,1);impossible.add_supply(0,2);impossible.add_demand(1,2);assert(!impossible.feasible_flow().has_value());m1une::flow::BFlow<longlong>alias(1);assert(alias.size()==1);}voidtest_bounded_min_cost_flow(){m1une::flow::BoundedMinCostFlow<longlong,longlong>st(3);st.reserve_edges(3);inte01=st.add_edge(0,1,1,3,2);inte12=st.add_edge(1,2,1,3,1);inte02=st.add_edge(0,2,0,3,10);autoexact=st.min_cost_st_flow(0,2,3);assert(exact.has_value());assert(exact->cost==9);assert(exact->flow(e01)==3);assert(exact->flow(e12)==3);assert(exact->flow(e02)==0);autoexact_polynomial=st.min_cost_st_flow_polynomial(0,2,3);assert(exact_polynomial.has_value());assert(exact_polynomial->cost==9);std::vector<longlong>balance(3,0);for(constauto&edge:exact->edges){assert(edge.lower<=edge.flow&&edge.flow<=edge.upper);balance[edge.from]+=edge.flow;balance[edge.to]-=edge.flow;}assert((balance==std::vector<longlong>{3,0,-3}));m1une::flow::BoundedMinCostFlow<longlong,longlong>bf(3);intp01=bf.add_edge(0,1,0,2,1);intp12=bf.add_edge(1,2,0,2,1);intp02=bf.add_edge(0,2,0,2,5);bf.add_supply(0,2);bf.add_demand(2,2);autobflow=bf.min_cost_flow();assert(bflow.has_value());assert(bflow->cost==4);assert(bflow->flow(p01)==2);assert(bflow->flow(p12)==2);assert(bflow->flow(p02)==0);m1une::flow::BoundedMinCostFlow<longlong,longlong>negative(2);intneg=negative.add_edge(0,1,-5,5,2);negative.add_demand(0,3);negative.add_supply(1,3);autonegative_flow=negative.min_cost_flow();assert(negative_flow.has_value());assert(negative_flow->flow(neg)==-3);assert(negative_flow->cost==-6);m1une::flow::BoundedMinCostFlow<longlong,longlong>cycle(2);intc01=cycle.add_edge(0,1,0,1,-5);intc10=cycle.add_edge(1,0,0,1,3);autocirculation=cycle.min_cost_flow();assert(circulation.has_value());assert(circulation->flow(c01)==1);assert(circulation->flow(c10)==1);assert(circulation->cost==-2);autopolynomial_circulation=cycle.min_cost_flow_polynomial();assert(polynomial_circulation.has_value());assert(polynomial_circulation->cost==-2);usingImmediateFallback=m1une::flow::BoundedMinCostFlow<longlong,longlong,longlong,0>;ImmediateFallbackfallback_cycle(2);fallback_cycle.add_edge(0,1,0,1,-5);fallback_cycle.add_edge(1,0,0,1,3);autofallback_circulation=fallback_cycle.min_cost_flow();assert(fallback_circulation.has_value());assert(fallback_circulation->cost==-2);m1une::flow::BoundedMinCostFlow<longlong,longlong>impossible(2);impossible.add_edge(0,1,0,1,0);impossible.add_supply(0,2);impossible.add_demand(1,2);assert(!impossible.min_cost_flow().has_value());usingWideCostFlow=m1une::flow::BoundedMinCostFlow<longlong,longlong,__int128_t>;WideCostFlowwide_cost(1);wide_cost.reserve_edges(1);constexprlonglongtrillion=1000000000000LL;wide_cost.add_edge(0,0,trillion,trillion,trillion);autowide_result=wide_cost.min_cost_flow();assert(wide_result.has_value());assert(wide_result->cost==__int128_t(trillion)*trillion);std::mt19937random(987654321);for(intiteration=0;iteration<500;iteration++){intn=1+int(random()%4);intm=int(random()%7);structSmallEdge{intfrom;intto;longlonglower;longlongupper;longlongcost;};std::vector<SmallEdge>edges;m1une::flow::BoundedMinCostFlow<longlong,longlong>solver(n);for(inti=0;i<m;i++){intfrom=int(random()%n);intto=int(random()%n);longlonglower=static_cast<longlong>(random()%5)-2;longlongupper=lower+static_cast<longlong>(random()%4);longlongcost=static_cast<longlong>(random()%9)-4;edges.push_back(SmallEdge{from,to,lower,upper,cost});solver.add_edge(from,to,lower,upper,cost);}std::vector<longlong>required_balance(n,0);longlongbalance_sum=0;for(intvertex=0;vertex+1<n;vertex++){required_balance[vertex]=static_cast<longlong>(random()%7)-3;balance_sum+=required_balance[vertex];}required_balance.back()=-balance_sum;boolfeasible=false;longlongbest_cost=0;std::vector<longlong>flow(m);autoenumerate=[&](auto&&self,intedge_id)->void{if(edge_id!=m){for(flow[edge_id]=edges[edge_id].lower;flow[edge_id]<=edges[edge_id].upper;flow[edge_id]++){self(self,edge_id+1);}return;}std::vector<longlong>actual_balance(n,0);longlongcost=0;for(inti=0;i<m;i++){actual_balance[edges[i].from]+=flow[i];actual_balance[edges[i].to]-=flow[i];cost+=flow[i]*edges[i].cost;}if(actual_balance!=required_balance)return;if(!feasible||cost<best_cost)best_cost=cost;feasible=true;};enumerate(enumerate,0);autovalidate_result=[&](constauto&result){assert(result.has_value()==feasible);if(!result.has_value())return;assert(result->cost==best_cost);std::vector<longlong>actual_balance(n,0);for(constauto&edge:result->edges){assert(edge.lower<=edge.flow&&edge.flow<=edge.upper);actual_balance[edge.from]+=edge.flow;actual_balance[edge.to]-=edge.flow;longlongreduced_cost=edge.cost+result->potential[edge.from]-result->potential[edge.to];if(edge.flow<edge.upper)assert(0<=reduced_cost);if(edge.lower<edge.flow)assert(reduced_cost<=0);}assert(actual_balance==required_balance);};validate_result(solver.min_cost_flow(required_balance));validate_result(solver.min_cost_flow_polynomial(required_balance));}m1une::flow::BMinCostFlow<longlong,longlong>alias(1);assert(alias.size()==1);m1une::flow::BMinCostFlow<longlong,longlong,__int128_t>wide_alias(1);assert(wide_alias.size()==1);m1une::flow::BMinCostFlow<longlong,longlong,longlong,0>fallback_alias(1);assert(fallback_alias.size()==1);}voidtest_min_cost_flow(){m1une::flow::MinCostFlow<longlong,longlong>mcf(4);mcf.add_edge(0,1,2,1);mcf.add_edge(0,2,1,2);mcf.add_edge(1,2,1,0);mcf.add_edge(1,3,1,3);mcf.add_edge(2,3,2,1);autoresult=mcf.flow(0,3,2);assert(result.first==2);assert(result.second==5);autoedges=mcf.edges();longlongtotal_source_flow=0;for(constauto&e:edges){if(e.from==0)total_source_flow+=e.flow;assert(0<=e.flow&&e.flow<=e.cap);}assert(total_source_flow==2);m1une::flow::MinCostFlow<longlong,longlong>negative(3);negative.add_edge(0,1,1,-5);negative.add_edge(1,2,1,2);negative.add_edge(0,2,1,10);autoslope=negative.slope(0,2,2);std::vector<std::pair<longlong,longlong>>expected_slope={std::pair<longlong,longlong>{0,0},std::pair<longlong,longlong>{1,-3},std::pair<longlong,longlong>{2,7},};assert(slope==expected_slope);std::mt19937random(31415926);for(intiteration=0;iteration<200;iteration++){intn=2+int(random()%7);intm=int(random()%30);longlongflow_limit=random()%16;m1une::flow::MinCostFlow<longlong,longlong>tested(n);m1une::flow::BoundedMinCostFlow<longlong,longlong>expected(n);m1une::flow::MaxFlow<longlong>maximum(n);tested.reserve_edges(m);expected.reserve_edges(m);maximum.reserve_edges(m);for(intedge=0;edge<m;edge++){intfrom=int(random()%(n-1));intto=from+1+int(random()%(n-from-1));longlongcap=random()%6;longlongcost=int(random()%21)-10;tested.add_edge(from,to,cap,cost);expected.add_edge(from,to,0,cap,cost);maximum.add_edge(from,to,cap);}longlongsent=maximum.max_flow(0,n-1,flow_limit);autoexpected_result=expected.min_cost_st_flow(0,n-1,sent);assert(expected_result.has_value());autotested_result=tested.flow(0,n-1,flow_limit);assert(tested_result.first==sent);assert(tested_result.second==expected_result->cost);}// At least eight terminal arcs are required, exercising the guarded// one-shot network-simplex path rather than successive shortest paths.m1une::flow::MinCostFlow<longlong,longlong>fast(10);m1une::flow::MinCostFlow<longlong,longlong>reference(10);fast.reserve_edges(64);reference.reserve_edges(64);for(intv=1;v<=8;v++){fast.add_edge(0,v,1,v);fast.add_edge(v,9,1,10-v);reference.add_edge(0,v,1,v);reference.add_edge(v,9,1,10-v);}while(fast.edge_count()<64){intfrom=1+fast.edge_count()%8;intto=1+(fast.edge_count()*3+1)%8;if(from==to)to=to==8?1:to+1;fast.add_edge(from,to,1,100);reference.add_edge(from,to,1,100);}autofast_result=fast.flow(0,9,10);autoreference_result=reference.slope(0,9,10).back();assert(fast_result==reference_result);std::pair<longlong,longlong>expected_fast_result(8,80);assert(fast_result==expected_fast_result);assert(fast.flow(0,9).first==0);// After one shortest-path augmentation, a limit below both remaining// terminal capacities exercises the exact residual s-t solve without// terminal contraction, including negative-cost reverse arcs.m1une::flow::MinCostFlow<longlong,longlong>partial(12);m1une::flow::MinCostFlow<longlong,longlong>partial_reference(12);for(intv=1;v<=10;v++){partial.add_edge(0,v,1,v);partial.add_edge(v,11,1,0);partial_reference.add_edge(0,v,1,v);partial_reference.add_edge(v,11,1,0);}while(partial.edge_count()<64){intfrom=1+partial.edge_count()%10;intto=1+(partial.edge_count()*7+1)%10;if(from==to)to=to==10?1:to+1;partial.add_edge(from,to,1,100);partial_reference.add_edge(from,to,1,100);}assert(partial.flow(0,11,1)==partial_reference.slope(0,11,1).back());assert(partial.flow(0,11,8)==partial_reference.slope(0,11,8).back());assert(partial.flow(0,11)==partial_reference.slope(0,11).back());// The terminal cut permits eight units but an internal bottleneck permits// only four, exercising the max-flow fallback before the exact-cost solve.m1une::flow::MinCostFlow<longlong,longlong>bottleneck(20);for(intv=1;v<=8;v++){bottleneck.add_edge(0,v,1,0);bottleneck.add_edge(v,9,1,0);bottleneck.add_edge(10,v+10,1,0);bottleneck.add_edge(v+10,19,1,0);}bottleneck.add_edge(9,10,4,1);while(bottleneck.edge_count()<64){intfrom=1+bottleneck.edge_count()%8;intto=1+(bottleneck.edge_count()*5+1)%8;if(from==to)to=to==8?1:to+1;bottleneck.add_edge(from,to,1,100);}std::pair<longlong,longlong>expected_bottleneck(4,4);assert(bottleneck.flow(0,19,8)==expected_bottleneck);// Randomized comparisons cover the terminal-contraction optimization on// feasible non-negative-cost instances.for(intiteration=0;iteration<100;iteration++){constexprintn=18;m1une::flow::MinCostFlow<longlong,longlong>adaptive(n);m1une::flow::MinCostFlow<longlong,longlong>shortest_paths(n);adaptive.reserve_edges(80);shortest_paths.reserve_edges(80);for(intv=1;v<=8;v++){longlongfirst_cost=random()%11;longlongmiddle_cost=random()%11;longlonglast_cost=random()%11;adaptive.add_edge(0,v,1,first_cost);adaptive.add_edge(v,v+8,1,middle_cost);adaptive.add_edge(v+8,17,1,last_cost);shortest_paths.add_edge(0,v,1,first_cost);shortest_paths.add_edge(v,v+8,1,middle_cost);shortest_paths.add_edge(v+8,17,1,last_cost);}while(adaptive.edge_count()<80){intfrom=1+int(random()%16);intto=1+int(random()%16);if(from==to)to=to==16?1:to+1;longlongcap=1+random()%3;longlongcost=random()%21;adaptive.add_edge(from,to,cap,cost);shortest_paths.add_edge(from,to,cap,cost);}assert(adaptive.flow(0,17,8)==shortest_paths.slope(0,17,8).back());}}intmain(){m1une::utilities::FastInputfast_input;m1une::utilities::FastOutputfast_output;test_max_flow();test_gomory_hu();test_bounded_flow();test_bounded_min_cost_flow();test_min_cost_flow();longlonga,b;fast_input>>a>>b;fast_output<<a+b<<'\n';}
verify/graph/flow/flow_algorithms.test.cpp