Overview
In this post I will explain with an example one way to code the Edmonds Karp max flow algorithm to solve a problem.
The Problem
The problem we study is Array and Operations http://mirror.codeforces.com/problemset/problem/498/C. The basic idea is you have a sequence of numbers and a bunch of special pairs taken from the sequence. For any of the pairs you can divide each number by a common factor greater than 1 and replace the elements in the pair with the result of this division. The question asks what is the maximum number of times you can do this.
Converting this to a max flow problem
This section is specific to the problem. You can skip it if you just want to know how to implement Edmonds Karp on a given flow network.
Before we convert this problem to a max flow problem we make some observations. First if we have a pair of numbers the maximum amount of times we can perform the operation described on this pair is equal to the number of prime divisors the pair have in common counting multiplicity (i.e if the same prime divides the pair multiple times we count it multiple times). Next we note this is the same as the number of prime divisors in the greatest common divisor of the pair counting multiplicity.
Next assume we have one number that is paired with multiple other numbers. Then to find the maximum number of times we can perform the operation on this number we first find the gcd between it and it first pairing. Then we count the number of prime divisors in this gcd. Next we replace the original number by the original number divided by the gcd found and repeat the process with the second pairing. We could instead find the gcd of the number and the product of all the numbers paired with it and then count the prime divisors of this gcd to get the same result. However this approach may lead to overflow issues.
Now for any given number we know how to count maximum the number of times we can perform the operation on a pair and the maximum number of times we can perform this operation on the number in total. If we think of the flow as the maximum number of times we can perform the operation we can construct a flow network as follows. We let each number that appears in special pair have two nodes associated with it. Call them the left node and the right node. The left node has a directed edge running from the source to it with a capacity equal to the maximum number of times the operation can be performed on this number in total. The right node has a directed edge running from it to the the sink with the same capacity. Next for each number in a special pair we connect its left node to the right node of the number with which it is paired with a directed edge from left to right with capacity equal to the number of times we can perform the operation on this pair. This gives us our flow network. We also mention it is different from the flow network given in the editorial but it also works.
Of course to do all this we need functions to calculate the number of prime divisors and greatest common divisors. For calculating the greatest common divisor the Euclidean algorithm is fairly easy to implement. Since the numbers given can go up to 109 calculating the number of prime divisor is harder. One is to calculate all primes up to 1000 by brute force. Then we can those as a sieve to get all primes up to 32000. Finally since 320002 > 109 with all these primes we can find the number of prime divisors for any number up to 109.
Set up three graphs
Now that we have our network and edge capacities we create three graphs. We store these in adjacency matrix form we can use a 2 dimensional array or a vector<vector> here. The first graph is the flow where all entries are set to 0. The second graph is the residual where all the edge values are set to the capacity. Finally the third graph is the network graph where we have a 1 if two edges are connected and 0 otherwise.
Find an augmenting path
Next we create a function to find an augmenting path that returns pair<int,vector> where the first term is the capacity and the second is the path. To do this we perform a breadth first search starting from the source. We only need to store the predecessor of each vertex and not the distance. With the predecessors it is easy to find a path by back tracking from the sink. If there is no path we just return an empty pair and deal with it later when we augment along the path. Next we find the capacity of the path by taking the minimum of all capacities of edges included in the path.
Augmenting along the path
Next we create a function that takes the path we found before and update both the flow and the residual. It is useful to let this function return a bool which if true if we made an update and false otherwise. First we call our previous function to find a path. If the path is an empty pair we return false. Otherwise take an edge from point x to point y in the path. For the residual we decrease the capacity of the edge from x to y by the capacity of the path and increase the capacity of the edgy from y to x by the capacity of the path. For the flow we if the edge from x to y is in the network graph we increase the flow from x to y by the capacity of the path otherwise we decrease the capacity of the flow from y to x by the capacity of the path. We do this for all edges in the path then return true.
Finding the maximal flow
To find the maximal flow we just keep calling the augment function until it returns false.
Find the answer
To find the answer we first find the total flow by adding all the flows leaving the source. We could also have stored total flow as we were augmenting. Finally to get the answer we divide the total flow by since our network counts the application of each operation twice.
Code
#include<algorithm>
#include<vector>
#include<list>
#include<iostream>
#include<utility>
#include<string>
#include<set>
#include<queue>
#include<stack>
using namespace std;
//ps stores all primes below 32000
vector<long long> ps;
//graph is the network graph, fl is the flow
//res is the residual
vector<vector<int>> graph,fl,res;
//find an augmenting path
pair<int, vector<int>> getPath(){
pair<int, vector<int>> ans;
int sz=fl.size();
//find predecessors using a breadth first search
vector<int> preds(sz,-1);
preds[0]=0;
queue<int> q;
q.push(0);
while(!q.empty()){
int cur=q.front();
q.pop();
for(int i=0; i<sz; i++){
if(res[cur][i] && preds[i]==-1){
preds[i]=cur;
q.push(i);
}
}
}
//if there is no path to the sink
//return an empty answer
if(preds[sz-1]==-1) return ans;
//create an augmenting path
//from predecessors
int pos=sz-1;
vector<int> vi;
while(pos){
vi.push_back(pos);
pos=preds[pos];
}
vi.push_back(0);
reverse(vi.begin(),vi.end());
int cap=res[vi[0]][vi[1]];
//find the capacity of the path
for(int i=1; i<vi.size()-1; i++){
cap=min(cap,res[vi[i]][vi[i+1]]);
}
return make_pair(cap,vi);
}
//udpate the residual and flow
bool augment(){
//find a path
auto path=getPath();
auto vi=path.second;
//return false if the path is empty
if(!vi.size()) return false;
//get the capacity of the path
int cap=path.first;
for(int i=0; i<vi.size()-1; i++){
//update the residual
res[vi[i]][vi[i+1]]-=cap;
res[vi[i+1]][vi[i]]+=cap;
//update the flow
if(graph[vi[i]][vi[i+1]]){
fl[vi[i]][vi[i+1]]+=cap;
}
else{
fl[vi[i+1]][vi[i]]-=cap;
}
}
return true;
}
//count the number of prime
//divisor in a number
int divCt(long long num){
if(num==1) return 0;
int ct=0;
//isP checks we get
//a prime greater than
//32 000
bool isP=false;
int cur=num;
//keep dividing by primes until we get
//1 or a prime greater than 32 000
while(cur>1 && !isP){
bool found=false;
for(int i=0; i<ps.size() && !found; i++){
if(!(cur%ps[i])){
cur=cur/ps[i];
found=true;
ct++;
}
}
if(!found && cur>1) isP=true;
}
return ct+isP;
}
//get the gcd of two numbers using the
//Euclidean algorithm
long long gcd(long long a, long long b){
if(b>a) return gcd(b,a);
while(b>1){
int r=a%b;
if(!r) return b;
a=b;
b=r;
}
return 1;
}
int main(){
vector<int> psj;
//find primes below 1000
//by brute force and store
//it in psj
for(int i=2; i<1000; i++){
int cnt=0;
for(int j=2; j<i; j++){
if(!(i%j)){
cnt++;
}
}
if(!cnt){
psj.push_back(i);
}
}
//find primes below 32 000
//using psj to sieve
vector<bool> psb(32000,true);
psb[0]=false;
psb[1]=false;
for(int i:psj){
for(int j=2; i*j<=32000; j++){
psb[i*j]=false;
}
}
for(int i=0; i<32000; i++){
if(psb[i]){
ps.push_back(i);
}
}
//store the sequence given in
//the question
vector<long long> sequence;
sequence.push_back(0);
int n,m;
cin>>n>>m;
for(int i=0; i<n; i++){
int num;
cin>>num;
sequence.push_back(num);
}
//store whether two points are in a pair or not
vector<vector<int>> basicGraph(n+1,vector<int>(n+1,0));
for(int i=0; i<m; i++){
int x,y;
cin>>x>>y;
basicGraph[x][y]=1;
basicGraph[y][x]=1;
}
//initialise flow to all 0s
vector<vector<int>> flow(2*n+2,vector<int>(2*n+2,0));
fl=flow;
//initialise the residual
vector<vector<int>> residual(2*n+2,vector<int>(2*n+2,0));
//fill in values for edges not involving the
//source or the sink
for(int i=1; i<=n; i++){
for(int j=1; j<=n; j++){
if(basicGraph[i][j]){
auto num=gcd(sequence[i],sequence[j]);
residual[i][j+n]=divCt(num);
}
}
}
//fill in values for edges involving the source
//or sink
for(int i=1; i<=n; i++){
int ct=0;
long long cur=sequence[i];
for(int j=1; j<=n; j++){
if(basicGraph[i][j]){
long long other=sequence[j];
long long num=gcd(cur,other);
ct+=divCt(num);
cur=cur/num;
}
}
residual[0][i]=ct;
residual[n+i][2*n+1]=ct;
}
res=residual;
//fill in edge values for the
//network graph
graph=residual;
for(int i=0; i<2*n+2; i++){
for(int j=0; j<2*n+2; j++){
graph[i][j]=min(1,res[i][j]);
}
}
int ans=0;
//keep augmenting until
//we cannot find an augmenting
//path
while(augment()){}
//set ans to the max flow
for(int i=0; i<2*n+2; i++){
ans+=fl[0][i];
}
cout<<ans/2;
return 0;
}
