Not 100% sure if this would work but one possible algorithm would be to get the bounding box of the shapefile, convert this to Spherical Mercartor (see here for some tips), then work out the tile extent. This might go something like this:
xmaxextent = 40075016.685578488 (max extent of spherical mercartor x)
ymaxextent = 40075016.685578488 (max extent of spherical mercartor y)
mintile_x = ((min x of bbox + (xmaxextent / 2)) / xmaxextent) * (2 ^ (zoomlayer + 1))
maxtile_x = ((max x of bbox + (xmaxextent / 2)) / xmaxextent) * (2 ^ (zoomlayer + 1))
mintile_y = ((min y of bbox + (ymaxextent / 2)) / ymaxextent) * (2 ^ (zoomlayer + 1))
maxtile_y = ((max y of bbox + (ymaxextent / 2)) / ymaxextent) * (2 ^ (zoomlayer + 1))
You will have to round the min numbers down and the max numbers up because they will be calculated as decimals, but essentially that should give you the range for tile extents for any given zoom layer, i.e. (@z2,5-6x,4-5y). You'll then probably need to write a script or something to convert those ranges into real tile numbers. In this hypothetical case that would be:
2/5/4
2/5/5
2/6/4
2/6/5
This link gives some useful background info and a python script to work out the tile number corresponding to any given lat/lon.