From 2de1c9c9b481d3b6b63fe42026f80045b93432d3 Mon Sep 17 00:00:00 2001
From: Moritz Moeller <virtualritz@protonmail.com>
Date: Sun, 7 May 2023 15:26:15 +0200
Subject: [PATCH] Added `smootherstep()`.

---
 src/lib.rs | 149 +++++++++++++++++++++++++++--------------------------
 src/num.rs |  53 +++++++++++--------
 src/vec.rs |  74 +++++++++++++++-----------
 3 files changed, 152 insertions(+), 124 deletions(-)

diff --git a/src/lib.rs b/src/lib.rs
index c54a5c1..2864b54 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -92,7 +92,7 @@ pub fn copysign<T: Float, V: FloatOps<T>>(a: V, sign: T) -> V {
     V::copysign(a, sign)
 }
 
-/// returns the absolute (positive) value of `a` 
+/// returns the absolute (positive) value of `a`
 pub fn abs<T: SignedNumber, V: SignedNumberOps<T>>(a: V) -> V {
     V::abs(a)
 }
@@ -202,22 +202,27 @@ pub fn slerp<T: Float + NumberOps<T> + FloatOps<T>, V: Slerp<T>>(e0: V, e1: V, t
     V::slerp(e0, e1, t)
 }
 
-/// returns the hermite interpolated value between edge `e0` and `e1` by percentage `t`
+/// returns the cubic hermite interpolated value between edge `e0` and `e1` by percentage `t`
 pub fn smoothstep<T: Float, V: FloatOps<T>>(e0: V, e1: V, t: T) -> V {
     V::smoothstep(e0, e1, t)
 }
 
+/// returns the quintic hermite interpolated value between edge `e0` and `e1` by percentage `t`
+pub fn smootherstep<T: Float, V: FloatOps<T>>(e0: V, e1: V, t: T) -> V {
+    V::smootherstep(e0, e1, t)
+}
+
 /// returns the saturated value of `x` into to the 0-1 range, this is the same as `clamp(x, 0.0, 1.0)`
 pub fn saturate<T: Float, V: FloatOps<T>>(x: V) -> V {
     V::saturate(x)
 }
 
-/// returns the vector cross product of `a x b`, makes sense only for 3 or 7 dimensional vectors 
+/// returns the vector cross product of `a x b`, makes sense only for 3 or 7 dimensional vectors
 pub fn cross<T: Number, V: Cross<T>>(a: V, b: V) -> V {
     V::cross(a, b)
 }
 
-/// returns the scalar triple product of `a x b x c`, makes sense only for 3 dimensional vectors 
+/// returns the scalar triple product of `a x b x c`, makes sense only for 3 dimensional vectors
 pub fn scalar_triple<T: Number + SignedNumber, V: Triple<T>>(a: V, b: V, c: V) -> T {
     V::scalar_triple(a, b, c)
 }
@@ -230,7 +235,7 @@ pub fn vector_triple<T: Number + SignedNumber, V: Triple<T>>(a: V, b: V, c: V) -
 /// returns the perpedicular vector of `a` performing anti-clockwise rotation by 90 degrees
 pub fn perp<T: SignedNumber>(a: Vec2<T>) -> Vec2<T> {
     Vec2 {
-        x: -a.y, 
+        x: -a.y,
         y: a.x
     }
 }
@@ -401,15 +406,15 @@ pub fn log<T: Float, V: FloatOps<T>>(v: V, base: T) -> V {
     V::log(v, base)
 }
 
-/// returns the vec2 rotated anti-clockwise rotation by radian `angle` 
+/// returns the vec2 rotated anti-clockwise rotation by radian `angle`
 pub fn rotate_2d<T: Float + FloatOps<T>>(v: Vec2<T>, angle: T) -> Vec2<T> {
     let c = cos(angle);
     let s = sin(angle);
     Vec2::new(c * v.x - s * v.y, s * v.x + c * v.y)
 }
-   
+
 /// returns a convex hull wound clockwise from point cloud `points`
-pub fn convex_hull_from_points<T: Float + SignedNumberOps<T> + NumberOps<T> + FloatOps<T>>(points: &Vec<Vec2<T>>) -> Vec<Vec2<T>> {    
+pub fn convex_hull_from_points<T: Float + SignedNumberOps<T> + NumberOps<T> + FloatOps<T>>(points: &Vec<Vec2<T>>) -> Vec<Vec2<T>> {
     //find right most point
     let mut cur = points[0];
     let mut curi = 0;
@@ -420,7 +425,7 @@ pub fn convex_hull_from_points<T: Float + SignedNumberOps<T> + NumberOps<T> + Fl
             curi = i;
         }
     }
-    
+
     // wind the hull clockwise by using cross product to test which side of an edge a point lies on
     // discarding points that do not form the perimeter
     let mut hull = vec![cur];
@@ -445,7 +450,7 @@ pub fn convex_hull_from_points<T: Float + SignedNumberOps<T> + NumberOps<T> + Fl
         if approx(x1, hull[0], T::small_epsilon()) {
             break;
         }
-            
+
         cur = x1;
         curi = rm;
         hull.push(x1);
@@ -702,11 +707,11 @@ pub fn closest_point_on_cone<T: Float, V: VecN<T> + VecFloatOps<T>>(p: V, cp: V,
     // find point onbase plane and clamp to the extent
     let cplane = closest_point_on_plane(p, l2, cv);
     let extent = l2 + normalize(cplane - l2) * r;
-    
+
     // test closest point on line with the axis along the side and bottom of the cone
     let e1 = closest_point_on_line_segment(p, cp, extent);
     let e2 = closest_point_on_line_segment(p, l2, extent);
-    
+
     if dist2(p, e1) < dist2(p, e2) {
         e1
     }
@@ -729,7 +734,7 @@ pub fn closest_point_on_polygon<T: Float + FloatOps<T>>(p: Vec2<T>, poly: &Vec<V
         let i2 = (i+1)%ncp;
         let cpp = closest_point_on_line_segment(p, poly[i], poly[i2]);
         let cppd2 = dist2(p, cpp);
-        if dist2(p, cpp) < cd2 { 
+        if dist2(p, cpp) < cd2 {
             cp = cpp;
             cd2 = cppd2;
         }
@@ -818,7 +823,7 @@ pub fn point_inside_frustum<T: Number>(p: Vec3<T>, planes: &[Vec4<T>; 6]) -> boo
     }
     true
 }
-    
+
 /// returns the classification of point `p` vs the plane defined by point on plane `x` and normal `n`
 pub fn point_vs_plane<T: SignedNumber + SignedNumberOps<T>>(p: Vec3<T>,  x: Vec3<T>, n: Vec3<T>) -> Classification {
     let pd = plane_distance(x, n);
@@ -903,7 +908,7 @@ pub fn capsule_vs_plane<T: Float + FloatOps<T> + SignedNumber + SignedNumberOps<
     }
 }
 
-/// return the classification of cone defined by position `cp`, direction `cv` with height `h` and radius at the base of `r` 
+/// return the classification of cone defined by position `cp`, direction `cv` with height `h` and radius at the base of `r`
 /// vs the plane defined by point `x` and normal `n`
 pub fn cone_vs_plane<T: Float + SignedNumberOps<T>, V: VecN<T> + Cross<T> + Dot<T> + SignedVecN<T> + VecFloatOps<T>>(cp: V, cv: V, h: T, r: T, x: V, n: V) -> Classification {
     let l2 = cp + cv * h;
@@ -989,7 +994,7 @@ pub fn aabb_vs_aabb<T: Number, V: VecN<T> + NumberOps<T>>(aabb_min1: V, aabb_max
     for i in 0..V::len() {
         if aabb_max1[i] < aabb_min2[i] || aabb_min1[i] > aabb_max2[i] {
             return false;
-        } 
+        }
     }
     true
 }
@@ -1007,7 +1012,7 @@ pub fn aabb_vs_obb<T: Number + Float + SignedNumber + SignedNumberOps<T> + Numbe
         Vec3::<T>::new( T::one(),  T::one(),  T::one()),
         Vec3::<T>::new(-T::one(),  T::one(),  T::one()),
     ];
-    
+
     // aabb
     let verts0 = vec![
         aabb_min,
@@ -1019,13 +1024,13 @@ pub fn aabb_vs_obb<T: Number + Float + SignedNumber + SignedNumberOps<T> + Numbe
         Vec3::<T>::new(aabb_max.x, aabb_max.y, aabb_min.z),
         aabb_max
     ];
-    
+
     // obb from corners
     let mut verts1 = Vec::new();
     for corner in corners {
         verts1.push(obb * corner);
     }
-    
+
     gjk_3d(verts0, verts1)
 }
 
@@ -1050,7 +1055,7 @@ pub fn obb_vs_obb<T: Number + Float + SignedNumber + SignedNumberOps<T> + Number
         Vec3::<T>::new( T::one(),  T::one(),  T::one()),
         Vec3::<T>::new(-T::one(),  T::one(),  T::one()),
     ];
-    
+
     // obb from corners
     let mut verts0 = Vec::new();
     for corner in corners {
@@ -1061,7 +1066,7 @@ pub fn obb_vs_obb<T: Number + Float + SignedNumber + SignedNumberOps<T> + Number
     for corner in corners {
         verts1.push(obb1 * corner);
     }
-    
+
     gjk_3d(verts0, verts1)
 }
 
@@ -1087,7 +1092,7 @@ pub fn capsule_vs_capsule<T: Float + FloatOps<T> + SignedNumberOps<T>, V: VecN<T
         if d < r2 {
             let l1 = normalize(cp3 - cp2);
             let t1 = dot(cp0 - cp2, l1);
-            
+
             // now check if the capsule axes overlap
             if t0 >= T::zero() && t0*t0 < dist2(cp1, cp0) {
                 true
@@ -1157,8 +1162,8 @@ pub fn ray_vs_aabb<T: Number + NumberOps<T>, V: VecN<T>>(r0: V, rv: V, aabb_min:
 }
 
 /// returns the intersection of the 3D ray with origin `r0` and direction `rv` with the obb defined by `mat`
-pub fn ray_vs_obb<T: Float + NumberOps<T>, 
-    V: VecFloatOps<T> + NumberOps<T> + SignedNumberOps<T> + VecN<T> + SignedVecN<T>, 
+pub fn ray_vs_obb<T: Float + NumberOps<T>,
+    V: VecFloatOps<T> + NumberOps<T> + SignedNumberOps<T> + VecN<T> + SignedVecN<T>,
     M: MatTranslate<V> + MatInverse<T> + MatRotate3D<T, V> + MatN<T, V>
     + Into<Mat3<T>> + Copy>
     (r0: V, rv: V, mat: M) -> Option<V> where Mat3<T> : MatN<T, V> {
@@ -1225,18 +1230,18 @@ pub fn ray_vs_capsule<T: Float + FloatOps<T> + NumberOps<T> + SignedNumberOps<T>
             let b = c1;
             let v = rv;
             let r = cr;
-            
+
             let ab = b - a;
             let ao = r0 - a;
             let aoxab = cross(ao, ab);
             let vxab = cross(v, ab);
             let ab2 = dot(ab, ab);
-            
+
             let aa = dot(vxab, vxab);
             let bb = T::two() * dot(vxab, aoxab);
             let cc = dot(aoxab, aoxab) - (r*r * ab2);
             let dd = bb * bb - T::four() * aa * cc;
-            
+
             if dd >= T::zero() {
                 let t = (-bb - sqrt(dd)) / (T::two() * aa);
                 if t >= T::zero() {
@@ -1306,18 +1311,18 @@ pub fn ray_vs_cylinder<T: Float + FloatOps<T> + NumberOps<T> + SignedNumberOps<T
     let b = c1;
     let v = rv;
     let r = cr;
-    
+
     let ab = b - a;
     let ao = r0 - a;
     let aoxab = cross(ao, ab);
     let vxab = cross(v, ab);
     let ab2 = dot(ab, ab);
-    
+
     let aa = dot(vxab, vxab);
     let bb = T::two() * dot(vxab, aoxab);
     let cc = dot(aoxab, aoxab) - (r*r * ab2);
     let dd = bb * bb - T::four() * aa * cc;
-    
+
     if dd >= T::zero() {
         let t = (-bb - sqrt(dd)) / (T::two() * aa);
         if t >= T::zero() {
@@ -1330,7 +1335,7 @@ pub fn ray_vs_cylinder<T: Float + FloatOps<T> + NumberOps<T> + SignedNumberOps<T
             }
         }
     }
-    
+
     // intersect with the top and bottom circles
     let ip_top = ray_vs_plane(r0, rv, c0, normalize(c0 - c1));
     let ip_bottom = ray_vs_plane(r0, rv, c1, normalize(c1 - c0));
@@ -1359,7 +1364,7 @@ pub fn ray_vs_cylinder<T: Float + FloatOps<T> + NumberOps<T> + SignedNumberOps<T
             }
         }
     }
-    
+
     if btop {
         if let Some(ip_top) = ip_top {
             return Some(ip_top);
@@ -1470,7 +1475,7 @@ pub fn ray_vs_line_segment<T: Float + SignedNumberOps<T> + FloatOps<T>>(r0: Vec3
 /// returns the shortest line segment between 2 line segments `p1-p2` and `p3-p4` as an option tuple where `.0` is the point on line segment 1 and `.1` is the point on line segment 2
 pub fn shortest_line_segment_between_line_segments<T: Float + SignedNumberOps<T> + FloatOps<T>, V: VecN<T> + SignedNumberOps<T> + FloatOps<T> + Dot<T>>(p1: V, p2: V, p3: V, p4: V) -> Option<(V, V)> {
     // https://web.archive.org/web/20120404121511/http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline3d/lineline.c
-    
+
     let p13 = p1 - p3;
     let p43 = p4 - p3;
 
@@ -1507,7 +1512,7 @@ pub fn shortest_line_segment_between_line_segments<T: Float + SignedNumberOps<T>
 /// returns the shortest line segment between 2 lines `p1-p2` and `p3-p4` as an option tuple where `.0` is the point on line segment 1 and `.1` is the point on line segment 2
 pub fn shortest_line_segment_between_lines<T: Float + SignedNumberOps<T> + FloatOps<T>, V: VecN<T> + SignedNumberOps<T> + FloatOps<T> + Dot<T>>(p1: V, p2: V, p3: V, p4: V) -> Option<(V, V)> {
     // https://web.archive.org/web/20120404121511/http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline3d/lineline.c
-    
+
     let p13 = p1 - p3;
     let p43 = p4 - p3;
 
@@ -1544,7 +1549,7 @@ pub fn shortest_line_segment_between_lines<T: Float + SignedNumberOps<T> + Float
 /// returns the shortest line segment between 2 line segments `p1-p2` and `p3-p4` as an option tuple where `.0` is the point on line segment 1 and `.1` is the point on line segment 2
 pub fn shortest_line_segment_between_line_and_line_segment<T: Float + SignedNumberOps<T> + FloatOps<T>, V: VecN<T> + SignedNumberOps<T> + FloatOps<T> + Dot<T>>(p1: V, p2: V, p3: V, p4: V) -> Option<(V, V)> {
     // https://web.archive.org/web/20120404121511/http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline3d/lineline.c
-    
+
     let p13 = p1 - p3;
     let p43 = p4 - p3;
 
@@ -1652,7 +1657,7 @@ pub fn exp_step<T: SignedNumber + Float, X: Base<T> + FloatOps<T> + SignedNumber
 /// returns gain (y position on a graph for `x`); remapping the unit interval into the unit interval by expanding the sides and compressing the center
 pub fn gain<T: SignedNumber + Float + FloatOps<T>>(x: T, k: T) -> T {
     // inigo quilez: https://iquilezles.org/articles/functions/
-    let y = if x < T::point_five() { 
+    let y = if x < T::point_five() {
         x
     }
     else {
@@ -1690,7 +1695,7 @@ pub fn sinc<T: SignedNumber + Float, X: Base<T> + FloatOps<T> + SignedNumberOps<
 /// returns a hsv value in 0-1 range converted from `rgb` in 0-1 range
 pub fn rgb_to_hsv<T: Float + SignedNumberOps<T> + Cast<T>>(rgb: Vec3<T>) -> Vec3<T> {
     // from Foley & van Dam p592
-    // optimized: http://lolengine.net/blog/2013/01/13/fast-rgb-to-hsv 
+    // optimized: http://lolengine.net/blog/2013/01/13/fast-rgb-to-hsv
     let mut r = rgb.x;
     let mut g = rgb.y;
     let mut b = rgb.z;
@@ -1714,14 +1719,14 @@ pub fn rgb_to_hsv<T: Float + SignedNumberOps<T> + Cast<T>>(rgb: Vec3<T>) -> Vec3
         z: r
     }
 }
- 
+
 /// returns an rgb value in 0-1 range converted from `hsv` in 0-1 range
 pub fn hsv_to_rgb<T: Float + FloatOps<T> + Cast<T>>(hsv: Vec3<T>) -> Vec3<T> {
     // from Foley & van Dam p593: http://en.wikipedia.org/wiki/HSL_and_HSV
     let h = hsv.x;
     let s = hsv.y;
     let v = hsv.z;
-        
+
     if s == T::zero() {
         // gray
         return Vec3 {
@@ -1940,7 +1945,7 @@ pub fn gjk_mesh_support_function<T: Float + FloatOps<T> + NumberOps<T> + SignedN
         }
         fv
     };
-    
+
     // selects the furthest points on the 2 meshes in opposite directions
     let fp0 = furthest_point(dir, convex0);
     let fp1 = furthest_point(-dir, convex1);
@@ -1955,23 +1960,23 @@ fn handle_simplex_2d<T: Float + FloatOps<T> + NumberOps<T> + SignedNumber + Sign
             let b = simplex[0];
             let ab = b - a;
             let ao = -a;
-            
+
             *dir = vector_triple(ab, ao, ab);
-            
+
             false
         },
         3 => {
             let a = simplex[2];
             let b = simplex[1];
             let c = simplex[0];
-            
+
             let ab = b - a;
             let ac = c - a;
             let ao = -a;
-            
+
             let abperp = vector_triple(ac, ab, ab);
             let acperp = vector_triple(ab, ac, ac);
-            
+
             if dot(abperp, ao) > T::zero() {
                 simplex.remove(0);
                 *dir = abperp;
@@ -1995,29 +2000,29 @@ fn handle_simplex_2d<T: Float + FloatOps<T> + NumberOps<T> + SignedNumber + Sign
 /// returns true if the 2d convex hull `convex0` overlaps with `convex1` using the gjk algorithm
 pub fn gjk_2d<T: Float + FloatOps<T> + NumberOps<T> + SignedNumber + SignedNumberOps<T>, V: VecN<T> + FloatOps<T> + SignedNumberOps<T> + VecFloatOps<T> + Triple<T>>(convex0: Vec<V>, convex1: Vec<V>) -> bool {
     // implemented following details in this insightful video: https://www.youtube.com/watch?v=ajv46BSqcK4
-    
+
     // start with arbitrary direction
     let mut dir = V::unit_x();
     let support = gjk_mesh_support_function(&convex0, &convex1, dir);
     dir = normalize(-support);
-    
+
     // iterative build and test simplex
     let mut simplex = vec![support];
-    
+
     let max_iters = 32;
     for _i in 0..max_iters {
         let a = gjk_mesh_support_function(&convex0, &convex1, dir);
-        
+
         if dot(a, dir) < T::zero() {
             return false;
         }
         simplex.push(a);
-        
+
         if handle_simplex_2d(&mut simplex, &mut dir) {
             return true;
         }
     }
-    
+
     // if we reach here we likely have got stuck in a simplex building loop, we assume the shapes are touching but not intersecting
     false
 }
@@ -2033,30 +2038,30 @@ fn handle_simplex_3d<T: Float + FloatOps<T> + NumberOps<T> + SignedNumber + Sign
         2 => {
             let a = simplex[1];
             let b = simplex[0];
-            
+
             let ab = b - a;
             let ao = -a;
-            
+
             *dir = vector_triple(ab, ao, ab);
-            
+
             false
         },
         3 => {
             let a = simplex[2];
             let b = simplex[1];
             let c = simplex[0];
-            
+
             let ab = b - a;
             let ac = c - a;
             let ao = -a;
-            
+
             *dir = cross(ac, ab);
-    
+
             // flip normal so it points toward the origin
             if dot(*dir, ao) < T::zero() {
                 *dir = -*dir;
             }
-    
+
             false
         },
         4 => {
@@ -2064,35 +2069,35 @@ fn handle_simplex_3d<T: Float + FloatOps<T> + NumberOps<T> + SignedNumber + Sign
             let b = simplex[2];
             let c = simplex[1];
             let d = simplex[0];
-            
+
             let centre = (a+b+c+d) / T::four();
-            
+
             let ab = b - a;
             let ac = c - a;
             let ad = d - a;
             let ao = -a;
-            
+
             let mut abac = cross(ab, ac);
             let mut acad = cross(ac, ad);
             let mut adab = cross(ad, ab);
-            
+
             // flip the normals so they always face outward
             let centre_abc = (a + b + c) / T::three();
             let centre_acd = (a + c + d) / T::three();
             let centre_adb = (a + d + b) / T::three();
-            
+
             if dot(centre - centre_abc, abac) > T::zero() {
                 abac = -abac;
             }
-            
+
             if dot(centre - centre_acd, acad) > T::zero() {
                 acad = -acad;
             }
-            
+
             if dot(centre - centre_adb, adab) > T::zero() {
                 adab = -adab;
             }
-            
+
             if dot(abac, ao) > T::zero() {
                 // erase c
                 simplex.remove(0);
@@ -2124,29 +2129,29 @@ fn handle_simplex_3d<T: Float + FloatOps<T> + NumberOps<T> + SignedNumber + Sign
 /// returns true if the 3D convex hull `convex0` overlaps with `convex1` using the gjk algorithm
 pub fn gjk_3d<T: Float + FloatOps<T> + NumberOps<T> + SignedNumber + SignedNumberOps<T>, V: VecN<T> + FloatOps<T> + SignedNumberOps<T> + VecFloatOps<T> + Triple<T> + Cross<T>>(convex0: Vec<V>, convex1: Vec<V>) -> bool {
     // implemented following details in this insightful video: https://www.youtube.com/watch?v=ajv46BSqcK4
-    
+
     // start with arbitrary direction
     let mut dir = V::unit_x();
     let support = gjk_mesh_support_function(&convex0, &convex1, dir);
     dir = normalize(-support);
-    
+
     // iterative build and test simplex
     let mut simplex = vec![support];
-    
+
     let max_iters = 32;
     for _i in 0..max_iters {
         let a = gjk_mesh_support_function(&convex0, &convex1, dir);
-        
+
         if dot(a, dir) < T::zero() {
             return false;
         }
         simplex.push(a);
-        
+
         if handle_simplex_3d(&mut simplex, &mut dir) {
             return true;
         }
     }
-    
+
     // if we reach here we likely have got stuck in a simplex building loop, we assume the shapes are touching but not intersecting
     false
 }
diff --git a/src/num.rs b/src/num.rs
index b4a04fa..cce54ff 100644
--- a/src/num.rs
+++ b/src/num.rs
@@ -25,12 +25,12 @@ use std::cmp::PartialOrd;
 
 use std::fmt::Display;
 
-/// base trait for scalar and vector numerical operations, arithmetic and generic constants 
+/// base trait for scalar and vector numerical operations, arithmetic and generic constants
 pub trait Base<T: Number>:
     Copy + Display +
     Add<Output=Self> + AddAssign +
-    Sub<Output=Self> + SubAssign + 
-    Mul<Output=Self> + MulAssign + 
+    Sub<Output=Self> + SubAssign +
+    Mul<Output=Self> + MulAssign +
     Div<Output=Self> + DivAssign +
     Rem<Output=Self> + RemAssign
     where Self: Sized {
@@ -135,8 +135,10 @@ pub trait FloatOps<T: Float>: Lerp<T> where Self: Sized {
     fn ceil(a: Self) -> Self;
     /// returns value `a` with the same sign as the second parameter `sign`
     fn copysign(a: Self, sign: T) -> Self;
-    /// returns hermite interpolation between `0-1` of `t` between edges `e0` and `e1`
+    /// returns cubic hermite interpolation between `0-1` of `t` between edges `e0` and `e1`
     fn smoothstep(e0: Self, e1: Self, t: T) -> Self;
+    /// returns quintic hermite interpolation between `0-1` of `t` between edges `e0` and `e1`
+    fn smootherstep(e0: Self, e1: Self, t: T) -> Self;
     /// returns `a` rounded component wise
     fn round(a: Self) -> Self;
     /// returns true if `a` is not a number (`nan`)
@@ -157,7 +159,7 @@ pub trait FloatOps<T: Float>: Lerp<T> where Self: Sized {
     fn frac(v: Self) -> Self;
     /// returns the integer part of float value `v` truncating the decimal part
     fn trunc(v: Self) -> Self;
-    /// returns a tuple containing `(frac(v), trunc(v))` breaking the float into 2 parts 
+    /// returns a tuple containing `(frac(v), trunc(v))` breaking the float into 2 parts
     fn modf(v: Self) -> (Self, Self);
     /// returns the cosine of `v` where the value `v` is in radians
     fn cos(v: Self) -> Self;
@@ -263,7 +265,7 @@ macro_rules! number_impl {
             fn step(a: Self, b: Self) -> Self {
                 if a >= b {
                     Self::one()
-                } 
+                }
                 else {
                     Self::zero()
                 }
@@ -272,31 +274,31 @@ macro_rules! number_impl {
 
         impl Cast<$t> for $t {
             fn from_f32(v: f32) -> Self {
-                v as Self 
+                v as Self
             }
 
             fn from_f64(v: f64) -> Self {
-                v as Self 
+                v as Self
             }
 
             fn from_u32(v: u32) -> Self {
-                v as Self 
+                v as Self
             }
 
             fn from_i32(v: i32) -> Self {
-                v as Self 
+                v as Self
             }
 
             fn from_u64(v: u64) -> Self {
-                v as Self 
+                v as Self
             }
 
             fn from_i64(v: i64) -> Self {
-                v as Self 
+                v as Self
             }
 
             fn from_usize(v: usize) -> Self {
-                v as Self 
+                v as Self
             }
 
             fn as_f32(&self) -> f32 {
@@ -318,7 +320,7 @@ macro_rules! number_impl {
             fn as_u64(&self) -> u64 {
                 *self as u64
             }
-            
+
             fn as_i64(&self) -> i64 {
                 *self as i64
             }
@@ -399,7 +401,7 @@ macro_rules! float_impl {
                 0.5 as Self
             }
 
-            #[allow(clippy::excessive_precision)] 
+            #[allow(clippy::excessive_precision)]
             fn pi() -> Self {
                 3.14159265358979323846264338327950288 as Self
             }
@@ -412,7 +414,7 @@ macro_rules! float_impl {
                 1.0 as Self / Self::pi() as Self
             }
 
-            #[allow(clippy::excessive_precision)] 
+            #[allow(clippy::excessive_precision)]
             fn phi() -> Self {
                 1.618033988749894 as Self
             }
@@ -421,7 +423,7 @@ macro_rules! float_impl {
                 1.0 as Self / Self::phi() as Self
             }
 
-            #[allow(clippy::excessive_precision)] 
+            #[allow(clippy::excessive_precision)]
             fn tau() -> Self {
                 6.2831853071795864 as Self
             }
@@ -461,6 +463,13 @@ macro_rules! float_impl {
                 x * x * (3 as Self - 2 as Self * x)
             }
 
+            fn smootherstep(e0: Self, e1: Self, t: Self) -> Self {
+                if t < e0 { return Self::zero(); }
+                if (t >= e1) { return Self::one(); }
+                let x = (t - e0) / (e1 - e0);
+                x * x * x * (3 as Self * x * (2 as Self * x - 5 as Self) + 10 as Self)
+            }
+
             fn saturate(v: Self) -> Self {
                 Self::max(Self::min(v, 1.0), 0.0)
             }
@@ -468,7 +477,7 @@ macro_rules! float_impl {
             fn powi(v: Self, exp: i32) -> Self {
                 v.powi(exp)
             }
-            
+
             fn powf(v: Self, exp: $t) -> Self {
                 v.powf(exp)
             }
@@ -515,12 +524,12 @@ macro_rules! float_impl {
 macro_rules! integer_trait_impl {
     ($($func:ident),*) => {
         /// integer point trait for various sized integers
-        pub trait Integer: Number + 
+        pub trait Integer: Number +
             Shl<Output=Self> + ShlAssign +
-            Shr<Output=Self> + ShrAssign + 
+            Shr<Output=Self> + ShrAssign +
             BitOr<Output=Self> + BitOrAssign +
-            BitAnd<Output=Self> + BitAndAssign + 
-            BitXor<Output=Self> + BitXorAssign { 
+            BitAnd<Output=Self> + BitAndAssign +
+            BitXor<Output=Self> + BitXorAssign {
         }
         integer_impl!(i8 { $($func),* });
         integer_impl!(u8 { $($func),* });
diff --git a/src/vec.rs b/src/vec.rs
index 300e75c..fd1b717 100644
--- a/src/vec.rs
+++ b/src/vec.rs
@@ -22,13 +22,13 @@ use std::fmt::Formatter;
 use crate::num::*;
 
 //
-// Vec Traits 
+// Vec Traits
 //
 
 /// generic vec trait to allow sized vectors to be treated generically
-pub trait VecN<T: Number>: 
-    Base<T> + Dot<T> + 
-    Index<usize, Output=T> + IndexMut<usize> + 
+pub trait VecN<T: Number>:
+    Base<T> + Dot<T> +
+    Index<usize, Output=T> + IndexMut<usize> +
     Add<T, Output=Self> + Sub<T, Output=Self> +
     Mul<T, Output=Self> + Div<T, Output=Self> {
     /// returns the count of elements in the vector type
@@ -105,8 +105,10 @@ pub trait VecFloatOps<T: Float>: SignedVecN<T> + Magnitude<T> {
     fn refract(i: Self, n: Self, eta: T) -> Self;
     /// returns linear interpolation between `e0` and `e1`, `t` specifies the ratio to interpolate between the values, with component-wise `t`
     fn vlerp(e0: Self, e1: Self, t: Self) -> Self;
-    /// returns vector with component wise hermite interpolation between `0-1`, with component-wise `t`
+    /// returns vector with component wise cubic hermite interpolation between `0-1`, with component-wise `t`
     fn vsmoothstep(e0: Self, e1: Self, t: Self) -> Self;
+    /// returns vector with component wise quintic hermite interpolation between `0-1`, with component-wise `t`
+    fn vsmootherstep(e0: Self, e1: Self, t: Self) -> Self;
     /// returns vector `a` raised to component wise power `exp`
     fn powfn(a: Self, exp: Self) -> Self;
 }
@@ -144,7 +146,7 @@ pub trait Cross<T> {
 impl<T> Cross<T> for Vec3<T> where T: Number {
     fn cross(a: Self, b: Self) -> Self {
         Vec3 {
-            x: (a.y * b.z) - (a.z * b.y), 
+            x: (a.y * b.z) - (a.z * b.y),
             y: (a.z * b.x) - (a.x * b.z),
             z: (a.x * b.y) - (a.y * b.x),
         }
@@ -178,7 +180,7 @@ impl<T> Triple<T> for Vec2<T> where T: Number + SignedNumber {
 
     fn vector_triple(a: Self, b: Self, c: Self) -> Self {
         let a3 = Vec3::<T>::new(a.x, a.y, T::zero());
-        let b3 = Vec3::<T>::new(b.x, b.y, T::zero()); 
+        let b3 = Vec3::<T>::new(b.x, b.y, T::zero());
         let c3 = Vec3::<T>::new(c.x, c.y, T::zero());
         let v3 = Vec3::<T>::cross(Vec3::<T>::cross(a3, b3), c3);
         Self {
@@ -200,14 +202,14 @@ pub trait Nlerp<T: Float> {
     fn nlerp(e0: Self, e1: Self, t: T) -> Self;
 }
 
-// 
+//
 // Macro Implementation
 //
 
 /// macro to stamp out various n-dimensional vectors, all of their ops and functions
 macro_rules! vec_impl {
     ($VecN:ident { $($field:ident, $field_index:expr),* }, $len:expr, $module:ident) => {
-        
+
         #[cfg_attr(feature="serde", derive(serde::Serialize, serde::Deserialize))]
         #[derive(Debug, Copy, Clone)]
         #[repr(C)]
@@ -478,6 +480,12 @@ macro_rules! vec_impl {
                 }
             }
 
+            fn vsmootherstep(e0: Self, e1: Self, t: Self) -> Self {
+                Self {
+                    $($field: T::smootherstep(e0.$field, e1.$field, t.$field),)+
+                }
+            }
+
             fn powfn(a: Self, exp: Self) -> Self {
                 Self {
                     $($field: T::powf(a.$field, exp.$field),)+
@@ -488,7 +496,7 @@ macro_rules! vec_impl {
         impl<T> Slerp<T> for $VecN<T> where T: Float + FloatOps<T> + NumberOps<T> {
             fn slerp(e0: Self, e1: Self, t: T) -> Self {
                 // https://blog.demofox.org/2016/02/19/normalized-vector-interpolation-tldr/
-                let dot = Self::dot(e0, e1);     
+                let dot = Self::dot(e0, e1);
                 let dot = T::clamp(dot, T::minus_one(), T::one());
                 let theta = T::acos(dot) * t;
                 let v = Self::normalize(e1 - e0 * dot);
@@ -664,12 +672,18 @@ macro_rules! vec_impl {
                 }
             }
 
+            fn smootherstep(e0: Self, e1: Self, t: T) -> Self {
+                Self {
+                    $($field: T::smootherstep(e0.$field, e1.$field, t),)+
+                }
+            }
+
             fn round(a: Self) -> Self {
                 Self {
                     $($field: T::round(a.$field),)+
                 }
             }
-        
+
             fn is_nan(a: Self) -> Self {
                 Self {
                     $($field: T::is_nan(a.$field),)+
@@ -923,13 +937,13 @@ macro_rules! vec_impl {
                 *self + other
             }
         }
-        
+
         impl<T> AddAssign<Self> for $VecN<T> where T: Number {
             fn add_assign(&mut self, other: Self) {
                 $(self.$field += other.$field;)+
             }
         }
-        
+
         impl<T> AddAssign<T> for $VecN<T> where T: Number {
             fn add_assign(&mut self, other: T) {
                 $(self.$field += other;)+
@@ -942,7 +956,7 @@ macro_rules! vec_impl {
                 $(self.$field += other.$field;)+
             }
         }
-        
+
         impl<T> Sub<Self> for $VecN<T> where T: Number {
             type Output = Self;
             fn sub(self, other: Self) -> Self {
@@ -951,7 +965,7 @@ macro_rules! vec_impl {
                 }
             }
         }
-        
+
         impl<T> Sub<T> for $VecN<T> where T: Number {
             type Output = Self;
             fn sub(self, other: T) -> Self {
@@ -982,13 +996,13 @@ macro_rules! vec_impl {
                 *self - other
             }
         }
-        
+
         impl<T> SubAssign<Self> for $VecN<T> where T: Number {
             fn sub_assign(&mut self, other: Self) {
                 $(self.$field -= other.$field;)+
             }
         }
-        
+
         impl<T> SubAssign<T> for $VecN<T> where T: Number {
             fn sub_assign(&mut self, other: T) {
                 $(self.$field -= other;)+
@@ -1001,7 +1015,7 @@ macro_rules! vec_impl {
                 $(self.$field -= other.$field;)+
             }
         }
-        
+
         impl<T> Mul<Self> for $VecN<T> where T: Number {
             type Output = Self;
             fn mul(self, other: Self) -> Self {
@@ -1010,7 +1024,7 @@ macro_rules! vec_impl {
                 }
             }
         }
-        
+
         impl<T> Mul<T> for $VecN<T> where T: Number {
             type Output = Self;
             fn mul(self, other: T) -> Self {
@@ -1041,13 +1055,13 @@ macro_rules! vec_impl {
                 *self * other
             }
         }
-        
+
         impl<T> MulAssign<Self> for $VecN<T> where T: Number {
             fn mul_assign(&mut self, other: Self) {
                 $(self.$field *= other.$field;)+
             }
         }
-        
+
         impl<T> MulAssign<T> for $VecN<T> where T: Number {
             fn mul_assign(&mut self, other: T) {
                 $(self.$field *= other;)+
@@ -1060,7 +1074,7 @@ macro_rules! vec_impl {
                 $(self.$field *= other.$field;)+
             }
         }
-        
+
         impl<T> Div<Self> for $VecN<T> where T: Number {
             type Output = Self;
             fn div(self, other: Self) -> Self {
@@ -1069,7 +1083,7 @@ macro_rules! vec_impl {
                 }
             }
         }
-        
+
         impl<T> Div<T> for $VecN<T> where T: Number {
             type Output = Self;
             fn div(self, other: T) -> Self {
@@ -1100,13 +1114,13 @@ macro_rules! vec_impl {
                 *self / other
             }
         }
-        
+
         impl<T> DivAssign<Self> for $VecN<T> where T: Number {
             fn div_assign(&mut self, other: Self) {
                 $(self.$field /= other.$field;)+
             }
         }
-        
+
         impl<T> DivAssign<T> for $VecN<T> where T: Number {
             fn div_assign(&mut self, other: T) {
                 $(self.$field /= other;)+
@@ -1128,7 +1142,7 @@ macro_rules! vec_impl {
                 }
             }
         }
-        
+
         impl<T> Rem<T> for $VecN<T> where T: Number {
             type Output = Self;
             fn rem(self, other: T) -> Self {
@@ -1137,19 +1151,19 @@ macro_rules! vec_impl {
                 }
             }
         }
-        
+
         impl<T> RemAssign<Self> for $VecN<T> where T: Number {
             fn rem_assign(&mut self, other: Self) {
                 $(self.$field %= other.$field;)+
             }
         }
-        
+
         impl<T> RemAssign<T> for $VecN<T> where T: Number {
             fn rem_assign(&mut self, other: T) {
                 $(self.$field %= other;)+
             }
         }
-        
+
         impl<T> Neg for $VecN<T> where T: SignedNumber {
             type Output = Self;
             fn neg(self) -> Self::Output {
@@ -1158,7 +1172,7 @@ macro_rules! vec_impl {
                 }
             }
         }
-        
+
         impl<T> Eq for $VecN<T> where T: Eq  {}
         impl<T> PartialEq for $VecN<T> where T: PartialEq  {
             fn eq(&self, other: &Self) -> bool {