Yet with formal verification we can prove that the implementation actually satisfies the specification for all of infinitely many possible inputs.
One complication: For union the spec would look a little bit different than for intersection. I defined the solid associated to a mesh in a way that does not include the surface of the mesh itself. That works for the intersection. For union we would have:
solidWithSurface (meshUnion M₁ M₂) = solidWithSurface M₁ ∪ solidWithSurface M₂
I think it would not be hard to reduce the implementation and the proof of correctness to the intersection: There are some tricks which you can do like flipping the orientation of a mesh and putting a large cube that contains both meshes around it and then intersect these. (Hopefully the implementation would not actually do this, but this construction could be used in the proof.)
I'm also curious if you've consider fuzz testing this system to verify the results.