If I were planning on using this regularly I'd probably write a wrapper function to turn a sage expression into a pygsl-compatible function automatically (N.B. I installed pygsl to make life a little easier):

```
## warning: untested!
import numpy as np
import pygsl.monte as monte
import pygsl.rng
def sage_expr_wrap(some_expr,varorder=None):
if varorder == None: varorder=sorted(some_expr.variables())
def numpy_fn(xx, params):
return np.float(some_expr.subs(dict(zip(varorder, xx))))
return monte.gsl_monte_function(numpy_fn, None, len(varorder))
def sage_vegas(some_expr, xl, xu, ncalls):
xl = np.array(xl,dtype=float)
xu = np.array(xu,dtype=float)
ncalls = int(ncalls)
grand = pygsl.rng.mt19937_1999()
G = sage_expr_wrap(some_expr) # warning: defaults to alphabetical order
s = monte.vegas(len(some_expr.variables()))
s.init()
res, err = s.integrate(G, xl, xu, ncalls, grand)
return res, err
sage: var("x y z")
(x, y, z)
sage: print sage_vegas(3*(x+y+z),[0,0,0],[1,1,1],100)
(4.5695212037401234, 0.0432533513425419)
sage: print sage_vegas(3*(x+y+z),[0,0,0],[1,1,1],1000)
(4.5050689799289199, 0.0051246768259351639)
sage: print sage_vegas(3*(x+y+z),[0,0,0],[1,1,1],10000)
(4.5002245924643995, 0.00031658038710181448)
```

but if speed's a concern obviously the cythonic version will win. :^)