diff --git a/lib/core/array.h b/lib/core/array.h
index 431de96..77f5f1e 100644
--- a/lib/core/array.h
+++ b/lib/core/array.h
@@ -34,6 +34,7 @@
 #include "macro.h"
 #include "vector_t.h"
 #include "exceptions.h"
+#include "random"
 
 #ifdef GCTL_EIGEN
 /*The followings bypass a bug that VsCode's IntelliSense reports incorrect errors when using eigen3 library*/
@@ -1070,6 +1071,503 @@ namespace gctl
 		}
 		return os;
 	}
+
+	template <typename T>
+    void sequence(array<T> &out_arr, const T &start, const T &inc)
+    {
+        static_assert(std::is_arithmetic<T>::value, "gctl::sequence(...) could only be used with an arithmetic type.");
+
+        for (size_t i = 0; i < out_arr.size(); i++)
+        {
+            out_arr[i] = start + i*inc;
+        }
+        return;
+    }
+
+	template <typename T>
+    void linespace(const T &start, const T &end, unsigned int size, array<T> &out_arr)
+    {
+        //static_assert(std::is_arithmetic<T>::value, "gctl::linespace(...) could only be used with an arithmetic type.");
+        if (size < 1) throw invalid_argument("Invalid array size. From gctl::linespace(...)");
+
+        out_arr.resize(size);
+
+        if (size == 1)
+        {
+            out_arr[0] = 0.5*(start + end);
+            return;
+        }
+
+        T space = 1.0/(size-1)*(end - start);
+        for (int i = 0; i < size; i++)
+        {
+            out_arr[i] = start + i*space;
+        }
+        return;
+    }
+
+	template <typename T>
+    void gridspace(const T &xs, const T &xe, const T &ys, const T &ye, unsigned int xn, 
+        unsigned int yn, array<T> &out_arr)
+    {
+        //static_assert(std::is_arithmetic<T>::value, "gctl::gridspace(...) could only be used with an arithmetic type.");
+        if (xn < 1 || yn < 1) throw invalid_argument("Invalid grid size. From gctl::gridspace(...)");
+
+        array<T> out_x, out_y;
+        linespace(xs, xe, xn, out_x);
+        linespace(ys, ye, yn, out_y);
+
+        out_arr.resize(xn*yn);
+        for (int i = 0; i < yn; ++i)
+        {
+            for (int j = 0; j < xn; ++j)
+            {
+                out_arr[j+i*xn] = out_x[j] + out_y[i];
+            }
+        }
+        return;
+    }
+
+    template <typename T>
+    void meshspace(const T &xs, const T &xe, const T &ys, const T &ye, const T &zs, const T &ze, 
+        unsigned int xn, unsigned int yn, unsigned int zn, array<T> &out_arr)
+    {
+        //static_assert(std::is_arithmetic<T>::value, "gctl::meshspace(...) could only be used with an arithmetic type.");
+        if (xn < 1 || yn < 1 || zn < 1) throw invalid_argument("Invalid grid size. From gctl::meshspace(...)");
+
+        array<T> out_x, out_y, out_z;
+        linespace(xs, xe, xn, out_x);
+        linespace(ys, ye, yn, out_y);
+        linespace(zs, ze, zn, out_z);
+
+        out_arr.resize(xn*yn*zn);
+        for (int i = 0; i < zn; ++i)
+        {
+            for (int j = 0; j < yn; ++j)
+            {
+                for (int k = 0; k < xn; ++k)
+                {
+                    out_arr[k+j*xn+i*xn*yn] = out_x[k] + out_y[j] + out_z[i];
+                }
+            }
+        }
+        return;
+    }
+
+    template <typename T>
+    void logspace(const T &base, const T &start, const T &end, size_t size, array<T> &out_arr)
+    {
+        static_assert(std::is_arithmetic<T>::value, "gctl::logspace(...) could only be used with an arithmetic type.");
+        if (size < 1) throw invalid_argument("Invalid array size. From gctl::logspace(...)");
+
+        out_arr.resize(size);
+
+        if (size == 1)
+        {
+            out_arr[0] = pow(base, 0.5*(start + end));
+            return;
+        }
+
+        T space = 1.0/(size-1)*(end - start);
+        for (size_t i = 0; i < size; i++)
+        {
+            out_arr[i] = pow(base, start + i*space);
+        }
+        return;
+    }
+
+    template <typename T>
+    void linear2log(const T &base, const array<T> &in_arr, array<T> &out_arr)
+    {
+        static_assert(std::is_arithmetic<T>::value, "gctl::linear2log(...) could only be used with an arithmetic type.");
+        if (out_arr.size() != in_arr.size()) out_arr.resize(in_arr.size());
+
+        for (size_t i = 0; i < in_arr.size(); ++i)
+        {
+            out_arr[i] = log(in_arr[i])/log(base);
+        }
+        return;
+    }
+
+    template <typename T>
+    T mean(const array<T> &val_arr, const T &zero = 0)
+    {
+        static_assert(std::is_arithmetic<T>::value, "gctl::mn(...) could only be used with an arithmetic type.");
+
+        int size = val_arr.size();
+
+        T mn = zero;
+        for (int i = 0; i < size; i++)
+        {
+            mn = mn + val_arr[i];
+        }
+        return mn/size;
+    }
+
+    template <typename T>
+    T variance(const array<T> &val_arr, const T &zero = 0)
+    {
+        static_assert(std::is_arithmetic<T>::value, "gctl::variance(...) could only be used with an arithmetic type.");
+
+        T mn = mean(val_arr);
+        int size = val_arr.size();
+
+        T d = zero;
+        for (int i = 0; i < size; i++)
+        {
+            d = d + (val_arr[i] - mn)*(val_arr[i] - mn);
+        }
+        return d/size;
+    }
+
+    template <typename T>
+    T std(const array<T> &val_arr, const T &zero = 0)
+    {
+        static_assert(std::is_arithmetic<T>::value, "gctl::std(...) could only be used with an arithmetic type.");
+
+        T mn = mean(val_arr);
+        int size = val_arr.size();
+
+        T d = zero;
+        for (int i = 0; i < size; i++)
+        {
+            d = d + (val_arr[i] - mn)*(val_arr[i] - mn);
+        }
+        return sqrt(d/size);
+    }
+
+    template <typename T>
+    T std_unbiased(const array<T> &val_arr, const T &zero = 0)
+    {
+        static_assert(std::is_arithmetic<T>::value, "gctl::std_unbiased(...) could only be used with an arithmetic type.");
+
+        T mn = mean(val_arr);
+        int size = val_arr.size();
+        if (size < 2) throw std::runtime_error("[gctl::std_unbiased] Invalid array size.");
+
+        T d = zero;
+        for (int i = 0; i < size; i++)
+        {
+            d = d + (val_arr[i] - mn)*(val_arr[i] - mn);
+        }
+        return sqrt(d/(size - 1));
+    }
+
+    template <typename T>
+    T rms(const array<T> &val_arr, const T &zero = 0)
+    {
+        static_assert(std::is_arithmetic<T>::value, "gctl::rms(...) could only be used with an arithmetic type.");
+
+        int size = val_arr.size();
+
+        T mn = zero;
+        for (int i = 0; i < size; i++)
+        {
+            mn = mn + val_arr[i]*val_arr[i];
+        }
+        return sqrt(mn/size);
+    }
+
+    template <typename T>
+	T max(const array<T> &val_arr)
+	{
+		static_assert(std::is_arithmetic<T>::value, "gctl::max(...) could only be used with an arithmetic type.");
+
+		T m = val_arr[0];
+		for (size_t i = 1; i < val_arr.size(); i++)
+		{
+			m = std::max(val_arr[i], m);
+		}
+		return m;
+	}
+
+	template <typename T>
+	T min(const array<T> &val_arr)
+	{
+		static_assert(std::is_arithmetic<T>::value, "gctl::min(...) could only be used with an arithmetic type.");
+
+		T m = val_arr[0];
+		for (size_t i = 1; i < val_arr.size(); i++)
+		{
+			m = std::min(val_arr[i], m);
+		}
+		return m;
+	}
+
+    template <typename T>
+	void scale(const array<T> &val_arr, T factor)
+	{
+		static_assert(std::is_arithmetic<T>::value, "gctl::scale(...) could only be used with an arithmetic type.");
+
+		for (size_t i = 0; i < val_arr.size(); i++)
+		{
+			val_arr[i] = factor*val_arr[i];
+		}
+		return;
+	}
+
+    /**
+     * @brief 范围约束类型
+     * 
+     */
+    enum range_type_e
+    {
+        HardScale, // 所有数据按比例映射至min和max范围内
+        SoftScale, // 所有数据按比例映射至max(min, min(arr))和min(max, max(arr))范围内
+        CutOff, // 超过min和max范围数据直接设置为边界值
+    };
+
+    template <typename T>
+	void set2range(const array<T> &arr, T min, T max, range_type_e rt = HardScale)
+	{
+		static_assert(std::is_arithmetic<T>::value, "gctl::set2range(...) could only be used with an arithmetic type.");
+
+        T amin = arr[0], amax = arr[0];
+        for (size_t i = 1; i < arr.size(); i++)
+        {
+            amin = amin<arr[i]?amin:arr[i];
+            amax = amax>arr[i]?amax:arr[i];
+        }
+
+        if (rt == HardScale)
+        {
+            for (size_t i = 0; i < arr.size(); i++)
+            {
+                arr[i] = (max - min)*(arr[i] - amin)/(amax - amin) + min;
+            }
+            return;
+        }
+
+        if (amin >= min && amax <= max) return; // aleardy in the range, nothing to be done
+
+        if (rt == CutOff)
+        {
+            for (size_t i = 0; i < arr.size(); i++)
+            {
+                if (arr[i] > max) arr[i] = max;
+                if (arr[i] < min) arr[i] = min;
+            }
+            return;
+        }
+
+        // Soft Scale
+        T min2 = amin>min?amin:min;
+        T max2 = amax<max?amax:max;
+        for (size_t i = 0; i < arr.size(); i++)
+        {
+            arr[i] = (max2 - min2)*(arr[i] - amin)/(amax - amin) + min2;
+        }
+        return;
+	}
+
+    template <typename T>
+    T normalize(array<T> &in_arr, T eps = 1e-8)
+    {
+        T norm_val = module(in_arr, L2);
+        if (norm_val < eps)
+            return 0.0;
+
+        for (int i = 0; i < in_arr.size(); i++)
+        {
+            in_arr[i] /= norm_val;
+        }
+        return norm_val;
+    }
+
+    template <typename T>
+	void normalize(array<T> &val_arr, T norm, norm_type_e n_type)
+	{
+        static_assert(std::is_arithmetic<T>::value, "gctl::normalize(...) could only be used with an arithmetic type.");
+
+		T norm_sum = 0;
+        if (n_type == L0)
+        {
+			int n = 0;
+            for (size_t i = 0; i < val_arr.size(); i++)
+            {
+                if (val_arr[i] != 0)
+                {
+                    n++;
+                }
+            }
+			norm_sum = (T) n;
+        }
+
+        if (n_type == L1)
+        {
+            for (size_t i = 0; i < val_arr.size(); i++)
+            {
+                norm_sum += GCTL_FABS(val_arr[i]);
+            }
+        }
+
+        if (n_type == L2)
+        {
+            for (size_t i = 0; i < val_arr.size(); i++)
+            {
+                norm_sum += val_arr[i] * val_arr[i];
+            }
+            norm_sum = sqrt(norm_sum);
+        }
+
+		if (norm_sum <= GCTL_ZERO)
+		{
+			return;
+		}
+
+		for (size_t i = 0; i < val_arr.size(); i++)
+		{
+			val_arr[i] = val_arr[i]*norm/norm_sum;
+		}
+		return;
+	}
+
+    /**
+     * @brief Return a random number in the range [low, hig]
+     * 
+     * @note Call srand(seed) to initiate the random squence before using this function.
+     * 
+     * @tparam T Value type
+     * @param low Lower bound
+     * @param hig Higher bound
+     * @return T Random value
+     */
+    template <typename T>
+    T random(T low, T hig)
+    {
+        T f = (T) rand()/RAND_MAX;
+        return f*(hig-low)+low;
+    }
+
+    template <typename T>
+    void random(array<T> &val_arr, T p1, T p2, random_type_e mode = RdNormal, unsigned int seed = 0)
+    {
+        static_assert(std::is_arithmetic<T>::value, "gctl::random(...) could only be used with an arithmetic type.");
+
+        size_t size = val_arr.size();
+
+        if (seed == 0) seed = std::chrono::system_clock::now().time_since_epoch().count();
+
+        std::default_random_engine generator(seed);
+        if (mode == RdNormal)
+        {
+            //添加高斯分布的随机值
+            std::normal_distribution<T> dist(p1, p2);
+            for (int i = 0; i < size; i++)
+            {
+                val_arr[i] = dist(generator);
+            }
+            return;
+        }
+
+        //添加均匀分布的随机值
+        std::uniform_real_distribution<T> dist(p1, p2);
+        for (int i = 0; i < size; i++)
+        {
+            val_arr[i] = dist(generator);
+        }
+        return;
+    }
+
+	/**
+     * @brief      计算数组的模长
+     *
+     * @param[in]  in_arr      输入数组
+     * @param[in]  n_type      模长的计算类型
+     * 
+     * @return     返回的模长
+     */
+    template <typename T>
+    T module(const array<T> &in_arr, norm_type_e n_type = L2)
+    {
+        static_assert(std::is_arithmetic<T>::value, "gctl::module(...) could only be used with an arithmetic type.");
+
+        T tmp = 0.0;
+        if (n_type == L0)
+        {
+            for (int i = 0; i < in_arr.size(); i++)
+            {
+                if (in_arr[i] != 0.0)
+                {
+                    tmp += 1.0;
+                }
+            }
+            return tmp;
+        }
+
+        if (n_type == L1)
+        {
+            for (int i = 0; i < in_arr.size(); i++)
+            {
+                    tmp += GCTL_FABS(in_arr[i]);
+            }
+            return tmp;
+        }
+
+        if (n_type == L2)
+        {
+            for (int i = 0; i < in_arr.size(); i++)
+            {
+                    tmp += in_arr[i] * in_arr[i];
+            }
+            return sqrt(tmp);
+        }
+
+        tmp = in_arr[0];
+        for (int i = 1; i < in_arr.size(); i++)
+        {
+                tmp = GCTL_MAX(tmp, in_arr[i]);
+        }
+        return sqrt(tmp);
+    }
+
+    /**
+	 * @brief      计算两个数组的点积
+	 *
+	 * @param[in]  a     输入的第一个数组
+	 * @param[in]  b     输入的第二个数组
+	 *
+	 * @return     数组的点积
+	 */
+    template <typename T>
+    T dot_product(const array<T> &a, const array<T> &b)
+    {
+        static_assert(std::is_arithmetic<T>::value, "gctl::dot_product(...) could only be used with an arithmetic type.");
+
+        if (a.size() != b.size())
+            throw runtime_error("Incompatible arrays' sizes. Thrown by gctl::dot_product(...)");
+
+        T p = (T) 0;
+        for (size_t i = 0; i < a.size(); i++)
+        {
+            p = p + a[i]*b[i];
+        }
+        return p;
+    }
+
+    /**
+     * @brief      计算一个向量相对于另一个向量的正交向量
+     *
+     * @param[in]  a     输入的第一个数组
+     * @param      b     输入的第二个数组，计算后的正交向量以引用的形式返回
+     */
+    template <typename T>
+    void orth(const array<T> &a, array<T> &b)
+    {
+        static_assert(std::is_arithmetic<T>::value, "gctl::orth(...) could only be used with an arithmetic type.");
+
+        if (a.size() != b.size())
+        {
+            throw runtime_error("The arrays have different sizes. Thrown by gctl::orth(...)");
+        }
+
+        T product = dot_product(a, b);
+        for (int i = 0; i < a.size(); i++)
+        {
+            b[i] -= product*a[i];
+        }
+        return;
+    }
 }
 
 #endif // _GCTL_ARRAY_H
diff --git a/lib/maths/mathfunc_t.h b/lib/maths/mathfunc_t.h
index 680deaf..e21db46 100644
--- a/lib/maths/mathfunc_t.h
+++ b/lib/maths/mathfunc_t.h
@@ -32,9 +32,6 @@
 #include "../core/array.h"
 #include "../core/matrix.h"
 #include "../core/exceptions.h"
-
-#include "vector"
-#include "random"
 #include "cmath"
 
 namespace gctl
@@ -153,404 +150,11 @@ namespace gctl
         }
         return ang;
     }
+}
 
-    template <typename T>
-    void sequence(array<T> &out_arr, const T &start, const T &inc)
-    {
-        static_assert(std::is_arithmetic<T>::value, "gctl::sequence(...) could only be used with an arithmetic type.");
-
-        for (size_t i = 0; i < out_arr.size(); i++)
-        {
-            out_arr[i] = start + i*inc;
-        }
-        return;
-    }
-
-    template <typename T>
-    void linespace(const T &start, const T &end, unsigned int size, array<T> &out_arr)
-    {
-        //static_assert(std::is_arithmetic<T>::value, "gctl::linespace(...) could only be used with an arithmetic type.");
-        if (size < 1) throw invalid_argument("Invalid array size. From gctl::linespace(...)");
-
-        out_arr.resize(size);
-
-        if (size == 1)
-        {
-            out_arr[0] = 0.5*(start + end);
-            return;
-        }
-
-        T space = 1.0/(size-1)*(end - start);
-        for (int i = 0; i < size; i++)
-        {
-            out_arr[i] = start + i*space;
-        }
-        return;
-    }
-
-    template <typename T>
-    void gridspace(const T &xs, const T &xe, const T &ys, const T &ye, unsigned int xn, 
-        unsigned int yn, array<T> &out_arr)
-    {
-        //static_assert(std::is_arithmetic<T>::value, "gctl::gridspace(...) could only be used with an arithmetic type.");
-        if (xn < 1 || yn < 1) throw invalid_argument("Invalid grid size. From gctl::gridspace(...)");
-
-        array<T> out_x, out_y;
-        linespace(xs, xe, xn, out_x);
-        linespace(ys, ye, yn, out_y);
-
-        out_arr.resize(xn*yn);
-        for (int i = 0; i < yn; ++i)
-        {
-            for (int j = 0; j < xn; ++j)
-            {
-                out_arr[j+i*xn] = out_x[j] + out_y[i];
-            }
-        }
-        return;
-    }
-
-    template <typename T>
-    void meshspace(const T &xs, const T &xe, const T &ys, const T &ye, const T &zs, const T &ze, 
-        unsigned int xn, unsigned int yn, unsigned int zn, array<T> &out_arr)
-    {
-        //static_assert(std::is_arithmetic<T>::value, "gctl::meshspace(...) could only be used with an arithmetic type.");
-        if (xn < 1 || yn < 1 || zn < 1) throw invalid_argument("Invalid grid size. From gctl::meshspace(...)");
-
-        array<T> out_x, out_y, out_z;
-        linespace(xs, xe, xn, out_x);
-        linespace(ys, ye, yn, out_y);
-        linespace(zs, ze, zn, out_z);
-
-        out_arr.resize(xn*yn*zn);
-        for (int i = 0; i < zn; ++i)
-        {
-            for (int j = 0; j < yn; ++j)
-            {
-                for (int k = 0; k < xn; ++k)
-                {
-                    out_arr[k+j*xn+i*xn*yn] = out_x[k] + out_y[j] + out_z[i];
-                }
-            }
-        }
-        return;
-    }
-
-    template <typename T>
-    void logspace(const T &base, const T &start, const T &end, size_t size, array<T> &out_arr)
-    {
-        static_assert(std::is_arithmetic<T>::value, "gctl::logspace(...) could only be used with an arithmetic type.");
-        if (size < 1) throw invalid_argument("Invalid array size. From gctl::logspace(...)");
-
-        out_arr.resize(size);
-
-        if (size == 1)
-        {
-            out_arr[0] = pow(base, 0.5*(start + end));
-            return;
-        }
-
-        T space = 1.0/(size-1)*(end - start);
-        for (size_t i = 0; i < size; i++)
-        {
-            out_arr[i] = pow(base, start + i*space);
-        }
-        return;
-    }
-
-    template <typename T>
-    void linear2log(const T &base, const array<T> &in_arr, array<T> &out_arr)
-    {
-        static_assert(std::is_arithmetic<T>::value, "gctl::linear2log(...) could only be used with an arithmetic type.");
-        if (out_arr.size() != in_arr.size()) out_arr.resize(in_arr.size());
-
-        for (size_t i = 0; i < in_arr.size(); ++i)
-        {
-            out_arr[i] = log(in_arr[i])/log(base);
-        }
-        return;
-    }
-
-    template <typename T>
-    T mean(const array<T> &val_arr, const T &zero = 0)
-    {
-        static_assert(std::is_arithmetic<T>::value, "gctl::mn(...) could only be used with an arithmetic type.");
-
-        int size = val_arr.size();
-
-        T mn = zero;
-        for (int i = 0; i < size; i++)
-        {
-            mn = mn + val_arr[i];
-        }
-        return mn/size;
-    }
-
-    template <typename T>
-    T variance(const array<T> &val_arr, const T &zero = 0)
-    {
-        static_assert(std::is_arithmetic<T>::value, "gctl::variance(...) could only be used with an arithmetic type.");
-
-        T mn = mean(val_arr);
-        int size = val_arr.size();
-
-        T d = zero;
-        for (int i = 0; i < size; i++)
-        {
-            d = d + (val_arr[i] - mn)*(val_arr[i] - mn);
-        }
-        return d/size;
-    }
-
-    template <typename T>
-    T std(const array<T> &val_arr, const T &zero = 0)
-    {
-        static_assert(std::is_arithmetic<T>::value, "gctl::std(...) could only be used with an arithmetic type.");
-
-        T mn = mean(val_arr);
-        int size = val_arr.size();
-
-        T d = zero;
-        for (int i = 0; i < size; i++)
-        {
-            d = d + (val_arr[i] - mn)*(val_arr[i] - mn);
-        }
-        return sqrt(d/size);
-    }
-
-    template <typename T>
-    T std_unbiased(const array<T> &val_arr, const T &zero = 0)
-    {
-        static_assert(std::is_arithmetic<T>::value, "gctl::std_unbiased(...) could only be used with an arithmetic type.");
-
-        T mn = mean(val_arr);
-        int size = val_arr.size();
-        if (size < 2) throw std::runtime_error("[gctl::std_unbiased] Invalid array size.");
-
-        T d = zero;
-        for (int i = 0; i < size; i++)
-        {
-            d = d + (val_arr[i] - mn)*(val_arr[i] - mn);
-        }
-        return sqrt(d/(size - 1));
-    }
-
-    template <typename T>
-    T rms(const array<T> &val_arr, const T &zero = 0)
-    {
-        static_assert(std::is_arithmetic<T>::value, "gctl::rms(...) could only be used with an arithmetic type.");
-
-        int size = val_arr.size();
-
-        T mn = zero;
-        for (int i = 0; i < size; i++)
-        {
-            mn = mn + val_arr[i]*val_arr[i];
-        }
-        return sqrt(mn/size);
-    }
-
-    template <typename T>
-	T max(const array<T> &val_arr)
-	{
-		static_assert(std::is_arithmetic<T>::value, "gctl::max(...) could only be used with an arithmetic type.");
-
-		T m = val_arr[0];
-		for (size_t i = 1; i < val_arr.size(); i++)
-		{
-			m = std::max(val_arr[i], m);
-		}
-		return m;
-	}
-
-	template <typename T>
-	T min(const array<T> &val_arr)
-	{
-		static_assert(std::is_arithmetic<T>::value, "gctl::min(...) could only be used with an arithmetic type.");
-
-		T m = val_arr[0];
-		for (size_t i = 1; i < val_arr.size(); i++)
-		{
-			m = std::min(val_arr[i], m);
-		}
-		return m;
-	}
-
-    template <typename T>
-	void scale(const array<T> &val_arr, T factor)
-	{
-		static_assert(std::is_arithmetic<T>::value, "gctl::scale(...) could only be used with an arithmetic type.");
-
-		for (size_t i = 0; i < val_arr.size(); i++)
-		{
-			val_arr[i] = factor*val_arr[i];
-		}
-		return;
-	}
-
-    /**
-     * @brief 范围约束类型
-     * 
-     */
-    enum range_type_e
-    {
-        HardScale, // 所有数据按比例映射至min和max范围内
-        SoftScale, // 所有数据按比例映射至max(min, min(arr))和min(max, max(arr))范围内
-        CutOff, // 超过min和max范围数据直接设置为边界值
-    };
-
-    template <typename T>
-	void set2range(const array<T> &arr, T min, T max, range_type_e rt = HardScale)
-	{
-		static_assert(std::is_arithmetic<T>::value, "gctl::set2range(...) could only be used with an arithmetic type.");
-
-        T amin = arr[0], amax = arr[0];
-        for (size_t i = 1; i < arr.size(); i++)
-        {
-            amin = amin<arr[i]?amin:arr[i];
-            amax = amax>arr[i]?amax:arr[i];
-        }
-
-        if (rt == HardScale)
-        {
-            for (size_t i = 0; i < arr.size(); i++)
-            {
-                arr[i] = (max - min)*(arr[i] - amin)/(amax - amin) + min;
-            }
-            return;
-        }
-
-        if (amin >= min && amax <= max) return; // aleardy in the range, nothing to be done
-
-        if (rt == CutOff)
-        {
-            for (size_t i = 0; i < arr.size(); i++)
-            {
-                if (arr[i] > max) arr[i] = max;
-                if (arr[i] < min) arr[i] = min;
-            }
-            return;
-        }
-
-        // Soft Scale
-        T min2 = amin>min?amin:min;
-        T max2 = amax<max?amax:max;
-        for (size_t i = 0; i < arr.size(); i++)
-        {
-            arr[i] = (max2 - min2)*(arr[i] - amin)/(amax - amin) + min2;
-        }
-        return;
-	}
-
-    template <typename T>
-    T normalize(array<T> &in_arr, T eps = 1e-8)
-    {
-        T norm_val = module(in_arr, L2);
-        if (norm_val < eps)
-            return 0.0;
-
-        for (int i = 0; i < in_arr.size(); i++)
-        {
-            in_arr[i] /= norm_val;
-        }
-        return norm_val;
-    }
-
-    template <typename T>
-	void normalize(array<T> &val_arr, T norm, norm_type_e n_type)
-	{
-        static_assert(std::is_arithmetic<T>::value, "gctl::normalize(...) could only be used with an arithmetic type.");
-
-		T norm_sum = 0;
-        if (n_type == L0)
-        {
-			int n = 0;
-            for (size_t i = 0; i < val_arr.size(); i++)
-            {
-                if (val_arr[i] != 0)
-                {
-                    n++;
-                }
-            }
-			norm_sum = (T) n;
-        }
-
-        if (n_type == L1)
-        {
-            for (size_t i = 0; i < val_arr.size(); i++)
-            {
-                norm_sum += GCTL_FABS(val_arr[i]);
-            }
-        }
-
-        if (n_type == L2)
-        {
-            for (size_t i = 0; i < val_arr.size(); i++)
-            {
-                norm_sum += val_arr[i] * val_arr[i];
-            }
-            norm_sum = sqrt(norm_sum);
-        }
-
-		if (norm_sum <= GCTL_ZERO)
-		{
-			return;
-		}
-
-		for (size_t i = 0; i < val_arr.size(); i++)
-		{
-			val_arr[i] = val_arr[i]*norm/norm_sum;
-		}
-		return;
-	}
-
-    /**
-     * @brief Return a random number in the range [low, hig]
-     * 
-     * @note Call srand(seed) to initiate the random squence before using this function.
-     * 
-     * @tparam T Value type
-     * @param low Lower bound
-     * @param hig Higher bound
-     * @return T Random value
-     */
-    template <typename T>
-    T random(T low, T hig)
-    {
-        T f = (T) rand()/RAND_MAX;
-        return f*(hig-low)+low;
-    }
-
-    template <typename T>
-    void random(array<T> &val_arr, T p1, T p2, random_type_e mode = RdNormal, unsigned int seed = 0)
-    {
-        static_assert(std::is_arithmetic<T>::value, "gctl::random(...) could only be used with an arithmetic type.");
-
-        size_t size = val_arr.size();
-
-        if (seed == 0) seed = std::chrono::system_clock::now().time_since_epoch().count();
-
-        std::default_random_engine generator(seed);
-        if (mode == RdNormal)
-        {
-            //添加高斯分布的随机值
-            std::normal_distribution<T> dist(p1, p2);
-            for (int i = 0; i < size; i++)
-            {
-                val_arr[i] = dist(generator);
-            }
-            return;
-        }
-
-        //添加均匀分布的随机值
-        std::uniform_real_distribution<T> dist(p1, p2);
-        for (int i = 0; i < size; i++)
-        {
-            val_arr[i] = dist(generator);
-        }
-        return;
-    }
+#endif //_GCTL_MATHFUNC_TEMPLATE_H
 
+/*
     template <typename T>
     void matrix_vector(const matrix<T> &mat, const array<T> &vec, array<T> &ret, T zero = 0, matrix_layout_e layout = NoTrans)
     {
@@ -619,110 +223,6 @@ namespace gctl
         return;
     }
 
-    /**
-     * @brief      计算数组的模长
-     *
-     * @param[in]  in_arr      输入数组
-     * @param[in]  n_type      模长的计算类型
-     * 
-     * @return     返回的模长
-     */
-    template <typename T>
-    T module(const array<T> &in_arr, norm_type_e n_type = L2)
-    {
-        static_assert(std::is_arithmetic<T>::value, "gctl::module(...) could only be used with an arithmetic type.");
-
-        T tmp = 0.0;
-        if (n_type == L0)
-        {
-            for (int i = 0; i < in_arr.size(); i++)
-            {
-                if (in_arr[i] != 0.0)
-                {
-                    tmp += 1.0;
-                }
-            }
-            return tmp;
-        }
-
-        if (n_type == L1)
-        {
-            for (int i = 0; i < in_arr.size(); i++)
-            {
-                    tmp += GCTL_FABS(in_arr[i]);
-            }
-            return tmp;
-        }
-
-        if (n_type == L2)
-        {
-            for (int i = 0; i < in_arr.size(); i++)
-            {
-                    tmp += in_arr[i] * in_arr[i];
-            }
-            return sqrt(tmp);
-        }
-
-        tmp = in_arr[0];
-        for (int i = 1; i < in_arr.size(); i++)
-        {
-                tmp = GCTL_MAX(tmp, in_arr[i]);
-        }
-        return sqrt(tmp);
-    }
-
-    /**
-	 * @brief      计算两个数组的点积
-	 *
-	 * @param[in]  a     输入的第一个数组
-	 * @param[in]  b     输入的第二个数组
-	 *
-	 * @return     数组的点积
-	 */
-    template <typename T>
-    T dot_product(const array<T> &a, const array<T> &b)
-    {
-        static_assert(std::is_arithmetic<T>::value, "gctl::dot_product(...) could only be used with an arithmetic type.");
-
-        if (a.size() != b.size())
-            throw runtime_error("Incompatible arrays' sizes. Thrown by gctl::dot_product(...)");
-
-        T p = (T) 0;
-        for (size_t i = 0; i < a.size(); i++)
-        {
-            p = p + a[i]*b[i];
-        }
-        return p;
-    }
-
-    /**
-     * @brief      计算一个向量相对于另一个向量的正交向量
-     *
-     * @param[in]  a     输入的第一个数组
-     * @param      b     输入的第二个数组，计算后的正交向量以引用的形式返回
-     */
-    template <typename T>
-    void orth(const array<T> &a, array<T> &b)
-    {
-        static_assert(std::is_arithmetic<T>::value, "gctl::orth(...) could only be used with an arithmetic type.");
-
-        if (a.size() != b.size())
-        {
-            throw runtime_error("The arrays have different sizes. Thrown by gctl::orth(...)");
-        }
-
-        T product = dot_product(a, b);
-        for (int i = 0; i < a.size(); i++)
-        {
-            b[i] -= product*a[i];
-        }
-        return;
-    }
-}
-
-#endif //_GCTL_MATHFUNC_TEMPLATE_H
-
-/*
     template <typename T>
     void linespace(const T &start, const T &end, unsigned int size, std::vector<T> &out_vec)
     {