diff --git a/include/cglm/bezier.h b/include/cglm/bezier.h
new file mode 100644
index 0000000..f36ac1d
--- /dev/null
+++ b/include/cglm/bezier.h
@@ -0,0 +1,56 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#ifndef cglm_bezier_h
+#define cglm_bezier_h
+
+#define GLM_BEZIER_MAT_INIT  {{-1.0f,  3.0f, -3.0f, 1.0f},                    \
+                              { 3.0f, -6.0f,  3.0f, 0.0f},                    \
+                              {-3.0f,  3.0f,  0.0f, 0.0f},                    \
+                              { 1.0f,  0.0f,  0.0f, 0.0f}}
+
+/* for C only */
+#define GLM_BEZIER_MAT ((mat4)GLM_BEZIER_MAT_INIT)
+
+#define CGLM_DECASTEL_EPS   1e-9
+#define CGLM_DECASTEL_MAX   1000
+#define CGLM_DECASTEL_SMALL 1e-20
+
+/*!
+ * @brief cubic bezier interpolation
+ *
+ * Formula:
+ *  B(s) = P0*(1-s)^3 + 3*C0*s*(1-s)^2 + 3*C1*s^2*(1-s) + P1*s^3
+ *
+ * similar result using matrix:
+ *  B(s) = glm_smc(t, GLM_BEZIER_MAT, (vec4){p0, c0, c1, p1})
+ *
+ * glm_eq(glm_smc(...), glm_bezier(...)) should return TRUE
+ *
+ * @param[in]  s    parameter between 0 and 1
+ * @param[in]  p0   begin point
+ * @param[in]  c0   control point 1
+ * @param[in]  c1   control point 2
+ * @param[in]  p1   end point
+ *
+ * @return B(s)
+ */
+CGLM_INLINE
+float
+glm_bezier(float s, float p0, float c0, float c1, float p1) {
+  float x, xx, ss, xs3, a;
+
+  x   = 1.0f - s;
+  xx  = x * x;
+  ss  = s * s;
+  xs3 = (s - ss) * 3.0f;
+  a   = p0 * xx + c0 * xs3;
+
+  return a + s * (c1 * xs3 + p1 * ss - a);
+}
+
+#endif /* cglm_bezier_h */
diff --git a/include/cglm/call/bezier.h b/include/cglm/call/bezier.h
new file mode 100644
index 0000000..531c15d
--- /dev/null
+++ b/include/cglm/call/bezier.h
@@ -0,0 +1,23 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#ifndef cglmc_curve_h
+#define cglmc_curve_h
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include "../cglm.h"
+
+CGLM_EXPORT
+float
+glmc_bezier(float s, float p0, float c0, float c1, float p1);
+
+#ifdef __cplusplus
+}
+#endif
+#endif /* cglmc_curve_h */
diff --git a/include/cglm/cglm.h b/include/cglm/cglm.h
index d79a88e..7c301bf 100644
--- a/include/cglm/cglm.h
+++ b/include/cglm/cglm.h
@@ -27,5 +27,6 @@
 #include "sphere.h"
 #include "ease.h"
 #include "curve.h"
+#include "bezier.h"
 
 #endif /* cglm_h */
diff --git a/src/bezier.c b/src/bezier.c
new file mode 100644
index 0000000..244e580
--- /dev/null
+++ b/src/bezier.c
@@ -0,0 +1,15 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#include "../include/cglm/cglm.h"
+#include "../include/cglm/call.h"
+
+CGLM_EXPORT
+float
+glmc_bezier(float s, float p0, float c0, float c1, float p1) {
+  return glm_bezier(s, p0, c0, c1, p1);
+}
diff --git a/test/src/test_bezier.c b/test/src/test_bezier.c
new file mode 100644
index 0000000..32d8c35
--- /dev/null
+++ b/test/src/test_bezier.c
@@ -0,0 +1,41 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#include "test_common.h"
+
+CGLM_INLINE
+float
+test_bezier_plain(float s, float p0, float c0, float c1, float p1) {
+  float x, xx, xxx, ss, sss;
+
+  x   = 1.0f - s;
+  xx  = x * x;
+  xxx = xx * x;
+  ss  = s * s;
+  sss = ss * s;
+
+  return p0 * xxx + 3.0f * (c0 * s * xx + c1 * ss * x) + p1 * sss;
+}
+
+void
+test_bezier(void **state) {
+  float s, p0, p1, c0, c1, smc, Bs, Bs_plain;
+
+  s        = test_rand();
+  p0       = test_rand();
+  p1       = test_rand();
+  c0       = test_rand();
+  c1       = test_rand();
+
+  smc      = glm_smc(s, GLM_BEZIER_MAT, (vec4){p0, c0, c1, p1});
+  Bs       = glm_bezier(s, p0, c0, c1, p1);
+  Bs_plain = test_bezier_plain(s, p0, c0, c1, p1);
+
+  assert_true(glm_eq(Bs,  Bs_plain));
+  assert_true(glm_eq(smc, Bs_plain));
+  assert_true(glm_eq(Bs,  smc));
+}
diff --git a/test/src/test_main.c b/test/src/test_main.c
index ff77b02..8ce1673 100644
--- a/test/src/test_main.c
+++ b/test/src/test_main.c
@@ -38,7 +38,10 @@ main(int argc, const char * argv[]) {
     cmocka_unit_test(test_vec3),
 
     /* affine */
-    cmocka_unit_test(test_affine)
+    cmocka_unit_test(test_affine),
+
+    /* bezier */
+    cmocka_unit_test(test_bezier)
   };
 
   return cmocka_run_group_tests(tests, NULL, NULL);
diff --git a/test/src/test_tests.h b/test/src/test_tests.h
index 7b9cf0a..618cc9f 100644
--- a/test/src/test_tests.h
+++ b/test/src/test_tests.h
@@ -40,4 +40,7 @@ test_vec3(void **state);
 void
 test_affine(void **state);
 
+void
+test_bezier(void **state);
+
 #endif /* test_tests_h */