matools.h

Go to the documentation of this file.

 
/*
Copyright (c) 2026 MrXie1109
 
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
 
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
 
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
 
#pragma once
#define _USE_MATH_DEFINES
#include <cmath>
#include <complex>
#include <iostream>
#include <numeric>
 
#include "../async/pool.h"
#include "../container/multiarray.h"
#include "../core/sfinae.h"
#include "./random.h"
 
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
 
namespace console {
 
    // ---------------------------- 统计 ----------------------------
    template <class T, size_t... Dims>
    double mean(const MultiArray<T, Dims...> &arr) {
        return double(sum(arr)) / arr.fsize();
    }

    template <class T, size_t... Dims>
    double variance(const MultiArray<T, Dims...> &arr, bool sample = true) {
        double m      = mean(arr);
        double sq_sum = 0;
        arr.for_each([&](const T &x) {
            double d = double(x) - m;
            sq_sum += d * d;
        });
        return sq_sum / (arr.fsize() - (sample ? 1 : 0));
    }

    template <class T, size_t... Dims>
    double stddev(const MultiArray<T, Dims...> &arr, bool sample = true) {
        return std::sqrt(variance(arr, sample));
    }
 
    // ---------------------------- 向量运算(一维) ----------------------------
    template <class T, size_t N>
    T dot(const MultiArray<T, N> &a, const MultiArray<T, N> &b) {
        return std::inner_product(a.fbegin(), a.fend(), b.fbegin(), T{});
    }

    template <class T, size_t N>
    double norm(const MultiArray<T, N> &a) {
        return std::sqrt(double(dot(a, a)));
    }

    template <class T, size_t N>
    double cosine(const MultiArray<T, N> &a, const MultiArray<T, N> &b) {
        return dot(a, b) / (norm(a) * norm(b));
    }

    template <class T, size_t N>
    MultiArray<T, N> normalize(const MultiArray<T, N> &a) {
        double len = norm(a);
        if (std::abs(len) < std::numeric_limits<double>::epsilon()) return a;
        MultiArray<T, N> result;
        auto             ai = a.fbegin();
        auto             ri = result.fbegin();
        while (ri != result.fend()) *ri++ = *ai++ / len;
        return result;
    }

    template <class T, size_t N>
    double euclidean(const MultiArray<T, N> &a, const MultiArray<T, N> &b) {
        double sum = 0;
        auto   ai  = a.fbegin();
        auto   bi  = b.fbegin();
        while (ai != a.fend()) {
            double d = *ai++ - *bi++;
            sum += d * d;
        }
        return std::sqrt(sum);
    }

    template <class T, size_t N>
    double manhattan(const MultiArray<T, N> &a, const MultiArray<T, N> &b) {
        double sum = 0;
        auto   ai  = a.fbegin();
        auto   bi  = b.fbegin();
        while (ai != a.fend()) sum += std::abs(*ai++ - *bi++);
        return sum;
    }
 
    // ---------------------------- 矩阵运算(二维) ----------------------------
#ifndef CONSOLE_MATMUL_BLOCK_SIZE
#define CONSOLE_MATMUL_BLOCK_SIZE 128
#endif

    template <class T, size_t M, size_t N, size_t K>
    MultiArray<T, M, K>
    matmul_blocked(const MultiArray<T, M, N> &A, const MultiArray<T, N, K> &B) {
        using std::min;
        MultiArray<T, M, K> C(T{});
        constexpr size_t    BLOCK = CONSOLE_MATMUL_BLOCK_SIZE;
        for (size_t i0 = 0; i0 < M; i0 += BLOCK)
            for (size_t j0 = 0; j0 < K; j0 += BLOCK)
                for (size_t k0 = 0; k0 < N; k0 += BLOCK)
                    for (size_t i = i0; i < min(i0 + BLOCK, M); ++i)
                        for (size_t j = j0; j < min(j0 + BLOCK, K); ++j)
                            for (size_t k = k0; k < min(k0 + BLOCK, N); ++k)
                                C[i][j] += A[i][k] * B[k][j];
        return C;
    }

    template <class T, size_t M, size_t N, size_t K>
    MultiArray<T, M, K>
    matmul_simple(const MultiArray<T, M, N> &A, const MultiArray<T, N, K> &B) {
        MultiArray<T, M, K> C(T{});
        for (size_t i = 0; i < M; ++i)
            for (size_t j = 0; j < K; ++j)
                for (size_t k = 0; k < N; ++k) C[i][j] += A[i][k] * B[k][j];
        return C;
    }

    template <class T, size_t M, size_t N, size_t K>
    inline MultiArray<T, M, K> matmul_impl(const MultiArray<T, M, N> &A,
        const MultiArray<T, N, K>                                    &B,
        std::true_type) {
        return matmul_blocked(A, B);
    }

    template <class T, size_t M, size_t N, size_t K>
    inline MultiArray<T, M, K> matmul_impl(const MultiArray<T, M, N> &A,
        const MultiArray<T, N, K>                                    &B,
        std::false_type) {
        return matmul_simple(A, B);
    }

    template <class T, size_t M, size_t N, size_t K>
    inline MultiArray<T, M, K>
    matmul(const MultiArray<T, M, N> &A, const MultiArray<T, N, K> &B) {
        constexpr size_t BLOCK       = CONSOLE_MATMUL_BLOCK_SIZE;
        constexpr bool   use_blocked = M >= BLOCK && N >= BLOCK && K >= BLOCK;
        return matmul_impl(A, B, std::integral_constant<bool, use_blocked>{});
    }

    template <class T, size_t M, size_t N, size_t K>
    MultiArray<T, M, K> matmul_parallel(
        const MultiArray<T, M, N> &A, const MultiArray<T, N, K> &B) {
        using std::min;
        MultiArray<T, M, K>            C(T{});
        constexpr size_t               BLOCK        = CONSOLE_MATMUL_BLOCK_SIZE;
        const size_t                   num_i_blocks = (M + BLOCK - 1) / BLOCK;
        const size_t                   num_j_blocks = (K + BLOCK - 1) / BLOCK;
        std::vector<std::future<void>> futures;
        futures.reserve(num_i_blocks * num_j_blocks);
        for (size_t i0_idx = 0; i0_idx < num_i_blocks; ++i0_idx)
            for (size_t j0_idx = 0; j0_idx < num_j_blocks; ++j0_idx) {
                size_t i0 = i0_idx * BLOCK;
                size_t j0 = j0_idx * BLOCK;
                futures.push_back(pool::submit([&, i0, j0]() {
                    for (size_t k0 = 0; k0 < N; k0 += BLOCK)
                        for (size_t i = i0; i < min(i0 + BLOCK, M); ++i)
                            for (size_t j = j0; j < min(j0 + BLOCK, K); ++j) {
                                T sum = C[i][j];
                                for (size_t k = k0; k < min(k0 + BLOCK, N); ++k)
                                    sum += A[i][k] * B[k][j];
                                C[i][j] = sum;
                            }
                }));
            }
        for (auto &f : futures) f.get();
        return C;
    }

    template <class T, size_t M, size_t N>
    MultiArray<T, N, M> transpose(const MultiArray<T, M, N> &A) {
        MultiArray<T, N, M> B;
        for (size_t i = 0; i < M; i++)
            for (size_t j = 0; j < N; j++) B[j][i] = A[i][j];
        return B;
    }

    template <class T, size_t N>
    MultiArray<T, N, N> identity() {
        MultiArray<T, N, N> I(T{});
        for (size_t i = 0; i < N; i++) I[i][i] = T{1};
        return I;
    }

    template <class T, size_t N>
    T trace(const MultiArray<T, N, N> &A) {
        T result{};
        for (size_t i = 0; i < N; i++) result += A[i][i];
        return result;
    }

    template <class T>
    MultiArray<T, 3>
    cross(const MultiArray<T, 3> &a, const MultiArray<T, 3> &b) {
        return MultiArray<T, 3>{a[1] * b[2] - a[2] * b[1],
            a[2] * b[0] - a[0] * b[2],
            a[0] * b[1] - a[1] * b[0]};
    }
 
    // ---------------------------- 元素级运算 ----------------------------
    template <class T, size_t... Dims>
    MultiArray<T, Dims...>
    clamp(const MultiArray<T, Dims...> &arr, T low, T high) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) {
            *ri = *ai < low ? low : (*ai > high ? high : *ai);
            ++ri;
            ++ai;
        }
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> abs(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::abs(*ai++);
        return result;
    }
 
    // ---------------------------- 随机初始化 ----------------------------
    template <class T, size_t... Dims>
    void randomize(MultiArray<T, Dims...> &arr, T min = 0, T max = 100) {
        uniform_distribution_t<T> dis(min, max);
        arr.for_each([&](T &x) { x = dis(thread_rng()); });
    }

    template <class T, size_t... Dims>
    void
    randomize_normal(MultiArray<T, Dims...> &arr, T mean = 0, T stddev = 1) {
        std::normal_distribution<T> dis(mean, stddev);
        arr.for_each([&](T &x) { x = dis(thread_rng()); });
    }

    template <class T, size_t... Dims>
    void linspace(MultiArray<T, Dims...> &arr, T start, T end) {
        size_t n = arr.fsize();
        for (size_t i = 0; i < n; i++)
            arr.fbegin()[i] = start + (end - start) * i / (n - 1);
    }
 
    // ---------------------------- 其他 ----------------------------
    template <class T, size_t... Dims>
    T product(const MultiArray<T, Dims...> &arr) {
        return std::accumulate(
            arr.fbegin(), arr.fend(), T{1}, std::multiplies<T>());
    }

    template <class T, size_t... Dims>
    T kth_smallest(MultiArray<T, Dims...> arr, size_t k) {
        if (k >= arr.fsize()) k = arr.fsize() - 1;
        std::nth_element(arr.fbegin(), arr.fbegin() + k, arr.fend());
        return arr.fbegin()[k];
    }

    template <class T, size_t N, size_t K>
    MultiArray<T, N + K - 1>
    convolve1d(const MultiArray<T, N> &signal, const MultiArray<T, K> &kernel) {
        MultiArray<T, N + K - 1> result(T{});
        for (size_t i = 0; i < N; i++) {
            T sum{};
            for (size_t j = 0; j < K; j++)
                result[i + j] += signal[i] * kernel[j];
        }
        return result;
    }
 
    // ---------------------------- 三角函数 ----------------------------
    template <class T, size_t... Dims>
    MultiArray<T, Dims...> sin(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::sin(*ai++);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> cos(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::cos(*ai++);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> tan(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::tan(*ai++);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> asin(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::asin(*ai++);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> acos(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::acos(*ai++);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> atan(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::atan(*ai++);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> sinh(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::sinh(*ai++);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> cosh(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::cosh(*ai++);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> tanh(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::tanh(*ai++);
        return result;
    }
 
    // ---------------------------- 指数对数 ----------------------------
    template <class T, size_t... Dims>
    MultiArray<T, Dims...> exp(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::exp(*ai++);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> log(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::log(*ai++);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> log10(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::log10(*ai++);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> pow(const MultiArray<T, Dims...> &arr, T exponent) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::pow(*ai++, exponent);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...>
    pow(const MultiArray<T, Dims...> &base, const MultiArray<T, Dims...> &exp) {
        MultiArray<T, Dims...> result;
        auto                   bi = base.fbegin();
        auto                   ei = exp.fbegin();
        auto                   ri = result.fbegin();
        while (ri != result.fend()) *ri++ = std::pow(*bi++, *ei++);
        return result;
    }
 
    // ---------------------------- 取整函数 ----------------------------
    template <class T, size_t... Dims>
    MultiArray<T, Dims...> floor(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::floor(*ai++);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> ceil(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::ceil(*ai++);
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> round(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) *ri++ = std::round(*ai++);
        return result;
    }

    template < //
        class T,
        size_t... Dims,
        class CharT,
        class Traits = std::char_traits<CharT>>
    void print_stats(std::basic_ostream<CharT, Traits> &os,
        const MultiArray<T, Dims...>                   &arr,
        const std::string                              &name = "") {
        if (!name.empty()) os << "=== " << name << " ===" << '\n';
        os << "  sum   : " << sum(arr) << '\n';
        os << "  mean  : " << mean(arr) << '\n';
        os << "  min   : " << min(arr) << '\n';
        os << "  max   : " << max(arr) << '\n';
        os << "  stddev: " << stddev(arr) << std::endl;
    }

    template <class T, size_t... Dims>
    double covariance(
        const MultiArray<T, Dims...> &a, const MultiArray<T, Dims...> &b) {
        double mean_a = mean(a);
        double mean_b = mean(b);
        double sum    = 0;
        auto   ai     = a.fbegin();
        auto   bi     = b.fbegin();
        while (ai != a.fend()) {
            sum += (double(*ai) - mean_a) * (double(*bi) - mean_b);
            ++ai;
            ++bi;
        }
        return sum / (a.fsize() - 1);
    }

    template <class T, size_t... Dims>
    double correlation(
        const MultiArray<T, Dims...> &a, const MultiArray<T, Dims...> &b) {
        double cov   = covariance(a, b);
        double std_a = stddev(a);
        double std_b = stddev(b);
        if (std_a == 0 || std_b == 0) return 0;
        return cov / (std_a * std_b);
    }

    template <class T, size_t... Dims>
    double median(MultiArray<T, Dims...> arr) {
        size_t n     = arr.fsize();
        auto   begin = arr.fbegin();
        auto   end   = arr.fend();
        std::nth_element(begin, begin + n / 2, end);
        if (n % 2 == 1) {
            return double(arr.fbegin()[n / 2]);
        } else {
            double a = double(arr.fbegin()[n / 2 - 1]);
            double b = double(arr.fbegin()[n / 2]);
            return (a + b) / 2;
        }
    }

    enum class QuantileMethod {
        Nearest,  
        Linear,   
        Midpoint, 
        Lower,    
        Upper,    
        ExcelInc, 
        ExcelExc  
    };

    template <class T, size_t... Dims>
    double quantile(MultiArray<T, Dims...> arr,
        double                             q,
        QuantileMethod                     method = QuantileMethod::Nearest) {
        size_t n = arr.fsize();
        if (n == 0) return 0;
        if (q < 0) q = 0;
        if (q > 1) q = 1;
        auto begin = arr.fbegin();
        auto end   = arr.fend();
        switch (method) {
        case QuantileMethod::Nearest: {
            size_t idx = size_t(q * (n - 1) + 0.5);
            if (idx >= n) idx = n - 1;
            std::nth_element(begin, begin + idx, end);
            return double(arr.fbegin()[idx]);
        }
        case QuantileMethod::Linear: {
            double pos  = q * (n - 1);
            size_t idx  = size_t(pos);
            double frac = pos - idx;
            if (idx >= n - 1) {
                std::nth_element(begin, begin + n - 1, end);
                return double(arr.fbegin()[n - 1]);
            }
            std::nth_element(begin, begin + idx + 1, end);
            double a = double(arr.fbegin()[idx]);
            double b = double(arr.fbegin()[idx + 1]);
            return a + (b - a) * frac;
        }
        case QuantileMethod::Midpoint: {
            double pos = q * (n - 1);
            size_t idx = size_t(pos);
            if (idx >= n - 1) {
                std::nth_element(begin, begin + n - 1, end);
                return double(arr.fbegin()[n - 1]);
            }
            std::nth_element(begin, begin + idx + 1, end);
            double a = double(arr.fbegin()[idx]);
            double b = double(arr.fbegin()[idx + 1]);
            return (a + b) / 2;
        }
        case QuantileMethod::Lower: {
            size_t idx = size_t(q * (n - 1));
            if (idx >= n) idx = n - 1;
            std::nth_element(begin, begin + idx, end);
            return double(arr.fbegin()[idx]);
        }
        case QuantileMethod::Upper: {
            size_t idx = size_t(q * (n - 1));
            if (idx >= n - 1) return double(arr.fbegin()[n - 1]);
            idx++;
            std::nth_element(begin, begin + idx, end);
            return double(arr.fbegin()[idx]);
        }
        case QuantileMethod::ExcelInc: {
            double pos  = (n - 1) * q;
            size_t idx  = size_t(pos);
            double frac = pos - idx;
            if (idx >= n - 1) {
                return double(arr.fbegin()[n - 1]);
            }
            std::nth_element(begin, begin + idx + 1, end);
            double a = double(arr.fbegin()[idx]);
            double b = double(arr.fbegin()[idx + 1]);
            return a + (b - a) * frac;
        }
        case QuantileMethod::ExcelExc: {
            double pos  = (n + 1) * q;
            size_t idx  = size_t(pos);
            double frac = pos - idx;
            if (idx < 1) idx = 1;
            if (idx >= n) {
                if (idx == n) {
                    std::nth_element(begin, begin + n - 1, end);
                    return double(arr.fbegin()[n - 1]);
                }
                return double(arr.fbegin()[n - 1]);
            }
            std::nth_element(begin, begin + idx, end);
            double a = double(arr.fbegin()[idx - 1]);
            double b = double(arr.fbegin()[idx]);
            return a + (b - a) * frac;
        }
        default:
            return quantile(arr, q, QuantileMethod::Nearest);
        }
    }

    template <class T>
    T det2x2(const MultiArray<T, 2, 2> &A) {
        return A[0][0] * A[1][1] - A[0][1] * A[1][0];
    }

    template <class T>
    T det3x3(const MultiArray<T, 3, 3> &A) {
        return A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1])
               - A[0][1] * (A[1][0] * A[2][2] - A[1][2] * A[2][0])
               + A[0][2] * (A[1][0] * A[2][1] - A[1][1] * A[2][0]);
    }

    template <class T>
    MultiArray<T, 2, 2> inv2x2(const MultiArray<T, 2, 2> &A) {
        T d = det2x2(A);
        if (d == 0) return MultiArray<T, 2, 2>{};
        return MultiArray<T, 2, 2>{
            A[1][1] / d, -A[0][1] / d, -A[1][0] / d, A[0][0] / d};
    }

    template <class T>
    MultiArray<T, 3, 3> inv3x3(const MultiArray<T, 3, 3> &A) {
        T d = det3x3(A);
        if (d == 0) return MultiArray<T, 3, 3>{};
        MultiArray<T, 3, 3> B;
        B[0][0] = (A[1][1] * A[2][2] - A[1][2] * A[2][1]) / d;
        B[0][1] = (A[0][2] * A[2][1] - A[0][1] * A[2][2]) / d;
        B[0][2] = (A[0][1] * A[1][2] - A[0][2] * A[1][1]) / d;
        B[1][0] = (A[1][2] * A[2][0] - A[1][0] * A[2][2]) / d;
        B[1][1] = (A[0][0] * A[2][2] - A[0][2] * A[2][0]) / d;
        B[1][2] = (A[0][2] * A[1][0] - A[0][0] * A[1][2]) / d;
        B[2][0] = (A[1][0] * A[2][1] - A[1][1] * A[2][0]) / d;
        B[2][1] = (A[0][1] * A[2][0] - A[0][0] * A[2][1]) / d;
        B[2][2] = (A[0][0] * A[1][1] - A[0][1] * A[1][0]) / d;
        return B;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> sigmoid(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result(T{});
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) {
            *ri++ = T(1) / (T(1) + std::exp(-*ai++));
        }
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> relu(const MultiArray<T, Dims...> &arr) {
        MultiArray<T, Dims...> result;
        auto                   ri = result.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != result.fend()) {
            *ri++ = *ai < 0 ? 0 : *ai;
            ++ai;
        }
        return result;
    }

    template <class T, size_t... Dims>
    MultiArray<T, Dims...> softmax(const MultiArray<T, Dims...> &arr) {
        auto                   max_val = max(arr);
        MultiArray<T, Dims...> exp_arr;
        auto                   ri = exp_arr.fbegin();
        auto                   ai = arr.fbegin();
        while (ri != exp_arr.fend()) {
            *ri++ = std::exp(*ai++ - max_val);
        }
        T                      sum_exp = sum(exp_arr);
        MultiArray<T, Dims...> result;
        ri = result.fbegin();
        ai = exp_arr.fbegin();
        while (ri != result.fend()) {
            *ri++ = *ai++ / sum_exp;
        }
        return result;
    }

    template <class T, size_t H, size_t W, size_t KH, size_t KW>
    MultiArray<T, H - KH + 1, W - KW + 1> convolve2d(
        const MultiArray<T, H, W> &image, const MultiArray<T, KH, KW> &kernel) {
        constexpr size_t      OH = H - KH + 1;
        constexpr size_t      OW = W - KW + 1;
        MultiArray<T, OH, OW> result(T{});
        for (size_t i = 0; i < OH; ++i) {
            for (size_t j = 0; j < OW; ++j) {
                T sum{};
                for (size_t ki = 0; ki < KH; ++ki) {
                    for (size_t kj = 0; kj < KW; ++kj) {
                        sum += image[i + ki][j + kj] * kernel[ki][kj];
                    }
                }
                result[i][j] = sum;
            }
        }
        return result;
    }

    template <class T, size_t M, size_t N>
    void print_matrix(std::ostream &os,
        const MultiArray<T, M, N>  &A,
        const char                 *name = nullptr) {
        if (name) os << name << " = \n";
        for (size_t i = 0; i < M; ++i) {
            os << "  [";
            for (size_t j = 0; j < N; ++j) {
                repr(A[i][j], os);
                if (j < N - 1) os << ", ";
            }
            os << "]\n";
        }
    }

    template <class T, size_t N>
    MultiArray<std::complex<T>, N>
    fft(const MultiArray<std::complex<T>, N> &x) {
        static_assert((N & (N - 1)) == 0, "FFT length must be a power of 2");
        MultiArray<std::complex<T>, N> result = x;
        for (size_t i = 0; i < N; ++i) {
            size_t j = 0;
            for (size_t k = 1; k < N; k <<= 1) {
                if (i & k)
                    j = (j << 1) | 1;
                else
                    j <<= 1;
            }
            if (i < j) std::swap(result[i], result[j]);
        }
        for (size_t len = 2; len <= N; len <<= 1) {
            T               angle = T(-2 * M_PI) / len;
            std::complex<T> wlen(std::cos(angle), std::sin(angle));
            for (size_t i = 0; i < N; i += len) {
                std::complex<T> w = 1;
                for (size_t j = 0; j < len / 2; ++j) {
                    auto u                  = result[i + j];
                    auto v                  = result[i + j + len / 2] * w;
                    result[i + j]           = u + v;
                    result[i + j + len / 2] = u - v;
                    w *= wlen;
                }
            }
        }
        return result;
    }

    template <class T, size_t N>
    MultiArray<std::complex<T>, N>
    ifft(const MultiArray<std::complex<T>, N> &X) {
        MultiArray<std::complex<T>, N> conj_X;
        for (size_t i = 0; i < N; ++i) conj_X[i] = std::conj(X[i]);
 
        auto result = fft(conj_X);
 
        for (size_t i = 0; i < N; ++i) {
            result[i] = std::conj(result[i]) / T(N);
        }
 
        return result;
    }

    template <class T, size_t N>
    MultiArray<std::complex<T>, N> rfft(const MultiArray<T, N> &x) {
        MultiArray<std::complex<T>, N> cx;
        for (size_t i = 0; i < N; ++i) cx[i] = std::complex<T>(x[i], 0);
        return fft(cx);
    }

    template <class T, size_t N>
    MultiArray<T, N> irfft(const MultiArray<std::complex<T>, N> &X) {
        auto             cx = ifft(X);
        MultiArray<T, N> result;
        for (size_t i = 0; i < N; ++i) result[i] = cx[i].real();
        return result;
    }

    template <class T, size_t N>
    MultiArray<T, N / 2 + 1> spectrum(const MultiArray<T, N> &x) {
        auto                     f = rfft(x);
        MultiArray<T, N / 2 + 1> mag(T{});
        const T                  norm = T(1) / T(N);
        for (size_t i = 0; i <= N / 2; ++i) mag[i] = std::abs(f[i]) * norm;
        return mag;
    }

    template <class T, size_t N>
    MultiArray<T, N / 2 + 1> power_spectrum(const MultiArray<T, N> &x) {
        auto mag = spectrum(x);
        for (size_t i = 0; i < mag.fsize(); ++i) {
            mag[i] = mag[i] * mag[i];
        }
        return mag;
    }

    template <class T, size_t N>
    MultiArray<T, N / 2 + 1> phase_spectrum(const MultiArray<T, N> &x) {
        auto                     f = rfft(x);
        MultiArray<T, N / 2 + 1> phase;
        for (size_t i = 0; i <= N / 2; ++i) {
            phase[i] = std::arg(f[i]);
        }
        return phase;
    }

    template <class T, size_t N>
    MultiArray<T, N / 2 + 1> log_spectrum(const MultiArray<T, N> &x) {
        auto                     mag = spectrum(x);
        MultiArray<T, N / 2 + 1> logmag;
        for (size_t i = 0; i < mag.fsize(); ++i) {
            T val     = mag[i];
            logmag[i] = (val > 0) ? T(20) * std::log10(val) : T(-100);
        }
        return logmag;
    }

    template <class T, size_t M, size_t N>
    MultiArray<T, M>
    matvec(const MultiArray<T, M, N> &A, const MultiArray<T, N> &x) {
        MultiArray<T, M> y(T{});
        for (size_t i = 0; i < M; ++i) {
            T sum{};
            for (size_t j = 0; j < N; ++j) {
                sum += A[i][j] * x[j];
            }
            y[i] = sum;
        }
        return y;
    }

    template <class T, size_t M, size_t N>
    MultiArray<T, M, N>
    outer(const MultiArray<T, M> &a, const MultiArray<T, N> &b) {
        MultiArray<T, M, N> result(T{});
        for (size_t i = 0; i < M; ++i) {
            for (size_t j = 0; j < N; ++j) {
                result[i][j] = a[i] * b[j];
            }
        }
        return result;
    }

    template <class T, size_t M, size_t N>
    double frobenius_norm(const MultiArray<T, M, N> &A) {
        double sum{};
        for (auto it = A.fbegin(); it != A.fend(); ++it) {
            double v = double(*it);
            sum += v * v;
        }
        return std::sqrt(sum);
    }

    template <class T, size_t M, size_t N>
    T column_norm(const MultiArray<T, M, N> &A) {
        T max_col{};
        for (size_t j = 0; j < N; ++j) {
            T col_sum{};
            for (size_t i = 0; i < M; ++i) {
                col_sum += std::abs(A[i][j]);
            }
            if (j == 0 || col_sum > max_col) max_col = col_sum;
        }
        return max_col;
    }

    template <class T, size_t M, size_t N>
    T row_norm(const MultiArray<T, M, N> &A) {
        T max_row{};
        for (size_t i = 0; i < M; ++i) {
            T row_sum{};
            for (size_t j = 0; j < N; ++j) {
                row_sum += std::abs(A[i][j]);
            }
            if (i == 0 || row_sum > max_row) max_row = row_sum;
        }
        return max_row;
    }

    template <class T, size_t N>
    MultiArray<T, N> hann_window() {
        MultiArray<T, N> win;
        for (size_t i = 0; i < N; ++i) {
            win[i] = T(0.5) * (T(1) - std::cos(T(2 * M_PI) * i / (N - 1)));
        }
        return win;
    }

    template <class T, size_t N>
    MultiArray<T, N> hamming_window() {
        MultiArray<T, N> win;
        for (size_t i = 0; i < N; ++i) {
            win[i] = T(0.54) - T(0.46) * std::cos(T(2 * M_PI) * i / (N - 1));
        }
        return win;
    }

    template <class T, size_t N>
    MultiArray<T, N> blackman_window() {
        MultiArray<T, N> win;
        for (size_t i = 0; i < N; ++i) {
            T a    = T(2 * M_PI) * i / (N - 1);
            win[i] = T(0.42) - T(0.5) * std::cos(a) + T(0.08) * std::cos(2 * a);
        }
        return win;
    }

    template <class T, size_t N>
    MultiArray<T, N> fftshift(const MultiArray<T, N> &x) {
        MultiArray<T, N> result;
        size_t           half = N / 2;
        for (size_t i = 0; i < half; ++i) {
            result[i]        = x[i + half];
            result[i + half] = x[i];
        }
        return result;
    }

    template <class T, size_t N, size_t M>
    MultiArray<T, N + M - 1>
    cross_correlation(const MultiArray<T, N> &a, const MultiArray<T, M> &b) {
        MultiArray<T, N + M - 1> result(T{});
        for (size_t i = 0; i < N; ++i) {
            for (size_t j = 0; j < M; ++j) {
                result[i + j] += a[i] * b[j];
            }
        }
        return result;
    }
 // end of matools group
}