Pursuing your own 'prune the unnecessary' idea led me to this idea:

Preparation: Find the overall envelope and divide it into equal rectangles creating a grid isomorphic to (0..root(max(x) - min(x))) X (0..root(max(y) - min(y))). This creates discrete grid co-ordinates which can be computed for each point. Store the polygons first as an array and then make a hash(x-grid-coord) of hash(y-grid-coord) of point(=small array) of hash of indexes into the polygon array. Wrap all of this into another hash to store both in a Storable. This has the effect of making a new-style hash (associative array) function as an old-style hash (i.e. indirect addressing that buckets neighbours together in the storage space, to drastically reduce lookup time).

Iteration: To find which polygons contain a test point then requires computing the grid reference of the point and looking up in the hash of hash for which polygons include it. After the one-off preparation is complete, although that is large, this reduces subsequent iterations to one per point, with the nested hash lookups substituting for what was previously a cross-product of iterations.

Update: Could also let the points themselves (ungridded) loose on a hash of hash and test to see which approach is quicker - I have some expectation that the overhead of computing co-ordinates (in the grid) (bucketed hashing) should win versus lookup with large hash element counts when using raw hashing, but by no means certainty.