diff --git a/Cargo.toml b/Cargo.toml
index f9495b48..9861577d 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "peroxide"
-version = "0.34.1"
+version = "0.34.2"
 authors = ["axect <axect@outlook.kr>"]
 edition = "2018"
 description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax"
@@ -27,6 +27,7 @@ order-stat = "0.1"
 puruspe = "0.2"
 matrixmultiply = { version = "0.3", features = ["threading"] }
 peroxide-ad = "0.3"
+peroxide-num = "0.1"
 #num-complex = "0.3"
 netcdf = { version = "0.7", optional = true, default_features = false }
 pyo3 = { version = "0.19", optional = true }
@@ -34,7 +35,7 @@ blas = { version = "0.22", optional = true }
 lapack = { version = "0.19", optional = true }
 serde = { version = "1.0", features = ["derive"], optional = true }
 json = { version = "0.12", optional = true }
-arrow2 = { version = "0.17", features = ["io_parquet", "io_parquet_compression"], optional = true }
+arrow2 = { version = "0.18", features = ["io_parquet", "io_parquet_compression"], optional = true }
 
 [package.metadata.docs.rs]
 rustdoc-args = [ "--html-in-header", "katex-header.html", "--cfg", "docsrs"]
diff --git a/RELEASES.md b/RELEASES.md
index b233812f..19cdc54c 100644
--- a/RELEASES.md
+++ b/RELEASES.md
@@ -1,3 +1,18 @@
+# Release 0.34.2 (2023-11-20)
+
+## New sub-crate : `peroxide-num`
+
+* Add new sub-crate : `peroxide-num`
+* Change all dependencies of `ExpLogOps, PowOps, TrigOps` to `peroxide-num`
+
+## Fix errata
+
+* R example in `structure/matrix` (#56) (Thanks to [@rdavis120](https://github.com/rdavis120))
+
+## Fix old syntax
+
+* Fix old syntax - e.g. explicit `into_iter`, `Vec::with_capacity` & `set_len`
+
 # Release 0.34.1 (2023-08-03)
 
 ## Modify Statistics of `WeightedUniform`
diff --git a/examples/ad.rs b/examples/ad.rs
index ae17b89e..08149163 100644
--- a/examples/ad.rs
+++ b/examples/ad.rs
@@ -32,7 +32,6 @@ fn main() {
     (b - 1f64).print();
     (b * 2f64).print();
     (b / 2f64).print();
-    assert_eq!(1f64 + b, b + 1f64);
     assert_eq!(-(1f64 - b), b - 1f64);
     assert_eq!(2f64 * b, b * 2f64);
     assert_eq!(1f64 / (2f64 / b), b / 2f64);
diff --git a/peroxide-num/.gitignore b/peroxide-num/.gitignore
new file mode 100644
index 00000000..2f7896d1
--- /dev/null
+++ b/peroxide-num/.gitignore
@@ -0,0 +1 @@
+target/
diff --git a/peroxide-num/Cargo.toml b/peroxide-num/Cargo.toml
new file mode 100644
index 00000000..46be56fc
--- /dev/null
+++ b/peroxide-num/Cargo.toml
@@ -0,0 +1,16 @@
+[package]
+name = "peroxide-num"
+version = "0.1.3"
+edition = "2021"
+authors = ["Axect <ax2ct@outlook.com>"]
+description = "Numerical traits for Peroxide"
+repository = "https://github.com/Axect/Peroxide"
+license = "MIT OR Apache-2.0"
+categories = ["science"]
+readme = "README.md"
+keywords = ["numerical", "mathematics"]
+exclude = ["examples/"]
+
+# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+
+[dependencies]
diff --git a/peroxide-num/README.md b/peroxide-num/README.md
new file mode 100644
index 00000000..62277328
--- /dev/null
+++ b/peroxide-num/README.md
@@ -0,0 +1,72 @@
+# Peroxide-num
+
+`Peroxide-num` is a Rust crate dedicated to providing comprehensive numeric structures and operations, extending generic programming capabilities with a focus on mathematical computations.
+
+## Overview
+
+This crate defines a set of traits for numeric operations, including basic arithmetic, power functions, trigonometric functions, exponential, and logarithmic functions. It is designed to be used with numeric types that implement these traits, allowing for a wide range of mathematical operations.
+
+## Traits
+
+- `PowOps`: Operations related to powers and roots.
+- `TrigOps`: Trigonometric functions.
+- `ExpLogOps`: Exponential and logarithmic functions.
+- `Float`: Define own floating point type (`f32` and `f64` are implemented as default).
+- `Numeric`: A comprehensive trait that encompasses all of the above along with standard arithmetic operations.
+
+## Example 1: Defining a Simple Numeric Type
+
+Below is an example of how you can define your own simple numeric type that implements the `Numeric` trait.
+
+```rust
+#[derive(Debug, Clone, Copy, PartialOrd)]
+struct SimpleNumber(f64);
+
+impl PowOps for SimpleNumber {
+    type Float = Self;
+
+    fn powi(&self, n: i32) -> Self {
+        SimpleNumber(self.0.powi(n))
+    }
+
+    fn powf(&self, f: Self::Float) -> Self {
+        SimpleNumber(self.0.powf(f.0))
+    }
+
+    fn pow(&self, f: Self) -> Self {
+        SimpleNumber(self.0.powf(f.0))
+    }
+
+    fn sqrt(&self) -> Self {
+        SimpleNumber(self.0.sqrt())
+    }
+}
+
+// Implement other required operations for SimpleNumber...
+// - Add, Sub, Mul, Div, Neg
+// - PowOps (implemented above)
+// - TrigOps
+// - ExpLogOps
+
+impl Numeric<f64> for SimpleNumber {}
+```
+
+This `SimpleNumber` struct wraps a `f64` and implements the `Numeric` trait, making it capable of all the operations defined in the `Peroxide-num` crate.
+
+## Example 2: `Vec3D` (3D Vector)
+
+[examples/vec3d.rs](examples/vec3d.rs)
+
+## Usage
+
+To use this type in your own computations:
+
+```rust
+let num = SimpleNumber(2.0);
+let result = num.sin(); // Compute the sine of 2.0
+println!("{:?}", result); // Should display the sine of 2.0
+```
+
+The `Peroxide-num` crate is designed to be flexible and extensible, allowing for easy integration with the larger `Peroxide` ecosystem for scientific computing in Rust.
+
+For more information and advanced usage, please refer to the documentation and examples provided with the crate.
diff --git a/peroxide-num/RELEASE.md b/peroxide-num/RELEASE.md
new file mode 100644
index 00000000..bbda1dac
--- /dev/null
+++ b/peroxide-num/RELEASE.md
@@ -0,0 +1,18 @@
+# Ver 0.1.3 (2023-11-09)
+
+- Add `Group` and `Ring`
+  - `Group`: require `Add<Self, Output=Self>` and has `zero()`
+  - `Ring`: require `Group + Mul<Self, Output=Self>` and has `one()`
+  - `f32` and `f64` are `Ring`
+
+# Ver 0.1.2 (2023-11-09)
+
+- Remove unnecessary constraints for `Numeric<T>`
+  - `T: Add<Self, Output=Self> + Sub<Self, Output=Self> + Mul<Self, Output=Self> + Div<Self, Output=Self>`
+- More specify `Neg` for `Numeric` & `Float`: `Neg` -> `Neg<Output=Self>` 
+
+# Ver 0.1.1 (2023-11-09)
+
+- Add constraints for `Numeric<T>`
+  - `T: Add<Self, Output=Self> + Sub<Self, Output=Self> + Mul<Self, Output=Self> + Div<Self, Output=Self>`
+- Add example : `Vec3D`
diff --git a/peroxide-num/examples/vec3d.rs b/peroxide-num/examples/vec3d.rs
new file mode 100644
index 00000000..851c5d1b
--- /dev/null
+++ b/peroxide-num/examples/vec3d.rs
@@ -0,0 +1,246 @@
+use peroxide_num::{PowOps, ExpLogOps, TrigOps, Numeric};
+use std::ops::{Neg, Add, Sub, Mul, Div};
+use std::f64::consts::PI;
+
+fn main() {
+    let a = Vec3D::new(0f64, 0.5 * PI, PI);
+    println!("{:?}", a);
+    println!("{:?}", -a);
+    println!("{:?}", a + a);
+    println!("{:?}", a - a);
+    println!("{:?}", a * 2f64);
+    println!("{:?}", 2f64 * a);
+    println!("{:?}", a / 2f64);
+    println!("{:?}", 2f64 / a);
+    println!("{:?}", a.sin());
+    println!("{:?}", a.cos());
+    println!("{:?}", a.tan());
+    println!("{:?}", a.asin());
+    println!("{:?}", a.acos());
+    println!("{:?}", a.atan());
+    println!("{:?}", a.sinh());
+    println!("{:?}", a.cosh());
+    println!("{:?}", a.tanh());
+    println!("{:?}", a.asinh());
+    println!("{:?}", a.acosh());
+    println!("{:?}", a.atanh());
+    println!("{:?}", a.exp());
+    println!("{:?}", a.ln());
+}
+
+#[derive(Debug, Copy, Clone)]
+struct Vec3D {
+    x: f64,
+    y: f64,
+    z: f64
+}
+
+impl Vec3D {
+    fn new(x: f64, y: f64, z: f64) -> Vec3D {
+        Vec3D { x, y, z }
+    }
+}
+
+impl Neg for Vec3D {
+    type Output = Vec3D;
+
+    fn neg(self) -> Vec3D {
+        Vec3D::new(-self.x, -self.y, -self.z)
+    }
+}
+
+impl Add for Vec3D {
+    type Output = Vec3D;
+
+    fn add(self, other: Vec3D) -> Vec3D {
+        Vec3D::new(self.x + other.x, self.y + other.y, self.z + other.z)
+    }
+}
+
+impl Sub for Vec3D {
+    type Output = Vec3D;
+
+    fn sub(self, other: Vec3D) -> Vec3D {
+        Vec3D::new(self.x - other.x, self.y - other.y, self.z - other.z)
+    }
+}
+
+impl Mul for Vec3D {
+    type Output = Vec3D;
+
+    fn mul(self, _other: Vec3D) -> Vec3D {
+        unimplemented!()
+    }
+}
+
+impl Div for Vec3D {
+    type Output = Vec3D;
+
+    fn div(self, _other: Vec3D) -> Vec3D {
+        unimplemented!()
+    }
+}
+
+impl Add<f64> for Vec3D {
+    type Output = Vec3D;
+
+    fn add(self, scalar: f64) -> Vec3D {
+        Vec3D::new(self.x + scalar, self.y + scalar, self.z + scalar)
+    }
+}
+
+impl Sub<f64> for Vec3D {
+    type Output = Vec3D;
+
+    fn sub(self, scalar: f64) -> Vec3D {
+        Vec3D::new(self.x - scalar, self.y - scalar, self.z - scalar)
+    }
+}
+
+impl Mul<f64> for Vec3D {
+    type Output = Vec3D;
+
+    fn mul(self, scalar: f64) -> Vec3D {
+        Vec3D::new(self.x * scalar, self.y * scalar, self.z * scalar)
+    }
+}
+
+impl Div<f64> for Vec3D {
+    type Output = Vec3D;
+
+    fn div(self, scalar: f64) -> Vec3D {
+        Vec3D::new(self.x / scalar, self.y / scalar, self.z / scalar)
+    }
+}
+
+impl Add<Vec3D> for f64 {
+    type Output = Vec3D;
+
+    fn add(self, other: Vec3D) -> Vec3D {
+        Vec3D::new(self + other.x, self + other.y, self + other.z)
+    }
+}
+
+impl Sub<Vec3D> for f64 {
+    type Output = Vec3D;
+
+    fn sub(self, other: Vec3D) -> Vec3D {
+        Vec3D::new(self - other.x, self - other.y, self - other.z)
+    }
+}
+
+impl Mul<Vec3D> for f64 {
+    type Output = Vec3D;
+
+    fn mul(self, other: Vec3D) -> Vec3D {
+        Vec3D::new(self * other.x, self * other.y, self * other.z)
+    }
+}
+
+impl Div<Vec3D> for f64 {
+    type Output = Vec3D;
+
+    fn div(self, other: Vec3D) -> Vec3D {
+        Vec3D::new(self / other.x, self / other.y, self / other.z)
+    }
+}
+
+impl PowOps for Vec3D {
+    type Float = f64;
+
+    fn pow(&self, _power: Self) -> Self {
+        unimplemented!()
+    }
+
+    fn powf(&self, power: Self::Float) -> Self {
+        Vec3D::new(self.x.powf(power), self.y.powf(power), self.z.powf(power))
+    }
+
+    fn powi(&self, power: i32) -> Self {
+        Vec3D::new(self.x.powi(power), self.y.powi(power), self.z.powi(power))
+    }
+
+    fn sqrt(&self) -> Self {
+        Vec3D::new(self.x.sqrt(), self.y.sqrt(), self.z.sqrt())
+    }
+}
+
+impl TrigOps for Vec3D {
+    fn sin(&self) -> Self {
+        Vec3D::new(self.x.sin(), self.y.sin(), self.z.sin())
+    }
+
+    fn cos(&self) -> Self {
+        Vec3D::new(self.x.cos(), self.y.cos(), self.z.cos())
+    }
+
+    fn tan(&self) -> Self {
+        Vec3D::new(self.x.tan(), self.y.tan(), self.z.tan())
+    }
+
+    fn asin(&self) -> Self {
+        Vec3D::new(self.x.asin(), self.y.asin(), self.z.asin())
+    }
+
+    fn acos(&self) -> Self {
+        Vec3D::new(self.x.acos(), self.y.acos(), self.z.acos())
+    }
+
+    fn atan(&self) -> Self {
+        Vec3D::new(self.x.atan(), self.y.atan(), self.z.atan())
+    }
+
+    fn sinh(&self) -> Self {
+        Vec3D::new(self.x.sinh(), self.y.sinh(), self.z.sinh())
+    }
+
+    fn cosh(&self) -> Self {
+        Vec3D::new(self.x.cosh(), self.y.cosh(), self.z.cosh())
+    }
+
+    fn tanh(&self) -> Self {
+        Vec3D::new(self.x.tanh(), self.y.tanh(), self.z.tanh())
+    }
+
+    fn asinh(&self) -> Self {
+        Vec3D::new(self.x.asinh(), self.y.asinh(), self.z.asinh())
+    }
+
+    fn acosh(&self) -> Self {
+        Vec3D::new(self.x.acosh(), self.y.acosh(), self.z.acosh())
+    }
+
+    fn atanh(&self) -> Self {
+        Vec3D::new(self.x.atanh(), self.y.atanh(), self.z.atanh())
+    }
+
+    fn sin_cos(&self) -> (Self, Self) {
+        (self.sin(), self.cos())
+    }
+}
+
+impl ExpLogOps for Vec3D {
+    type Float = f64;
+
+    fn exp(&self) -> Self {
+        Vec3D::new(self.x.exp(), self.y.exp(), self.z.exp())
+    }
+
+    fn ln(&self) -> Self {
+        Vec3D::new(self.x.ln(), self.y.ln(), self.z.ln())
+    }
+
+    fn log(&self, base: Self::Float) -> Self {
+        Vec3D::new(self.x.log(base), self.y.log(base), self.z.log(base))
+    }
+
+    fn log2(&self) -> Self {
+        Vec3D::new(self.x.log2(), self.y.log2(), self.z.log2())
+    }
+
+    fn log10(&self) -> Self {
+        Vec3D::new(self.x.log10(), self.y.log10(), self.z.log10())
+    }
+}
+
+impl Numeric<f64> for Vec3D {}
diff --git a/peroxide-num/src/lib.rs b/peroxide-num/src/lib.rs
new file mode 100644
index 00000000..7415bd8c
--- /dev/null
+++ b/peroxide-num/src/lib.rs
@@ -0,0 +1,231 @@
+//! Missing operations & comprehensive number structures
+//!
+//! ## `Numeric` trait
+//!
+//! * `Numeric` requires `PowOps, TrigOps, ExpLogOps` & `std::Ops<Self>` & `std::Ops<f64>`
+
+use std::ops::{Add, Div, Mul, Neg, Sub};
+
+pub trait PowOps: Sized {
+    type Float;
+    fn powi(&self, n: i32) -> Self;
+    fn powf(&self, f: Self::Float) -> Self;
+    fn pow(&self, f: Self) -> Self;
+    fn sqrt(&self) -> Self;
+}
+
+pub trait TrigOps: Sized + Div<Output = Self> {
+    fn sin_cos(&self) -> (Self, Self);
+    fn sin(&self) -> Self {
+        let (s, _) = self.sin_cos();
+        s
+    }
+    fn cos(&self) -> Self {
+        let (_, c) = self.sin_cos();
+        c
+    }
+    fn tan(&self) -> Self {
+        let (s, c) = self.sin_cos();
+        s / c
+    }
+    fn sinh(&self) -> Self;
+    fn cosh(&self) -> Self;
+    fn tanh(&self) -> Self;
+    fn asin(&self) -> Self;
+    fn acos(&self) -> Self;
+    fn atan(&self) -> Self;
+    fn asin_acos(&self) -> (Self, Self) {
+        (self.asin(), self.acos())
+    }
+    fn asinh(&self) -> Self;
+    fn acosh(&self) -> Self;
+    fn atanh(&self) -> Self;
+    fn asinh_acosh(&self) -> (Self, Self) {
+        (self.asinh(), self.acosh())
+    }
+}
+
+pub trait ExpLogOps: Sized {
+    type Float;
+    fn exp(&self) -> Self;
+    fn ln(&self) -> Self;
+    fn log(&self, base: Self::Float) -> Self;
+    fn log2(&self) -> Self;
+    fn log10(&self) -> Self;
+}
+
+pub trait Float:
+    PowOps<Float = Self>
+    + TrigOps
+    + ExpLogOps<Float = Self>
+    + Neg<Output = Self>
+    + Add<Output = Self>
+    + Mul<Output = Self>
+    + Div<Output = Self>
+    + Sub<Output = Self>
+    + PartialOrd
+    + Copy
+    + Clone
+{
+}
+
+pub trait Numeric<T: Float>:
+    PowOps<Float = T>
+    + TrigOps
+    + ExpLogOps<Float = T>
+    + Neg<Output = Self>
+    + Add<Output = Self>
+    + Mul<Output = Self>
+    + Div<Output = Self>
+    + Sub<Output = Self>
+    + Add<T, Output = Self>
+    + Mul<T, Output = Self>
+    + Div<T, Output = Self>
+    + Sub<T, Output = Self>
+    + Clone
+{
+}
+
+// =============================================================================
+// Numeric Traits for f32, f64
+// =============================================================================
+
+macro_rules! impl_float {
+    ($($t:ty),*) => {
+        $(
+            impl PowOps for $t {
+                type Float = $t;
+                fn powi(&self, n: i32) -> Self {
+                    (*self).powi(n)
+                }
+
+                fn powf(&self, f: Self::Float) -> Self {
+                    (*self).powf(f)
+                }
+
+                fn pow(&self, f: Self) -> Self {
+                    (*self).powf(f)
+                }
+
+                fn sqrt(&self) -> Self {
+                    (*self).sqrt()
+                }
+            }
+
+            impl TrigOps for $t {
+                fn sin_cos(&self) -> (Self, Self) {
+                    (*self).sin_cos()
+                }
+
+                fn sin(&self) -> Self {
+                    (*self).sin()
+                }
+
+                fn cos(&self) -> Self {
+                    (*self).cos()
+                }
+
+                fn tan(&self) -> Self {
+                    (*self).tan()
+                }
+
+                fn sinh(&self) -> Self {
+                    (*self).sinh()
+                }
+
+                fn cosh(&self) -> Self {
+                    (*self).cosh()
+                }
+
+                fn tanh(&self) -> Self {
+                    (*self).tanh()
+                }
+
+                fn asin(&self) -> Self {
+                    (*self).asin()
+                }
+
+                fn acos(&self) -> Self {
+                    (*self).acos()
+                }
+
+                fn atan(&self) -> Self {
+                    (*self).atan()
+                }
+
+                fn asinh(&self) -> Self {
+                    (*self).asinh()
+                }
+
+                fn acosh(&self) -> Self {
+                    (*self).acosh()
+                }
+
+                fn atanh(&self) -> Self {
+                    (*self).atanh()
+                }
+            }
+
+            impl ExpLogOps for $t {
+                type Float = $t;
+                fn exp(&self) -> Self {
+                    (*self).exp()
+                }
+                fn ln(&self) -> Self {
+                    (*self).ln()
+                }
+                fn log(&self, base: Self::Float) -> Self {
+                    (*self).log(base)
+                }
+                fn log2(&self) -> Self {
+                    (*self).log2()
+                }
+                fn log10(&self) -> Self {
+                    (*self).log10()
+                }
+            }
+
+            impl Float for $t {}
+        )*
+    };
+}
+
+impl_float!(f32, f64);
+
+impl Numeric<f32> for f32 {}
+impl Numeric<f64> for f64 {}
+
+// ┌──────────────────────────────────────────────────────────┐
+//  Fundamental Traits
+// └──────────────────────────────────────────────────────────┘
+pub trait Group: Sized + Add<Self, Output = Self> {
+    fn zero() -> Self;
+}
+
+pub trait Ring: Group + Mul<Self, Output = Self> {
+    fn one() -> Self;
+}
+
+impl Group for f32 {
+    fn zero() -> Self {
+        0.0
+    }
+}
+
+impl Ring for f32 {
+    fn one() -> Self {
+        1.0
+    }
+}
+
+impl Group for f64 {
+    fn zero() -> Self {
+        0.0
+    }
+}
+
+impl Ring for f64 {
+    fn one() -> Self {
+        1.0
+    }
+}
diff --git a/src/fuga/mod.rs b/src/fuga/mod.rs
index c8a497ba..f91274b6 100644
--- a/src/fuga/mod.rs
+++ b/src/fuga/mod.rs
@@ -151,12 +151,14 @@ pub use crate::macros::{julia_macro::*, matlab_macro::*, r_macro::*};
 
 pub use peroxide_ad::{ad_function, ad_closure};
 
+pub use peroxide_num::{ExpLogOps, PowOps, TrigOps};
+
 pub use crate::traits::{
     fp::{FPMatrix, FPVector},
     general::Algorithm,
     math::{InnerProduct, LinearOp, MatrixProduct, Norm, Normed, Vector, VectorProduct},
     mutable::{MutFP, MutMatrix},
-    num::{ExpLogOps, PowOps, Real, TrigOps},
+    num::Real,
     pointer::{MatrixPtr, Oxide, Redox, RedoxCommon},
     stable::StableFn,
     sugar::{Scalable, ScalableMut, VecOps, ConvToMat},
diff --git a/src/numerical/spline.rs b/src/numerical/spline.rs
index d504aa5f..6b70900b 100644
--- a/src/numerical/spline.rs
+++ b/src/numerical/spline.rs
@@ -198,7 +198,7 @@ use crate::structure::matrix::*;
 use crate::structure::polynomial::*;
 #[allow(unused_imports)]
 use crate::structure::vector::*;
-use crate::traits::num::PowOps;
+use peroxide_num::PowOps;
 #[allow(unused_imports)]
 use crate::util::non_macro::*;
 use crate::util::useful::zip_range;
diff --git a/src/prelude/mod.rs b/src/prelude/mod.rs
index 4905251f..ea576bb7 100644
--- a/src/prelude/mod.rs
+++ b/src/prelude/mod.rs
@@ -154,11 +154,13 @@ pub use crate::traits::{
     general::Algorithm,
     math::{InnerProduct, LinearOp, MatrixProduct, Vector, VectorProduct},
     mutable::{MutFP, MutMatrix},
-    num::{ExpLogOps, PowOps, Real, TrigOps},
+    num::Real,
     pointer::{MatrixPtr, Oxide, Redox, RedoxCommon},
     sugar::{Scalable, ScalableMut, VecOps, ConvToMat},
 };
 
+pub use peroxide_num::{ExpLogOps, TrigOps, PowOps};
+
 pub use simpler::SimpleNorm;
 
 #[allow(unused_imports)]
diff --git a/src/structure/ad.rs b/src/structure/ad.rs
index 15f2c85e..55d0b007 100644
--- a/src/structure/ad.rs
+++ b/src/structure/ad.rs
@@ -17,7 +17,7 @@
 //! * `FromIterator<&f64>`
 //! * `Index`, `IndexMut`
 //! * `std::ops::{Neg, Add, Sub, Mul, Div}`
-//! * `peroxide::traits::num::{PowOps, ExpLogOps, TrigOps}`
+//! * `peroxide_num::{PowOps, ExpLogOps, TrigOps}`
 //! * `peroxide::util::print::Printable`
 //!
 //! ## Iterator of `AD`
@@ -105,9 +105,9 @@
 //! ```
 //!
 
+use peroxide_num::{ExpLogOps, PowOps, TrigOps};
 use crate::statistics::ops::C;
 use crate::traits::{
-    num::{ExpLogOps, PowOps, TrigOps},
     stable::StableFn,
     fp::FPVector,
     math::Vector,
@@ -173,6 +173,10 @@ impl AD {
         }
     }
 
+    pub fn is_empty(&self) -> bool {
+        self.len() == 0
+    }
+
     pub fn iter(&self) -> ADIter {
         self.into_iter()
     }
@@ -309,9 +313,9 @@ impl AD {
     #[allow(dead_code)]
     unsafe fn x_ptr(&self) -> Option<*const f64> {
         match self {
-            AD0(x) => Some(&*x),
-            AD1(x,_) => Some(&*x),
-            AD2(x,_,_) => Some(&*x),
+            AD0(x) => Some(x),
+            AD1(x,_) => Some(x),
+            AD2(x,_,_) => Some(x),
         }
     }
 
@@ -319,8 +323,8 @@ impl AD {
     unsafe fn dx_ptr(&self) -> Option<*const f64> {
         match self {
             AD0(_) => None,
-            AD1(_,dx) => Some(&*dx),
-            AD2(_,dx,_) => Some(&*dx),
+            AD1(_,dx) => Some(dx),
+            AD2(_,dx,_) => Some(dx),
         }
     }
 
@@ -329,7 +333,7 @@ impl AD {
         match self {
             AD0(_) => None,
             AD1(_,_) => None,
-            AD2(_,_,ddx) => Some(&*ddx),
+            AD2(_,_,ddx) => Some(ddx),
         }
     }
 
@@ -619,7 +623,7 @@ impl Add<AD> for AD {
         let (a, b) = (self.to_order(ord), rhs.to_order(ord));
 
         a.into_iter()
-            .zip(b.into_iter())
+            .zip(b)
             .map(|(x, y)| x + y)
             .collect()
     }
@@ -633,7 +637,7 @@ impl Sub<AD> for AD {
         let (a, b) = (self.to_order(ord), rhs.to_order(ord));
 
         a.into_iter()
-            .zip(b.into_iter())
+            .zip(b)
             .map(|(x, y)| x - y)
             .collect()
     }
@@ -646,7 +650,7 @@ impl Mul<AD> for AD {
         let ord = self.order().max(rhs.order());
         let (a, b) = (self.to_order(ord), rhs.to_order(ord));
 
-        let mut z = a.clone();
+        let mut z = a;
         for t in 0..z.len() {
             z[t] = a
                 .into_iter()
@@ -666,7 +670,7 @@ impl Div<AD> for AD {
         let ord = self.order().max(rhs.order());
         let (a, b) = (self.to_order(ord), rhs.to_order(ord));
 
-        let mut z = a.clone();
+        let mut z = a;
         z[0] = a[0] / b[0];
         let y0 = 1f64 / b[0];
         for i in 1 .. z.len() {
@@ -681,6 +685,8 @@ impl Div<AD> for AD {
 }
 
 impl ExpLogOps for AD {
+    type Float = f64;
+
     fn exp(&self) -> Self {
         let mut z = self.empty();
         z[0] = self[0].exp();
@@ -711,11 +717,21 @@ impl ExpLogOps for AD {
     fn log(&self, base: f64) -> Self {
         self.ln().iter().map(|x| x / base.ln()).collect()
     }
+
+    fn log2(&self) -> Self {
+        self.log(2f64)
+    }
+
+    fn log10(&self) -> Self {
+        self.log(10f64)
+    }
 }
 
 impl PowOps for AD {
+    type Float = f64;
+
     fn powi(&self, n: i32) -> Self {
-        let mut z = self.clone();
+        let mut z = *self;
         for _i in 1 .. n {
             z = z * *self;
         }
@@ -750,6 +766,10 @@ impl PowOps for AD {
         }
         z
     }
+
+    fn sqrt(&self) -> Self {
+        self.powf(0.5f64)
+    }
 }
 
 impl TrigOps for AD {
@@ -773,7 +793,7 @@ impl TrigOps for AD {
         (u, v)
     }
 
-    fn sinh_cosh(&self) -> (Self, Self) {
+    fn sinh(&self) -> Self {
         let mut u = self.empty();
         let mut v = self.empty();
         u[0] = self[0].sinh();
@@ -790,7 +810,47 @@ impl TrigOps for AD {
                 .enumerate()
                 .fold(0f64, |x, (k, (&u1, &x1))| x + (C(i-1, k) as f64) * x1 * u1);
         }
-        (u, v)
+        u
+    }
+
+    fn cosh(&self) -> Self {
+        let mut u = self.empty();
+        let mut v = self.empty();
+        u[0] = self[0].sinh();
+        v[0] = self[0].cosh();
+        for i in 1 .. u.len() {
+            u[i] = v.iter()
+                .take(i)
+                .zip(self.iter().skip(1).take(i).rev())
+                .enumerate()
+                .fold(0f64, |x, (k, (&v1, &x1))| x + (C(i-1, k) as f64) * x1 * v1);
+            v[i] = u.iter()
+                .take(i)
+                .zip(self.iter().skip(1).take(i).rev())
+                .enumerate()
+                .fold(0f64, |x, (k, (&u1, &x1))| x + (C(i-1, k) as f64) * x1 * u1);
+        }
+        v
+    }
+
+    fn tanh(&self) -> Self {
+        let mut u = self.empty();
+        let mut v = self.empty();
+        u[0] = self[0].sinh();
+        v[0] = self[0].cosh();
+        for i in 1 .. u.len() {
+            u[i] = v.iter()
+                .take(i)
+                .zip(self.iter().skip(1).take(i).rev())
+                .enumerate()
+                .fold(0f64, |x, (k, (&v1, &x1))| x + (C(i-1, k) as f64) * x1 * v1);
+            v[i] = u.iter()
+                .take(i)
+                .zip(self.iter().skip(1).take(i).rev())
+                .enumerate()
+                .fold(0f64, |x, (k, (&u1, &x1))| x + (C(i-1, k) as f64) * x1 * u1);
+        }
+        u / v
     }
 
     fn asin(&self) -> Self {
@@ -1216,13 +1276,13 @@ impl ADVec for Vec<f64> {
     }
 
     fn to_f64_vec(&self) -> Vec<f64> {
-        self.iter().map(|&t| t).collect()
+        self.clone()
     }
 }
 
 impl ADVec for Vec<AD> {
     fn to_ad_vec(&self) -> Vec<AD> {
-        self.iter().map(|&t| t).collect()
+        self.clone()
     }
 
     fn to_f64_vec(&self) -> Vec<f64> {
@@ -1249,7 +1309,7 @@ impl FPVector for Vec<AD> {
     fn filter<F>(&self, f: F) -> Self
     where
             F: Fn(Self::Scalar) -> bool {
-        self.iter().filter(|&x| f(*x)).map(|&t| t).collect()
+        self.iter().filter(|&x| f(*x)).cloned().collect()
     }
 
     fn zip_with<F>(&self, f: F, other: &Self) -> Self
@@ -1259,11 +1319,11 @@ impl FPVector for Vec<AD> {
     }
 
     fn take(&self, n: usize) -> Self {
-        self.iter().take(n).map(|&t| t).collect()
+        self.iter().take(n).cloned().collect()
     }
 
     fn skip(&self, n: usize) -> Self {
-        self.iter().skip(n).map(|&t| t).collect()
+        self.iter().skip(n).cloned().collect()
     }
 
     fn sum(&self) -> Self::Scalar {
@@ -1280,11 +1340,11 @@ impl FPVector for Vec<AD> {
 impl Vector for Vec<AD> {
     type Scalar = AD;
 
-    fn add_vec<'a, 'b>(&'a self, rhs: &'b Self) -> Self {
+    fn add_vec(&self, rhs: &Self) -> Self {
         self.add_v(rhs)
     }
 
-    fn sub_vec<'a, 'b>(&'a self, rhs: &'b Self) -> Self {
+    fn sub_vec(&self, rhs: &Self) -> Self {
         self.sub_v(rhs)
     }
 
diff --git a/src/structure/dataframe.rs b/src/structure/dataframe.rs
index a11bfd6b..2d8feec8 100644
--- a/src/structure/dataframe.rs
+++ b/src/structure/dataframe.rs
@@ -914,12 +914,12 @@ fn vtype_to_dtype(dv: netcdf::types::BasicType) -> DType {
 }
 
 #[cfg(feature= "nc")]
-fn nc_put_value<'f, T: Numeric>(var: &mut VariableMut<'f>, v: Vec<T>) -> Result<(), netcdf::error::Error> {
+fn nc_put_value<T: Numeric>(var: &mut VariableMut, v: Vec<T>) -> Result<(), netcdf::error::Error> {
     var.put_values(&v, None, None)
 }
 
 #[cfg(feature= "nc")]
-fn nc_read_value<'f, T: Numeric + Default + Clone>(val: &Variable<'f>, v: Vec<T>) -> Result<Series, netcdf::error::Error> where Series: TypedVector<T> {
+fn nc_read_value<T: Numeric + Default + Clone>(val: &Variable, v: Vec<T>) -> Result<Series, netcdf::error::Error> where Series: TypedVector<T> {
     let mut v = v;
     v.resize_with(val.len(), Default::default);
     val.values_to(&mut v, None, None)?;
@@ -1011,7 +1011,7 @@ macro_rules! dtype_match_to_arrow {
 }
 
 #[cfg(feature= "parquet")]
-fn parquet_read_value<'f, T: Default + Clone + NativeType>(arr: &Box<dyn Array>, v: Vec<T>) -> Result<Series, arrow2::error::Error> where Series: TypedVector<T> {
+fn parquet_read_value<T: Default + Clone + NativeType>(arr: &Box<dyn Array>, _v: Vec<T>) -> Result<Series, arrow2::error::Error> where Series: TypedVector<T> {
     let x = arr.as_any().downcast_ref::<PrimitiveArray<T>>().unwrap();
     let x = x.values_iter().cloned().collect::<Vec<_>>();
 
@@ -1020,12 +1020,12 @@ fn parquet_read_value<'f, T: Default + Clone + NativeType>(arr: &Box<dyn Array>,
 
 fn add_vec<T: std::ops::Add<T, Output=T> + Clone>(v: Vec<T>, w: Vec<T>) -> Series 
 where Series: TypedVector<T> {
-    Series::new(v.into_iter().zip(w.into_iter()).map(|(x, y)| x + y).collect::<Vec<T>>())
+    Series::new(v.into_iter().zip(w).map(|(x, y)| x + y).collect::<Vec<T>>())
 }
 
 fn sub_vec<T: std::ops::Sub<T, Output=T> + Clone>(v: Vec<T>, w: Vec<T>) -> Series 
 where Series: TypedVector<T> {
-    Series::new(v.into_iter().zip(w.into_iter()).map(|(x, y)| x - y).collect::<Vec<T>>())
+    Series::new(v.into_iter().zip(w).map(|(x, y)| x - y).collect::<Vec<T>>())
 }
 
 fn mul_scalar<T: std::ops::Mul<T, Output=T> + Clone + Copy>(v: Vec<T>, s: T) -> Series
diff --git a/src/structure/matrix.rs b/src/structure/matrix.rs
index d5919b75..fe37ed64 100644
--- a/src/structure/matrix.rs
+++ b/src/structure/matrix.rs
@@ -77,7 +77,7 @@
 //!
 //!     ```R
 //!     # R
-//!     a = matrix(1:4:1, 2, 2, Row)
+//!     a = matrix(1:4, nrow = 2, ncol = 2, byrow = T)
 //!     print(a)
 //!     #      [,1] [,2]
 //!     # [1,]    1    2
@@ -3655,8 +3655,7 @@ impl MutMatrix for Matrix {
         assert!(idx < self.col, "Index out of range");
         match self.shape {
             Shape::Col => {
-                let mut v: Vec<*mut f64> = Vec::with_capacity(self.row);
-                v.set_len(self.row);
+                let mut v: Vec<*mut f64> = vec![&mut 0f64; self.row];
                 let start_idx = idx * self.row;
                 let p = self.mut_ptr();
                 for (i, j) in (start_idx..start_idx + v.len()).enumerate() {
@@ -3665,8 +3664,7 @@ impl MutMatrix for Matrix {
                 v
             }
             Shape::Row => {
-                let mut v: Vec<*mut f64> = Vec::with_capacity(self.row);
-                v.set_len(self.row);
+                let mut v: Vec<*mut f64> = vec![&mut 0f64; self.row];
                 let p = self.mut_ptr();
                 for i in 0..v.len() {
                     v[i] = p.add(idx + i * self.col);
@@ -3680,8 +3678,7 @@ impl MutMatrix for Matrix {
         assert!(idx < self.row, "Index out of range");
         match self.shape {
             Shape::Row => {
-                let mut v: Vec<*mut f64> = Vec::with_capacity(self.col);
-                v.set_len(self.col);
+                let mut v: Vec<*mut f64> = vec![&mut 0f64; self.col];
                 let start_idx = idx * self.col;
                 let p = self.mut_ptr();
                 for (i, j) in (start_idx..start_idx + v.len()).enumerate() {
@@ -3690,8 +3687,7 @@ impl MutMatrix for Matrix {
                 v
             }
             Shape::Col => {
-                let mut v: Vec<*mut f64> = Vec::with_capacity(self.col);
-                v.set_len(self.col);
+                let mut v: Vec<*mut f64> = vec![&mut 0f64; self.col];
                 let p = self.mut_ptr();
                 for i in 0..v.len() {
                     v[i] = p.add(idx + i * self.row);
diff --git a/src/structure/polynomial.rs b/src/structure/polynomial.rs
index 1a01fb53..0487f181 100644
--- a/src/structure/polynomial.rs
+++ b/src/structure/polynomial.rs
@@ -6,7 +6,8 @@ use crate::util::useful::*;
 #[cfg(feature = "serde")]
 use serde::{Deserialize, Serialize};
 
-use crate::traits::{fp::FPVector, num::PowOps};
+use peroxide_num::PowOps;
+use crate::traits::fp::FPVector;
 use std::cmp::{max, min};
 use std::convert;
 use std::fmt;
@@ -461,6 +462,7 @@ impl Mul<Polynomial> for f64 {
 // Extra operations for Polynomial
 // =============================================================================
 impl PowOps for Polynomial {
+    type Float = f64;
     fn powi(&self, n: i32) -> Self {
         let mut result = self.clone();
         for _i in 0..n - 1 {
diff --git a/src/traits/num.rs b/src/traits/num.rs
index a281812a..77d3fa55 100644
--- a/src/traits/num.rs
+++ b/src/traits/num.rs
@@ -26,68 +26,7 @@
 
 use std::ops::{Neg, Add, Sub, Mul, Div};
 use crate::structure::ad::AD;
-
-pub trait PowOps: Sized {
-    fn powi(&self, n: i32) -> Self;
-    fn powf(&self, f: f64) -> Self;
-    fn pow(&self, f: Self) -> Self;
-    fn sqrt(&self) -> Self {
-        self.powf(0.5)
-    }
-}
-
-pub trait TrigOps: Sized + Div<Output = Self> {
-    fn sin_cos(&self) -> (Self, Self);
-    fn sin(&self) -> Self {
-        let (s, _) = self.sin_cos();
-        s
-    }
-    fn cos(&self) -> Self {
-        let (_, c) = self.sin_cos();
-        c
-    }
-    fn tan(&self) -> Self {
-        let (s, c) = self.sin_cos();
-        s / c
-    }
-    fn sinh_cosh(&self) -> (Self, Self);
-    fn sinh(&self) -> Self {
-        let (s, _) = self.sinh_cosh();
-        s
-    }
-    fn cosh(&self) -> Self {
-        let (_, c) = self.sinh_cosh();
-        c
-    }
-    fn tanh(&self) -> Self {
-        let (s, c) = self.sinh_cosh();
-        s / c
-    }
-    fn asin(&self) -> Self;
-    fn acos(&self) -> Self;
-    fn atan(&self) -> Self;
-    fn asin_acos(&self) -> (Self, Self) {
-        (self.asin(), self.acos())
-    }
-    fn asinh(&self) -> Self;
-    fn acosh(&self) -> Self;
-    fn atanh(&self) -> Self;
-    fn asinh_acosh(&self) -> (Self, Self) {
-        (self.asinh(), self.acosh())
-    }
-}
-
-pub trait ExpLogOps: Sized {
-    fn exp(&self) -> Self;
-    fn ln(&self) -> Self;
-    fn log(&self, base: f64) -> Self;
-    fn log2(&self) -> Self {
-        self.log(2f64)
-    }
-    fn log10(&self) -> Self {
-        self.log(10f64)
-    }
-}
+use peroxide_num::{PowOps, TrigOps, ExpLogOps};
 
 pub trait Real:
     PowOps
@@ -111,107 +50,6 @@ pub trait Real:
     fn to_ad(&self) -> AD;
 }
 
-// =============================================================================
-// Real Traits for f64
-// =============================================================================
-impl PowOps for f64 {
-    fn powi(&self, n: i32) -> Self {
-        (*self).powi(n)
-    }
-
-    fn powf(&self, f: f64) -> Self {
-        (*self).powf(f)
-    }
-
-    fn pow(&self, f: f64) -> Self {
-        (*self).powf(f)
-    }
-
-    fn sqrt(&self) -> Self {
-        (*self).sqrt()
-    }
-}
-
-impl TrigOps for f64 {
-    fn sin(&self) -> Self {
-        (*self).sin()
-    }
-
-    fn cos(&self) -> Self {
-        (*self).cos()
-    }
-
-    fn tan(&self) -> Self {
-        (*self).tan()
-    }
-
-    fn asin(&self) -> Self {
-        (*self).asin()
-    }
-
-    fn acos(&self) -> Self {
-        (*self).acos()
-    }
-
-    fn atan(&self) -> Self {
-        (*self).atan()
-    }
-
-    fn sinh(&self) -> Self {
-        (*self).sinh()
-    }
-
-    fn cosh(&self) -> Self {
-        (*self).cosh()
-    }
-
-    fn sinh_cosh(&self) -> (Self, Self) {
-        ((*self).sinh(), (*self).cosh())
-    }
-
-    fn tanh(&self) -> Self {
-        (*self).tanh()
-    }
-
-    fn asinh(&self) -> Self {
-        (*self).asinh()
-    }
-
-    fn acosh(&self) -> Self {
-        (*self).acosh()
-    }
-
-    fn atanh(&self) -> Self {
-        (*self).atanh()
-    }
-
-    fn sin_cos(&self) -> (Self, Self) {
-        (*self).sin_cos()
-    }
-}
-
-impl ExpLogOps for f64 {
-    fn exp(&self) -> Self {
-        (*self).exp()
-    }
-
-    fn ln(&self) -> Self {
-        (*self).ln()
-    }
-
-    fn log(&self, base: f64) -> Self {
-        (*self).log(base)
-    }
-
-    fn log2(&self) -> Self {
-        (*self).log2()
-    }
-
-    fn log10(&self) -> Self {
-        (*self).log10()
-    }
-}
-
 impl Real for f64 {
     fn to_f64(&self) -> f64 {
         *self
diff --git a/src/traits/pointer.rs b/src/traits/pointer.rs
index 03eb3ef4..65fab349 100644
--- a/src/traits/pointer.rs
+++ b/src/traits/pointer.rs
@@ -284,8 +284,7 @@ impl MatrixPtr for Matrix {
         assert!(idx < self.col, "Index out of range");
         match self.shape {
             Shape::Row => {
-                let mut v: Vec<*const f64> = Vec::with_capacity(self.col);
-                v.set_len(self.col);
+                let mut v: Vec<*const f64> = vec![&0f64; self.col];
                 let start_idx = idx * self.col;
                 let p = self.ptr();
                 for (i, j) in (start_idx..start_idx + v.len()).enumerate() {
@@ -294,8 +293,7 @@ impl MatrixPtr for Matrix {
                 v
             }
             Shape::Col => {
-                let mut v: Vec<*const f64> = Vec::with_capacity(self.col);
-                v.set_len(self.col);
+                let mut v: Vec<*const f64> = vec![&0f64; self.col];
                 let p = self.ptr();
                 for (i, elem) in v.iter_mut().enumerate() {
                     *elem = p.add(idx + i * self.row);
@@ -309,8 +307,7 @@ impl MatrixPtr for Matrix {
         assert!(idx < self.col, "Index out of range");
         match self.shape {
             Shape::Col => {
-                let mut v: Vec<*const f64> = Vec::with_capacity(self.row);
-                v.set_len(self.row);
+                let mut v: Vec<*const f64> = vec![&0f64; self.row];
                 let start_idx = idx * self.row;
                 let p = self.ptr();
                 for (i, j) in (start_idx..start_idx + v.len()).enumerate() {
@@ -319,8 +316,7 @@ impl MatrixPtr for Matrix {
                 v
             }
             Shape::Row => {
-                let mut v: Vec<*const f64> = Vec::with_capacity(self.row);
-                v.set_len(self.row);
+                let mut v: Vec<*const f64> = vec![&0f64; self.row];
                 let p = self.ptr();
                 for (i, elem) in v.iter_mut().enumerate() {
                     *elem = p.add(idx + i * self.col);
diff --git a/src/util/useful.rs b/src/util/useful.rs
index faa80cf3..c288debf 100644
--- a/src/util/useful.rs
+++ b/src/util/useful.rs
@@ -26,9 +26,7 @@ where
     let e = 1e-7;
     let p: f64 = x.into().abs();
     let q: f64 = y.into().abs();
-    if (p - q).abs() < e {
-        b = true;
-    } else if (p - q).abs() / (p + q).min(MAX) < e {
+    if (p - q).abs() < e || (p - q).abs() / (p + q).min(MAX) < e {
         b = true;
     }
     b
@@ -40,7 +38,7 @@ pub fn tab(s: &str, space: usize) -> String {
     let mut m: String = String::new();
     let fs = format!("{}{}", " ".repeat(space - l), s);
     m.push_str(&fs);
-    return m;
+    m
 }
 
 pub fn quot_rem(x: usize, y: usize) -> (i32, i32) {
@@ -120,7 +118,7 @@ pub fn sgn(x: usize) -> f64 {
 }
 
 /// Vector compare
-pub fn eq_vec(x: &Vec<f64>, y: &Vec<f64>, tol: f64) -> bool {
+pub fn eq_vec(x: &[f64], y: &[f64], tol: f64) -> bool {
     x.iter().zip(y.iter()).all(|(x, y)| (x - y).abs() <= tol)
 }
 
@@ -141,7 +139,7 @@ pub fn eq_vec(x: &Vec<f64>, y: &Vec<f64>, tol: f64) -> bool {
 pub fn auto_zip<T: Clone>(x: &Vec<T>) -> Vec<(T, T)> {
     let x_head = x[0 .. x.len() - 1].to_vec();
     let x_tail = x[1 .. x.len()].to_vec();
-    x_head.into_iter().zip(x_tail.into_iter()).collect()
+    x_head.into_iter().zip(x_tail).collect()
 }
 
 /// Find the index of interval of x
@@ -180,7 +178,7 @@ pub fn find_interval<T: PartialOrd + PartialEq>(sorted_intervals: &Vec<(T, T)>,
             return mid;
         }
     }
-    return i;
+    i
 }
 
 /// Generate Range of Intervals
@@ -234,11 +232,11 @@ pub fn gen_range<T: Clone + PartialOrd>(x: &[T]) -> Vec<Range<T>> {
 /// assert_eq!(r, answer);
 /// ```
 pub fn zip_range<T: Clone + PartialOrd, U: Clone>(x: &[T], y: &[U]) -> Vec<(Range<T>, U)> {
-    y[0 .. x.len() - 1].into_iter()
+    y[0 .. x.len() - 1].iter()
         .enumerate()
         .map(|(i, yi)| (Range {
             start: x[i].clone(),
             end: x[i + 1].clone(),
         }, yi.clone()))
         .collect()
-}
\ No newline at end of file
+}