#include<math.h>#include<iostream>#include<fstream>#include<sstream>#include<utility>#include<vector>#include<iomanip>#include<queue>#include<map>#include<set>usingnamespacestd;/** DEBUG ONLY: If compiled -DDEBUG, each data set will print a C*C table* showing the distances in kilometers between every pair of cities in the data* set. It also prints the distances that were selected for the minimum* spanning tree along with the total minimum cost of the tree. Finally, it* also produces a series of files named "planetsN.dot" (one for each data set)* which are Graphviz representations of the minimum spanning tree that was* computed. These can be fed to the "neato" program to produce a visual graph* layout for verification purposes.*//* SANITY CHECK: assertion macro for verifying input data and internal state */#defineASSERT(e) {if(!(e)) { cerr << #e << endl;throw; } }/** DEBUG ONLY: The "neato" output file currently being written out*/#ifdefDEBUG ofstream graphout;#endif/** DEBUG ONLY: Data type for decimal numbers. This can be changed from double* precision to single precision float in order to evaluate the effect of* round off errors on the final answer.*/typedefdouble decimal_t;/** This data structure holds the distances (in km) between every pair of cities* in the data set. It maps a pair of city numbers to the scalar distance* between them. Even though distances are commutative (i.e. distnace (N. M) is* the same as (M, N)), this data structure holds the full C*C matrix simply* because it's easier to look up these values when running the minimum spanning* tree algorithm.*/typedefpair<int, int> edge_t;typedefmap<edge_t, decimal_t> distmap_t;/** Priority queue for selecting the next lowest cost edge in Prim's minimum* spanning tree algorithm. The priority queue is ordered in ascending order by* edge cost (i.e. distance) and each queue entry also contains the (from, to)* city numbers that define this edge. The TO city number is needed to find the* next set of outgoing edges from this city on the next iteration of Prim's* algorithm. The FROM city number is there purely for debug output purposes.*/typedefpair<decimal_t, edge_t> queuedata_t;typedefpriority_queue<queuedata_t, vector<queuedata_t>, greater<queuedata_t> > queue_t;/* A (latitude, longitude) coordinate pair */typedefpair<decimal_t, decimal_t> coord_t;/* A 3D vector in rectangular coordinates */typedefstruct{ decimal_t x, y, z; } rect_t;/* Multiplying 3D vectors together computes dot product */inlinedecimal_toperator*(rect_tconst&a, rect_tconst&b) {return(a.x * b.x) + (a.y * b.y) + (a.z * b.z); }/** The (latitude, longitude) coordinate pair can be treated as spherical* coordinates in 3D space, where latitude is theta (i.e. zenith) and longitude* is phi (i.e. azimuth). This allows us to convert the coordinates into a* rectangular form. For this conversion we assume a unit sphere (radius 1).** Since the rectangular coordinates are only used for computing angles between* pairs of vectors, the sphere radius does not actually matter at this point.* In fact, using a unit sphere means that we don't have to find the magnitude* of each vector later on since all the vectors we return are of unit length.** Finally, the angular measurements in the dataset are given as degrees while* the math.h trigonometric functions expect radians so we have to do the* conversion before calling sin() and cos().*/rect_ttorect(coord_t coord) {/* Convert angle measurements from degrees to radians */coord.first = coord.first * 2 * M_PI / 360; coord.second = coord.second * 2 * M_PI / 360;/* Convert spherical coordinates to rectangular */rect_t rect; rect.x =cos(coord.first) *sin(coord.second); rect.y =cos(coord.first) *cos(coord.second); rect.z =sin(coord.first);returnrect; }/** Given two points on the surface of the sphere (specified a two 3D vectors)* and given the diameter of the planet, compute the distance between the two* points along the surface of the sphere.** The two vectors define a plane that passes through the origin of the sphere,* and the intersection of this plane with the surface of the sphere forms* a great circle (i.e. a circle that passes through the sphere's origin).* The shortest distance between the two points will be an arc along the* perimiter of this great circle.** The dot product formula states: x * y = cos(A) |x| |y|* In other words, the dot product of two vectors equals the product of the* magnitudes multiplied by the cosine of the angle between the two vectors.* By re-arranging the formula we can calculate the angle A with acos().* Also since we know that torect() only returns unit length vectors, we can* simplify the magnitude multiplication out of the dot product equation.** Once we know angle A, we can calculate the arc length by first computing* the perimeter of the great circle and then multiplying it by the fraction* A/2pi. Also note that the longest distance you can possibly have is half* the perimeter (i.e. opposite ends of the sphere like the North and South* poles).*/decimal_tdistance(rect_tconst&a, rect_tconst&b, decimal_t diameter) { decimal_t perimeter = M_PI * diameter; decimal_t fraction =acos(a * b) / (2 * M_PI);returnfraction * perimeter; }/** Add "city" to the "visited" set and add to the priority "queue" any edges* from "city" to all other unvisitied cities.*/voidprim_visit(int city, int city_num, distmap_t distmap, set<int> &visited, queue_t &queue) { visited.insert(city);for(int i = 0; i < city_num; i++) {if(visited.find(i) == visited.end()) { edge_tedge(city, i); queue.push(queuedata_t(distmap[edge], edge)); } } }/** Use Prim's algorithm to compute a minimum spanning tree that connects all of* the cities together. We start with a complete graph that connects every* pair of cities together. We then arbitrarily pick city number 0 as the* starting point of the algorithm.** DEBUG ONLY: Print out the edge costs of the select minimum spanning tree* to stderr, and write out the complete minimum spanning tree in Graphviz* format to "planetsX.neato" files.*/decimal_tprim(distmap_t &distmap, vector<coord_t> city) { decimal_t total = 0;/* Accumulated total length of visited edges */set<int> visited;/* Track which cities have been visited already */queue_t queue;/* Priority queue to select next minimum cost edge */#ifdefDEBUG cerr << "MSP: "; graphout << "graph G {" << endl;#endif/* Arbitrarily begin with the 0th city. */prim_visit(0, city.size(), distmap, visited, queue);/* Keep running until all cities are visited */while(visited.size() != city.size()) { queuedata_t data;/* Get minimum cost edge to an unvisited city from any visited one */do{/* SANITY CHECK: If unvisited cities remain then so must an edge */ASSERT(queue.size()); data = queue.top(); queue.pop();/* Skip over edges pointing TO already visited cities */}while(visited.find(data.second.second) != visited.end());#ifdefDEBUG coord_t from = city[data.second.first]; coord_t to = city[data.second.second]; graphout << " \"" << from.first << " " << from.second << "\""; graphout << " -- "; graphout << "\"" << to.first << " " << to.second << "\""; graphout << " [label = " << distmap[data.second] << "]" << endl; cerr << data.first << " ";#endif/* Accumulate edge cost and add the next set of edges to queue */total += data.first;prim_visit(data.second.second, city.size(), distmap, visited, queue); }#ifdefDEBUG graphout << "}" << endl;#endifreturntotal; }/* Main body of program */voidprocess(void) { int data_num, data_idx;/* Read how many data sets to process */cin >> data_num;/* Process each data set separately */for(data_idx = 0; data_idx < data_num; data_idx++) { decimal_t diameter, length; vector<coord_t> city; distmap_t distmap; int city_num;/* Read in diameter, length, and city count */cin >> diameter >> length >> city_num;ASSERT(1 <= diameter && diameter <= 1000000);ASSERT(1 <= length && length <= 1000000);ASSERT(1 <= city_num && city_num <= 100);/* Read in the list of city coordinates */for(int i = 0; i < city_num; i++) { decimal_t latitude, longitude; cin >> latitude >> longitude;ASSERT(-90 <= latitude && latitude <= 90);ASSERT(-180 <= longitude && longitude <= 180); city.push_back(coord_t(latitude, longitude)); }/** Precompute the C*C matrix of distances (in km) between all cities* and store the results in the "distance map" for later use by the* minimum spanning tree algorithm.** DEBUG ONLY: Print a C*C table showing the computed distances to* verify that the calculations are correct and to evaluate the effect* due to round off errors with floats vs doubles.*/for(int i = 0; i < city_num; i++) {for(int j = 0; j < city_num; j++) { decimal_t dist =distance(torect(city[i]),torect(city[j]), diameter);#ifdefDEBUG cerr <<setw(4) << dist << " ";#endifdistmap[edge_t(i, j)] = dist; }#ifdefDEBUG cerr << endl;#endif}#ifdefDEBUG/* DEBUG ONLY: Open "neato" graph file for output */ostringstream filename; filename << "planets" << data_idx + 1 << ".neato"; graphout.open(filename.str().c_str(), ios::trunc);#endif/* Compute the minimum cost spanning tree */decimal_t mincost =prim(distmap, city);/* Print the answer for the dataset */#ifdefDEBUG cerr << "= " << mincost << endl; cerr << (mincost <= length ? "IS POSSIBLE" : "IS NOT POSSIBLE") << endl; graphout.close();#endifcout << (mincost <= length ? "IS POSSIBLE" : "IS NOT POSSIBLE") << endl; } }/* Run program and print out any exceptions that occur */intmain(void) {/* Throw exceptions on EOF or failed data extraction in >> operator */cin.exceptions(ios::eofbit | ios::failbit);#ifdefDEBUG/* DEBUG ONLY: Setup floating point format for debug output on stderr */cerr.precision(0); cerr.setf(ios::fixed); cerr.fill('0');#endif/* Run main body of code */try{process(); }/* Catch unexpected EOF or bad input data */catch(ios::failureconst&e) { cerr << "Unexpected EOF or data type mismatch on input" << endl; }return0; }