This program solves TSP instances, input as a series of points. The points are represented as x, y floating point co-ordinates. The edge costs are calculated using the euclidian distance between the points. There are no heuristics used in the algorithm, so it should correctly solve general instances in the same amount of time. To modify it to wo…

tsp.cpp

#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <bitset>
#include <limits>

using namespace std;
using dist = float;
const dist INF = numeric_limits<dist>::max();

struct Point {
  dist x;
  dist y;
};

vector<vector<dist>> distcache;
vector<vector<int>> nckcache;

istream& operator>> (istream& is, Point& p) {
  return is >> p.x >> p.y;
};

// how many unordered subsets of cardinality k can we get from a set of cardinality n
long long nchoosek(int n, int k) {
  long long perms = 1, orderings = 1;

  if (k > n) return 0;

  if (k > (n/2)) {
    // cancel terms
    k = n-k;
  }

  if (k == n)
    return 1;
  if (k > n)
    return 0;

  for (int i = n; i > (n-k); i--) {
    perms *= i;
  }

  for (int i = 2; i <= k; i++) {
    orderings *= i;
  }

  return perms / orderings;
}

// first line of input file is N, the number of cities
// the subsequent N lines represent points in the form <Float> <Float>
vector<Point> loadGraph(istream& is) {
  vector<Point> points;
  int n;
  is >> n;
  cout << "Number of cities: " << n << endl;
  points.reserve(n);
  Point p;
  while (is >> p && is) {
    points.push_back(p);
  }

  return points;
}

dist operator- (const Point& p1, const Point& p2) {
  dist xdist = p2.x - p1.x;
  dist ydist = p2.y - p1.y;
  return sqrt(xdist*xdist + ydist*ydist);
}

ostream& operator<< (ostream& os, const Point &p) {
  return os << "["<<p.x<<","<<p.y<<"]";
}

// gosper's hack to generate the next k-combination
bool gosper(unsigned int& mask) {
  unsigned int u = mask & -mask;
  unsigned int v = u + mask;
  if (v == 0) {
    return false;
  }
  mask = v + (((v^mask)/u)>>2);
  return true;
}

// return the index of this k-combination in the combinatorial number system
// no longer used due to subidxcache
int combinadic_index (const unsigned int combo, const int k) {
  int idx = 0;
  int i = 0;
  unsigned int mask = 1;
  for (int kcount = 1; kcount <= k;) {
    if (i > sizeof(unsigned int)*8) {
      // we've checked all possible bits.  the subset must have not had k elements. error.
      cerr << "Error: subset " << bitset<32>(combo) << " didn't have " << k << " bits set." << endl;
      exit(1);
    }

    if (combo & mask) {
      idx += nckcache[i][kcount++];
    }

    i++;

    mask <<= 1;
  }

  return idx;
}

dist tsp(vector<Point> points) {
  int n = points.size();
  // allocate space for our memo
  // - we will be searching all subsets of n-1 cities (since all tours must have the first city in them)
  // - for each subset S of cardinality k, we want to store solutions to the following problem:
  //    * for all j in S, what is the length of the shortest 1-j path which visits all vertices in S
  //      and visits all elements in S exactly once each.
  // - subproblem size is defined by the cardinality of the set.  each iteration will work over all subsets
  //   of cardinality m.  so, in one iteration, we will only examine at most nchoosek(n-1, m) subsets.
  //   the largest such iteration will be the middle one, where m is ceil((n-1)/2.0), so allocate
  //   that much memory for our memo
  int disttoallocate = nckcache[n - 1][ceil((n-1)/2.0)] * ceil((n-1) / 2.0);
  cout << "Allocate array for " << disttoallocate << " dist vars."<<endl;
  dist *cur = new dist[disttoallocate];
  dist *last = new dist[disttoallocate];
  // an index of subproblem locations to speed things up
  int* subidxcache = new int[(int)ceil(pow(2.0,n-1))];

  // fill out base-cases - edges from 1 to each node
  for (int i = 1; i < n; i++) {
    dist d = distcache[0][i];
    last[i-1] = d;
  }

  for (int m = 2; m < n; m++) {
    // generate first subset of cardinality m
    int num_subsets = nckcache[n-1][m];
    // generate initial subset of index 0
    unsigned int subset = (1 << m) - 1;
    // solve all m-combinations
    for (int i = 0; i < num_subsets; i++) {
      //cout << "Subset " << i << " is " << bitset<24>(subset) << endl;
      subidxcache[subset] = i;
      // for each subset we need to solve a subproblem for all choices of the finishing node, j
      unsigned int jmask = 1;
      // jidx means we are solving the subproblem with node j being the jidx'th element of the subset
      int jidx = 0;
      for (int j = 1; j < n; j++) {
        if (subset & jmask) {
          // node j is in this subset.. so we have a subproblem to solve.
          // work out the index of this subproblem as an (m-1)-combination
          unsigned int submask = subset & ~jmask;
          //cout << "(submask = " << bitset<24>(submask) << endl;
          //int subidx = combinadic_index(submask, m-1);
          int subidx = subidxcache[subset];
          // we want to find node k in the submask which minimises memo[submask][k] + dist[k][j]
          dist bestdist = INF;
          int bestk = 0;
          for (int bit = 1, k = 0; bit <= sizeof(unsigned int)*8 && k < m-1; bit++) {
            if (1 << (bit-1) & submask) {
              if (last[(subidx*(m-1))+k] + distcache[bit][j] < bestdist) {
                bestdist = last[(subidx*(m-1))+k] + distcache[bit][j];
                bestk = k;
              }
              k++;
            }
          }
          cur[(i*m)+jidx] = bestdist;
          jidx++;
        }
        jmask <<= 1;
      }

      gosper(subset); // no need to check if it succeeds since we are bounded by num_subsets
    }

    // swap the cur and last memos for next iteration
    dist *tmp = last; last = cur; cur = tmp;
  }

  // now we have one final check to do.. choose which of the n biggest subproblems minimises the trip back to 1.
  dist bestdist = INF;
  for (int i = 0; i < n-1; i++) {
    if (last[i] + distcache[i+1][0] < bestdist) {
      bestdist = last[i] + distcache[i+1][0];
    }
  }

  delete cur;
  delete last;

  return bestdist;
}

int main (int argc, char *argv[]) {
  if (argc < 2) {
    cerr << "Usage: " << argv[0] << " <infile>" << endl;
    exit(1);
  }

  ifstream is(argv[1]);
  if (!is) {
    cerr << "Couldn't open input file for reading: " << argv[1] << endl;
  }

  vector<Point> points = loadGraph(is);
  int n = points.size();
  if (n > 32) {
    cout << "Cannot solve for more than 32 cities, due to storing subsets as 32-bit ints." << endl;
    cout << "Code must be modified to use larger types (eg. long long) for bigger instances." << endl;
  }

  // now build a cache of edge differences so we only calculate them once.
  distcache.resize(n, vector<dist>(n));
  for (int i = 0; i < n; i++)
    for (int j = 0; j < n; j++)
      distcache[i][j] = abs(points[i]-points[j]);

  // now build a table of nchoosek's we will use
  nckcache.resize(n, vector<int>(n));
  for (int i = 0; i < n; i++)
    for (int j = 0; j < n; j++)
      nckcache[i][j] = nchoosek(i, j);

  // run TSP
  dist shortest = tsp(points);
  cout << "Shortest TSP tour: " << shortest << endl;
}