tests for multi-scale and direct basis conversion

algoim/bernstein.hpp

@@ -23,7 +23,8 @@
 #elif __has_include(<mkl_lapacke.h>)
 #include <mkl_lapacke.h>
 #else
-#error "Algoim requires a LAPACKE implementation to compute eigenvalues and SVD factorisations, but a suitable lapacke.h include file was not found; did you forget to specify its include path?"
+#error \
+    "Algoim requires a LAPACKE implementation to compute eigenvalues and SVD factorisations, but a suitable lapacke.h include file was not found; did you forget to specify its include path?"
 #endif
 
 namespace algoim::bernstein
@@ -36,14 +37,12 @@ namespace algoim::bernstein
     assert(P >= 1);
     const real* binom = Binomial::row(P - 1);
     T           p     = 1.0;
-        for (int i = 0; i < P; ++i)
-        {
+    for (int i = 0; i < P; ++i) {
         out[i]  = p * binom[i];
         p      *= x;
     }
     p = 1.0;
-        for (int i = P - 1; i >= 0; --i)
-        {
+    for (int i = P - 1; i >= 0; --i) {
         out[i] *= p;
         p      *= 1.0 - x;
     }
@@ -55,14 +54,11 @@ namespace algoim::bernstein
 {
     uvector<real*, N> basis;
     algoim_spark_alloc_vec(real, basis, beta.ext());
-        for (int i = 0; i < N; ++i)
-            evalBernsteinBasis(x(i), beta.ext(i), basis(i));
+    for (int i = 0; i < N; ++i) evalBernsteinBasis(x(i), beta.ext(i), basis(i));
     real r = 0.0;
-        for (auto i = beta.loop(); ~i; ++i)
-        {
+    for (auto i = beta.loop(); ~i; ++i) {
         real s = beta.l(i);
-            for (int dim = 0; dim < N; ++dim)
-                s *= basis(dim)[i(dim)];
+        for (int dim = 0; dim < N; ++dim) s *= basis(dim)[i(dim)];
         r += s;
     }
     return r;
@@ -75,11 +71,12 @@ namespace algoim::bernstein
     assert(P > 1);
     real *a, *b;
     algoim_spark_alloc(real, &a, P, &b, P);
-        a[0] = 1; for (int i = 1; i < P; ++i) a[i] = a[i-1] * x;
-        b[0] = 1; for (int i = 1; i < P; ++i) b[i] = b[i-1] * (1 - x);
+    a[0] = 1;
+    for (int i = 1; i < P; ++i) a[i] = a[i - 1] * x;
+    b[0] = 1;
+    for (int i = 1; i < P; ++i) b[i] = b[i - 1] * (1 - x);
     value = alpha[0] * b[P - 1] + alpha[P - 1] * a[P - 1];
-        for (int i = 1; i < P - 1; ++i)
-            value += alpha[i] * binom[i] * a[i] * b[P-1-i];
+    for (int i = 1; i < P - 1; ++i) value += alpha[i] * binom[i] * a[i] * b[P - 1 - i];
     deriv = (alpha[P - 1] * a[P - 2] - alpha[0] * b[P - 2]) * (P - 1);
     for (int i = 1; i < P - 1; ++i)
         deriv += alpha[i] * binom[i] * (a[i - 1] * b[P - 1 - i] * i - a[i] * b[P - 2 - i] * (P - 1 - i));
@@ -92,29 +89,23 @@ namespace algoim::bernstein
     uvector<real*, N> basis, prime;
     algoim_spark_alloc_vec(real, basis, beta.ext());
     algoim_spark_alloc_vec(real, prime, beta.ext());
-        for (int i = 0; i < N; ++i)
-        {
+    for (int i = 0; i < N; ++i) {
         int P = beta.ext(i);
         assert(P >= 1);
         evalBernsteinBasis(x(i), P, basis(i));
-            if (P > 1)
-            {
+        if (P > 1) {
             real* buff;
             algoim_spark_alloc(real, &buff, P - 1);
             evalBernsteinBasis(x(i), P - 1, buff);
             prime(i)[0]     = (P - 1) * (-buff[0]);
             prime(i)[P - 1] = (P - 1) * (buff[P - 2]);
-                for (int j = 1; j < P - 1; ++j)
-                    prime(i)[j] = (P - 1) * (buff[j-1] - buff[j]);
-            }
-            else
+            for (int j = 1; j < P - 1; ++j) prime(i)[j] = (P - 1) * (buff[j - 1] - buff[j]);
+        } else
             prime(i)[0] = 0.0;
     }
     uvector<real, N> g = real(0.0);
-        for (auto i = beta.loop(); ~i; ++i)
-        {
-            for (int j = 0; j < N; ++j)
-            {
+    for (auto i = beta.loop(); ~i; ++i) {
+        for (int j = 0; j < N; ++j) {
             real s = beta.l(i);
             for (int dim = 0; dim < N; ++dim)
                 if (dim == j)
@@ -132,13 +123,10 @@ namespace algoim::bernstein
 void reverseScaledCoeff(xarray<real, N>& p)
 {
     uvector<const real*, N> binom;
-        for (int i = 0; i < N; ++i)
-            binom(i) = Binomial::row(p.ext(i) - 1);
-        for (auto i = p.loop(); i; ++i)
-        {
+    for (int i = 0; i < N; ++i) binom(i) = Binomial::row(p.ext(i) - 1);
+    for (auto i = p.loop(); i; ++i) {
         real alpha = 1;
-            for (int dim = 0; dim < N; ++dim)
-                alpha *= binom(dim)[i(dim)];
+        for (int dim = 0; dim < N; ++dim) alpha *= binom(dim)[i(dim)];
         p.l(i) /= alpha;
     }
 }
@@ -149,21 +137,18 @@ namespace algoim::bernstein
 real squaredL2norm(const xarray<real, N>& p)
 {
     uvector<const real*, N> b1, b2;
-        for (int dim = 0; dim < N; ++dim)
-        {
+    for (int dim = 0; dim < N; ++dim) {
         b1(dim) = Binomial::row(p.ext(dim) - 1);
         b2(dim) = Binomial::row(2 * p.ext(dim) - 2);
     }
     real delta = 0;
-        for (auto i = p.loop(); ~i; ++i) for (auto j = p.loop(); ~j; ++j)
-        {
+    for (auto i = p.loop(); ~i; ++i)
+        for (auto j = p.loop(); ~j; ++j) {
             real g = 1;
-            for (int dim = 0; dim < N; ++dim)
-                g *= (b1(dim)[i(dim)] / b2(dim)[i(dim) + j(dim)]) * b1(dim)[j(dim)];
+            for (int dim = 0; dim < N; ++dim) g *= (b1(dim)[i(dim)] / b2(dim)[i(dim) + j(dim)]) * b1(dim)[j(dim)];
             delta += p.l(i) * p.l(j) * g;
         }
-        for (int dim = 0; dim < N; ++dim)
-            delta /= 2*p.ext(dim) - 1;
+    for (int dim = 0; dim < N; ++dim) delta /= 2 * p.ext(dim) - 1;
     return delta;
 }
 
@@ -173,27 +158,20 @@ namespace algoim::bernstein
 template <int N>
 void collapseAlongAxis(const xarray<real, N>& beta, const uvector<real, N - 1>& x0, int dim, real* out)
 {
-        if constexpr (N == 1)
-        {
+    if constexpr (N == 1) {
         assert(dim == 0);
-            for (int i = 0; i < beta.ext(0); ++i)
-                out[i] = beta[i];
-        }
-        else
-        {
+        for (int i = 0; i < beta.ext(0); ++i) out[i] = beta[i];
+    } else {
         assert(0 <= dim && dim < N);
         uvector<real*, N - 1> basis;
         algoim_spark_alloc_vec(real, basis, remove_component(beta.ext(), dim));
-            for (int i = 0; i < N - 1; ++i)
-            {
+        for (int i = 0; i < N - 1; ++i) {
             int P = beta.ext(i < dim ? i : i + 1);
             evalBernsteinBasis(x0(i), P, basis(i));
         }
         int P = beta.ext(dim);
-            for (int i = 0; i < P; ++i)
-                out[i] = 0.0;
-            for (auto i = beta.loop(); ~i; ++i)
-            {
+        for (int i = 0; i < P; ++i) out[i] = 0.0;
+        for (auto i = beta.loop(); ~i; ++i) {
             real s = beta.l(i);
             for (int j = 0; j < N; ++j)
                 if (j < dim)
@@ -218,8 +196,7 @@ namespace algoim::bernstein
     algoim_spark_alloc(real, &basis, P);
     evalBernsteinBasis(x, P, basis);
     out = 0.0;
-        for (auto i = beta.loop(); i; ++i)
-            out.m(remove_component(i(), dim)) += beta.l(i) * basis[i(dim)];
+    for (auto i = beta.loop(); i; ++i) out.m(remove_component(i(), dim)) += beta.l(i) * basis[i(dim)];
 }
 
 // Normalise polynomial by its largest (in absolute value) coefficient
@@ -227,8 +204,7 @@ namespace algoim::bernstein
 void normalise(xarray<real, N>& alpha)
 {
     real x = alpha.maxNorm();
-        if (x > 0)
-            alpha *= real(1) / x;
+    if (x > 0) alpha *= real(1) / x;
 }
 
 // Applying a simple examination of coefficient signs, returns +1 if the
@@ -239,8 +215,7 @@ namespace algoim::bernstein
 {
     int s = util::sign(beta[0]);
     for (int i = 1; i < beta.size(); ++i)
-            if (util::sign(beta[i]) != s)
-                return 0;
+        if (util::sign(beta[i]) != s) return 0;
     return s;
 }
 
@@ -251,8 +226,7 @@ namespace algoim::bernstein
 void bernsteinDerivative(const T* alpha, int P, T* out)
 {
     assert(P >= 2);
-        for (int i = 0; i < P - 1; ++i)
-        {
+    for (int i = 0; i < P - 1; ++i) {
         out[i]  = alpha[i + 1];
         out[i] -= alpha[i];
         out[i] *= P - 1;
@@ -266,8 +240,7 @@ namespace algoim::bernstein
     assert(all(out.ext() == inc_component(a.ext(), dim, -1)));
     int P = a.ext(dim);
     assert(P >= 2);
-        for (auto i = out.loop(); ~i; ++i)
-            out.l(i) = a.m(i.shifted(dim, 1)) - a.m(i());
+    for (auto i = out.loop(); ~i; ++i) out.l(i) = a.m(i.shifted(dim, 1)) - a.m(i());
     out *= P - 1;
 }
 
@@ -278,14 +251,14 @@ namespace algoim::bernstein
 {
     assert(all(out.ext() == a.ext()) && 0 <= dim && dim < N);
     int P = a.ext(dim);
-        for (auto i = a.loop(); ~i; ++i)
-        {
+    for (auto i = a.loop(); ~i; ++i) {
         if (i(dim) == 0)
             out.l(i) = (a.m(i.shifted(dim, 1)) - a.l(i)) * (P - 1);
         else if (i(dim) == P - 1)
             out.l(i) = (a.l(i) - a.m(i.shifted(dim, -1))) * (P - 1);
         else
-                out.l(i) = a.m(i.shifted(dim, -1)) * (-i(dim)) + a.l(i) * (2*i(dim) - P + 1) + a.m(i.shifted(dim, 1)) * (P - 1 - i(dim));
+            out.l(i) =
+                a.m(i.shifted(dim, -1)) * (-i(dim)) + a.l(i) * (2 * i(dim) - P + 1) + a.m(i.shifted(dim, 1)) * (P - 1 - i(dim));
     }
 }
 
@@ -295,8 +268,7 @@ namespace algoim::bernstein
 {
     int P = alpha.ext(0);
     for (int i = 1; i < P; ++i)
-            for (int j = P - 1; j >= i; --j)
-            {
+        for (int j = P - 1; j >= i; --j) {
             alpha.a(j) *= tau; // here the ptr to the array is included in the return
             alpha.a(j) += alpha.a(j - 1) * (1.0 - tau);
         }
@@ -308,8 +280,7 @@ namespace algoim::bernstein
 {
     int P = alpha.ext(0);
     for (int i = 1; i < P; ++i)
-            for (int j = 0; j < P - i; ++j)
-            {
+        for (int j = 0; j < P - i; ++j) {
             alpha.a(j) *= (1.0 - tau);
             alpha.a(j) += alpha.a(j + 1) * tau;
         }
@@ -322,34 +293,25 @@ namespace algoim::bernstein
 template <int N, bool B = false>
 void deCasteljau(xarray<real, N>& alpha, const real* a, const real* b)
 {
-        using std::swap;
     using std::abs;
-        if constexpr (N == 1 || B)
-        {
+    using std::swap;
+    if constexpr (N == 1 || B) {
         int P = alpha.ext(0);
-            if (*b < *a)
-            {
+        if (*b < *a) {
             deCasteljau<N, B>(alpha, b, a);
-                for (int i = 0; i < P / 2; ++i)
-                    swap(alpha.a(i), alpha.a(P - 1 - i));
+            for (int i = 0; i < P / 2; ++i) swap(alpha.a(i), alpha.a(P - 1 - i));
             return;
         }
-            if (abs(*b) >= abs(*a - 1))
-            {
+        if (abs(*b) >= abs(*a - 1)) {
             deCasteljauLeft(alpha, *b);
             deCasteljauRight(alpha, *a / *b);
-            }
-            else
-            {
+        } else {
             deCasteljauRight(alpha, *a);
             deCasteljauLeft(alpha, (*b - *a) / (real(1) - *a));
         }
-        }
-        else
-        {
+    } else {
         deCasteljau<2, true>(alpha.flatten().ref(), a, b);
-            for (int i = 0; i < alpha.ext(0); ++i)
-                deCasteljau(alpha.slice(i).ref(), a + 1, b + 1);
+        for (int i = 0; i < alpha.ext(0); ++i) deCasteljau(alpha.slice(i).ref(), a + 1, b + 1);
     }
 }
 
@@ -369,24 +331,17 @@ namespace algoim::bernstein
 void bernsteinElevate(const xarray<real, N>& alpha, xarray<real, N>& beta)
 {
     assert(all(beta.ext() >= alpha.ext()));
-        if constexpr (N == 1 || B)
-        {
+    if constexpr (N == 1 || B) {
         int P = alpha.ext(0), Q = beta.ext(0);
-            if (P == Q)
-            {
-                for (int k = 0; k < P; ++k)
-                    beta.a(k) = alpha.a(k);
-            }
-            else
-            {
+        if (P == Q) {
+            for (int k = 0; k < P; ++k) beta.a(k) = alpha.a(k);
+        } else {
             int n = P - 1;
             int r = Q - 1 - n;
-                if (r == 1)
-                {
+            if (r == 1) {
                 beta.a(0)     = alpha.a(0);
                 beta.a(n + 1) = alpha.a(n);
-                    for (int k = 1; k <= n; ++k)
-                    {
+                for (int k = 1; k <= n; ++k) {
                     beta.a(k)  = alpha.a(k - 1) * (real(k) / real(n + 1));
                     beta.a(k) += alpha.a(k) * (real(1) - real(k) / real(n + 1));
                 }
@@ -395,21 +350,17 @@ namespace algoim::bernstein
             const real* bn  = Binomial::row(n);
             const real* br  = Binomial::row(r);
             const real* bnr = Binomial::row(n + r);
-                for (int k = 0; k <= n + r; ++k)
-                {
+            for (int k = 0; k <= n + r; ++k) {
                 beta.a(k) = 0.0;
                 for (int j = std::max(0, k - r); j <= std::min(n, k); ++j)
                     beta.a(k) += alpha.a(j) * ((br[k - j] * bn[j]) / bnr[k]);
             }
         }
-        }
-        else
-        {
+    } else {
         xarray<real, N> gamma(nullptr, set_component(alpha.ext(), 0, beta.ext(0)));
         algoim_spark_alloc(real, gamma);
         bernsteinElevate<2, true>(alpha.flatten(), gamma.flatten().ref());
-            for (int i = 0; i < beta.ext(0); ++i)
-                bernsteinElevate(gamma.slice(i), beta.slice(i).ref());
+        for (int i = 0; i < beta.ext(0); ++i) bernsteinElevate(gamma.slice(i), beta.slice(i).ref());
     }
 }
 
@@ -428,18 +379,15 @@ namespace algoim::bernstein
     assert(b.ext(0) == P + 1 && b.ext(1) > 0);
     real* gamma;
     algoim_spark_alloc_def(real, 0, &gamma, P);
-            for (int i = 0; i < P; ++i)
-            {
+    for (int i = 0; i < P; ++i) {
         real c, s;
         util::givens_get(alpha[i], beta[i], c, s);
         util::givens_rotate(alpha[i], beta[i], c, s);
         if (i < P - 1) util::givens_rotate(gamma[i + 1], alpha[i + 1], c, s);
-                for (int k = 0; k < b.ext(1); ++k)
-                    util::givens_rotate(b(i,k), b(i+1,k), c, s);
+        for (int k = 0; k < b.ext(1); ++k) util::givens_rotate(b(i, k), b(i + 1, k), c, s);
     }
     b.a(P - 1) /= alpha[P - 1];
-            for (int i = P - 2; i >= 0; --i)
-            {
+    for (int i = P - 2; i >= 0; --i) {
         b.a(i) -= b.a(i + 1) * gamma[i + 1];
         b.a(i) /= alpha[i];
     }
@@ -453,35 +401,28 @@ namespace algoim::bernstein
 void bernsteinReduction(xarray<real, N>& alpha, int dim)
 {
     assert(all(alpha.ext() >= 1) && 0 <= dim && dim < N && alpha.ext(dim) >= 2);
-        if (dim == 0)
-        {
+    if (dim == 0) {
         int   P = alpha.ext(0) - 1;
         real *a, *b;
         algoim_spark_alloc(real, &a, P, &b, P);
         a[0]     = 1;
         b[P - 1] = 1;
-            for (int k = 1; k < P; ++k)
-            {
+        for (int k = 1; k < P; ++k) {
             a[k]     = real(1) - real(k) / real(P);
             b[k - 1] = real(k) / real(P);
         }
         xarray<real, 2> view(alpha.data(), uvector<int, 2>{P + 1, prod(alpha.ext(), 0)});
         detail::lsqr_bidiagonal(a, b, P, view);
-        }
-        else if constexpr (N > 1)
-        {
-            for (int i = 0; i < alpha.ext(0); ++i)
-                bernsteinReduction<N-1,true>(alpha.slice(i).ref(), dim - 1);
+    } else if constexpr (N > 1) {
+        for (int i = 0; i < alpha.ext(0); ++i) bernsteinReduction<N - 1, true>(alpha.slice(i).ref(), dim - 1);
     }
 
-        if (!B)
-        {
+    if (!B) {
         xarray<real, N> beta(nullptr, alpha.ext());
         algoim_spark_alloc(real, beta);
         beta = alpha;
         alpha.alterExtent(inc_component(alpha.ext(), dim, -1));
-            for (auto i = alpha.loop(); ~i; ++i)
-                alpha.l(i) = beta.m(i());
+        for (auto i = alpha.loop(); ~i; ++i) alpha.l(i) = beta.m(i());
     }
 }
 
@@ -489,13 +430,11 @@ namespace algoim::bernstein
 template <int N>
 bool autoReduction(xarray<real, N>& alpha, real tol = 1.0e3 * std::numeric_limits<real>::epsilon(), int dim = 0)
 {
-        using std::sqrt;
     using std::abs;
-        if (dim < 0 || dim >= N || tol <= 0)
-            return false;
+    using std::sqrt;
+    if (dim < 0 || dim >= N || tol <= 0) return false;
     bool stay = false;
-        if (alpha.ext(dim) >= 2)
-        {
+    if (alpha.ext(dim) >= 2) {
         xarray<real, N> beta(nullptr, alpha.ext()), gamma(nullptr, alpha.ext());
         algoim_spark_alloc(real, beta, gamma);
         beta = alpha;
@@ -504,19 +443,16 @@ namespace algoim::bernstein
         gamma      -= alpha;
         real delta  = sqrt(abs(squaredL2norm(gamma)));
         real norm   = sqrt(abs(squaredL2norm(alpha)));
-            if (delta < tol * norm)
-            {
+        if (delta < tol * norm) {
             alpha.alterExtent(beta.ext());
             alpha = beta;
             stay  = true;
         }
     }
-        if (stay)
-        {
+    if (stay) {
         autoReduction<N>(alpha, tol, dim);
         return true;
-        }
-        else
+    } else
         return autoReduction<N>(alpha, tol, dim + 1);
 }
 
@@ -527,27 +463,22 @@ namespace algoim::bernstein
 {
     assert(sign == 0 || sign == -1 || sign == 1);
     assert(all(x.ext() == y.ext()));
-        using std::min;
-        using std::max;
-        using std::isinf;
     using std::abs;
-        if (sign == 0)
-            return orthantTestBase(x, y, -1) || orthantTestBase(x, y, 1);
+    using std::isinf;
+    using std::max;
+    using std::min;
+    if (sign == 0) return orthantTestBase(x, y, -1) || orthantTestBase(x, y, 1);
     real alphaMax = std::numeric_limits<real>::infinity();
     real alphaMin = -std::numeric_limits<real>::infinity();
-        for (int i = 0; i < x.size(); ++i)
-        {
-            if (y[i] == 0.0 && x[i] * sign <= 0.0)
-                return false;
+    for (int i = 0; i < x.size(); ++i) {
+        if (y[i] == 0.0 && x[i] * sign <= 0.0) return false;
         if (y[i] > 0.0)
             alphaMin = max(alphaMin, -x[i] / y[i] * sign);
         else if (y[i] < 0.0)
             alphaMax = min(alphaMax, -x[i] / y[i] * sign);
     }
-        if (isinf(alphaMin) || isinf(alphaMax))
-            return true;
-        if (alphaMax - alphaMin > 1.0e5 * std::numeric_limits<real>::epsilon() * max(abs(alphaMin), abs(alphaMax)))
-            return true;
+    if (isinf(alphaMin) || isinf(alphaMax)) return true;
+    if (alphaMax - alphaMin > 1.0e5 * std::numeric_limits<real>::epsilon() * max(abs(alphaMin), abs(alphaMax))) return true;
     return false;
 }
 
@@ -559,8 +490,7 @@ namespace algoim::bernstein
 {
     if (all(f.ext() == g.ext()))
         return orthantTestBase(f, g);
-        else
-        {
+    else {
         uvector<int, N> ext = max(f.ext(), g.ext());
         xarray<real, N> fe(nullptr, ext), ge(nullptr, ext);
         algoim_spark_alloc(real, fe, ge);
@@ -579,8 +509,7 @@ namespace algoim::bernstein
 }
 
 // Methods to compute, and cache, the SVD for Bernstein interpolation based on modified Chebysev nodes
-    struct BernsteinVandermondeSVD
-    {
+struct BernsteinVandermondeSVD {
     struct SVD // 将矩阵A分解为U * diag(sigma) * Vt，其中U和Vt是正交矩阵，sigma是对角矩阵
     {
         real* U;
@@ -592,26 +521,25 @@ namespace algoim::bernstein
     {
         assert(P >= 1);
         static thread_local std::unordered_map<int, std::vector<real>> cache;
-            if (cache.count(P) == 1)
-            {
+        if (cache.count(P) == 1) {
             real* base = cache.at(P).data();
             return SVD{base, base + P * P, base + 2 * P * P};
         }
 
         real *A, *superb, *basis;
         algoim_spark_alloc(real, &A, P * P, &superb, P, &basis, P);
-            for (int i = 0; i < P; ++i)
-            {
+        for (int i = 0; i < P; ++i) {
             evalBernsteinBasis(modifiedChebyshevNode(i, P), P, basis);
-                for (int j = 0; j < P; ++j)
-                    A[i*P + j] = basis[j];
+            for (int j = 0; j < P; ++j) A[i * P + j] = basis[j];
         }
 
         cache[P].resize(P * P + P * P + P);
         real* base = cache[P].data();
         SVD   result{base, base + P * P, base + 2 * P * P};
 
-            static_assert(std::is_same_v<real, double>, "Algoim's default LAPACK code assumes real == double; a custom SVD solver is required when real != double");
+        static_assert(
+            std::is_same_v<real, double>,
+            "Algoim's default LAPACK code assumes real == double; a custom SVD solver is required when real != double");
         int info = LAPACKE_dgesvd(LAPACK_ROW_MAJOR, 'A', 'A', P, P, A, P, result.sigma, result.U, P, result.Vt, P, superb);
         assert(info == 0 && "LAPACKE_dgesvd call failed (algoim::bernstein::BernsteinVandermondeSVD::get)");
         return result;
@@ -625,8 +553,7 @@ namespace algoim::bernstein
 void bernsteinInterpolate(const xarray<real, N>& f, real tol, xarray<real, N>& out)
 {
     assert(all(out.ext() == f.ext()));
-        if constexpr (N == 1 || B)
-        {
+    if constexpr (N == 1 || B) {
         int P = f.ext(0);
         int O = prod(f.ext(), 0);
         assert(P >= 1 && O >= 1);
@@ -636,34 +563,26 @@ namespace algoim::bernstein
 
         auto svd = BernsteinVandermondeSVD::get(P);
 
-            for (int i = 0; i < P * O; ++i)
-                tmp[i] = 0.0;
+        for (int i = 0; i < P * O; ++i) tmp[i] = 0.0;
         for (int i = 0; i < P; ++i)
             for (int j = 0; j < P; ++j)
-                    for (int k = 0; k < O; ++k)
-                        tmp[i*O + k] += svd.U[j*P + i] * f[j*O + k];
+                for (int k = 0; k < O; ++k) tmp[i * O + k] += svd.U[j * P + i] * f[j * O + k];
 
         real minsigma = tol * svd.sigma[0];
-            for (int i = 0; i < P; ++i)
-            {
+        for (int i = 0; i < P; ++i) {
             real alpha = (svd.sigma[i] >= minsigma) ? (real(1) / svd.sigma[i]) : 0.0;
-                for (int k = 0; k < O; ++k)
-                    tmp[i*O + k] *= alpha;
+            for (int k = 0; k < O; ++k) tmp[i * O + k] *= alpha;
         }
 
         out = 0;
         for (int i = 0; i < P; ++i)
             for (int j = 0; j < P; ++j)
-                    for (int k = 0; k < O; ++k)
-                        out[i*O + k] += svd.Vt[j*P + i] * tmp[j*O + k];
-        }
-        else
-        {
+                for (int k = 0; k < O; ++k) out[i * O + k] += svd.Vt[j * P + i] * tmp[j * O + k];
+    } else {
         xarray<real, N> gamma(nullptr, f.ext());
         algoim_spark_alloc(real, gamma);
         bernsteinInterpolate<2, true>(f.flatten(), tol, gamma.flatten().ref());
-            for (int i = 0; i < f.ext(0); ++i)
-                bernsteinInterpolate(gamma.slice(i), tol, out.slice(i).ref());
+        for (int i = 0; i < f.ext(0); ++i) bernsteinInterpolate(gamma.slice(i), tol, out.slice(i).ref());
     }
 }
 
@@ -673,11 +592,9 @@ namespace algoim::bernstein
 {
     xarray<real, N> ff(nullptr, out.ext());
     algoim_spark_alloc(real, ff);
-        for (auto i = ff.loop(); ~i; ++i)
-        {
+    for (auto i = ff.loop(); ~i; ++i) {
         uvector<real, N> x;
-            for (int dim = 0; dim < N; ++dim)
-                x(dim) = modifiedChebyshevNode(i(dim), out.ext(dim));
+        for (int dim = 0; dim < N; ++dim) x(dim) = modifiedChebyshevNode(i(dim), out.ext(dim));
         ff.l(i) = f(x);
     }
     bernsteinInterpolate(ff, std::pow(100.0 * std::numeric_limits<real>::epsilon(), 1.0 / N), out);
@@ -696,14 +613,38 @@ namespace algoim::bernstein
     algoim_spark_alloc(real, &alphar, N, &alphai, N, &beta, N, &lscale, N, &rscale, N);
     real abnrm, bbnrm;
     int  ilo, ihi;
-            static_assert(std::is_same_v<real, double>, "Algoim's default LAPACK code assumes real == double; a custom generalised eigenvalue solver is required when real != double");
-            int info = LAPACKE_dggevx(LAPACK_ROW_MAJOR, 'B', 'N', 'N', 'N', N, A.data(), N, B.data(), N, alphar, alphai, beta, nullptr, N, nullptr, N, &ilo, &ihi, lscale, rscale, &abnrm, &bbnrm, nullptr, nullptr);
+    static_assert(std::is_same_v<real, double>,
+                  "Algoim's default LAPACK code assumes real == double; a custom generalised eigenvalue solver is required "
+                  "when real != double");
+    int info = LAPACKE_dggevx(LAPACK_ROW_MAJOR,
+                              'B',
+                              'N',
+                              'N',
+                              'N',
+                              N,
+                              A.data(),
+                              N,
+                              B.data(),
+                              N,
+                              alphar,
+                              alphai,
+                              beta,
+                              nullptr,
+                              N,
+                              nullptr,
+                              N,
+                              &ilo,
+                              &ihi,
+                              lscale,
+                              rscale,
+                              &abnrm,
+                              &bbnrm,
+                              nullptr,
+                              nullptr);
     assert(info == 0 && "LAPACKE_dggevx call failed (algoim::bernstein::detail::generalisedEigenvalues)");
-            for (int i = 0; i < N; ++i)
-            {
+    for (int i = 0; i < N; ++i) {
         if (beta[i] != 0.0)
-                    out(i,0) = alphar[i] / beta[i],
-                    out(i,1) = alphai[i] / beta[i];
+            out(i, 0) = alphar[i] / beta[i], out(i, 1) = alphai[i] / beta[i];
         else
             out(i, 0) = out(i, 1) = std::numeric_limits<real>::infinity();
     }
@@ -716,17 +657,15 @@ namespace algoim::bernstein
 void rootsBernsteinPoly(const real* alpha, int P, xarray<real, 2>& out)
 {
     assert(P >= 2 && out.ext(0) == P - 1 && out.ext(1) == 2);
-        using std::max;
     using std::abs;
+    using std::max;
 
     real* beta;
     algoim_spark_alloc(real, &beta, P);
     real tol = 0.0;
-        for (int i = 0; i < P; ++i)
-            tol = max(tol, abs(alpha[i]));
+    for (int i = 0; i < P; ++i) tol = max(tol, abs(alpha[i]));
     tol *= util::sqr(std::numeric_limits<real>::epsilon());
-        for (int i = 0; i < P; ++i)
-            beta[i] = (abs(alpha[i]) > tol) ? alpha[i] : 0;
+    for (int i = 0; i < P; ++i) beta[i] = (abs(alpha[i]) > tol) ? alpha[i] : 0;
 
     int             N = P - 1;
     xarray<real, 2> A(nullptr, uvector<int, 2>{N, N});
@@ -734,13 +673,10 @@ namespace algoim::bernstein
     algoim_spark_alloc(real, A, B);
     A = 0;
     B = 0;
-        for (int i = 0; i < N - 1; ++i)
-            A(i, i + 1) = B(i, i + 1) = 1.0;
-        for (int i = 0; i < N; ++i)
-            A(N - 1, i) = B(N - 1, i) = -beta[i];
+    for (int i = 0; i < N - 1; ++i) A(i, i + 1) = B(i, i + 1) = 1.0;
+    for (int i = 0; i < N; ++i) A(N - 1, i) = B(N - 1, i) = -beta[i];
     B(N - 1, N - 1) += beta[N] / N;
-        for (int i = 0; i < N - 1; ++i)
-            B(i, i) = real(N - i) / real(i + 1);
+    for (int i = 0; i < N - 1; ++i) B(i, i) = real(N - i) / real(i + 1);
 
     detail::generalisedEigenvalues(A, B, out);
 }
@@ -756,55 +692,44 @@ namespace algoim::bernstein
     real f0, f1, dummy;
     f(x0, f0, dummy);
     f(x1, f1, dummy);
-            if ((f0 > 0.0 && f1 > 0.0) || (f0 < 0.0 && f1 < 0.0))
-                return false;
-            if (f0 == real(0.0))
-            {
+    if ((f0 > 0.0 && f1 > 0.0) || (f0 < 0.0 && f1 < 0.0)) return false;
+    if (f0 == real(0.0)) {
         root = x0;
         return true;
     }
-            if (f1 == real(0.0))
-            {
+    if (f1 == real(0.0)) {
         root = x1;
         return true;
     }
 
     // x0 and x1 define the bracket; x0 always corresponds to negative value of f; x1 positive value of f
-            if (f1 < 0.0)
-                std::swap(x0, x1);
+    if (f1 < 0.0) std::swap(x0, x1);
 
     // Initial guess is midpoint
     real x = (x0 + x1) * 0.5;
     real fx, fpx;
     f(x, fx, fpx);
     real dx = x1 - x0;
-            for (int step = 0; step < maxsteps; ++step)
-            {
-                if ((fpx*(x - x0) - fx)*(fpx*(x - x1) - fx) < 0.0 && abs(fx) < abs(dx*fpx)*0.5)
-                {
+    for (int step = 0; step < maxsteps; ++step) {
+        if ((fpx * (x - x0) - fx) * (fpx * (x - x1) - fx) < 0.0 && abs(fx) < abs(dx * fpx) * 0.5) {
             // Step in Newton's method falls within bracket and is less than half the previous step size
             dx         = -fx / fpx;
             real xold  = x;
             x         += dx;
-                    if (xold == x)
-                    {
+            if (xold == x) {
                 root = x;
                 return true;
             }
-                }
-                else
-                {
+        } else {
             // Revert to bisection
             dx = (x1 - x0) * 0.5;
             x  = x0 + dx;
-                    if (x == x0)
-                    {
+            if (x == x0) {
                 root = x;
                 return true;
             }
         }
-                if (abs(dx) < tol)
-                {
+        if (abs(dx) < tol) {
             root = x;
             return true;
         }
@@ -821,7 +746,7 @@ namespace algoim::bernstein
     }
     return false;
 }
-    }
+} // namespace detail
 
 // Compute, if possible, a simple real root in [0,1] of a Bernstein polynomial using
 // Descartes' rule of signs:
@@ -836,22 +761,21 @@ namespace algoim::bernstein
     assert(P >= 2);
     using std::abs;
     for (int i = 0; i < P; ++i)
-            if (abs(alpha[i]) < tol)
-                return -1;
+        if (abs(alpha[i]) < tol) return -1;
     int count = 0;
     for (int i = 1; i < P; ++i)
-            if (alpha[i-1] < 0 && alpha[i] >= 0 || alpha[i-1] >= 0 && alpha[i] < 0)
-                ++count;
-        if (count == 0)
-            return 0;
-        if (count > 1)
-            return -1;
+        if (alpha[i - 1] < 0 && alpha[i] >= 0 || alpha[i - 1] >= 0 && alpha[i] < 0) ++count;
+    if (count == 0) return 0;
+    if (count > 1) return -1;
     real        newton_tol = 10.0 * std::numeric_limits<real>::epsilon();
     const real* binom      = Binomial::row(P - 1);
-        bool b = detail::newtonBisectionSearch([=](real x, real& value, real& prime)
-        {
-            bernsteinValueAndDerivative(alpha, P, binom, x, value, prime);
-        }, 0, 1, newton_tol, 12, root);
+    bool        b          = detail::newtonBisectionSearch(
+        [=](real x, real& value, real& prime) { bernsteinValueAndDerivative(alpha, P, binom, x, value, prime); },
+        0,
+        1,
+        newton_tol,
+        12,
+        root);
     return b ? 1 : -1;
 }
 
@@ -864,32 +788,27 @@ namespace algoim::bernstein
     real root;
     int  res = bernsteinSimpleRoot(alpha.data(), alpha.ext(0), tol, root);
     // If it worked with a guarantee of no roots, return
-        if (res == 0)
-            return 0;
+    if (res == 0) return 0;
     // If it worked with a guarantee of just one root computed accurately,
     // transform that root to the [a,b] interval, record it, and return
-        if (res == 1)
-        {
+    if (res == 1) {
         *out = a + (b - a) * root;
         return 1;
     }
     // Otherwise, the simple root method failed. Apply bisection, provided not already too deep
-        if (depth >= 4)
-            return -1;
+    if (depth >= 4) return -1;
     xarray<real, 1> beta(nullptr, alpha.ext());
     algoim_spark_alloc(real, beta);
     // Apply to left half
     beta = alpha;
     deCasteljauLeft(beta, 0.5);
     int r1 = rootsBernsteinPolyFast(beta, a, a + (b - a) * 0.5, depth + 1, tol, out);
-        if (r1 < 0)
-            return -1;
+    if (r1 < 0) return -1;
     // Apply to right half, shifting buffer by r1
     beta = alpha;
     deCasteljauRight(beta, 0.5);
     int r2 = rootsBernsteinPolyFast(beta, a + (b - a) * 0.5, b, depth + 1, tol, out + r1);
-        if (r2 < 0)
-            return -1;
+    if (r2 < 0) return -1;
     return r1 + r2;
 }
 
@@ -903,10 +822,8 @@ namespace algoim::bernstein
     rootsBernsteinPoly(alpha, P, roots);
     real tol   = 1.0e4 * std::numeric_limits<real>::epsilon(); // nearly-real-root tolerance
     int  count = 0;
-        for (int j = 0; j < P - 1; ++j)
-        {
-            if (0 <= roots(j,0) && roots(j,0) <= 1 && abs(roots(j,1)) < tol)
-            {
+    for (int j = 0; j < P - 1; ++j) {
+        if (0 <= roots(j, 0) && roots(j, 0) <= 1 && abs(roots(j, 1)) < tol) {
             *(out + count) = roots(j, 0);
             ++count;
         }
@@ -919,14 +836,13 @@ namespace algoim::bernstein
 // if failed, returns -1
 int bernsteinUnitIntervalRealRoots_fast(const real* alpha, int P, real* out)
 {
-        using std::max;
     using std::abs;
+    using std::max;
     // Compute a tolerance by which to declare a nearly-zero coefficient as being
     // too close to zero (for Descartes' rule of signs and to avoid problems where
     // a root lies close to a subinterval endpoint which can confuse bisection)
     real tol = 0;
-        for (int i = 0; i < P; ++i)
-            tol = max(tol, abs(alpha[i]));
+    for (int i = 0; i < P; ++i) tol = max(tol, abs(alpha[i]));
     tol *= 1.0e4 * std::numeric_limits<real>::epsilon(); // nearly zero coeff tolerance, can be loose
     return rootsBernsteinPolyFast(xarray<real, 1>(const_cast<real*>(alpha), P), 0, 1, 0, tol, out);
 }
@@ -939,12 +855,10 @@ namespace algoim::bernstein
 int bernsteinUnitIntervalRealRoots(const real* alpha, int P, real* out)
 {
     using std::sqrt;
-        if (P == 1)
-            return 0;
+    if (P == 1) return 0;
 
     // Direct method for linear polynomials
-        if (P == 2)
-        {
+    if (P == 2) {
         if (alpha[0] == alpha[1]) return 0;
         real x = alpha[0] / (alpha[0] - alpha[1]);
         if (x < 0 || x > 1) return 0;
@@ -953,8 +867,7 @@ namespace algoim::bernstein
     }
 
     // Direct method for quadratic polynomials, using numerically-stable quadratic formula
-        if (P == 3)
-        {
+    if (P == 3) {
         real a     = alpha[0] - alpha[1] * 2 + alpha[2];
         real b     = (alpha[1] - alpha[0]) * 2;
         real c     = alpha[0];
@@ -964,15 +877,20 @@ namespace algoim::bernstein
         real r1    = q / a;
         real r2    = c / q;
         int  count = 0;
-            if (0 <= r1 && r1 <= 1) { *out = r1; ++count; }
-            if (0 <= r2 && r2 <= 1) { *(out + count) = r2; ++count; }
+        if (0 <= r1 && r1 <= 1) {
+            *out = r1;
+            ++count;
+        }
+        if (0 <= r2 && r2 <= 1) {
+            *(out + count) = r2;
+            ++count;
+        }
         return count;
     }
 
     // Apply fast method, if possible, and resort to eigenvalue method if it fails
     int count = bernsteinUnitIntervalRealRoots_fast(alpha, P, out);
-        if (count >= 0)
-            return count;
+    if (count >= 0) return count;
     return bernsteinUnitIntervalRealRoots_eigenvalue(alpha, P, out);
 }
 
@@ -986,11 +904,9 @@ namespace algoim::bernstein
     const real* bPQ = Binomial::row(P + Q - 3);
     out             = 0;
     for (int i = 0; i < Q - 1; ++i)
-            for (int j = 0; j < P; ++j)
-                out(i,j+i) = a[j] * (bP[j] / bPQ[j + i]);
+        for (int j = 0; j < P; ++j) out(i, j + i) = a[j] * (bP[j] / bPQ[j + i]);
     for (int i = 0; i < P - 1; ++i)
-            for (int j = 0; j < Q; ++j)
-                out(i+Q-1,j+i) = b[j] * (bQ[j] / bPQ[j + i]);
+        for (int j = 0; j < Q; ++j) out(i + Q - 1, j + i) = b[j] * (bQ[j] / bPQ[j + i]);
 }
 
 // Build Bezout matrix for Bernstein polynomials of equal degree P-1
@@ -1000,16 +916,14 @@ namespace algoim::bernstein
     assert(P >= 2 && out.ext(0) == P - 1 && out.ext(1) == P - 1);
     const int n = P - 1;
     out         = 0;
-        for (int i = 1; i <= n; ++i)
-            out(i-1,0) = (a[i] * b[0] - a[0] * b[i]) * real(n) / real(i);
-        for (int j = 1; j <= n - 1; ++j)
-            out(n-1,j) = (a[n] * b[j] - a[j] * b[n]) * real(n) / real(n - j);
+    for (int i = 1; i <= n; ++i) out(i - 1, 0) = (a[i] * b[0] - a[0] * b[i]) * real(n) / real(i);
+    for (int j = 1; j <= n - 1; ++j) out(n - 1, j) = (a[n] * b[j] - a[j] * b[n]) * real(n) / real(n - j);
     for (int i = n - 1; i >= 1; --i)
         for (int j = 1; j <= i - 1; ++j)
-                out(i-1,j) = (a[i] * b[j] - a[j] * b[i]) * real(n*n) / real(i*(n - j)) + out(i,j-1) * real(j*(n-i)) / real(i*(n-j));
+            out(i - 1, j) = (a[i] * b[j] - a[j] * b[i]) * real(n * n) / real(i * (n - j))
+                            + out(i, j - 1) * real(j * (n - i)) / real(i * (n - j));
     for (int i = 0; i < n; ++i)
-            for (int j = i + 1; j < n; ++j)
-                out(i,j) = out(j,i);
+        for (int j = i + 1; j < n; ++j) out(i, j) = out(j, i);
 }
 } // namespace algoim::bernstein
 

examples/examples_quad_multipoly.cpp

@@ -118,6 +118,52 @@ void qConv(const Phi&      phi,
     }
 }
 
+template <typename Phi, typename F>
+void debug3D(const Phi&      phi,
+             real            xmin,
+             real            xmax,
+             uvector<int, 3> P,
+             const F&        integrand,
+             int             q,
+             real            volume_exact,
+             real            surf_exact)
+{
+    // Construct Bernstein polynomial by mapping [0,1] onto bounding box [xmin,xmax]
+    xarray<real, 3> phipoly(nullptr, P);
+    algoim_spark_alloc(real, phipoly);
+    bernstein::bernsteinInterpolate<3>([&](const uvector<real, 3>& x) { return phi(xmin + x * (xmax - xmin)); }, phipoly);
+    DebugXArray<3>()(phipoly, [&](const uvector<real, 3>& x) { return phi(xmin + x * (xmax - xmin)); });
+
+    uvector<real, 3> testX(0., 0., 0.5);
+    real             testEvalBernstein = bernstein::evalBernsteinPoly(phipoly, testX);
+    std::cout << "eval bernstein using interpolation:" << testEvalBernstein << std::endl;
+
+    // Build quadrature hierarchy
+    ImplicitPolyQuadrature<3> ipquad(phipoly);
+
+    // Functional to evaluate volume and surface integrals of given integrand
+    real volume, surf;
+    auto compute = [&](int q) {
+        volume = 0.0;
+        surf   = 0.0;
+        // compute volume integral over {phi < 0} using AutoMixed strategy
+        ipquad.integrate(AutoMixed, q, [&](const uvector<real, 3>& x, real w) {
+            if (bernstein::evalBernsteinPoly(phipoly, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
+        });
+        // compute surface integral over {phi == 0} using AutoMixed strategy
+        ipquad.integrate_surf(AutoMixed, q, [&](const uvector<real, 3>& x, real w, const uvector<real, 3>& wn) {
+            surf += w * integrand(xmin + x * (xmax - xmin));
+        });
+        // scale appropriately
+        volume *= pow(xmax - xmin, 3);
+        surf   *= pow(xmax - xmin, 3 - 1);
+    };
+
+    // Compute results for all q and output in a convergence table
+    compute(q);
+    std::cout << "volume: " << volume << std::endl;
+}
+
 // Given a set of quadrature points and weights, output them to an VTP XML file for visualisation
 // purposes, e.g., using ParaView
 template <int N>
@@ -252,7 +298,7 @@ int main(int argc, char* argv[])
     std::cout << std::scientific << std::setprecision(10);
 
     // q-convergence study for a 2D ellipse
-    if (true) {
+    if (false) {
         auto ellipse      = [](const uvector<real, 2>& x) { return x(0) * x(0) + x(1) * x(1) * 4 - 1; };
         auto integrand    = [](const uvector<real, 2>& x) { return 1.0; };
         real volume_exact = algoim::util::pi / 2;
@@ -338,6 +384,17 @@ int main(int argc, char* argv[])
         std::cout << "(XML VTP file format).\n";
     }
 
+    // a 3D sphere
+    if (true) {
+        auto ellipsoid    = [](const uvector<real, 3>& x) { return x(0) * x(0) + x(1) * x(1) + x(2) * x(2) - 1; };
+        auto integrand    = [](const uvector<real, 3>& x) { return 1.0; };
+        real volume_exact = (algoim::util::pi * 2) / 9;
+        real surf_exact   = 4.400809564664970341600200389229705943483674323377145800356686868037845;
+        std::cout << "\n\nEllipsoid q-convergence test\n";
+        std::cout << "q      volume(q)         surf(q)        vol error       surf error\n";
+        debug3D(ellipsoid, -1., 1., 3, integrand, 10, volume_exact, surf_exact);
+    }
+
     // module_test();
     testMain();
 

gjj/analyticSolution.py

@@ -0,0 +1,35 @@
+import math
+
+class Sphere:
+  def __init__(self, radius, center):
+    self.radius = radius
+    self.center = center
+
+
+# specify class for the input parameters
+def caseSpheresBOnAFace(rA, rB):
+  #intersection:
+  rAsq = rA**2
+  rBsq = rB**2
+  crownHeightA = rBsq / (2 * rA)
+  crownHeightB = rB - crownHeightA
+  crownVolumeA = math.pi * crownHeightA**2 * (3 * rA - crownHeightA) / 3
+  crownVolumeB = math.pi * crownHeightB**2 * (3 * rB - crownHeightB) / 3
+  volumeInter = crownVolumeA + crownVolumeB
+
+  #union:
+  volumeA = 4/3 * math.pi * rAsq * rA
+  volumeB = 4/3 * math.pi * rBsq * rB
+  volumeUnion = volumeA + volumeB - volumeInter
+
+  #difference:
+  volumeDiff = volumeA - volumeInter
+
+  return volumeInter, volumeUnion, volumeDiff
+
+
+if __name__=="__main__":
+    volumeInter, volumeUnion, volumeDiff = caseSpheresBOnAFace(800000, 1000)
+    print("volumeInter = ", volumeInter)
+    print("volumeUnion = ", volumeUnion)
+    print("volumeDiff = ", volumeDiff)

gjj/myDebug.hpp

@@ -1,3 +1,4 @@
+#include <array>
 #include <bitset>
 #include <iostream>
 #include <booluarray.hpp>
@@ -10,6 +11,7 @@
 #include <vector>
 #include "bernstein.hpp"
 #include "quadrature_multipoly.hpp"
+#include "binomial.hpp"
 
 #include "real.hpp"
 #include "uvector.hpp"
@@ -68,6 +70,8 @@ void qConv1(const Phi&      phi,
     }
 }
 
+enum BoolOp { Union, Intersection, Difference };
+
 // Driver method which takes two phi functors defining two polynomials in the reference
 // rectangle [xmin, xmax]^N, each of of Bernstein degree P, builds a quadrature scheme with the
 // given q, and outputs it for visualisation in a set of VTP XML files
@@ -79,13 +83,13 @@ void qConvMultiPloy(const F1&              fphi1,
                     const uvector<int, N>& P,
                     const F&               integrand,
                     int                    q,
-                    std::string            qfile)
+                    BoolOp                 op)
 {
     // Construct phi by mapping [0,1] onto bounding box [xmin,xmax]
     xarray<real, N> phi1(nullptr, P), phi2(nullptr, P);
     algoim_spark_alloc(real, phi1, phi2);
     bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi1(xmin + x * (xmax - xmin)); }, phi1);
-    // bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi2(xmin + x * (xmax - xmin)); }, phi2);
+    bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi2(xmin + x * (xmax - xmin)); }, phi2);
 
     // Build quadrature hierarchy
     ImplicitPolyQuadrature<N> ipquad(phi1, phi2);
@@ -98,7 +102,16 @@ void qConvMultiPloy(const F1&              fphi1,
         surf   = 0.0;
         // compute volume integral over {phi < 0} using AutoMixed strategy
         ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
-            if (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
+            if (op == Union) {
+                if (bernstein::evalBernsteinPoly(phi1, x) < 0 || bernstein::evalBernsteinPoly(phi2, x) < 0)
+                    volume += w * integrand(xmin + x * (xmax - xmin));
+            } else if (op == Intersection) {
+                if (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) < 0)
+                    volume += w * integrand(xmin + x * (xmax - xmin));
+            } else if (op == Difference) {
+                if (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) > 0)
+                    volume += w * integrand(xmin + x * (xmax - xmin));
+            }
             // if (bernstein::evalBernsteinPoly(phi1, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
         });
         // compute surface integral over {phi == 0} using AutoMixed strategy
@@ -114,13 +127,7 @@ void qConvMultiPloy(const F1&              fphi1,
 }
 
 template <int N, typename F>
-void qConvMultiPloy2(
-                    real                   xmin,
-                    real                   xmax,
-                    const uvector<int, N>& P,
-                    const F&               integrand,
-                    int                    q,
-                    std::string            qfile)
+void qConvMultiPloy2(real xmin, real xmax, const uvector<int, N>& P, const F& integrand, int q, std::string qfile)
 {
     // Construct phi by mapping [0,1] onto bounding box [xmin,xmax]
 
@@ -135,7 +142,7 @@ void qConvMultiPloy2(
             real x = xx(0) - centerX;
             real y = xx(1) - centerY;
             real z = xx(2) - centerZ;
-				return x * x + y * y + z * z - 1;
+            return x * x + y * y + z * z - r * r;
         };
         xarray<real, N> phi(nullptr, P);
         algoim_spark_alloc(real, phi);
@@ -146,7 +153,6 @@ void qConvMultiPloy2(
     ImplicitPolyQuadrature<N> ipquad(phis);
 
 
-
     // Build quadrature hierarchy
     // ImplicitPolyQuadrature<N> ipquad(phi1, phi2);
     // ImplicitPolyQuadrature<N> ipquad(phi1);
@@ -158,8 +164,9 @@ void qConvMultiPloy2(
         surf   = 0.0;
         // compute volume integral over {phi < 0} using AutoMixed strategy
         ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
-            // if (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
-            // if (bernstein::evalBernsteinPoly(phi1, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
+            // if (bernstein::evalBernsteinPoly(phi1, x) < 0 && bernstein::evalBernsteinPoly(phi2, x) < 0) volume += w *
+            // integrand(xmin + x * (xmax - xmin)); if (bernstein::evalBernsteinPoly(phi1, x) < 0) volume += w * integrand(xmin
+            // + x * (xmax - xmin));
             volume += w * integrand(xmin + x * (xmax - xmin));
         });
         // compute surface integral over {phi == 0} using AutoMixed strategy
@@ -178,9 +185,7 @@ void qConvMultiPloy2(
     std::cout << "q volume: " << volume << std::endl;
 }
 
-void fromCommon2Bernstein() {
-
-}
+void fromCommon2Bernstein() {}
 
 void testMultiPolys()
 {
@@ -219,12 +224,233 @@ void testMultiPolys()
         // };
         auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
         // qConvMultiPloy<3>(phi0, phi1, -1.0, 1.0, 3, integrand, 3, "exampleC");
-        qConvMultiPloy2<3>(-1.0, 1.0, 3, integrand, 20, "exampleC");
-        std::cout << "\n\nQuadrature visualisation of a 2D implicitly-defined domain involving the\n";
-        std::cout << "intersection of two polynomials, corresponding to the top-left example of Figure 15,\n";
-        std::cout << "https://doi.org/10.1016/j.jcp.2021.110720, written to exampleC-xxxx.vtp files\n";
-        std::cout << "(XML VTP file format).\n";
+        qConvMultiPloy2<3>(-1.0, 1.0, 3, integrand, 10, "exampleC");
+    }
+}
+
+template <int N>
+void power2BernsteinTensorProduct(const xarray<real, N>& phiPower, xarray<real, N>& phiBernsetin)
+{
+    for (auto i = phiPower.loop(); ~i; ++i) {
+        // phi.l(i) = powerFactors.l(i);
+        real factorBase = phiPower.l(i);
+        if (factorBase == 0) continue;
+        auto                           traverseRange = phiPower.ext() - i();
+        std::vector<std::vector<real>> decompFactors(N, std::vector<real>(max(traverseRange), 0.));
+        for (int dim = 0; dim < N; ++dim) {
+            // Sigma
+            size_t      nDim      = phiPower.ext()(dim) - 1;
+            const real* binomNDim = Binomial::row(nDim);
+            for (int j = i(dim); j <= nDim; ++j) {
+                const real* binomJ             = Binomial::row(j);
+                decompFactors[dim][j - i(dim)] = binomJ[i(dim)] / binomNDim[i(dim)];
+            }
+        }
+        xarray<real, N> subgrid(nullptr, traverseRange);
+        // algoim_spark_alloc(real, subgrid);
+
+        for (auto ii = subgrid.loop(); ~ii; ++ii) {
+            real factor = factorBase;
+            for (int dim = 0; dim < N; ++dim) { factor *= decompFactors[dim][ii(dim)]; }
+            phiBernsetin.m(i() + ii()) += factor;
+        }
+    }
+}
+
+template <int N>
+real powerEvaluation(const xarray<real, N>& phi, const uvector<real, N>& x)
+{
+    real res = 0;
+    for (auto i = phi.loop(); ~i; ++i) {
+        real item = phi.l(i);
+        auto exps = i();
+        for (int dim = 0; dim < N; ++dim) { item *= pow(x(dim), exps(dim)); }
+        res += item;
+    }
+    return res;
+}
+
+auto xarray2StdVector(const xarray<real, 3>& array)
+{
+    auto                                        ext = array.ext();
+    std::vector<std::vector<std::vector<real>>> v(ext(0), std::vector<std::vector<real>>(ext(1), std::vector<real>(ext(2), 0)));
+    for (int i = 0; i < ext(0); ++i) {
+        for (int j = 0; j < ext(1); ++j) {
+            for (int k = 0; k < ext(2); ++k) { v[i][j][k] = array.m(uvector<int, 3>(i, j, k)); }
+        }
+    }
+    return v;
+}
+
+template <int N, typename F>
+void qConvBernstein(const xarray<real, N>&                  phi,
+                    real                                    xmin,
+                    real                                    xmax,
+                    const F&                                integrand,
+                    int                                     q,
+                    const std::vector<std::array<real, 4>>& halfFaces)
+{
+    uvector<real, 3> testX(0., 0., 0.5);
+    real             testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX);
+    auto             vec1              = xarray2StdVector(phi);
+    std::cout << "eval bernstein without interpolation:" << testEvalBernstein << std::endl;
+
+    // Build quadrature hierarchy
+    ImplicitPolyQuadrature<N> ipquad(phi);
+    // ImplicitPolyQuadrature<N> ipquad(phi1);
+
+    // Functional to evaluate volume and surface integrals of given integrand
+    real volume;
+    auto compute = [&](int q) {
+        volume = 0.0;
+        // compute volume integral over {phi < 0} using AutoMixed strategy
+        if (!halfFaces.empty()) {
+            ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
+                if (bernstein::evalBernsteinPoly(phi, x) >= 0) return;
+                bool in = true;
+                for (int i = 0; i < halfFaces.size(); ++i) {
+                    uvector<real, 3> factors(halfFaces[i][0], halfFaces[i][1], halfFaces[i][2]);
+                    if (prod(factors * x) + halfFaces[i][3] < 0) {
+                        in = false;
+                        break;
+                    }
+                }
+                if (in) volume += w * integrand(xmin + x * (xmax - xmin));
+            });
+
+        } else
+            ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
+                if (bernstein::evalBernsteinPoly(phi, x) < 0) volume += w * integrand(xmin + x * (xmax - xmin));
+            });
+        // scale appropriately
+        volume *= pow(xmax - xmin, N);
+    };
+    compute(q);
+    std::cout << "q volume: " << volume << std::endl;
+}
+
+template <typename T, int N>
+void xarrayInit(xarray<T, N>& arr)
+{
+    for (auto i = arr.loop(); ~i; ++i) { arr.l(i) = 0; }
+}
+
+void testSpherePowerDirectly()
+{
+    // a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
+    uvector<real, 3> xmin = -1, xmax = 1;
+    uvector<real, 3> range = xmax - xmin;
+    assert(all(range != 0));
+    // uvector<real, 3> k    = 1 / range;
+    // uvector<real, 3> bias = -k * xmin;
+    uvector<real, 3> k    = range;
+    uvector<real, 3> bias = xmin;
+    real             a[3] = {1, 1, 1};
+    real             c[3] = {0, 0, 0};
+    real             r    = 1;
+    uvector<int, 3>  ext  = 3;
+    xarray<real, 3>  phiPower(nullptr, ext), phiBernstein(nullptr, ext);
+    algoim_spark_alloc(real, phiPower, phiBernstein);
+    xarrayInit(phiPower);
+    xarrayInit(phiBernstein);
+    auto v  = xarray2StdVector(phiPower);
+    auto vv = xarray2StdVector(phiBernstein);
+    for (int dim = 0; dim < 3; ++dim) {
+        uvector<int, 3> idx  = 0;
+        phiPower.m(idx)     += a[dim] * (bias(dim) - c[dim]) * (bias(dim) - c[dim]);
+        idx(dim)             = 1;
+        phiPower.m(idx)      = 2 * a[dim] * k(dim) * (bias(dim) - c[dim]);
+        idx(dim)             = 2;
+        phiPower.m(idx)      = a[dim] * k(dim) * k(dim);
+    }
+    phiPower.m(0)  -= r * r;
+    auto integrand  = [](const uvector<real, 3>& x) { return 1.0; };
+    power2BernsteinTensorProduct(phiPower, phiBernstein);
+
+    v = xarray2StdVector(phiPower);
+    uvector<real, 3> testX(0., 0., 0.5);
+    real             testEval = powerEvaluation(phiPower, testX);
+    std::cout << "eval power:" << testEval << std::endl;
+
+    real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
+    auto vec               = xarray2StdVector(phiBernstein);
+    std::cout << "eval bernstein:" << testEval << std::endl;
+    qConvBernstein(phiBernstein, -1, 1, integrand, 10, {});
+}
+
+void testCylinderPowerDirectly()
+{
+    // a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
+    uvector<real, 3> xmin = -1, xmax = 1;
+    uvector<real, 3> range = xmax - xmin;
+    assert(all(range != 0));
+    // uvector<real, 3> k    = 1 / range;
+    // uvector<real, 3> bias = -k * xmin;
+    uvector<real, 3> k    = range;
+    uvector<real, 3> bias = xmin;
+    real             c[3] = {0, 0, 0};
+    real             r = 1, h = 1;
+    uvector<int, 3>  ext = 3;
+    xarray<real, 3>  phiPower(nullptr, ext), phiBernstein(nullptr, ext);
+    algoim_spark_alloc(real, phiPower, phiBernstein);
+    xarrayInit(phiPower);
+    xarrayInit(phiBernstein);
+    auto v   = xarray2StdVector(phiPower);
+    auto vv  = xarray2StdVector(phiBernstein);
+    real top = c[2] + h * 0.5, bottom = c[2] - h * 0.5;
+    phiPower.m(uvector<int, 3>(2, 0, 2)) = -1;
+    phiPower.m(uvector<int, 3>(0, 2, 2)) = -1;
+    phiPower.m(uvector<int, 3>(0, 0, 2)) = r * r;
+    phiPower.m(uvector<int, 3>(2, 0, 1)) = top + bottom;
+    phiPower.m(uvector<int, 3>(0, 2, 1)) = top + bottom;
+    phiPower.m(uvector<int, 3>(0, 0, 1)) = -(bottom + top) * r * r;
+    phiPower.m(uvector<int, 3>(1, 0, 0)) = -bottom * top;
+    phiPower.m(uvector<int, 3>(0, 1, 0)) = -bottom * top;
+    phiPower.m(uvector<int, 3>(0, 0, 0)) = bottom * top * r * r;
+    auto integrand                       = [](const uvector<real, 3>& x) { return 1.0; };
+    power2BernsteinTensorProduct(phiPower, phiBernstein);
+
+    v = xarray2StdVector(phiPower);
+    uvector<real, 3> testX(0., 0., 0.5);
+    real             testEval = powerEvaluation(phiPower, testX);
+    std::cout << "eval power:" << testEval << std::endl;
+
+    real testEvalBernstein = bernstein::evalBernsteinPoly(phiBernstein, testX);
+    auto vec               = xarray2StdVector(phiBernstein);
+    std::cout << "eval bernstein:" << testEval << std::endl;
+    qConvBernstein(phiBernstein,
+                   -1,
+                   1,
+                   integrand,
+                   50,
+                   {
+                       {0, 0, 1,  -bottom},
+                       {0, 0, -1, top    }
+    });
 }
+
+void testMultiScale()
+{
+    auto phi0 = [](const uvector<real, 3>& xx) {
+        real x = xx(0) + 100000.;
+        real y = xx(1);
+        real z = xx(2) - 100000.;
+        real r = 800000.;
+        return x * x + y * y + z * z - r * r;
+    };
+    auto phi1 = [](const uvector<real, 3>& xx) {
+        // real x = xx(0);
+        // real y = xx(1);
+        // real z = xx(2);
+        // return x * x + y * y + z * z - 1;
+        real x = xx(0) + 100000.;
+        real y = xx(1) - 800000.;
+        real z = xx(2) - 100000.;
+        real r = 1000.;
+        return x * x + y * y + z * z - r * r;
+    };
+    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
+    qConvMultiPloy<3>(phi0, phi1, -1000000., 1000000., 3, integrand, 20, Difference);
 }
 
 void testBitSet()
@@ -247,4 +473,7 @@ void testMain()
 {
     testBooluarray();
     testMultiPolys();
+    testMultiScale();
+    testSpherePowerDirectly();
+    // testCylinderPowerDirectly();
 }