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://codeforces.net/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 10⁹ calculating the number of prime divisor is harder. One approach is to calculate all primes up to 1000 by brute force. Then we can use those as a sieve to get all primes up to 32000. Finally since 32000² > 10⁹ with all these primes we can find the number of prime divisors for any number up to 10⁹.

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
//along an augmenting path
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;
}