diff --git a/README.md b/README.md
index 107f751..55ae01e 100644
--- a/README.md
+++ b/README.md
@@ -977,6 +977,9 @@ T fract(T x) noexcept;
 template < floating_point T >
 T fmod(T x, T y) noexcept;
 
+template < floating_point T >
+T modf(T x, T* y) noexcept;
+
 template < floating_point T >
 T copysign(T x, T s) noexcept;
 
@@ -1044,6 +1047,9 @@ vec<T, Size> fmod(const vec<T, Size>& xs, T y);
 template < typename T, size_t Size >
 vec<T, Size> fmod(const vec<T, Size>& xs, const vec<T, Size>& ys);
 
+template < typename T, size_t Size >
+vec<T, Size> modf(const vec<T, Size>& xs, vec<T, Size>* is);
+
 template < typename T, size_t Size >
 vec<T, Size> copysign(const vec<T, Size>& xs, T s);
 
diff --git a/headers/vmath.hpp/vmath_fun.hpp b/headers/vmath.hpp/vmath_fun.hpp
index d0500e1..5417869 100644
--- a/headers/vmath.hpp/vmath_fun.hpp
+++ b/headers/vmath.hpp/vmath_fun.hpp
@@ -78,6 +78,12 @@ namespace vmath_hpp
         return std::fmod(x, y);
     }
 
+    template < typename T >
+    [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
+    modf(T x, T* y) noexcept {
+        return std::modf(x, y);
+    }
+
     template < typename T >
     [[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, T>
     copysign(T x, T s) noexcept {
diff --git a/headers/vmath.hpp/vmath_vec_fun.hpp b/headers/vmath.hpp/vmath_vec_fun.hpp
index 8b66c28..4404f37 100644
--- a/headers/vmath.hpp/vmath_vec_fun.hpp
+++ b/headers/vmath.hpp/vmath_vec_fun.hpp
@@ -505,6 +505,20 @@ namespace vmath_hpp
         return map_join([](T x, T y) { return fmod(x, y); }, xs, ys);
     }
 
+    namespace impl
+    {
+        template < typename T, std::size_t Size, std::size_t... Is >
+        VMATH_HPP_FORCE_INLINE
+        vec<T, Size> modf_impl(const vec<T, Size>& xs, vec<T, Size>* is, std::index_sequence<Is...>) {
+            return { modf(xs[Is], &(*is)[Is])... };
+        }
+    }
+
+    template < typename T, std::size_t Size >
+    vec<T, Size> modf(const vec<T, Size>& xs, vec<T, Size>* is) {
+        return impl::modf_impl(xs, is, std::make_index_sequence<Size>{});
+    }
+
     template < typename T, std::size_t Size >
     [[nodiscard]] vec<T, Size> copysign(const vec<T, Size>& xs, T s) {
         return map_join([s](T x) { return copysign(x, s); }, xs);
diff --git a/untests/vmath_fix_tests.cpp b/untests/vmath_fix_tests.cpp
index 12f03df..0a9e21e 100644
--- a/untests/vmath_fix_tests.cpp
+++ b/untests/vmath_fix_tests.cpp
@@ -34,6 +34,7 @@ namespace
         constexpr friend fix fract(const fix& l) { using vmath_hpp::fract; return { fix{fract(l.underlying())} }; }
 
         constexpr friend fix fmod(const fix& l, const fix& r) { using vmath_hpp::fmod; return { fix{fmod(l.underlying(), r.underlying())} }; }
+        constexpr friend fix modf(const fix& l, fix* r) { using vmath_hpp::modf; return { fix{modf(l.underlying(), &r->underlying_)} }; }
         constexpr friend fix copysign(const fix& l, const fix& r) { using vmath_hpp::copysign; return { fix{copysign(l.underlying(), r.underlying())} }; }
 
         constexpr friend fix min(const fix& l, const fix& r) { using vmath_hpp::min; return fix{min(l.underlying(), r.underlying())}; }
@@ -67,31 +68,31 @@ namespace
 
         constexpr friend fix radians(const fix& l) { using vmath_hpp::radians; return fix{radians(l.underlying())}; }
         constexpr friend fix degrees(const fix& l) { using vmath_hpp::degrees; return fix{degrees(l.underlying())}; }
-        friend fix sin(const fix& l) { using vmath_hpp::sin; return fix{sin(l.underlying())}; }
-        friend fix cos(const fix& l) { using vmath_hpp::cos; return fix{cos(l.underlying())}; }
-        friend fix tan(const fix& l) { using vmath_hpp::tan; return fix{tan(l.underlying())}; }
-        friend fix asin(const fix& l) { using vmath_hpp::asin; return fix{asin(l.underlying())}; }
-        friend fix acos(const fix& l) { using vmath_hpp::acos; return fix{acos(l.underlying())}; }
-        friend fix atan(const fix& l) { using vmath_hpp::atan; return fix{atan(l.underlying())}; }
-        friend fix atan2(const fix& l, const fix& r) { using vmath_hpp::atan2; return fix{atan2(l.underlying(), r.underlying())}; }
-        friend fix sinh(const fix& l) { using vmath_hpp::sinh; return fix{sinh(l.underlying())}; }
-        friend fix cosh(const fix& l) { using vmath_hpp::cosh; return fix{cosh(l.underlying())}; }
-        friend fix tanh(const fix& l) { using vmath_hpp::tanh; return fix{tanh(l.underlying())}; }
-        friend fix asinh(const fix& l) { using vmath_hpp::asinh; return fix{asinh(l.underlying())}; }
-        friend fix acosh(const fix& l) { using vmath_hpp::acosh; return fix{acosh(l.underlying())}; }
-        friend fix atanh(const fix& l) { using vmath_hpp::atanh; return fix{atanh(l.underlying())}; }
-        friend std::pair<fix,fix> sincos(const fix& l) { return { sin(l), cos(l) }; }
-        friend void sincos(const fix& l, fix* s, fix* c) { *s = sin(l); *c = cos(l); }
+        constexpr friend fix sin(const fix& l) { using vmath_hpp::sin; return fix{sin(l.underlying())}; }
+        constexpr friend fix cos(const fix& l) { using vmath_hpp::cos; return fix{cos(l.underlying())}; }
+        constexpr friend fix tan(const fix& l) { using vmath_hpp::tan; return fix{tan(l.underlying())}; }
+        constexpr friend fix asin(const fix& l) { using vmath_hpp::asin; return fix{asin(l.underlying())}; }
+        constexpr friend fix acos(const fix& l) { using vmath_hpp::acos; return fix{acos(l.underlying())}; }
+        constexpr friend fix atan(const fix& l) { using vmath_hpp::atan; return fix{atan(l.underlying())}; }
+        constexpr friend fix atan2(const fix& l, const fix& r) { using vmath_hpp::atan2; return fix{atan2(l.underlying(), r.underlying())}; }
+        constexpr friend fix sinh(const fix& l) { using vmath_hpp::sinh; return fix{sinh(l.underlying())}; }
+        constexpr friend fix cosh(const fix& l) { using vmath_hpp::cosh; return fix{cosh(l.underlying())}; }
+        constexpr friend fix tanh(const fix& l) { using vmath_hpp::tanh; return fix{tanh(l.underlying())}; }
+        constexpr friend fix asinh(const fix& l) { using vmath_hpp::asinh; return fix{asinh(l.underlying())}; }
+        constexpr friend fix acosh(const fix& l) { using vmath_hpp::acosh; return fix{acosh(l.underlying())}; }
+        constexpr friend fix atanh(const fix& l) { using vmath_hpp::atanh; return fix{atanh(l.underlying())}; }
+        constexpr friend std::pair<fix,fix> sincos(const fix& l) { return { sin(l), cos(l) }; }
+        constexpr friend void sincos(const fix& l, fix* s, fix* c) { *s = sin(l); *c = cos(l); }
 
         //
 
-        friend fix pow(const fix& l, const fix& r) { using vmath_hpp::pow; return fix{pow(l.underlying(), r.underlying())}; }
-        friend fix exp(const fix& l) { using vmath_hpp::exp; return fix{exp(l.underlying())}; }
-        friend fix log(const fix& l) { using vmath_hpp::log; return fix{log(l.underlying())}; }
-        friend fix exp2(const fix& l) { using vmath_hpp::exp2; return fix{exp2(l.underlying())}; }
-        friend fix log2(const fix& l) { using vmath_hpp::log2; return fix{log2(l.underlying())}; }
-        friend fix sqrt(const fix& l) { using vmath_hpp::sqrt; return fix{sqrt(l.underlying())}; }
-        friend fix rsqrt(const fix& l) { using vmath_hpp::rsqrt; return fix{rsqrt(l.underlying())}; }
+        constexpr friend fix pow(const fix& l, const fix& r) { using vmath_hpp::pow; return fix{pow(l.underlying(), r.underlying())}; }
+        constexpr friend fix exp(const fix& l) { using vmath_hpp::exp; return fix{exp(l.underlying())}; }
+        constexpr friend fix log(const fix& l) { using vmath_hpp::log; return fix{log(l.underlying())}; }
+        constexpr friend fix exp2(const fix& l) { using vmath_hpp::exp2; return fix{exp2(l.underlying())}; }
+        constexpr friend fix log2(const fix& l) { using vmath_hpp::log2; return fix{log2(l.underlying())}; }
+        constexpr friend fix sqrt(const fix& l) { using vmath_hpp::sqrt; return fix{sqrt(l.underlying())}; }
+        constexpr friend fix rsqrt(const fix& l) { using vmath_hpp::rsqrt; return fix{rsqrt(l.underlying())}; }
 
         //
 
@@ -169,6 +170,8 @@ namespace vmath_hpp
     template fix2f fmod(const fix2f&, fix<float>);
     template fix2f fmod(const fix2f&, const fix2f&);
 
+    template fix2f modf(const fix2f&, fix2f*);
+
     template fix2f copysign(const fix2f&, fix<float>);
     template fix2f copysign(const fix2f&, const fix2f&);
 
diff --git a/untests/vmath_fun_tests.cpp b/untests/vmath_fun_tests.cpp
index 484afa7..a681dc0 100644
--- a/untests/vmath_fun_tests.cpp
+++ b/untests/vmath_fun_tests.cpp
@@ -77,6 +77,12 @@ TEST_CASE("vmath/fun") {
 
         CHECK(fmod(1.7f, 1.2f) == uapprox(0.5f));
 
+        {
+            float out_i{};
+            CHECK(modf(1.7f, &out_i) == uapprox(0.7f));
+            CHECK(out_i == uapprox(1.f));
+        }
+
         STATIC_CHECK(min(0.f, 1.f) == uapprox(0.f));
         STATIC_CHECK(max(0.f, 1.f) == uapprox(1.f));
 
diff --git a/untests/vmath_vec_fun_tests.cpp b/untests/vmath_vec_fun_tests.cpp
index 2919ed8..ea893fe 100644
--- a/untests/vmath_vec_fun_tests.cpp
+++ b/untests/vmath_vec_fun_tests.cpp
@@ -205,6 +205,12 @@ TEST_CASE("vmath/vec_fun") {
         CHECK(fmod(float2(1.7f), 1.2f) == uapprox2(0.5f));
         CHECK(fmod(float2(1.7f), float2(1.2f)) == uapprox2(0.5f));
 
+        {
+            float2 out_i{};
+            CHECK(modf(float2(1.7f), &out_i) == uapprox2(0.7f));
+            CHECK(out_i.x == uapprox(1.f));
+        }
+
         STATIC_CHECK(min(int2(1,2)) == 1);
         STATIC_CHECK(min(int2(1,2), 1) == int2(1,1));
         STATIC_CHECK(min(1, int2(1,2)) == int2(1,1));