Oops. Yes, of course.
Is there a reference I can look at for this stuff?
Hmm. I would myself use google and wikipedia mainly. “Computational Homology” by Kaczynski, Mischaikow, Mrozek (available on ligben) contains some context, and I know that Mischaikow’s group is also in the business of producing code, not just theorems. That being said, I failed to find any references at all that explicitly talks about the connection between Lefschetz index and interval-math / rectangle based root finding. So this might turn out to become a research project yet 
The way I thought about what you are doing:
You decompose the domain into rectangles which intersect at their boundaries (so you leave no “gap” between x and nextfloat(x)). For each rectangle, and each coordinate function, write down a +1 if the coordinate is positive in the rectangle, a -1 if it is negative, and a 0 if it might be zero.
The rectangles that might contain roots are labeled with zeros(dim).
The intersection of R,R’ gets label zero if both R,R’ are labeled zero, plus if one is labeled plus, minus if one is labeled minus (for each coordinate separately). The plus/minus case cannot happen, due to the fact that the rectangles intersect at their boundaries.
Now, every rectangle R’s boundary decomposes into a union of intersections between R and some R’. The same holds for each subset M of rectangles; the boundary of the union decomposes into a union of intersections between one M-rectangle and one M-complement-rectangle.
Hence, you obtain (without extra calculation) a map from the boundary of the union of maybe-root-containing rectangles to the space of non-zero labels. Even better, this map is continuous, in the sense that the labels are non-contradictory at the transitions, where you move from one bounding rectangle to the next; non-contradictory means no jump from plus to minus (must go via zero).
This map contains the Lefshetz-index! And, by finding connected components of your set of rectangles, you can compute the Lefshetz indices separately (so: each connected maybe-root-containing region gets its own personal Lefshetz index, and each connected component of the boundary of such a region also gets its own index, e.g. if you are left with a maybe-zero-containing double-annulus, like a disc with two disjoint smaller discs punched out).
Is this picture of your procedure roughly correct? If this really is the information we can go on, then I’ll think a little about how to get a proper homology solver for this (with the goal of almost no additional f-evaluations, up to the difference between x and nextfloat(x)).