GIS (as a field) hasn't done too hot when it comes to really grappling with the surface of the globe.
For example, your problem isn't fully defined. Unlike in 2D, where we know the edges of a polygon are composed of straight lines, what are they on the globe? Arcs of great circles, minimizing distance between vertices, are a good choice but not the only one. Straight lines (that thus travel below the surface of the globe) are another choice, for example. (Rhumb lines are a choice as well but potentially a dumb one, especially near the poles...) Here's a great link WRT edge interpretation issues: http://blog.opengeo.org/2010/08/10/shape-of-a-polygon/
Another big issue with the project-to-2d-then-intersect approach, besides the singularities others have mentioned, is that any newly created intersection points (where polygon edges cross) will be out-of-place and dependent on the projection used. (I believe this is where the recommendation for densification comes from: By adding a ton of intermediate vertices, you get a reduced error from projection of the middle of polygon edges.)
Assuming you don't want to make all the compromises and workarounds of projection to 2D and are looking to think through and code something up yourself, I've done a little bit of prototype (and prototype-only!) code for a client dealing with this.
Here's a sketch of the approach. You'll need to know what a vector is as well as the meaning of the dot and cross products. (Caveat: wikipedia is convenient for quick links, but sufficient if it's your first introduction to the topics. A good 3D graphics tutorial will help.)
- Represent a point on the sphere by a unit Cartesian 3D vector. The point on the surface of the earth is where, if you extended the vector into an infinitely-long ray, it would intersect the surface of the earth.
- Represent great circles by a plane through the origin. (In 3D a single unit vector suffices to define a plane through the origin; it is the normal to the plane.) The great circle is the intersection of the entire plane with the surface of the earth.
- You can find the intersection points of two great circles by intersecting their two planes.
- Define an arc of a great circle by two points. The vector that defines the great circle normal is the cross product of the start and end points' vectors.
- To determine if a point that we know is on the great circle lies within the arc, create two planes such that: they are perpendicular to the plane of the great circle, one contains the start point and the other contains the end point, and they are oriented facing each other. Then the point is on the arc iff it lies to the 'inside' of both planes. (To help visualize: you have created pac-man's smile as he chomps down upon the arc. If the test point is between his jaws, then it lies on the arc, as we already know it is upon the great circle.)
- To determine if two arcs cross: find the two points of intersection of their corresponding great circles, then test each point to see if it lies within both arcs.
- A less-ambiguous definition: a polygon is a collection of points, each of which are connected by edges which consist of arcs of great circles. The points are ordered such that if you walk along the surface of the earth along the polygon's edges, the 'inside' of the polygon will lie to your left. Let's leave complex polygons (islands, holes, and self-intersections) out of it for now.
- You can tell whether a point is on the right or left side of a plane via the sign of the dot product of their corresponding vectors. (This is equivalent to, as you travel around the great circle, whether the point on the surface of the sphere lies to your left or your right.)
- An accurate test to determine if a point is inside a polygon: Does it lie on the left side of all the edges?
- Now we have the ability to test if a point is inside a polygon and to determine edge crossings: the ingredients for polygon-polygon intersections! The margin of this comment is too small to write out a full algorithm but the basic steps are to: (a) find all intersections, then (b) walk edges, alternating which polygon you're walking on as you encounter intersection points.
- Once all of the above is working, start thinking about indexing strategies to make it faster, since the point-in-polygon I outlines is O(n) in the number of edges, and polygon intersection O(m * n) in number of edges.
There are some great advantages to this approach: All of the operations above boil down to only multiplications and additions. (After converting data to this representation: e.g. lat/long coordinated require some trig to get a Cartesian XYZ vector.) There are no singularities at the poles or in the coordinate system, and not too many special cases to worry about (example special case: parallel overlapping edges).
Peeking at the code it looks like the Spheres package someone else linked follows some of this approach, although it also looks a bit half baked.
PostGIS may also use a similar approach over its geography datatype, but I haven't looked under the hood. I do know that for spatial indexing, at least, they use an R-tree over 3D cartesian.
(Note: this answer became long enough that I'm probably going to edit up into a blog post... Feedback/comments very welcome!)