Let’s say that I want to integrate some arbitrarily “nice” function (uniformly $ C^{\infty}$ -smooth, for example) over an ImplicitRegion in more than three dimensions. For example, let’s consider the integral $ \int_{D\times D} 1\,\mathrm{d}\mu$ , where $ D$ is the unit disk in two dimensions. We can approximate this in Mathematica as follows:

```
disk2 = ImplicitRegion[x1^2 + y1^2 <= 1 && x2^2 + y2^2 <= 1, {x1, y1, x2, y2}]
val[n_] := NIntegrate[1, Element[{x1, y1, x2, y2}, disk2], WorkingPrecision -> n]
```

The exact value of this is $ \pi^2$ , but the numerical value returned by Mathematica 11.1 doesn’t converge to $ \pi^2$ for increasing $ n$ . Namely, val[50]-Pi^2 and val[100]-Pi^2 both give me the same error, which is roughly $ 1.38\times 10^{-18}$ . With slightly more complicated functions (that are still $ C^1$ , at least), the error we get can be much worse; consider the following integral:

```
NIntegrate[Sqrt[(x1 - x2)^2 + (y1 - y2)^2], Element[{x1, y1, x2, y2}, disk2], WorkingPrecision -> n]
```

We should expect this integral to have a value of $ \frac{128\pi}{45}$ , but the output differs from that value by a little more than $ 4\times 10^{-5}$ , regardless of the setting for WorkingPrecision. Additionally, it doesn’t seem that DiscretizeRegion can really be applied to regions that can’t be expressed as the product of $ n$ -dimensional regions for $ n\leq 3$ , which in the case of a more general region rules out the strategy here. Are there any general strategies I can use to get with arbitrary precision and accuracy the values of numerical integrals of $ C^1$ functions on ImplicitRegions in more than three dimensions? Thanks for your time.