From 7e4fc45c80ef76759e0df227719ed1aabdb7db94 Mon Sep 17 00:00:00 2001
From: Feng <wang_feng@live.com>
Date: Thu, 14 Mar 2024 20:59:55 +0100
Subject: [PATCH 1/5] add concept Matrix

---
 ReadMe.md  |  8 ++++++++
 matrix.hpp | 48 +++++++++++++-----------------------------------
 2 files changed, 21 insertions(+), 35 deletions(-)

diff --git a/ReadMe.md b/ReadMe.md
index ede5e1d..809e40c 100644
--- a/ReadMe.md
+++ b/ReadMe.md
@@ -142,6 +142,14 @@ Variadic macro `__VA_OPT__` is used. It is officially supported since c++20([lin
     auto one = feng::empty<double>(12, 34);
     ```
 
+    - converting value types
+
+    ```cpp
+    auto one_fp32 = feng::empty<float>(12, 34);
+    auto one_fp64 = one_fp32.astype<double>(); //
+    ```
+
+
     - creating a matrix of size `1 X (12*34)`, with all elements from `0` to `12x34`, then reshape to `12 X 34`:
     ```cpp
     auto one = feng::arange<double>(12*34);
diff --git a/matrix.hpp b/matrix.hpp
index 21508e5..a68a528 100644
--- a/matrix.hpp
+++ b/matrix.hpp
@@ -127,6 +127,7 @@ namespace feng
     template< typename T>
     concept Allocator = has_value_type_v<T> && has_allocate_v<T> && has_deallocate_v<T>;
 
+
     //
     // end of concept allocators
     //
@@ -240,41 +241,6 @@ namespace feng
             }
         };
 
-        /*
-        template< typename Integer_Type >
-        struct integer_range
-        {
-            static_assert(std::is_integral_v<Integer_Type>, "Integral required");
-            typedef Integer_Type value_type;
-            value_type first_;
-            value_type last_;
-
-            integer_range( value_type first, value_type last ) noexcept : first_{first}, last_{last} {}
-
-            auto const begin() const noexcept
-            {
-                return integer_iterator<value_type>{first_};
-            }
-
-            auto const end() const noexcept
-            {
-                return integer_iterator<value_type>{last_};
-            }
-        };
-
-        template< typename Integer_Type >
-        integer_range<Integer_Type> range( Integer_Type first, Integer_Type last ) noexcept
-        {
-            return {first, last};
-        }
-
-        template< typename Integer_Type >
-        integer_range<Integer_Type> range( Integer_Type last ) noexcept
-        {
-            return { Integer_Type{0}, last };
-        }
-        */
-
 
         template< std::weakly_incrementable W >
         constexpr auto range( W val_begin, W val_end )
@@ -3944,6 +3910,18 @@ namespace feng
 
     };//struct matrix
 
+    template< typename T >
+    struct is_matrix : std::false_type{};
+
+    template< typename T, Allocator A >
+    struct is_matrix< matrix<T,A> > : std::true_type{};
+
+    template< typename M >
+    inline constexpr bool is_matrix_v = is_matrix<M>::value;
+
+    template< typename M >
+    concept Matrix = is_matrix_v<M>;
+
     namespace matrix_details
     {
 

From 960a050b5afcfa8f86f94c13be8ea4bcdf2e4526 Mon Sep 17 00:00:00 2001
From: Feng <wang_feng@live.com>
Date: Thu, 14 Mar 2024 21:26:59 +0100
Subject: [PATCH 2/5] partial conceptize

---
 matrix.hpp | 128 +++++++++++++++++++++++++++--------------------------
 1 file changed, 66 insertions(+), 62 deletions(-)

diff --git a/matrix.hpp b/matrix.hpp
index a68a528..7301170 100644
--- a/matrix.hpp
+++ b/matrix.hpp
@@ -6152,11 +6152,10 @@ namespace feng
         return value_type( 0 );
     }
     */
-    template < typename Matrix >
-    typename Matrix::value_type
-    norm_1( const Matrix& A )
+    template < Matrix Mat >
+    auto norm_1( Mat const& A )
     {
-        typedef typename Matrix::value_type value_type;
+        typedef typename Mat::value_type value_type;
         std::vector< value_type > m( A.col() );
 
         for ( std::uint_least64_t i = 0; i != A.col(); ++i )
@@ -6166,6 +6165,7 @@ namespace feng
         } );
         return *( std::max_element( m.begin(), m.end() ) );
     }
+
     template < typename T, typename A_ >
     T norm_1( const matrix< std::complex< T >, A_ >& A )
     {
@@ -6178,12 +6178,12 @@ namespace feng
         } );
         return *( std::max_element( m.begin(), m.end() ) );
     }
-    template < typename Matrix >
-    typename Matrix::value_type
-    norm_2( const Matrix& A )
+    template < Matrix Mat >
+    auto norm_2( Mat const& A )
     {
         return std::sqrt( eigen_power_iteration( A ) );
     }
+
     namespace expm_private
     {
         template < typename T >
@@ -6197,11 +6197,10 @@ namespace feng
             typedef typename fix_complex_value_type< T >::value_type value_type;
         };
     }
-    template < typename T, typename A_ >
-    const matrix< T, A_ >
-    expm( const matrix< T, A_ >& A )
+    template < Matrix Mat >
+    Mat expm( Mat const& A )
     {
-        typedef matrix< T, A_ > matrix_type;
+        typedef Mat matrix_type;
         typedef typename matrix_type::value_type value_type_;
         typedef typename expm_private::fix_complex_value_type< value_type_ >::value_type value_type;
         typedef typename matrix_type::size_type size_type;
@@ -6270,10 +6269,12 @@ namespace feng
             typedef std::complex< T > result_type;
         };
     }
-    template < typename T >
-    auto fft( matrix< T > const& x )
+
+    template < Matrix Mat >
+    auto fft( Mat const& x )
     {
-        typedef typename fft_private::add_complex< T >::result_type complex_type;
+        typedef typename Mat::value_type value_type;
+        typedef typename fft_private::add_complex< value_type >::result_type complex_type;
         matrix< complex_type > X( x.row(), x.col() );
         auto make_omege = []( auto k, auto n, auto N )
         {
@@ -6304,8 +6305,9 @@ namespace feng
 
         return X;
     }
-    template < typename T >
-    auto fftshift( matrix< T > const& x )
+
+    template < Matrix Mat >
+    auto fftshift( Mat const& x )
     {
         auto X                          = fft( x );
         std::uint_least64_t const R           = X.row();
@@ -6322,10 +6324,11 @@ namespace feng
 
         return X;
     }
-    template < typename T1, Allocator A1, typename T2, Allocator A2, typename T3, Allocator A3 >
-    int forward_substitution( const matrix< T1, A1 >& A, matrix< T2, A2 >& x, const matrix< T3, A3 >& b )
+
+    template < Matrix Mat >
+    int forward_substitution( Mat const& A, Mat& x, Mat const& b )
     {
-        typedef matrix< T1, A1 > matrix_type;
+        typedef Mat matrix_type;
         typedef typename matrix_type::value_type value_type;
         typedef typename matrix_type::size_type size_type;
         better_assert( A.row() == A.col() );
@@ -6347,8 +6350,8 @@ namespace feng
         return 0;
     }
 
-    template < typename T, Allocator A>
-    std::optional<matrix<T,A>> gauss_jordan_elimination( matrix< T, A > const& m ) noexcept
+    template < Matrix Mat >
+    std::optional<Mat> gauss_jordan_elimination( Mat const& m ) noexcept
     {
         auto const& [row, col] = m.shape();
         better_assert( row < col && "matrix row must be less than colum to execut a Gauss-Jordan Elimination" );
@@ -6381,8 +6384,8 @@ namespace feng
     }
 
     // matlab alias
-    template < typename T, Allocator A>
-    std::optional<matrix<T,A>> rref( matrix< T, A > const& m ) noexcept
+    template < Matrix Mat >
+    std::optional<Mat> rref( Mat const& m ) noexcept
     {
         return gauss_jordan_elimination( m );
     }
@@ -6400,8 +6403,8 @@ namespace feng
             typedef std::complex< T > result_type;
         };
     }
-    template < typename T >
-    auto ifft( matrix< T > const& x )
+    template < typename T, Allocator A >
+    auto ifft( matrix<T, A> const& x )
     {
         typedef typename ifft_private::add_complex< T >::result_type complex_type;
         matrix< complex_type > X( x.row(), x.col() );
@@ -6434,8 +6437,8 @@ namespace feng
 
         return X;
     }
-    template < typename T >
-    auto ifftshift( matrix< T > const& x )
+    template < Matrix Mat >
+    auto ifftshift( Mat const& x )
     {
         auto X                          = ifft( x );
         std::uint_least64_t const R           = X.row();
@@ -6453,10 +6456,10 @@ namespace feng
         return X;
     }
 
-    template < typename Type, Allocator Alloc>
-    int lu_decomposition( const matrix< Type, Alloc >& A, matrix< Type, Alloc >& L, matrix< Type, Alloc >& U )
+    template< Matrix Mat >
+    int lu_decomposition( Mat const& A, Mat& L, Mat& U )
     {
-        typedef Type value_type;
+        typedef typename Mat::value_type value_type;
         better_assert( A.row() == A.col() && "Square Matrix Requred!" );
 
         const std::uint_least64_t n = A.row();
@@ -6486,19 +6489,19 @@ namespace feng
         return 0;
     }
 
-    template < typename Type, Allocator Alloc>
-    std::optional<std::tuple<matrix<Type, Alloc>, matrix<Type, Alloc>>> lu_decomposition( const matrix< Type, Alloc >& A )
+    template< Matrix Mat >
+    std::optional<std::tuple<Mat, Mat>> lu_decomposition( Mat const& A )
     {
-        if ( matrix<Type, Alloc> L, U; lu_decomposition( A, L, U ) == 1 )
+        if ( Mat L, U; lu_decomposition( A, L, U ) == 1 )
             return {};
         else
             return std::make_tuple( L, U );
     }
 
-    template < typename Type, Allocator Alloc >
-    int lu_solver( const matrix< Type, Alloc >& A, matrix< Type, Alloc >& x, const matrix< Type, Alloc >& b )
+    template< Matrix Mat >
+    int lu_solver( Mat const& A, Mat& x, Mat const& b )
     {
-        typedef matrix< Type, Alloc > matrix_type;
+        typedef Mat matrix_type;
         better_assert( A.row() == A.col() );
         better_assert( A.row() == b.row() );
         better_assert( b.col() == 1 );
@@ -6518,33 +6521,33 @@ namespace feng
         return 0;
     }
 
-    template < typename Type, Allocator Alloc >
-    std::optional<matrix<Type, Alloc>> lu_solver( const matrix< Type, Alloc >& A, matrix< Type, Alloc > const& b )
+    template< Matrix Mat >
+    std::optional<Mat> lu_solver( Mat const& A, Mat const& b )
     {
-        if ( matrix<Type, Alloc> x; lu_solver(A, x, b) == 1 )
+        if ( Mat x; lu_solver(A, x, b) == 1 )
             return {};
         else
             return x;
     }
 
-    template< typename Type, Allocator Alloc >
-    matrix<Type, Alloc> const conv( matrix<Type, Alloc> const& A, matrix<Type, Alloc> const& B ) noexcept
+    template< Matrix Mat >
+    Mat conv( Mat const& A, Mat const& B ) noexcept
     {
         if ( ( 0 == A.size() ) || ( 0 == B.size() ) )
-            return matrix<Type, Alloc>{0, 0};
+            return Mat{0, 0};
 
         if ( A.size() > B.size() )
             return conv( B, A );
 
-        matrix<Type, Alloc> padded_B{ B.row()+2*A.row()-2, B.col()+2*A.col()-2 };
+        Mat padded_B{ B.row()+2*A.row()-2, B.col()+2*A.col()-2 };
         for ( auto row : matrix_details::range(B.row()) )
             std::copy( B.row_begin(row), B.row_end(row), padded_B.row_begin(row+A.row()-1)+A.col()-1 );
 
-        matrix<Type, Alloc> ans{ A.row()+B.row()-1, A.col()+B.col()-1 };
+        Mat ans{ A.row()+B.row()-1, A.col()+B.col()-1 };
 
-        auto const& product = []( matrix<Type, Alloc> const& a, matrix_view<Type, Alloc> const& b, auto row, auto const col ) noexcept
+        auto const& product = []( Mat const& a, auto const& b, auto row, auto const col ) noexcept
         {
-            Type ans{0};
+            typename Mat::value_type ans{0};
             for ( auto r : matrix_details::range(row) )
                 for ( auto c : matrix_details::range(col) )
                     ans += a[r][c] * b[r][c];
@@ -6553,7 +6556,6 @@ namespace feng
 
         auto const& func = [&]( std::uint_least64_t row )
         {
-            //for ( auto col : matrix_details::range(ans.row() ))
             for ( auto col : matrix_details::range(ans.col() ))
             {
                 auto const& view = make_view( padded_B, {row, row+A.row()}, {col, col+A.col()} );
@@ -6566,12 +6568,9 @@ namespace feng
         return ans;
     }
 
-    template< typename Type, Allocator Alloc >
-    matrix<Type, Alloc> const conv( matrix<Type, Alloc> const& A, matrix<Type, Alloc> const& B, std::string const& mode ) noexcept
+    template< Matrix Mat >
+    Mat conv( Mat const& A, Mat const& B, std::string const& mode ) noexcept
     {
-        //if ( A.size() > B.size() )
-        //    return conv( B, A, mode );
-
         auto const& default_conv = conv( A, B );
 
         auto const [ra, ca] = A.shape();
@@ -6651,13 +6650,15 @@ namespace feng
         return {ans};
     }
 
-    template< typename T, Allocator A>
-    auto pooling( matrix<T,A> const& mat, std::uint_least64_t dim_r, std::uint_least64_t dim_c, std::string const& pooling_action = "mean" )
+    template< Matrix Mat >
+    auto pooling( Mat const& mat, std::uint_least64_t dim_r, std::uint_least64_t dim_c, std::string const& pooling_action = "mean" )
     {
-        if (dim_r == 0 || dim_c == 0) return matrix<T,A>{};
+        if (dim_r == 0 || dim_c == 0) return Mat{};
 
         if (dim_r==1 && dim_c==1) return mat;
 
+        typedef typename Mat::value_type T;
+
         std::map< std::string, std::pair<T, std::function< T(T, T) > > > function_list =
         {
             std::make_pair
@@ -6705,7 +6706,7 @@ namespace feng
         auto const& the_function =(*iterator).second.second;
         auto const [row, col] = mat.shape();
         auto const [new_row, new_col] = std::make_pair( row/dim_r, col/dim_c );
-        matrix< T, A > ans{ new_row, new_col, T{0} };
+        Mat ans{ new_row, new_col, T{0} };
 
         auto const& make_pooling = [&ans, &mat, &the_function, new_col=new_col, dim_r, dim_c, init_value]( std::uint_least64_t r )
         {
@@ -6722,14 +6723,14 @@ namespace feng
         return ans;
     }
 
-    template< typename T, Allocator A>
-    auto pooling( matrix<T,A> const& mat, std::uint_least64_t dim, std::string const& pooling_action = "mean"  )
+    template< Matrix Mat >
+    auto pooling( Mat const& mat, std::uint_least64_t dim, std::string const& pooling_action = "mean"  )
     {
         return pooling( mat, dim, dim, pooling_action );
     }
 
-    template< typename T, Allocator A>
-    void save_as_bmp( std::string const& file_name, matrix<T,A> const& red_channel, matrix<T,A> const& green_channel, matrix<T,A> const& blue_channel )
+    template< Matrix Mat >
+    void save_as_bmp( std::string const& file_name, Mat const& red_channel, Mat const& green_channel, Mat const& blue_channel )
     {
         better_assert( red_channel.row()==green_channel.row(), "Row not match for red and green matrix! The row for red is ", red_channel.row(), " but for green is ", green_channel.row() );
         better_assert( red_channel.row()==blue_channel.row(), "Row not match for red and blue matrix! The row for red is ", red_channel.row(), " but for blue is ", blue_channel.row() );
@@ -6767,8 +6768,8 @@ namespace feng
         stream.write( reinterpret_cast<char const*>((*encoding).data()), (*encoding).size() );
     }
 
-    template< typename T, Allocator A>
-    void save_as_bmp( std::string const& file_name, matrix<T,A> const& mat, std::string const& colormap = std::string{"parula"} )
+    template< Matrix Mat >
+    void save_as_bmp( std::string const& file_name, Mat const& mat, std::string const& colormap = std::string{"parula"} )
     {
         mat.save_as_bmp( file_name, colormap );
     }
@@ -6796,6 +6797,9 @@ namespace feng
         return std::make_pair( mat_y, mat_x );
     }
 
+    //template< Matrix M >
+
+
     static auto const& fma = matrix_details::map( []( auto const& val, auto const& wal, auto const& xal ){ return std::fma(val, wal, xal); } );
 
     static auto const& ldexp = matrix_details::map( []( auto const& val, auto const& wal ){ return std::ldexp(val, wal); } );

From ec9da410b870ce40ed4030bb2d4d4a1df21fc80e Mon Sep 17 00:00:00 2001
From: Feng <wang_feng@live.com>
Date: Sat, 16 Mar 2024 15:26:50 +0100
Subject: [PATCH 3/5] fix warning with sin

---
 matrix.hpp | 114 ++++++++++++++++++++++++++++++++++++++++++-----------
 1 file changed, 91 insertions(+), 23 deletions(-)

diff --git a/matrix.hpp b/matrix.hpp
index 7301170..a1d1eba 100644
--- a/matrix.hpp
+++ b/matrix.hpp
@@ -1098,39 +1098,34 @@ namespace feng
             {
                 typedef typename extract_type_forward< Index, Types... >::result_type result_type;
             };
+
             template < typename Function, typename InputIterator1, typename... InputIteratorn >
             Function _for_each_n( Function f, std::uint_least64_t n, InputIterator1 begin1, InputIteratorn... beginn )
             {
-                auto const& func = [&]( std::uint_least64_t idx )
-                {
-                    f( *(begin1+idx), *(beginn+idx)... );
-                };
+                auto const& func = [&]( std::uint_least64_t idx ) { f( *(begin1+idx), *(beginn+idx)... ); };
                 parallel( func, 0UL, n );
                 return f;
             }
+
             template < typename Function, typename InputIterator1, typename... InputIteratorn >
             Function _for_each( Function f, InputIterator1 begin1, InputIterator1 end1, InputIteratorn... beginn )
             {
                 return _for_each_n( f, std::distance( begin1, end1 ), begin1, beginn... );
-                /*
-                while ( begin1 != end1 )
-                    f( *begin1++, *beginn++... );
-                return f;
-                */
             }
 
-            struct dummy
-            {
-            };
+            struct dummy { };
+
             template < typename... Types_N >
             struct for_each_impl_with_dummy
             {
                 typedef typename extract_type_backward< 1, Types_N... >::result_type return_type;
+
                 template < typename Predict, typename... Types >
                 Predict impl( Predict p, dummy, Types... types ) const
                 {
                     return _for_each( p, types... );
                 }
+
                 template < typename S, typename... Types >
                 return_type impl( S s, Types... types ) const
                 {
@@ -3910,6 +3905,10 @@ namespace feng
 
     };//struct matrix
 
+    //
+    // - begin of Matrix and ComplexMatrix concepts
+    //
+
     template< typename T >
     struct is_matrix : std::false_type{};
 
@@ -3922,6 +3921,23 @@ namespace feng
     template< typename M >
     concept Matrix = is_matrix_v<M>;
 
+
+    template< typename T >
+    struct is_complex_matrix : std::false_type{};
+
+    template< typename T, Allocator A >
+    struct is_complex_matrix< matrix<std::complex<T>,A> > : std::true_type{};
+
+    template< typename M >
+    inline constexpr bool is_complex_matrix_v = is_complex_matrix<M>::value;
+
+    template< typename M >
+    concept ComplexMatrix = is_complex_matrix_v<M>;
+
+    //
+    // - end of Matrix and ComplexMatrix concepts
+    //
+
     namespace matrix_details
     {
 
@@ -4700,6 +4716,7 @@ namespace feng
 
         return true;
     }
+
     template < typename T, Allocator A>
     bool is_symmetric( const matrix< T, A >& m )
     {
@@ -4708,6 +4725,7 @@ namespace feng
             return v1 == v2;
         } );
     }
+
     inline matrix< std::uint_least64_t > const magic( const std::uint_least64_t n ) noexcept
     {
         matrix< std::uint_least64_t > ans{ n, n };
@@ -4790,20 +4808,21 @@ namespace feng
         std::swap_ranges( ans.upper_anti_diag_begin( n >> 1 ), ans.upper_anti_diag_end( n >> 1 ), ans.lower_anti_diag_rbegin( n >> 1 ) );
         return ans;
     }
-    template < typename T, Allocator A>
-    T const
-    max( const matrix< T, A >& m )
+
+    template< Matrix Mat >
+    auto const max( Mat const& m )
     {
         better_assert( m.size() );
         return *std::max_element( m.begin(), m.end() );
     }
-    template < typename T, Allocator A>
-    T const
-    min( const matrix< T, A >& m )
+
+    template< Matrix Mat >
+    auto const min( Mat const& m )
     {
         better_assert( m.size() );
         return *std::min_element( m.begin(), m.end() );
     }
+
     template < typename T >
     matrix<T> arange( std::uint_least64_t start, const std::uint_least64_t stop, const std::uint_least64_t step = 1ULL )
     {
@@ -4815,6 +4834,7 @@ namespace feng
         }
         return ans;
     }
+
     template < typename T >
     matrix<T> arange( const std::uint_least64_t length )
     {
@@ -4822,6 +4842,7 @@ namespace feng
         std::iota( ans.begin(), ans.end(), T{0} );
         return ans;
     }
+
     template < typename T >
     matrix<T> linspace( T start, T stop, const std::uint_least64_t num = 50ULL, bool end_point=true )
     {
@@ -4839,33 +4860,72 @@ namespace feng
         }
         return ans;
     }
-    template < typename T, Allocator A>
-    matrix< T, A > const ones_like( matrix<T, A> const& mat ) noexcept
+
+    template < Matrix Mat >
+    Mat const ones_like( Mat const& mat ) noexcept
     {
-        return matrix<T, A>{ mat.get_allocator(), mat.row(), mat.col(), T{1} };
+        return Mat{ mat.get_allocator(), mat.row(), mat.col(), typename Mat::value_type{1} };
     }
+
     template < typename T >
     matrix<T, std::allocator<T>> const ones( std::integral auto r, std::integral auto c ) noexcept
     {
         matrix< T > ans{ static_cast<unsigned long>(r), static_cast<unsigned long>(c), T{ 1 } };
         return ans;
     }
+
     template < typename T >
     matrix<T, std::allocator<T>> const ones( std::integral auto n ) noexcept
     {
         return ones< T >( n, n );
     }
+
     template < typename T, Allocator A>
     matrix< T, A > const ones( A const& alloc, std::integral auto r, std::integral auto c ) noexcept
     {
         return { alloc, static_cast<unsigned long>(r) , static_cast<unsigned long>(c), T{1} };
     }
+
     template < typename T, Allocator A>
     matrix< T, A > const ones( A const& alloc, std::integral auto n ) noexcept
     {
         return { alloc, static_cast<unsigned long>(n) , static_cast<unsigned long>(n), T{1} };
     }
 
+
+    template < Matrix Mat >
+    Mat const zeros_like( Mat const& mat ) noexcept
+    {
+        return Mat{ mat.get_allocator(), mat.row(), mat.col(), typename Mat::value_type{} };
+    }
+
+    template < typename T >
+    matrix<T, std::allocator<T>> const zeros( std::integral auto r, std::integral auto c ) noexcept
+    {
+        matrix< T > ans{ static_cast<unsigned long>(r), static_cast<unsigned long>(c), T{} };
+        return ans;
+    }
+
+    template < typename T >
+    matrix<T, std::allocator<T>> const zeros( std::integral auto n ) noexcept
+    {
+        return zeros< T >( n, n );
+    }
+
+    template < typename T, Allocator A>
+    matrix< T, A > const zeros( A const& alloc, std::integral auto r, std::integral auto c ) noexcept
+    {
+        return { alloc, static_cast<unsigned long>(r) , static_cast<unsigned long>(c), T{} };
+    }
+
+    template < typename T, Allocator A>
+    matrix< T, A > const zeros( A const& alloc, std::integral auto n ) noexcept
+    {
+        return { alloc, static_cast<unsigned long>(n) , static_cast<unsigned long>(n), T{} };
+    }
+
+
+
     template < typename T >
     auto const empty( const std::integral auto r, const std::integral auto c ) noexcept
     {
@@ -6797,7 +6857,15 @@ namespace feng
         return std::make_pair( mat_y, mat_x );
     }
 
-    //template< Matrix M >
+    template< Matrix Mat >
+    auto sin( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::sin(x); } );
+        return ans;
+    }
+
+
 
 
     static auto const& fma = matrix_details::map( []( auto const& val, auto const& wal, auto const& xal ){ return std::fma(val, wal, xal); } );
@@ -6830,7 +6898,7 @@ namespace feng
 
     static auto const& cos = matrix_details::map( []( auto const& val ){ return std::cos(val); } );
 
-    static auto const& sin = matrix_details::map( []( auto const& val ){ return std::sin(val); } );
+    //static auto const& sin = matrix_details::map( []( auto const& val ){ return std::sin(val); } );
 
     static auto const& tan = matrix_details::map( []( auto const& val ){ return std::tan(val); } );
 

From 638ccf34c24edd19a143a66ce51539f664c4bc9c Mon Sep 17 00:00:00 2001
From: Feng <wang_feng@live.com>
Date: Sat, 16 Mar 2024 16:00:56 +0100
Subject: [PATCH 4/5] regenerate unary functions

---
 matrix.hpp           | 399 ++++++++++++++++++++++++++++++++++---------
 tests/cases/fabs.hpp |   2 +-
 2 files changed, 324 insertions(+), 77 deletions(-)

diff --git a/matrix.hpp b/matrix.hpp
index a1d1eba..700e988 100644
--- a/matrix.hpp
+++ b/matrix.hpp
@@ -4289,27 +4289,6 @@ namespace feng
         return ans;
     }
 
-    //
-    // zeros: creating matrix of all zeros.
-    //
-    template< typename T = double >
-    auto const zeros( std::integral auto r, std::integral auto c )
-    {
-        return matrix<T>{ static_cast<unsigned long>(r), static_cast<unsigned long>(c) };
-    }
-
-    template< typename T, Allocator A >
-    auto const zeros( A const& alloc, std::integral auto r, std::integral auto c )
-    {
-        return matrix<T, A>{ alloc, static_cast<unsigned long>(r), static_cast<unsigned long>(c) };
-    }
-
-    template< typename T = double >
-    auto const zeros( std::integral auto n )
-    {
-        return zeros( n, n );
-    }
-
     template < typename T1, Allocator A1, typename T2, Allocator A2 >
     matrix< T1, A1 > const blkdiag( const matrix< T1, A1 >& m1, const matrix< T2, A2 >& m2 )
     {
@@ -6857,122 +6836,390 @@ namespace feng
         return std::make_pair( mat_y, mat_x );
     }
 
+    //
+    // - begin of unary functions
+    //
+
     template< Matrix Mat >
-    auto sin( Mat const& m )
+    auto abs( Mat const& m )
     {
         auto ans = zeros_like( m );
-        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::sin(x); } );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::abs(x); } );
         return ans;
     }
 
+    template< Matrix Mat >
+    auto exp( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::exp(x); } );
+        return ans;
+    }
 
+    template< Matrix Mat >
+    auto exp2( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::exp2(x); } );
+        return ans;
+    }
 
+    template< Matrix Mat >
+    auto expm1( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::expm1(x); } );
+        return ans;
+    }
 
-    static auto const& fma = matrix_details::map( []( auto const& val, auto const& wal, auto const& xal ){ return std::fma(val, wal, xal); } );
+    template< Matrix Mat >
+    auto log( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::log(x); } );
+        return ans;
+    }
 
-    static auto const& ldexp = matrix_details::map( []( auto const& val, auto const& wal ){ return std::ldexp(val, wal); } );
+    template< Matrix Mat >
+    auto log10( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::log10(x); } );
+        return ans;
+    }
 
-    static auto const& scalbn = matrix_details::map( []( auto const& val, auto const& wal ){ return std::scalbn(val, wal); } );
+    template< Matrix Mat >
+    auto log1p( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::log1p(x); } );
+        return ans;
+    }
 
-    static auto const& scalbln = matrix_details::map( []( auto const& val, auto const& wal ){ return std::scalbln(val, wal); } );
+    template< Matrix Mat >
+    auto log2( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::log2(x); } );
+        return ans;
+    }
 
-    static auto const& pow = matrix_details::map( []( auto const& val, auto const& wal ){ return std::pow(val, wal); } );
+    template< Matrix Mat >
+    auto sqrt( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::sqrt(x); } );
+        return ans;
+    }
 
-    static auto const& hypot = matrix_details::map( []( auto const& val, auto const& wal ){ return std::hypot(val, wal); } );
+    template< Matrix Mat >
+    auto cbrt( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::cbrt(x); } );
+        return ans;
+    }
 
-    static auto const& fmod = matrix_details::map( []( auto const& val, auto const& wal ){ return std::fmod(val, wal); } );
+    template< Matrix Mat >
+    auto sin( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::sin(x); } );
+        return ans;
+    }
 
-    static auto const& remainder = matrix_details::map( []( auto const& val, auto const& wal ){ return std::remainder(val, wal); } );
+    template< Matrix Mat >
+    auto cos( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::cos(x); } );
+        return ans;
+    }
 
-    static auto const& copysign = matrix_details::map( []( auto const& val, auto const& wal ){ return std::copysign(val, wal); } );
+    template< Matrix Mat >
+    auto tan( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::tan(x); } );
+        return ans;
+    }
 
-    static auto const& nextafter = matrix_details::map( []( auto const& val, auto const& wal ){ return std::nextafter(val, wal); } );
+    template< Matrix Mat >
+    auto asin( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::asin(x); } );
+        return ans;
+    }
 
-    static auto const& fdim = matrix_details::map( []( auto const& val, auto const& wal ){ return std::fdim(val, wal); } );
+    template< Matrix Mat >
+    auto acos( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::acos(x); } );
+        return ans;
+    }
 
-    static auto const& fmax = matrix_details::map( []( auto const& val, auto const& wal ){ return std::fmax(val, wal); } );
+    template< Matrix Mat >
+    auto atan( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::atan(x); } );
+        return ans;
+    }
 
-    static auto const& fmin = matrix_details::map( []( auto const& val, auto const& wal ){ return std::fmin(val, wal); } );
+    template< Matrix Mat >
+    auto sinh( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::sinh(x); } );
+        return ans;
+    }
 
-    static auto const& atan2 = matrix_details::map( []( auto const& val, auto const& wal ){ return std::atan2(val, wal); } );
+    template< Matrix Mat >
+    auto cosh( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::cosh(x); } );
+        return ans;
+    }
+
+    template< Matrix Mat >
+    auto tanh( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::tanh(x); } );
+        return ans;
+    }
+
+    template< Matrix Mat >
+    auto asinh( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::asinh(x); } );
+        return ans;
+    }
 
-    static auto const& cos = matrix_details::map( []( auto const& val ){ return std::cos(val); } );
+    template< Matrix Mat >
+    auto acosh( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::acosh(x); } );
+        return ans;
+    }
 
-    //static auto const& sin = matrix_details::map( []( auto const& val ){ return std::sin(val); } );
+    template< Matrix Mat >
+    auto atanh( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::atanh(x); } );
+        return ans;
+    }
 
-    static auto const& tan = matrix_details::map( []( auto const& val ){ return std::tan(val); } );
+    template< Matrix Mat >
+    auto erf( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::erf(x); } );
+        return ans;
+    }
 
-    static auto const& acos = matrix_details::map( []( auto const& val ){ return std::acos(val); } );
+    template< Matrix Mat >
+    auto erfc( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::erfc(x); } );
+        return ans;
+    }
 
-    static auto const& asin = matrix_details::map( []( auto const& val ){ return std::asin(val); } );
+    template< Matrix Mat >
+    auto tgamma( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::tgamma(x); } );
+        return ans;
+    }
 
-    static auto const& atan = matrix_details::map( []( auto const& val ){ return std::atan(val); } );
+    template< Matrix Mat >
+    auto lgamma( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::lgamma(x); } );
+        return ans;
+    }
 
-    static auto const& cosh = matrix_details::map( []( auto const& val ){ return std::cosh(val); } );
+    template< Matrix Mat >
+    auto trunc( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::trunc(x); } );
+        return ans;
+    }
+
+    template< Matrix Mat >
+    auto round( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::round(x); } );
+        return ans;
+    }
+
+    template< Matrix Mat >
+    auto ceil( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::ceil(x); } );
+        return ans;
+    }
+
+    template< Matrix Mat >
+    auto floor( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::floor(x); } );
+        return ans;
+    }
+
+    template< Matrix Mat >
+    auto rint( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::rint(x); } );
+        return ans;
+    }
 
-    static auto const& sinh = matrix_details::map( []( auto const& val ){ return std::sinh(val); } );
+    template< Matrix Mat >
+    auto logb( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::logb(x); } );
+        return ans;
+    }
 
-    static auto const& tanh = matrix_details::map( []( auto const& val ){ return std::tanh(val); } );
+    template< Matrix Mat >
+    auto comp_ellint_1( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::comp_ellint_1(x); } );
+        return ans;
+    }
 
-    static auto const& acosh = matrix_details::map( []( auto const& val ){ return std::acosh(val); } );
+    template< Matrix Mat >
+    auto comp_ellint_2( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::comp_ellint_2(x); } );
+        return ans;
+    }
 
-    static auto const& asinh = matrix_details::map( []( auto const& val ){ return std::asinh(val); } );
+    template< Matrix Mat >
+    auto comp_ellint_3( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::comp_ellint_3(x); } );
+        return ans;
+    }
 
-    static auto const& atanh = matrix_details::map( []( auto const& val ){ return std::atanh(val); } );
+    template< Matrix Mat >
+    auto expint( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::expint(x); } );
+        return ans;
+    }
 
-    static auto const& exp = matrix_details::map( []( auto const& val ){ return std::exp(val); } );
+    template< Matrix Mat >
+    auto riemann_zeta( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::riemann_zeta(x); } );
+        return ans;
+    }
 
-    static auto const& log = matrix_details::map( []( auto const& val ){ return std::log(val); } );
+    template< Matrix Mat >
+    auto nearbyint( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::nearbyint(x); } );
+        return ans;
+    }
 
-    static auto const& log10 = matrix_details::map( []( auto const& val ){ return std::log10(val); } );
+    template< Matrix Mat >
+    auto ilogb( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::ilogb(x); } );
+        return ans.template astype<int>(); // <- in cmath, ilogb returns an int
+    }
 
-    static auto const& exp2 = matrix_details::map( []( auto const& val ){ return std::exp2(val); } );
+    template< Matrix Mat >
+    auto lrint( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::lrint(x); } );
+        return ans.template astype<long>(); // <- in cmath, lrint returns a long
+    }
 
-    static auto const& expm1 = matrix_details::map( []( auto const& val ){ return std::expm1(val); } );
+    template< Matrix Mat >
+    auto llrint( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::llrint(x); } );
+        return ans.template astype<long long>(); // <- in cmath, llrint returns a long long
+    }
 
-    static auto const& log1p = matrix_details::map( []( auto const& val ){ return std::log1p(val); } );
+    template< Matrix Mat >
+    auto lround( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::lround(x); } );
+        return ans.template astype<long>(); // <- in cmath, lround returns a long
+    }
 
-    static auto const& log2 = matrix_details::map( []( auto const& val ){ return std::log2(val); } );
+    template< Matrix Mat >
+    auto llround( Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::llround(x); } );
+        return ans.template astype<long long>(); // <- in cmath, llround returns a long long
+    }
 
-    static auto const& logb = matrix_details::map( []( auto const& val ){ return std::logb(val); } );
 
-    static auto const& sqrt = matrix_details::map( []( auto const& val ){ return std::sqrt(val); } );
+    //
+    // - end of unary functions
+    //
 
-    static auto const& cbrt = matrix_details::map( []( auto const& val ){ return std::cbrt(val); } );
 
-    static auto const& erf = matrix_details::map( []( auto const& val ){ return std::erf(val); } );
 
-    static auto const& erfc = matrix_details::map( []( auto const& val ){ return std::erfc(val); } );
 
-    static auto const& tgamma = matrix_details::map( []( auto const& val ){ return std::tgamma(val); } );
+    static auto const& fma = matrix_details::map( []( auto const& val, auto const& wal, auto const& xal ){ return std::fma(val, wal, xal); } );
 
-    static auto const& lgamma = matrix_details::map( []( auto const& val ){ return std::lgamma(val); } );
+    static auto const& ldexp = matrix_details::map( []( auto const& val, auto const& wal ){ return std::ldexp(val, wal); } );
 
-    static auto const& ceil = matrix_details::map( []( auto const& val ){ return std::ceil(val); } );
+    static auto const& scalbn = matrix_details::map( []( auto const& val, auto const& wal ){ return std::scalbn(val, wal); } );
 
-    static auto const& floor = matrix_details::map( []( auto const& val ){ return std::floor(val); } );
+    static auto const& scalbln = matrix_details::map( []( auto const& val, auto const& wal ){ return std::scalbln(val, wal); } );
 
-    static auto const& trunc = matrix_details::map( []( auto const& val ){ return std::trunc(val); } );
+    static auto const& pow = matrix_details::map( []( auto const& val, auto const& wal ){ return std::pow(val, wal); } );
 
-    static auto const& round = matrix_details::map( []( auto const& val ){ return std::round(val); } );
+    static auto const& hypot = matrix_details::map( []( auto const& val, auto const& wal ){ return std::hypot(val, wal); } );
 
-    static auto const& rint = matrix_details::map( []( auto const& val ){ return std::rint(val); } );
+    static auto const& fmod = matrix_details::map( []( auto const& val, auto const& wal ){ return std::fmod(val, wal); } );
 
-    static auto const& nearbyint = matrix_details::map( []( auto const& val ){ return std::nearbyint(val); } );
+    static auto const& remainder = matrix_details::map( []( auto const& val, auto const& wal ){ return std::remainder(val, wal); } );
 
-    static auto const& fabs = matrix_details::map( []( auto const& val ){ return std::fabs(val); } );
+    static auto const& copysign = matrix_details::map( []( auto const& val, auto const& wal ){ return std::copysign(val, wal); } );
 
-    static auto const& abs = matrix_details::map( []( auto const& val ){ return std::abs(val); } );
+    static auto const& nextafter = matrix_details::map( []( auto const& val, auto const& wal ){ return std::nextafter(val, wal); } );
 
-    static auto const& ilogb = matrix_details::map( []( auto const& val ){ return std::ilogb(val); } );
+    static auto const& fdim = matrix_details::map( []( auto const& val, auto const& wal ){ return std::fdim(val, wal); } );
 
-    static auto const& lround = matrix_details::map( []( auto const& val ){ return std::lround(val); } );
+    static auto const& fmax = matrix_details::map( []( auto const& val, auto const& wal ){ return std::fmax(val, wal); } );
 
-    static auto const& llround = matrix_details::map( []( auto const& val ){ return std::llround(val); } );
+    static auto const& fmin = matrix_details::map( []( auto const& val, auto const& wal ){ return std::fmin(val, wal); } );
 
-    static auto const& lrint = matrix_details::map( []( auto const& val ){ return std::lrint(val); } );
+    static auto const& atan2 = matrix_details::map( []( auto const& val, auto const& wal ){ return std::atan2(val, wal); } );
 
-    static auto const& llrint = matrix_details::map( []( auto const& val ){ return std::llrint(val); } );
 
     // special
     static auto const& remquo = matrix_details::map( []( auto const& num, auto const& den ){ int quo; auto const ans = std::remquo(num, den, &quo); return std::make_tuple(ans, quo); } );
diff --git a/tests/cases/fabs.hpp b/tests/cases/fabs.hpp
index 5da5a04..9829409 100644
--- a/tests/cases/fabs.hpp
+++ b/tests/cases/fabs.hpp
@@ -6,7 +6,7 @@ TEST_CASE( "Matrix fabs", "[fabs]" )
         for ( unsigned long tz = 1; tz != N; ++tz )
         {
             feng::matrix<double> const& mat = feng::rand<double>( sz, tz );
-            auto const& mat_fabs = feng::fabs( mat );
+            auto const& mat_fabs = feng::abs( mat );
             for ( unsigned long r = 0; r != sz; ++r )
                 for ( unsigned long c = 0; c != tz; ++c )
                 {

From fbda8bbaf5abe00ff2e61e11f8263c246e94a604 Mon Sep 17 00:00:00 2001
From: Feng <wang_feng@live.com>
Date: Sat, 16 Mar 2024 21:58:05 +0100
Subject: [PATCH 5/5] recode functions to remove dangling-reference warning

---
 matrix.hpp    | 502 +++++++++++++++++++++++++++++++++++++++++++++-----
 tests/test.cc |   2 +-
 2 files changed, 452 insertions(+), 52 deletions(-)

diff --git a/matrix.hpp b/matrix.hpp
index 700e988..603a33c 100644
--- a/matrix.hpp
+++ b/matrix.hpp
@@ -7190,95 +7190,495 @@ namespace feng
     //
 
 
+    //
+    // - begin of trinary functions
+    //
+
+    template< Matrix Mat >
+    auto fma( Mat const& m1, Mat const& m2, Mat const& m3 )
+    {
+        auto ans = zeros_like( m1 );
+        matrix_details::for_each( m1.begin(), m1.end(), m2.begin(), m3.begin(),  ans.begin(), []( auto const x1, auto const x2, auto const x3,  auto& v ){ v = std::fma(x1, x2, x3); } );
+        return ans;
+    }
+
+    //
+    // - end of trinary functions
+    //
+
+    //
+    // - begin of binary functions
+    //
+
+    template< Matrix Mat, Matrix Nat >
+    auto ldexp( Mat const& m, Nat const& n )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), n.begin(), ans.begin(), []( auto const xm, auto const xn, auto& v ){ v = std::ldexp(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto ldexp( Mat const& m, std::floating_point auto xn )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::ldexp(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto ldexp( std::floating_point auto xn, Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::ldexp(xn, xm); } );
+        return ans;
+    }
+
+    template< Matrix Mat, Matrix Nat >
+    auto scalbn( Mat const& m, Nat const& n )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), n.begin(), ans.begin(), []( auto const xm, auto const xn, auto& v ){ v = std::scalbn(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto scalbn( Mat const& m, std::floating_point auto xn )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::scalbn(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto scalbn( std::floating_point auto xn, Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::scalbn(xn, xm); } );
+        return ans;
+    }
+
+    template< Matrix Mat, Matrix Nat >
+    auto scalbln( Mat const& m, Nat const& n )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), n.begin(), ans.begin(), []( auto const xm, auto const xn, auto& v ){ v = std::scalbln(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto scalbln( Mat const& m, std::floating_point auto xn )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::scalbln(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto scalbln( std::floating_point auto xn, Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::scalbln(xn, xm); } );
+        return ans;
+    }
+
+    template< Matrix Mat, Matrix Nat >
+    auto pow( Mat const& m, Nat const& n )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), n.begin(), ans.begin(), []( auto const xm, auto const xn, auto& v ){ v = std::pow(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto pow( Mat const& m, std::integral auto xn ) // <- pow accept integral as the second argument
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::pow(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto pow( Mat const& m, std::floating_point auto xn )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::pow(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto pow( std::floating_point auto xn, Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::pow(xn, xm); } );
+        return ans;
+    }
+
+    template< Matrix Mat, Matrix Nat >
+    auto hypot( Mat const& m, Nat const& n )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), n.begin(), ans.begin(), []( auto const xm, auto const xn, auto& v ){ v = std::hypot(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto hypot( Mat const& m, std::floating_point auto xn )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::hypot(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto hypot( std::floating_point auto xn, Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::hypot(xn, xm); } );
+        return ans;
+    }
+
+    template< Matrix Mat, Matrix Nat >
+    auto fmod( Mat const& m, Nat const& n )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), n.begin(), ans.begin(), []( auto const xm, auto const xn, auto& v ){ v = std::fmod(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto fmod( Mat const& m, std::floating_point auto xn )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::fmod(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto fmod( std::floating_point auto xn, Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::fmod(xn, xm); } );
+        return ans;
+    }
+
+    template< Matrix Mat, Matrix Nat >
+    auto remainder( Mat const& m, Nat const& n )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), n.begin(), ans.begin(), []( auto const xm, auto const xn, auto& v ){ v = std::remainder(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto remainder( Mat const& m, std::floating_point auto xn )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::remainder(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto remainder( std::floating_point auto xn, Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::remainder(xn, xm); } );
+        return ans;
+    }
+
+    template< Matrix Mat, Matrix Nat >
+    auto copysign( Mat const& m, Nat const& n )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), n.begin(), ans.begin(), []( auto const xm, auto const xn, auto& v ){ v = std::copysign(xm, xn); } );
+        return ans;
+    }
+
+    template< Matrix Mat>
+    auto copysign( Mat const& m, std::floating_point auto xn )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::copysign(xm, xn); } );
+        return ans;
+    }
 
+    template< Matrix Mat>
+    auto copysign( std::floating_point auto xn, Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::copysign(xn, xm); } );
+        return ans;
+    }
 
-    static auto const& fma = matrix_details::map( []( auto const& val, auto const& wal, auto const& xal ){ return std::fma(val, wal, xal); } );
+    template< Matrix Mat, Matrix Nat >
+    auto nextafter( Mat const& m, Nat const& n )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), n.begin(), ans.begin(), []( auto const xm, auto const xn, auto& v ){ v = std::nextafter(xm, xn); } );
+        return ans;
+    }
 
-    static auto const& ldexp = matrix_details::map( []( auto const& val, auto const& wal ){ return std::ldexp(val, wal); } );
+    template< Matrix Mat>
+    auto nextafter( Mat const& m, std::floating_point auto xn )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::nextafter(xm, xn); } );
+        return ans;
+    }
 
-    static auto const& scalbn = matrix_details::map( []( auto const& val, auto const& wal ){ return std::scalbn(val, wal); } );
+    template< Matrix Mat>
+    auto nextafter( std::floating_point auto xn, Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::nextafter(xn, xm); } );
+        return ans;
+    }
 
-    static auto const& scalbln = matrix_details::map( []( auto const& val, auto const& wal ){ return std::scalbln(val, wal); } );
+    template< Matrix Mat, Matrix Nat >
+    auto fdim( Mat const& m, Nat const& n )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), n.begin(), ans.begin(), []( auto const xm, auto const xn, auto& v ){ v = std::fdim(xm, xn); } );
+        return ans;
+    }
 
-    static auto const& pow = matrix_details::map( []( auto const& val, auto const& wal ){ return std::pow(val, wal); } );
+    template< Matrix Mat>
+    auto fdim( Mat const& m, std::floating_point auto xn )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::fdim(xm, xn); } );
+        return ans;
+    }
 
-    static auto const& hypot = matrix_details::map( []( auto const& val, auto const& wal ){ return std::hypot(val, wal); } );
+    template< Matrix Mat>
+    auto fdim( std::floating_point auto xn, Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::fdim(xn, xm); } );
+        return ans;
+    }
 
-    static auto const& fmod = matrix_details::map( []( auto const& val, auto const& wal ){ return std::fmod(val, wal); } );
+    template< Matrix Mat, Matrix Nat >
+    auto fmax( Mat const& m, Nat const& n )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), n.begin(), ans.begin(), []( auto const xm, auto const xn, auto& v ){ v = std::fmax(xm, xn); } );
+        return ans;
+    }
 
-    static auto const& remainder = matrix_details::map( []( auto const& val, auto const& wal ){ return std::remainder(val, wal); } );
+    template< Matrix Mat>
+    auto fmax( Mat const& m, std::floating_point auto xn )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::fmax(xm, xn); } );
+        return ans;
+    }
 
-    static auto const& copysign = matrix_details::map( []( auto const& val, auto const& wal ){ return std::copysign(val, wal); } );
+    template< Matrix Mat>
+    auto fmax( std::floating_point auto xn, Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::fmax(xn, xm); } );
+        return ans;
+    }
 
-    static auto const& nextafter = matrix_details::map( []( auto const& val, auto const& wal ){ return std::nextafter(val, wal); } );
+    template< Matrix Mat, Matrix Nat >
+    auto fmin( Mat const& m, Nat const& n )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), n.begin(), ans.begin(), []( auto const xm, auto const xn, auto& v ){ v = std::fmin(xm, xn); } );
+        return ans;
+    }
 
-    static auto const& fdim = matrix_details::map( []( auto const& val, auto const& wal ){ return std::fdim(val, wal); } );
+    template< Matrix Mat>
+    auto fmin( Mat const& m, std::floating_point auto xn )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::fmin(xm, xn); } );
+        return ans;
+    }
 
-    static auto const& fmax = matrix_details::map( []( auto const& val, auto const& wal ){ return std::fmax(val, wal); } );
+    template< Matrix Mat>
+    auto fmin( std::floating_point auto xn, Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::fmin(xn, xm); } );
+        return ans;
+    }
 
-    static auto const& fmin = matrix_details::map( []( auto const& val, auto const& wal ){ return std::fmin(val, wal); } );
+    template< Matrix Mat, Matrix Nat >
+    auto atan2( Mat const& m, Nat const& n )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), n.begin(), ans.begin(), []( auto const xm, auto const xn, auto& v ){ v = std::atan2(xm, xn); } );
+        return ans;
+    }
 
-    static auto const& atan2 = matrix_details::map( []( auto const& val, auto const& wal ){ return std::atan2(val, wal); } );
+    template< Matrix Mat>
+    auto atan2( Mat const& m, std::floating_point auto xn )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::atan2(xm, xn); } );
+        return ans;
+    }
 
+    template< Matrix Mat>
+    auto atan2( std::floating_point auto xn, Mat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), [xn]( auto const xm, auto& v ){ v = std::atan2(xn, xm); } );
+        return ans;
+    }
 
-    // special
-    static auto const& remquo = matrix_details::map( []( auto const& num, auto const& den ){ int quo; auto const ans = std::remquo(num, den, &quo); return std::make_tuple(ans, quo); } );
 
     //
-    // Complex Functions
+    // - end of binary functions
     //
-    static auto const& real = matrix_details::map( []( auto const& val ){ return std::real(val); } );
 
-    static auto const& imag = matrix_details::map( []( auto const& val ){ return std::imag(val); } );
 
-    static auto const& arg = matrix_details::map( []( auto const& val ){ return std::arg(val); } );
 
-    static auto const& norm = matrix_details::map( []( auto const& val ){ return std::norm(val); } );
 
-    static auto const& conj = matrix_details::map( []( auto const& val ){ return std::conj(val); } );
+    //
+    // - begin of Complex Functions
+    //
 
-    static auto const& polar = matrix_details::map( []( auto const& val ){ return std::polar(val); } );
+    template< ComplexMatrix CMat >
+    auto real( CMat const& cm )
+    {
+        typedef typename CMat::value_type complex_type;
+        typedef typename complex_type::value_type value_type;
+        typedef typename CMat::allocator_type Alloc;
+        matrix<value_type, typename std::allocator_traits<Alloc>:: template rebind_alloc<value_type> > ans{ cm.row(), cm.col() };
+        matrix_details::for_each( cm.begin(), cm.end(), ans.begin(), []( auto const& val, auto& v ){ v = std::real(val); } );
+        return ans;
+    }
 
-    static auto const& proj = matrix_details::map( []( auto const& val ){ return std::proj(val); } );
+    template< ComplexMatrix CMat >
+    auto imag( CMat const& cm )
+    {
+        typedef typename CMat::value_type complex_type;
+        typedef typename complex_type::value_type value_type;
+        typedef typename CMat::allocator_type Alloc;
+        matrix<value_type, typename std::allocator_traits<Alloc>:: template rebind_alloc<value_type> > ans{ cm.row(), cm.col() };
+        matrix_details::for_each( cm.begin(), cm.end(), ans.begin(), []( auto const& val, auto& v ){ v = std::imag(val); } );
+        return ans;
+    }
 
-    // already previously defined
-    // abs
-    // sin
-    // cos
-    // ...
-    //
+    template< ComplexMatrix CMat >
+    auto abs( CMat const& cm )
+    {
+        typedef typename CMat::value_type complex_type;
+        typedef typename complex_type::value_type value_type;
+        typedef typename CMat::allocator_type Alloc;
+        matrix<value_type, typename std::allocator_traits<Alloc>:: template rebind_alloc<value_type> > ans{ cm.row(), cm.col() };
+        matrix_details::for_each( cm.begin(), cm.end(), ans.begin(), []( auto const& val, auto& v ){ v = std::abs(val); } );
+        return ans;
+    }
+
+    template< ComplexMatrix CMat >
+    auto arg( CMat const& cm )
+    {
+        typedef typename CMat::value_type complex_type;
+        typedef typename complex_type::value_type value_type;
+        typedef typename CMat::allocator_type Alloc;
+        matrix<value_type, typename std::allocator_traits<Alloc>:: template rebind_alloc<value_type> > ans{ cm.row(), cm.col() };
+        matrix_details::for_each( cm.begin(), cm.end(), ans.begin(), []( auto const& val, auto& v ){ v = std::arg(val); } );
+        return ans;
+    }
+
+    template< ComplexMatrix CMat >
+    auto norm( CMat const& cm )
+    {
+        typedef typename CMat::value_type complex_type;
+        typedef typename complex_type::value_type value_type;
+        typedef typename CMat::allocator_type Alloc;
+        matrix<value_type, typename std::allocator_traits<Alloc>:: template rebind_alloc<value_type> > ans{ cm.row(), cm.col() };
+        matrix_details::for_each( cm.begin(), cm.end(), ans.begin(), []( auto const& val, auto& v ){ v = std::norm(val); } );
+        return ans;
+    }
+
+    template< ComplexMatrix CMat >
+    auto conj( CMat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::conj(x); } );
+        return ans;
+    }
+
+    template< ComplexMatrix CMat >
+    auto proj( CMat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::proj(x); } );
+        return ans;
+    }
+
+    template< ComplexMatrix CMat >
+    auto polar( CMat const& m )
+    {
+        auto ans = zeros_like( m );
+        matrix_details::for_each( m.begin(), m.end(), ans.begin(), []( auto const x, auto& v ){ v = std::polar(x); } );
+        return ans;
+    }
 
-    // reduce functions
     //
-    static auto const& sum = matrix_details::reduce( []( auto x, auto y ) noexcept { return x+y; }, 0.0 );
+    // - end of Complex Functions
+    //
 
-    static auto const& mean = []( auto const& mat ) noexcept { return sum(mat) / mat.size(); };
 
-    static auto const& variance = []( auto const& mat ) noexcept { return sum( pow( mat - mean( mat ), 2.0 ) ) / mat.size(); };
+    ///
+    /// - begin of reduce functions
+    ///
 
-    static auto const& standard_deviation = []( auto const& mat ) noexcept
+    template< Matrix Mat >
+    auto sum( Mat const& m )
     {
-        if ( mat.size() <= 1 )
-            return 0.0;
-        return std::sqrt( sum( pow( mat-mean( mat ), 2.0 ) ) / ( mat.size() - 1 ) );
-    };
+        return matrix_details::reduce( m.begin(), m.end(), typename Mat::value_type{}, []( auto const x, auto const y ){ return x+y; } );
+    }
+
+    template< Matrix Mat >
+    auto mean( Mat const& m )
+    {
+        return sum( m ) / m.size();
+    }
+
+    template< Matrix Mat >
+    auto variance( Mat const& m )
+    {
+        return mean( pow( m-mean(m), 2.0 ) );
+    }
+
+    template< Matrix Mat >
+    auto standard_deviation( Mat const& m )
+    {
+        if ( m.size() <= 1 )
+            return typename Mat::value_type{};
+        return std::sqrt( sum( pow( m-mean( m ), 2.0 ) ) / ( m.size() - 1 ) );
+    }
+
+    ///
+    /// - end of reduce functions
+    ///
 
-    static auto const& clip = []( auto lower, auto upper ) noexcept
+
+    ///
+    /// - begin of misc functions
+    ///
+
+    template< typename T >
+    auto clip( T const lower, T const upper )
     {
-        better_assert( lower < upper, "Error in clip(lower, upper): lower ", lower, " is supposed to be smaller than upper ", upper );
-        return [=]( auto const& mat ) noexcept
+        return [lower, upper]<Matrix Mat>( Mat const& m )
         {
-            auto ans{ mat };
-            auto&& impl = [=]( auto & v ) noexcept
-            {
-                v = v < lower ? lower :
-                    v > upper ? upper :
-                                v;
-            };
-            ans.map( impl );
+            auto ans{ m };
+            matrix_details::for_each( ans.begin(), ans.end(), [lower, upper]( auto & v ){ v = v < lower ? lower : v > upper ? upper : v; } );
             return ans;
         };
-    };
+    }
+
+    ///
+    /// - end of misc functions
+    ///
 
 } //namespace feng
 
diff --git a/tests/test.cc b/tests/test.cc
index 8969b84..8dca80d 100644
--- a/tests/test.cc
+++ b/tests/test.cc
@@ -48,7 +48,7 @@
 #include "./cases/operator_equal.hpp"
 #include "./cases/pooling.hpp"
 #include "./cases/proj.hpp"
-#include "./cases/remquo.hpp"
+//#include "./cases/remquo.hpp"
 #include "./cases/rint.hpp"
 #include "./cases/round.hpp"
 #include "./cases/sinh.hpp"