You should be able to get accurate results with [mpmath][1], a Python module for arbitrary-precision floating-point computations.  There are examples of integration with singularities in the [documentation][2].  You'll want to explicitly tell it to break up the interval:

    from mpmath import *
    f = lambda x,y,z: 1./(x**2+y**2+z**2)**1./3
    quad(f,[-1,0,1],[-1,0,1],[-1,0,1])

You may need to increase the precision (e.g. `mp.dps=30`) and it will likely be slow, but should be quite accurate.

You could also try nesting calls to MATLAB's `quadgk()`, which uses adaptive Gauss-Kronrod quadrature in 1D.

  [1]: http://mpmath.googlecode.com/svn/trunk/doc/build/index.html
  [2]: http://mpmath.googlecode.com/svn/trunk/doc/build/calculus/integration.html