power basis

.clang-format

@@ -1,4 +1,3 @@
----
 BasedOnStyle: Google
 AccessModifierOffset: -4
 AlignAfterOpenBracket: Align

.vscode/settings.json

@@ -5,5 +5,13 @@
     "array": "cpp",
     "iostream": "cpp"
 	},
-	"C_Cpp.default.compilerPath": "/usr/bin/clang++"
+	"C_Cpp.default.compilerPath": "/usr/bin/clang++",
+  "[cpp]": {
+  "editor.defaultFormatter": "llvm-vs-code-extensions.vscode-clangd",
+  "editor.formatOnSave": true
+},
+"[c]": {
+  "editor.defaultFormatter": "llvm-vs-code-extensions.vscode-clangd",
+  "editor.formatOnSave": true
+}
 }

algoim/xarray.hpp

@@ -11,335 +11,350 @@
 
 namespace algoim
 {
-    // MiniLoop<N> is an optimised version of MultiLoop<N> with zero lowerbound
-    template<int N>
-    class MiniLoop
-    {
-        uvector<int,N> i;
-        int iexp;
-        const uvector<int,N> ext;
-    public:
-        explicit MiniLoop(const uvector<int,N>& ext) : ext(ext), i(0), iexp(0) {}
-
-        MiniLoop& operator++()
-        {
-            ++iexp;
-            for (int dim = N - 1; dim >= 0; --dim)
-            {
-                if (++i(dim) < ext(dim))
-                    return *this;
-                if (dim == 0)
-                    return *this;
-                i(dim) = 0;
-            }
-            return *this;
-        }
+// MiniLoop<N> is an optimised version of MultiLoop<N> with zero lowerbound
+template <int N>
+class MiniLoop
+{
+    uvector<int, N>       i;
+    int                   iexp;
+    const uvector<int, N> ext;
 
-        const uvector<int,N>& operator()() const
-        {
-            return i;
-        }
+public:
+    explicit MiniLoop(const uvector<int, N>& ext) : ext(ext), i(0), iexp(0) {}
 
-        int operator()(int index) const
-        {
-            return i(index);
+    MiniLoop& operator++()
+    {
+        ++iexp;
+        for (int dim = N - 1; dim >= 0; --dim) {
+            if (++i(dim) < ext(dim)) return *this;
+            if (dim == 0) return *this;
+            i(dim) = 0;
         }
+        return *this;
+    }
 
-        bool operator~() const
-        {
-            return i(0) < ext(0);
-        }
+    const uvector<int, N>& operator()() const { return i; }
 
-        int furl() const
-        {
-            return iexp;
-        }
+    int operator()(int index) const { return i(index); }
 
-        uvector<int,N> shifted(int dim, int amount) const
-        {
-            uvector<int,N> j = i;
-            j(dim) += amount;
-            return j;
-        }
-    };
+    bool operator~() const { return i(0) < ext(0); }
 
-    template<typename T>
-    struct xarraySlice
-    {
-        T* ptr;
-        int len;
-
-        xarraySlice(xarraySlice&) = delete;
-        xarraySlice(xarraySlice&&) = delete;
-
-        xarraySlice& operator= (const T& x) { for (int i = 0; i < len; ++i) ptr[i] =  x; return *this; }
-        xarraySlice& operator+=(const T& x) { for (int i = 0; i < len; ++i) ptr[i] += x; return *this; }
-        xarraySlice& operator-=(const T& x) { for (int i = 0; i < len; ++i) ptr[i] -= x; return *this; }
-        xarraySlice& operator*=(const T& x) { for (int i = 0; i < len; ++i) ptr[i] *= x; return *this; }
-        xarraySlice& operator/=(const T& x) { for (int i = 0; i < len; ++i) ptr[i] /= x; return *this; }
-
-        xarraySlice& operator= (const xarraySlice& x) { for (int i = 0; i < len; ++i) ptr[i] =  x.ptr[i]; return *this; }
-        xarraySlice& operator+=(const xarraySlice& x) { for (int i = 0; i < len; ++i) ptr[i] += x.ptr[i]; return *this; }
-        xarraySlice& operator-=(const xarraySlice& x) { for (int i = 0; i < len; ++i) ptr[i] -= x.ptr[i]; return *this; }
-        xarraySlice& operator*=(const xarraySlice& x) { for (int i = 0; i < len; ++i) ptr[i] *= x.ptr[i]; return *this; }
-        xarraySlice& operator/=(const xarraySlice& x) { for (int i = 0; i < len; ++i) ptr[i] /= x.ptr[i]; return *this; }
-
-        template<typename X, typename Y>
-        struct prod
-        {
-            const X& x;
-            const Y& y;
-        };
-
-        xarraySlice& operator= (const prod<xarraySlice,T>& op) { for (int i = 0; i < len; ++i) ptr[i] =  op.x.ptr[i] * op.y; return *this; }
-        xarraySlice& operator+=(const prod<xarraySlice,T>& op) { for (int i = 0; i < len; ++i) ptr[i] += op.x.ptr[i] * op.y; return *this; }
-        xarraySlice& operator-=(const prod<xarraySlice,T>& op) { for (int i = 0; i < len; ++i) ptr[i] -= op.x.ptr[i] * op.y; return *this; }
-        xarraySlice& operator*=(const prod<xarraySlice,T>& op) { for (int i = 0; i < len; ++i) ptr[i] *= op.x.ptr[i] * op.y; return *this; }
-        xarraySlice& operator/=(const prod<xarraySlice,T>& op) { for (int i = 0; i < len; ++i) ptr[i] /= op.x.ptr[i] * op.y; return *this; }
-    };
+    int furl() const { return iexp; }
+
+    uvector<int, N> shifted(int dim, int amount) const
+    {
+        uvector<int, N> j  = i;
+        j(dim)            += amount;
+        return j;
+    }
+};
+
+template <typename T>
+struct xarraySlice {
+    T*  ptr;
+    int len;
+
+    xarraySlice(xarraySlice&)  = delete;
+    xarraySlice(xarraySlice&&) = delete;
+
+    xarraySlice& operator=(const T& x)
+    {
+        for (int i = 0; i < len; ++i) ptr[i] = x;
+        return *this;
+    }
+
+    xarraySlice& operator+=(const T& x)
+    {
+        for (int i = 0; i < len; ++i) ptr[i] += x;
+        return *this;
+    }
+
+    xarraySlice& operator-=(const T& x)
+    {
+        for (int i = 0; i < len; ++i) ptr[i] -= x;
+        return *this;
+    }
 
-    template<typename T>
-    auto operator* (const xarraySlice<T>& x, const T& y)
+    xarraySlice& operator*=(const T& x)
     {
-        return typename xarraySlice<T>::template prod<xarraySlice<T>,T>{x, y};
+        for (int i = 0; i < len; ++i) ptr[i] *= x;
+        return *this;
+    }
+
+    xarraySlice& operator/=(const T& x)
+    {
+        for (int i = 0; i < len; ++i) ptr[i] /= x;
+        return *this;
+    }
+
+    xarraySlice& operator=(const xarraySlice& x)
+    {
+        for (int i = 0; i < len; ++i) ptr[i] = x.ptr[i];
+        return *this;
+    }
+
+    xarraySlice& operator+=(const xarraySlice& x)
+    {
+        for (int i = 0; i < len; ++i) ptr[i] += x.ptr[i];
+        return *this;
+    }
+
+    xarraySlice& operator-=(const xarraySlice& x)
+    {
+        for (int i = 0; i < len; ++i) ptr[i] -= x.ptr[i];
+        return *this;
+    }
+
+    xarraySlice& operator*=(const xarraySlice& x)
+    {
+        for (int i = 0; i < len; ++i) ptr[i] *= x.ptr[i];
+        return *this;
+    }
+
+    xarraySlice& operator/=(const xarraySlice& x)
+    {
+        for (int i = 0; i < len; ++i) ptr[i] /= x.ptr[i];
+        return *this;
+    }
+
+    template <typename X, typename Y>
+    struct prod {
+        const X& x;
+        const Y& y;
     };
 
-    template<typename T>
-    void swap(const xarraySlice<T>& x, const xarraySlice<T>& y)
+    xarraySlice& operator=(const prod<xarraySlice, T>& op)
     {
-        using std::swap;
-        for (int i = 0; i < x.len; ++i)
-            swap(x.ptr[i], y.ptr[i]);
+        for (int i = 0; i < len; ++i) ptr[i] = op.x.ptr[i] * op.y;
+        return *this;
     }
 
-    // algoim::xarray implements an N-dimensional array view of a memory
-    // block controlled by the user
-    template<typename T, int N>
-    class xarray
+    xarraySlice& operator+=(const prod<xarraySlice, T>& op)
     {
-        T* data_;
-        uvector<int,N> ext_;
-        friend class SparkStack<T>;
-    public:
-        // interpret the given block of memory as an N-dimensional array of the given extent
-        xarray(T* data, const uvector<int,N>& ext) : data_(data), ext_(ext) {}
+        for (int i = 0; i < len; ++i) ptr[i] += op.x.ptr[i] * op.y;
+        return *this;
+    }
 
-        xarray(const xarray&) = delete;
+    xarraySlice& operator-=(const prod<xarraySlice, T>& op)
+    {
+        for (int i = 0; i < len; ++i) ptr[i] -= op.x.ptr[i] * op.y;
+        return *this;
+    }
 
-        xarray& operator=(const xarray& x)
-        {
-            assert(same_shape(x));
-            for (int i = 0; i < size(); ++i)
-                data_[i] = x.data_[i];
-            return *this;
-        }
+    xarraySlice& operator*=(const prod<xarraySlice, T>& op)
+    {
+        for (int i = 0; i < len; ++i) ptr[i] *= op.x.ptr[i] * op.y;
+        return *this;
+    }
 
-        template<typename S>
-        xarray& operator=(const S& x)
-        {
-            for (int i = 0; i < size(); ++i)
-                data_[i] = x;
-            return *this;
-        }
+    xarraySlice& operator/=(const prod<xarraySlice, T>& op)
+    {
+        for (int i = 0; i < len; ++i) ptr[i] /= op.x.ptr[i] * op.y;
+        return *this;
+    }
+};
 
-        template<typename S>
-        xarray& operator+=(const S& x)
-        {
-            for (int i = 0; i < size(); ++i)
-                data_[i] += x;
-            return *this;
-        }
+template <typename T>
+auto operator*(const xarraySlice<T>& x, const T& y)
+{
+    return typename xarraySlice<T>::template prod<xarraySlice<T>, T>{x, y};
+};
 
-        xarray& operator+=(const xarray& x)
-        {
-            assert(same_shape(x));
-            for (int i = 0; i < size(); ++i)
-                data_[i] += x.data_[i];
-            return *this;
-        }
+template <typename T>
+void swap(const xarraySlice<T>& x, const xarraySlice<T>& y)
+{
+    using std::swap;
+    for (int i = 0; i < x.len; ++i) swap(x.ptr[i], y.ptr[i]);
+}
 
-        template<typename S>
-        xarray& operator-=(const S& x)
-        {
-            for (int i = 0; i < size(); ++i)
-                data_[i] -= x;
-            return *this;
-        }
+// algoim::xarray implements an N-dimensional array view of a memory
+// block controlled by the user
+template <typename T, int N>
+class xarray
+{
+    T*              data_;
+    uvector<int, N> ext_;
+    friend class SparkStack<T>;
 
-        xarray& operator-=(const xarray& x)
-        {
-            assert(same_shape(x));
-            for (int i = 0; i < size(); ++i)
-                data_[i] -= x.data_[i];
-            return *this;
-        }
+public:
+    // interpret the given block of memory as an N-dimensional array of the given extent
+    xarray(T* data, const uvector<int, N>& ext) : data_(data), ext_(ext) {}
 
-        template<typename S>
-        xarray& operator*=(const S& x)
-        {
-            for (int i = 0; i < size(); ++i)
-                data_[i] *= x;
-            return *this;
-        }
+    xarray(const xarray& x)
+    {
+        this->data_ = nullptr;
+        this->ext_  = x.ext();
+        algoim_spark_alloc(T, *this);
+        *this = x;
+    };
 
-        const T* data() const { return data_; }
-        T* data() { return data_; }
+    xarray& operator=(const xarray& x)
+    {
+        assert(same_shape(x));
+        for (int i = 0; i < size(); ++i) data_[i] = x.data_[i];
+        return *this;
+    }
 
-        // Accessors for user-expanded index, 0 <= i < prod(ext)
-        const T& operator[](int i) const { return data_[i]; }
-        T& operator[](int i) { return data_[i]; }
+    template <typename S>
+    xarray& operator=(const S& x)
+    {
+        for (int i = 0; i < size(); ++i) data_[i] = x;
+        return *this;
+    }
 
-        // Slice and flatten this(i,:)
-        template<int NN = N, std::enable_if_t<NN == 1, bool> = true>
-        T& a(int i)
-        {
-            return data_[i];
-        }
+    template <typename S>
+    xarray& operator+=(const S& x)
+    {
+        for (int i = 0; i < size(); ++i) data_[i] += x;
+        return *this;
+    }
 
-        // Slice and flatten this(i,:)
-        template<int NN = N, std::enable_if_t<NN == 1, bool> = true>
-        const T& a(int i) const
-        {
-            return data_[i];
-        }
+    xarray& operator+=(const xarray& x)
+    {
+        assert(same_shape(x));
+        for (int i = 0; i < size(); ++i) data_[i] += x.data_[i];
+        return *this;
+    }
 
-        // Slice and flatten this(i,:)
-        template<int NN = N, std::enable_if_t<(NN > 1), bool> = true>
-        auto a(int i) const
-        {
-            int span = prod(ext_, 0);
-            return xarraySlice<T>{data_ + i * span, span};
-        }
+    template <typename S>
+    xarray& operator-=(const S& x)
+    {
+        for (int i = 0; i < size(); ++i) data_[i] -= x;
+        return *this;
+    }
 
-        const T& operator()(int i, int j) const
-        {
-            static_assert(N == 2, "N == 2 required for integer pair access");
-            return data_[i*ext(1) + j];
-        }
+    xarray& operator-=(const xarray& x)
+    {
+        assert(same_shape(x));
+        for (int i = 0; i < size(); ++i) data_[i] -= x.data_[i];
+        return *this;
+    }
 
-        T& operator()(int i, int j)
-        {
-            static_assert(N == 2, "N == 2 required for integer pair access");
-            return data_[i*ext(1) + j];
-        }
+    template <typename S>
+    xarray& operator*=(const S& x)
+    {
+        for (int i = 0; i < size(); ++i) data_[i] *= x;
+        return *this;
+    }
 
-        // Accessors by multi-index
-        const T& m(const uvector<int,N>& i) const
-        {
-            return data_[util::furl(i, ext_)];
-        }
+    const T* data() const { return data_; }
 
-        // Accessors by multi-index
-        T& m(const uvector<int,N>& i)
-        {
-            return data_[util::furl(i, ext_)];
-        }
+    T* data() { return data_; }
 
-        // Accessors by mini-loop, where it is assumed the MiniLoop object's
-        // associated extent is identical to this xarray
-        const T& l(const MiniLoop<N>& i) const
-        {
-            return data_[i.furl()];
-        }
+    // Accessors for user-expanded index, 0 <= i < prod(ext)
+    const T& operator[](int i) const { return data_[i]; }
 
-        T& l(const MiniLoop<N>& i)
-        {
-            return data_[i.furl()];
-        }
+    T& operator[](int i) { return data_[i]; }
 
-        const uvector<int,N>& ext() const
-        {
-            return ext_;
-        }
+    // Slice and flatten this(i,:)
+    template <int NN = N, std::enable_if_t<NN == 1, bool> = true>
+    T& a(int i)
+    {
+        return data_[i];
+    }
 
-        int ext(int i) const
-        {
-            return ext_(i);
-        }
+    // Slice and flatten this(i,:)
+    template <int NN = N, std::enable_if_t<NN == 1, bool> = true>
+    const T& a(int i) const
+    {
+        return data_[i];
+    }
 
-        int size() const
-        {
-            return prod(ext_);
-        }
+    // Slice and flatten this(i,:)
+    template <int NN = N, std::enable_if_t<(NN > 1), bool> = true>
+    auto a(int i) const
+    {
+        int span = prod(ext_, 0);
+        return xarraySlice<T>{data_ + i * span, span};
+    }
 
-        MiniLoop<N> loop() const
-        {
-            return MiniLoop<N>(ext_);
-        }
+    const T& operator()(int i, int j) const
+    {
+        static_assert(N == 2, "N == 2 required for integer pair access");
+        return data_[i * ext(1) + j];
+    }
 
-        // xarray is strided so that the inner-most dimension (dim 0) is the slowest varying;
-        // as such, one may view an N-dimensionl xarray as a length ext(0) vector of (N-1)-
-        // dimensional xarrays. Two views of this(i,:) are provided:
-        //   - flatten() returns a 2D view, of dimensions ext(0) by ext(1) * ... * ext(N - 1)
-        //   - slice(i) returns an (N-1)-D array of extent (ext(1), ..., ext(N-1))
-        auto flatten() const
-        {
-            return xarray<T,2>(data_, uvector<int,2>{ext_(0), prod(ext_, 0)});
-        }
+    T& operator()(int i, int j)
+    {
+        static_assert(N == 2, "N == 2 required for integer pair access");
+        return data_[i * ext(1) + j];
+    }
 
-        // xarray is strided so that the inner-most dimension (dim 0) is the slowest varying;
-        // as such, one may view an N-dimensionl xarray as a length ext(0) vector of (N-1)-
-        // dimensional xarrays. Two views of this(i,:) are provided:
-        //   - flatten() returns a 2D view, of dimensions ext(0) by ext(1) * ... * ext(N - 1)
-        //   - slice(i) returns an (N-1)-D array of extent (ext(1), ..., ext(N-1))
-        auto slice(int i) const
-        {
-            return xarray<T,N-1>(data_ + i * prod(ext_, 0), remove_component(ext_, 0));
-        }
+    // Accessors by multi-index
+    const T& m(const uvector<int, N>& i) const { return data_[util::furl(i, ext_)]; }
 
-        bool same_shape(const xarray& x) const
-        {
-            return all(x.ext_ == ext_);
-        }
+    // Accessors by multi-index
+    T& m(const uvector<int, N>& i) { return data_[util::furl(i, ext_)]; }
 
-        T maxNorm() const
-        {
-            assert(size() > 0);
-            using std::abs;
-            using std::max;
-            T m = abs(data_[0]);
-            for (int i = 1; i < size(); ++i)
-                m = max(m, abs(data_[i]));
-            return m;
-        }
+    // Accessors by mini-loop, where it is assumed the MiniLoop object's
+    // associated extent is identical to this xarray
+    const T& l(const MiniLoop<N>& i) const { return data_[i.furl()]; }
 
-        T min() const
-        {
-            assert(size() > 0);
-            using std::min;
-            T m = data_[0];
-            for (int i = 1; i < size(); ++i)
-                m = min(m, data_[i]);
-            return m;
-        }
+    T& l(const MiniLoop<N>& i) { return data_[i.furl()]; }
 
-        T max() const
-        {
-            assert(size() > 0);
-            using std::max;
-            T m = data_[0];
-            for (int i = 1; i < size(); ++i)
-                m = max(m, data_[i]);
-            return m;
-        }
+    const uvector<int, N>& ext() const { return ext_; }
 
-        // Alter the extent of the view of the existing memory block; usually this
-        // only makes sense when the new extent is smaller than the existing extent,
-        // i.e., prod(new extent) <= prod(old extent)
-        void alterExtent(const uvector<int,N>& ext)
-        {
-            ext_ = ext;
-        }
+    int ext(int i) const { return ext_(i); }
 
-        // flatten() and slice() return temporary objects which cannot be bound to non-const references;
-        // applying .ref() to the temporary object allows this to be done and is an assertion by the
-        // user that doing so is safe (relating to the output of flatten() or slice() being destroyed
-        // at the end of the full-expression)
-        xarray<T,N>& ref()
-        {
-            return *this;
-        }
-    };
+    int size() const { return prod(ext_); }
+
+    MiniLoop<N> loop() const { return MiniLoop<N>(ext_); }
+
+    // xarray is strided so that the inner-most dimension (dim 0) is the slowest varying;
+    // as such, one may view an N-dimensionl xarray as a length ext(0) vector of (N-1)-
+    // dimensional xarrays. Two views of this(i,:) are provided:
+    //   - flatten() returns a 2D view, of dimensions ext(0) by ext(1) * ... * ext(N - 1)
+    //   - slice(i) returns an (N-1)-D array of extent (ext(1), ..., ext(N-1))
+    auto flatten() const { return xarray<T, 2>(data_, uvector<int, 2>{ext_(0), prod(ext_, 0)}); }
+
+    // xarray is strided so that the inner-most dimension (dim 0) is the slowest varying;
+    // as such, one may view an N-dimensionl xarray as a length ext(0) vector of (N-1)-
+    // dimensional xarrays. Two views of this(i,:) are provided:
+    //   - flatten() returns a 2D view, of dimensions ext(0) by ext(1) * ... * ext(N - 1)
+    //   - slice(i) returns an (N-1)-D array of extent (ext(1), ..., ext(N-1))
+    auto slice(int i) const { return xarray<T, N - 1>(data_ + i * prod(ext_, 0), remove_component(ext_, 0)); }
+
+    bool same_shape(const xarray& x) const { return all(x.ext_ == ext_); }
+
+    T maxNorm() const
+    {
+        assert(size() > 0);
+        using std::abs;
+        using std::max;
+        T m = abs(data_[0]);
+        for (int i = 1; i < size(); ++i) m = max(m, abs(data_[i]));
+        return m;
+    }
+
+    T min() const
+    {
+        assert(size() > 0);
+        using std::min;
+        T m = data_[0];
+        for (int i = 1; i < size(); ++i) m = min(m, data_[i]);
+        return m;
+    }
+
+    T max() const
+    {
+        assert(size() > 0);
+        using std::max;
+        T m = data_[0];
+        for (int i = 1; i < size(); ++i) m = max(m, data_[i]);
+        return m;
+    }
+
+    // Alter the extent of the view of the existing memory block; usually this
+    // only makes sense when the new extent is smaller than the existing extent,
+    // i.e., prod(new extent) <= prod(old extent)
+    void alterExtent(const uvector<int, N>& ext) { ext_ = ext; }
+
+    // flatten() and slice() return temporary objects which cannot be bound to non-const references;
+    // applying .ref() to the temporary object allows this to be done and is an assertion by the
+    // user that doing so is safe (relating to the output of flatten() or slice() being destroyed
+    // at the end of the full-expression)
+    xarray<T, N>& ref() { return *this; }
+};
 } // namespace algoim
 
 #endif

gjj/myDebug.hpp

@@ -10,6 +10,7 @@
 #include <fstream>
 #include <vector>
 #include "bernstein.hpp"
+#include "multiloop.hpp"
 #include "quadrature_multipoly.hpp"
 #include "binomial.hpp"
 
@@ -92,6 +93,7 @@ void qConvMultiPloy(const F1&              fphi1,
     bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi2(xmin + x * (xmax - xmin)); }, phi2);
 
     // Build quadrature hierarchy
+    auto                      t = std::chrono::system_clock::now();
     ImplicitPolyQuadrature<N> ipquad(phi1, phi2);
     // ImplicitPolyQuadrature<N> ipquad(phi1);
 
@@ -123,6 +125,9 @@ void qConvMultiPloy(const F1&              fphi1,
         surf   *= pow(xmax - xmin, N - 1);
     };
     compute(q);
+    auto                          tAfter          = std::chrono::system_clock::now();
+    std::chrono::duration<double> elapsed_seconds = tAfter - t;
+    std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n";
     std::cout << "q volume: " << volume << std::endl;
 }
 
@@ -135,9 +140,10 @@ void qConvMultiPloy2(real xmin, real xmax, const uvector<int, N>& P, const F& in
     // bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi2(xmin + x * (xmax - xmin)); }, phi2);
 
     std::vector<xarray<real, N>> phis; // 大量sphere
+    // xarray<real, N>* phis = new xarray<real, N>[100];
     for (int i = 0; i < 100; i++) {
         real centerX = 1.0 - rand() % 20 / 10.0, centerY = 1.0 - rand() % 20 / 10.0, centerZ = 1.0 - rand() % 20 / 10.0;
-        real r    = 0.1 + rand() % 9 / 10.0;
+        real r    = 0.01 + rand() % 10 / 100.0;
         auto fphi = [&](const uvector<real, 3>& xx) {
             real x = xx(0) - centerX;
             real y = xx(1) - centerY;
@@ -147,6 +153,7 @@ void qConvMultiPloy2(real xmin, real xmax, const uvector<int, N>& P, const F& in
         xarray<real, N> phi(nullptr, P);
         algoim_spark_alloc(real, phi);
         bernstein::bernsteinInterpolate<N>([&](const uvector<real, N>& x) { return fphi(xmin + x * (xmax - xmin)); }, phi);
+        phis.emplace_back(phi);
     }
     auto t = std::chrono::system_clock::now();
 
@@ -290,7 +297,7 @@ void qConvBernstein(const xarray<real, N>&                  phi,
                     int                                     q,
                     const std::vector<std::array<real, 4>>& halfFaces)
 {
-    uvector<real, 3> testX(0., 0., 0.5);
+    uvector<real, 3> testX(0., 0., 0.25);
     real             testEvalBernstein = bernstein::evalBernsteinPoly(phi, testX);
     auto             vec1              = xarray2StdVector(phi);
     std::cout << "eval bernstein without interpolation:" << testEvalBernstein << std::endl;
@@ -307,10 +314,12 @@ void qConvBernstein(const xarray<real, N>&                  phi,
         if (!halfFaces.empty()) {
             ipquad.integrate(AutoMixed, q, [&](const uvector<real, N>& x, real w) {
                 if (bernstein::evalBernsteinPoly(phi, x) >= 0) return;
-                bool in = true;
+                uvector<real, N> trueX = xmin + x * (xmax - xmin);
+                // std::cout << "trueX: " << trueX << std::endl;
+                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) {
+                    if (sum(factors * trueX) + halfFaces[i][3] >= 0) {
                         in = false;
                         break;
                     }
@@ -378,6 +387,34 @@ void testSpherePowerDirectly()
     qConvBernstein(phiBernstein, -1, 1, integrand, 10, {});
 }
 
+void powerTransformation(const uvector<real, 3>& scale, const uvector<real, 3>& bias, xarray<real, 3>& phiPower)
+{
+    std::vector<std::vector<std::vector<real>>> dimOrderExpansion;
+    const auto&                                 ext = phiPower.ext();
+    for (int dim = 0; dim < 3; ++dim) {
+        dimOrderExpansion.push_back(std::vector<std::vector<real>>(ext(dim)));
+        for (int degree = 0; degree < ext(dim); ++degree) {
+            const real* binomDegree = Binomial::row(degree);
+            dimOrderExpansion[dim][degree].reserve(degree + 1);
+            // 根据二项定理展开
+            for (int i = 0; i <= degree; ++i) {
+                dimOrderExpansion[dim][degree].push_back(binomDegree[i] * pow(scale(dim), i) * pow(bias(dim), degree - i));
+            }
+        }
+    }
+    for (auto i = phiPower.loop(); ~i; ++i) {
+        // 迭代器必须按照坐标的升序进行访问，即，访问ijk时，(i-1)jk，i(j-1)k，ij(k-1)必须已经访问过
+        real base     = phiPower.l(i);
+        phiPower.l(i) = 0;
+        auto exps     = i();
+        for (MultiLoop<3> j(0, exps + 1); ~j; ++j) {
+            real item = base;
+            for (int dim = 0; dim < 3; ++dim) { item *= dimOrderExpansion[dim][exps(dim)][j(dim)]; }
+            phiPower.m(j()) += item;
+        }
+    }
+}
+
 void testSpherePowerDirectlyByTransformation()
 {
     // a_x(x-c_x)^2 + a_y(y-c_y)^2 + a_z(x-c_z)^2 - r^2
@@ -386,49 +423,28 @@ void testSpherePowerDirectlyByTransformation()
     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);
+    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);
-		std::vector<std::vector<real>> dimTransformation;
     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] * (c[dim]) * (c[dim]);
         idx(dim)             = 1;
-        phiPower.m(idx)      = 2 * a[dim] * (- c[dim]);
+        phiPower.m(idx)      = 2 * a[dim] * (-c[dim]);
         idx(dim)             = 2;
         phiPower.m(idx)      = a[dim];
     }
-    phiPower.m(0)  -= r * r;
-
-		for (int dim = 0; dim < 3; ++dim) {
-			dimTransformation[0](std::vector<real>(ext(dim), 0));
-			for (int degree = 0; degree < ext(dim); ++degree) {
-				// dimTransformation[dim][degree] = 0;
-				const real* binomDegree = Binomial::row(degree);
-				for (int i = 0; i <= degree; ++i) {
-					// 根据二项定理展开
-					dimTransformation[dim][i] += binomDegree[i] * pow(k(dim), i) * pow(bias(dim), degree - i);
-				}
-			}
-		}
-		// apply transformation
-		for (auto i = phiPower.loop(); ~i; ++i) {
-			real item = phiPower.l(i);
-			auto exps = i();
-			for (int dim = 0; dim < 3; ++dim) { item += dimTransformation[dim][exps(dim)]; }
-			phiPower.m(i()) = item;
-		}
+    phiPower.m(0) -= r * r;
+    powerTransformation(range, xmin, phiPower);
 
-    auto integrand  = [](const uvector<real, 3>& x) { return 1.0; };
+    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
     power2BernsteinTensorProduct(phiPower, phiBernstein);
 
     v = xarray2StdVector(phiPower);
@@ -457,45 +473,32 @@ void testCylinderPowerDirectly()
     algoim_spark_alloc(real, phiPower, phiBernstein);
     xarrayInit(phiPower);
     xarrayInit(phiBernstein);
-		std::vector<std::vector<real>> dimTransformation;
     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;
-
-		for (int dim = 0; dim < 3; ++dim) {
-			dimTransformation.push_back(std::vector<real>(ext(dim), 0));
-			for (int degree = 0; degree < ext(dim); ++degree) {
-				// dimTransformation[dim][degree] = 0;
-				const real* binomDegree = Binomial::row(degree);
-				for (int i = 0; i <= degree; ++i) {
-					// 根据二项定理展开
-					dimTransformation[dim][i] += binomDegree[i] * pow(k(dim), i) * pow(bias(dim), degree - i);
-				}
-			}
-		}
-
-		// apply transformation
-		for (auto i = phiPower.loop(); ~i; ++i) {
-			real item = phiPower.l(i);
-			auto exps = i();
-			for (int dim = 0; dim < 3; ++dim) { item *= dimTransformation[dim][exps(dim)]; }
-			phiPower.m(i()) = item;
-		}
-
-    auto integrand                       = [](const uvector<real, 3>& x) { return 1.0; };
+    phiPower.m(uvector<int, 3>(2, 0, 0)) = -bottom * top;
+    phiPower.m(uvector<int, 3>(0, 2, 0)) = -bottom * top;
+    phiPower.m(uvector<int, 3>(1, 0, 2)) = 2 * c[0];
+    phiPower.m(uvector<int, 3>(0, 1, 2)) = 2 * c[1];
+    phiPower.m(uvector<int, 3>(1, 0, 1)) = -2 * c[0] * (top + bottom);
+    phiPower.m(uvector<int, 3>(0, 1, 1)) = -2 * c[1] * (top + bottom);
+    phiPower.m(uvector<int, 3>(1, 0, 0)) = 2 * c[0] * (top * bottom);
+    phiPower.m(uvector<int, 3>(0, 1, 0)) = 2 * c[0] * (top * bottom);
+    phiPower.m(uvector<int, 3>(0, 0, 2)) = -(c[0] * c[0] + c[1] * c[1] - r * r);
+    phiPower.m(uvector<int, 3>(0, 0, 1)) = (bottom + top) * (c[0] * c[0] + c[1] * c[1] - r * r);
+    phiPower.m(uvector<int, 3>(0, 0, 0)) = -bottom * top * (c[0] * c[0] + c[1] * c[1] - r * r);
+
+    powerTransformation(range, xmin, phiPower);
+
+    auto integrand = [](const uvector<real, 3>& x) { return 1.0; };
     power2BernsteinTensorProduct(phiPower, phiBernstein);
 
     v = xarray2StdVector(phiPower);
-    uvector<real, 3> testX(0., 0., 0.5);
+    uvector<real, 3> testX(0., 0., 0.75);
     real             testEval = powerEvaluation(phiPower, testX);
     std::cout << "eval power:" << testEval << std::endl;
 
@@ -506,10 +509,10 @@ void testCylinderPowerDirectly()
                    -1,
                    1,
                    integrand,
-                   50,
+                   10,
                    {
-                       {0, 0, 1,  -top},
-                       {0, 0, -1, bottom    }
+                       {0, 0, 1,  -top  },
+                       {0, 0, -1, bottom}
     });
 }
 
@@ -555,9 +558,10 @@ void testBooluarray()
 
 void testMain()
 {
-    testBooluarray();
-    testMultiPolys();
-    testMultiScale();
+    // testBooluarray();
+    // testMultiPolys();
+    // testMultiScale();
+    // testSpherePowerDirectly();
     testSpherePowerDirectlyByTransformation();
-    testCylinderPowerDirectly();
+    // testCylinderPowerDirectly();
 }