diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..914110b
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,3 @@
+[build-system]
+requires = ["setuptools", "wheel", "cython>=3.0.3"]
+build-backend = "setuptools.build_meta"
\ No newline at end of file
diff --git a/setup.py b/setup.py
index a8f7da2..c2d98f0 100644
--- a/setup.py
+++ b/setup.py
@@ -2,17 +2,32 @@
 import io
 from distutils.core import setup
 
+# deciding to use cython or not based on its successful import
+USE_CYTHON : bool = True
+try:
+    from Cython.Build import cythonize
+except ModuleNotFoundError:
+    print('unable to retrieve cython, will resort to pure python installation')
+    USE_CYTHON = False
+
 # load the version from version.py
 version = {}
 with open("similaritymeasures/version.py") as fp:
     exec(fp.read(), version)
 
+# based on the option of using cython or not the ext_modules variable is set
+if USE_CYTHON:
+    ext_modules = cythonize('similaritymeasures/similaritymeasures.py', language='c++')
+else:
+    ext_modules = []
+
 setup(
     name='similaritymeasures',
     version=version["__version__"],
     author='Charles Jekel',
     author_email='cjekel@gmail.com',
     packages=['similaritymeasures'],
+    ext_modules=ext_modules,
     url='https://github.com/cjekel/similarity_measures',
     license='MIT License',
     description='Quantify the difference between two arbitrary curves in space',  # noqa E501
@@ -22,5 +37,6 @@
     install_requires=[
         "numpy >= 1.14.0",
         "scipy >= 0.19.0",
+        "cython >= 3.0.3"
     ],
 )
diff --git a/similaritymeasures/similaritymeasures.py b/similaritymeasures/similaritymeasures.py
index b03ff7f..f9e0794 100644
--- a/similaritymeasures/similaritymeasures.py
+++ b/similaritymeasures/similaritymeasures.py
@@ -1,6 +1,8 @@
 from __future__ import division
 import numpy as np
 from scipy.spatial import distance
+
+import cython
 # MIT License
 #
 # Copyright (c) 2018,2019 Charles Jekel
@@ -24,6 +26,8 @@
 # SOFTWARE.
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def poly_area(x, y):
     r"""
     A function that computes the polygonal area via the shoelace formula.
@@ -57,6 +61,8 @@ def poly_area(x, y):
     return 0.5*np.abs(np.dot(x, np.roll(y, 1))-np.dot(y, np.roll(x, 1)))
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def is_simple_quad(ab, bc, cd, da):
     r"""
     Returns True if a quadrilateral is simple
@@ -98,6 +104,8 @@ def is_simple_quad(ab, bc, cd, da):
     return sum(crossTF) > 2
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def makeQuad(x, y):
     r"""
     Calculate the area from the x and y locations of a quadrilateral
@@ -165,6 +173,8 @@ def makeQuad(x, y):
     return area
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def get_arc_length(dataset):
     r"""
     Obtain arc length distances between every point in 2-D space
@@ -203,6 +213,8 @@ def get_arc_length(dataset):
     return arcLength, arcLengths
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def area_between_two_curves(exp_data, num_data):
     r"""
     Calculates the area between two curves.
@@ -252,6 +264,9 @@ def area_between_two_curves(exp_data, num_data):
     #
     # then you can calculate the area as
     # area = area_between_two_curves(exp_data, num_data)
+    i = cython.declare(cython.int)
+    n_exp = cython.declare(cython.int)
+    n_num = cython.declare(cython.int)
 
     n_exp = len(exp_data)
     n_num = len(num_data)
@@ -300,6 +315,8 @@ def area_between_two_curves(exp_data, num_data):
     return np.sum(area)
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def get_length(x, y, norm_seg_length=True):
     r"""
     Compute arc lengths of an x y curve.
@@ -331,6 +348,8 @@ def get_length(x, y, norm_seg_length=True):
     >>> le, le_total, le_cum = get_length(x, y)
 
     """
+    i = cython.declare(cython.int)
+    n = cython.declare(cython.int)
     n = len(x)
 
     if norm_seg_length:
@@ -357,6 +376,8 @@ def get_length(x, y, norm_seg_length=True):
     return le, np.sum(le), l_sum
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def curve_length_measure(exp_data, num_data):
     r"""
     Compute the curve length based distance between two curves.
@@ -430,6 +451,8 @@ def curve_length_measure(exp_data, num_data):
     return np.sqrt(np.sum(r_sq))
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def frechet_dist(exp_data, num_data, p=2):
     r"""
     Compute the discrete Frechet distance
@@ -490,6 +513,10 @@ def frechet_dist(exp_data, num_data, p=2):
     >>> df = frechet_dist(exp_data, num_data)
 
     """
+    i = cython.declare(cython.int)
+    j = cython.declare(cython.int)
+    n = cython.declare(cython.int)
+    m = cython.declare(cython.int)
     n = len(exp_data)
     m = len(num_data)
     c = distance.cdist(exp_data, num_data, metric='minkowski', p=p)
@@ -507,6 +534,8 @@ def frechet_dist(exp_data, num_data, p=2):
     return ca[n-1, m-1]
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def normalizeTwoCurves(x, y, w, z):
     """
     Normalize two curves for PCM method.
@@ -561,6 +590,8 @@ def normalizeTwoCurves(x, y, w, z):
     return xi, eta, xiP, etaP
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def pcm(exp_data, num_data, norm_seg_length=False):
     """
     Compute the Partial Curve Mapping area.
@@ -677,6 +708,8 @@ def pcm(exp_data, num_data, norm_seg_length=False):
     return np.min(pcm_dists)
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def dtw(exp_data, num_data, metric='euclidean', **kwargs):
     r"""
     Compute the Dynamic Time Warping distance.
@@ -782,6 +815,8 @@ def dtw(exp_data, num_data, metric='euclidean', **kwargs):
     >>> r, d = dtw(exp_data, num_data, metric='cityblock')
 
     """
+    i = cython.declare(cython.int)
+    j = cython.declare(cython.int)
     c = distance.cdist(exp_data, num_data, metric=metric, **kwargs)
 
     d = np.zeros(c.shape)
@@ -797,6 +832,8 @@ def dtw(exp_data, num_data, metric='euclidean', **kwargs):
     return d[-1, -1], d
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def dtw_path(d):
     r"""
     Calculates the optimal DTW path from a given DTW cumulative distance
@@ -861,6 +898,8 @@ def dtw_path(d):
     >>> plt.show()
 
     """
+    i = cython.declare(cython.int)
+    j = cython.declare(cython.int)
     path = []
     i, j = d.shape
     i = i - 1
@@ -888,6 +927,8 @@ def dtw_path(d):
     return path[::-1]
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def mae(exp_data, num_data):
     """
     Compute the Mean Absolute Error (MAE).
@@ -913,6 +954,8 @@ def mae(exp_data, num_data):
     return np.mean(c)
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
 def mse(exp_data, num_data):
     """
     Compute the Mean Squared Error (MAE).