diff --git a/curvepy/bezier_curve.py b/curvepy/bezier_curve.py
index f28f97295930e97784d860b68651a2aa68ac3888..f8ea432eb2264fcc28601d0a160e01ec077dcd79 100644
--- a/curvepy/bezier_curve.py
+++ b/curvepy/bezier_curve.py
@@ -1,115 +1,19 @@
+from __future__ import annotations  # Needed until Py3.10, see PEP 563
 import numpy as np
 import sympy as sy
+import math
+import scipy.special as scs
 import matplotlib.pyplot as plt
-import threading as th
 import shapely.geometry as sg
 from abc import ABC, abstractmethod
-from typing import Tuple, Callable, Union
+from typing import List, Tuple, Callable, Union
+from functools import partial, cached_property
+import concurrent.futures
+from multiprocessing import cpu_count
 
-from curvepy.utilities import csv_read
-
-
-class Tholder:
-    """
-    Class holds Array with equidistant ts in [0,1] of length n
-
-    Parameters
-    ----------
-    n: int
-        Numbers of ts to calculate
-
-    Attributes
-    -------
-    _tArray : np.ndarray
-        array in which all ts are stored
-    _pointer : int
-        pointer pointing on next t to get
-    lockMe: th.Lock
-        lock for multithreading
-    """
-
-    def __init__(self, ts: np.ndarray = None) -> None:
-        self._tArray = ts
-        self._pointer = 0
-        self.lockMe = th.Lock()  # variable used to control access of threads
-
-    def get_next_t(self) -> float:
-        """
-        Method for threads to get next t
-
-        Returns
-        -------
-        float:
-            t to calculate the next De Casteljau step
-        """
-        if self._pointer == len(self._tArray):
-            return -1
-        res = self._tArray[self._pointer]
-        self._pointer += 1
-        return res
-
-
-class CasteljauThread(th.Thread):
-    """
-    Thread class computing the De Casteljau algorithm
-
-    Parameters
-    ----------
-    ts_holder: Tholder
-        Class which yields all ts
-    c: np.ndarray
-        Array with control points
-    f: Function
-         Function to transform t if necessary
-
-    Attributes
-    -------
-    _ts_holder : Tholder
-        instance of class Tholder so thread can get the ts for calculating the de Casteljau algorithm
-    _coords: np.ndarray
-        original control points
-    res: list
-        actual points on curve
-    _func: Function
-        function for transforming t
-    """
-
-    def __init__(self, ts_holder: Tholder, c: np.ndarray, f: Callable[[float], float] = lambda x: x) -> None:
-        th.Thread.__init__(self)
-        self._ts_holder = ts_holder
-        self._coords = c
-        self.res = []
-        self._func = f
-
-    def de_caes(self, t: float, n: int) -> None:
-        """
-        Method implementing the the De Casteljau algorithm
-
-        Parameters
-        ----------
-        t: float
-            value at which to be evaluated at
-        n: int
-            number of iterations TODO find more describing phrasing together
-        """
-        m = self._coords.copy()
-        t = self._func(t)
-        for r in range(n):
-            m[:, :(n - r - 1)] = (1 - t) * m[:, :(n - r - 1)] + t * m[:, 1:(n - r)]
-        self.res.append(m[:, 0])
-
-    def run(self) -> None:
-        """
-        Method calculates points until depot is empty
-        """
-        _, n = self._coords.shape
-        while True:
-            self._ts_holder.lockMe.acquire()
-            t = self._ts_holder.get_next_t()
-            self._ts_holder.lockMe.release()
-            if t == -1:
-                break
-            self.de_caes(t, n)
+from curvepy.de_caes import de_caes, subdivision
+from curvepy.utilities import min_max_box
+from curvepy.types import bernstein_polynomial
 
 
 class AbstractBezierCurve(ABC):
@@ -134,30 +38,23 @@ class AbstractBezierCurve(ABC):
         numbers of equidistant ts to calculate
     func: Callable
         function computing the Bezier Curve for a single point
-    _curve: list
-        list containing points belonging to actual curve
-    box: list
-        points describing minmax box of curve
     """
 
-    def __init__(self, m: np.ndarray, cnt_ts: int = 1000) -> None:
+    def __init__(self, m: np.ndarray, cnt_ts: int = 1000, use_parallel: bool = False,
+                 interval: Tuple[int, int] = (0, 1)) -> None:
         self._bezier_points = m
         self._dimension = self._bezier_points.shape[0]
         self._cnt_ts = cnt_ts
-        self.func = self.init_func(m)
-        self._curve = None
-        self.box = []
+        self.interval = interval
+        self.func = self.init_func()
+        self._use_parallel = use_parallel
 
     @abstractmethod
-    def init_func(self, m: np.ndarray) -> Callable:
+    def init_func(self) -> Callable[[float], np.ndarray]:
         """
+        # TODO: Dies ist eine absolut dreiste Luege!
         Method returns the function to calculate all values at once.
 
-        Parameters
-        ----------
-        m: np.ndarray:
-            array containing the Bezier Points
-
         Returns
         -------
         Callable:
@@ -165,23 +62,15 @@ class AbstractBezierCurve(ABC):
         """
         ...
 
-    def _get_or_create_values(self) -> np.ndarray:
-        """
-        Method returns minmax box of calculated curve
+    def parallel_execution(self, ts: np.ndarray):
+        with concurrent.futures.ThreadPoolExecutor(max_workers=cpu_count() * 2) as executor:
+            return executor.map(self.func, ts)
 
-        Returns
-        -------
-        np.ndarray:
-            calculating no. of points defined by variable cnt_ts
-        """
-        if self._curve is None:
-            ts = np.linspace(0, 1, self._cnt_ts)
-            self._curve = self.func(ts)
-            self.min_max_box()
-        return self._curve
+    def serial_execution(self, ts: np.ndarray):
+        return np.frompyfunc(self.func, 1, 1)(ts)
 
-    @property
-    def curve(self) -> Union[Tuple[list, list], Tuple[list, list, list]]:
+    @cached_property
+    def curve(self) -> Union[Tuple[List[float], List[float]], Tuple[List[float], List[float], List[float]]]:
         """
         Method returning coordinates of all calculated points
 
@@ -190,13 +79,16 @@ class AbstractBezierCurve(ABC):
         Union[tuple[list, list], tuple[list, list, list]]:
             first list for x coords, second for y coords and third for z if existing
         """
-        tmp = np.ravel([*self._get_or_create_values()])
+
+        ts = np.linspace(0, 1, self._cnt_ts)
+        curve = self.parallel_execution(ts) if self._use_parallel else self.serial_execution(ts)
+        tmp = np.ravel([*curve])
         if self._dimension == 2:
             return tmp[0::2], tmp[1::2]
         return tmp[0::3], tmp[1::3], tmp[2::3]
 
     @staticmethod
-    def intersect(t1: tuple, t2: tuple) -> bool:
+    def intersect(t1: np.ndarray, t2: np.ndarray) -> bool:
         """
         Method checking intersection of two given tuples TODO find more describing phrasing together
 
@@ -213,23 +105,18 @@ class AbstractBezierCurve(ABC):
             true if intersect otherwise false
         """
         return t2[0] <= t1[0] <= t2[1] \
-            or t2[0] <= t1[1] <= t2[1] \
-            or t1[0] <= t2[0] <= t1[1] \
-            or t1[0] <= t2[1] <= t1[1]
+               or t2[0] <= t1[1] <= t2[1] \
+               or t1[0] <= t2[0] <= t1[1] \
+               or t1[0] <= t2[1] <= t1[1]
 
-    def min_max_box(self) -> None:
+    @cached_property
+    def min_max_box(self) -> List[np.ndarray]:
         """
         Method creates minmax box for the corresponding curve
         """
-        xs = sorted([*self._bezier_points[0, :]])
-        ys = sorted([*self._bezier_points[1, :]])
-        if self._dimension == 2:
-            self.box = [(xs[0], xs[-1]), (ys[0], ys[-1])]
-            return
-        zs = sorted([*self._bezier_points[2, :]])
-        self.box.append((zs[0], zs[-1]))
+        return min_max_box(self._bezier_points)
 
-    def collision_check(self, other_curve) -> bool:
+    def collision_check(self, other_curve: BezierCurve) -> bool:
         """
         Method checking collision with given curve.
         Starts with a box check, if this didn't intersect it is checking the actual curves
@@ -249,7 +136,7 @@ class AbstractBezierCurve(ABC):
 
         return self.curve_collision_check(other_curve)
 
-    def box_collision_check(self, other_curve) -> bool:
+    def box_collision_check(self, other_curve: BezierCurve) -> bool:
         """
         Method checking box collision with the given curve.
 
@@ -263,15 +150,15 @@ class AbstractBezierCurve(ABC):
         bool:
             true if boxes collide otherwise false
         """
-        o_box = other_curve.box
-        box = self.box
+        o_box = other_curve.min_max_box
+        box = self.min_max_box
         for t1, t2 in zip(box, o_box):
             if not self.intersect(t1, t2):
                 return False
 
         return True
 
-    def curve_collision_check(self, other_curve) -> bool:
+    def curve_collision_check(self, other_curve: BezierCurve) -> bool:
         """
         Method checking curve collision with the given curve.
 
@@ -327,6 +214,130 @@ class AbstractBezierCurve(ABC):
             c.plot()
         plt.show()
 
+    def single_forward_difference(self, i: int = 0, r: int = 0) -> np.ndarray:
+        """
+        Method using equation 5.23 to calculate forward difference of degree r for specific point i
+
+        Parameters
+        ----------
+        i: int:
+            point i for which forward all_forward_differences are calculated
+
+        r: int:
+            degree of forward difference
+
+        Returns
+        -------
+        np.ndarray:
+                forward difference of degree r for point i
+
+        Notes
+        -----
+        Equation used for computing all_forward_differences:
+        math:: \\Delta^r b_i = \\sum_{j=0}^r \\binom{r}{j} (-1)^{r-j} b_{i+j}
+        """
+        return np.sum([scs.binom(r, j) * (-1) ** (r - j) * self._bezier_points[:, i + j] for j in range(0, r + 1)],
+                      axis=0)
+
+    def all_forward_differences(self, i: int = 0) -> np.ndarray:
+        """
+        Method using equation 5.23 to calculate all forward all_forward_differences for a given point i.
+        First entry is first difference, second entry is second difference and so on.
+
+        Parameters
+        ----------
+        i: int:
+            point i for which forward all_forward_differences are calculated
+
+        Returns
+        -------
+        np.ndarray:
+             array holds all forward all_forward_differences for given point i
+
+        Notes
+        -----
+        Equation used for computing all_forward_differences:
+        math:: \\Delta^r b_i = \\sum_{j=0}^r \\binom{r}{j} (-1)^{r-j} b_{i+j}
+        """
+        _, n = self._bezier_points.shape
+        diff = [self.single_forward_difference(i, r) for r in range(0, n)]
+        return np.array(diff).T
+
+    def derivative_bezier_curve(self, t: float = 1, r: int = 1) -> np.ndarray:
+        """
+        Method using equation 5.24 to calculate rth derivative of bezier curve at value t
+
+        Parameters
+        ----------
+        t: float:
+            value for which Bezier curves are calculated
+
+        r: int:
+            rth derivative
+
+        Returns
+        -------
+        np.ndarray:
+             point of the rth derivative at value t
+
+        Notes
+        -----
+        Equation used for computing:
+        math:: \\frac{d^r}{dt^r} b^n(t) = \\frac{n!}{(n-r)!} \\cdot \\sum_{j=0}^{n-r} \\Delta^r b_j \\cdot B_j^{n-r}(t)
+        """
+        _, n = self._bezier_points.shape
+        factor = scs.factorial(n) / scs.factorial(n - r)
+        tmp = [factor * self.single_forward_difference(j, r) * bernstein_polynomial(n - r, j, t) for j in range(n - r)]
+        return np.sum(tmp, axis=0)
+
+    @staticmethod
+    def barycentric_combination_bezier(m: AbstractBezierCurve, c: AbstractBezierCurve, alpha: float = 0,
+                                       beta: float = 1) -> AbstractBezierCurve:
+        """
+        TODO Docstrings
+        Method using 5.13 to calculate the barycentric combination of two given bezier curves
+
+        Parameters
+        ----------
+        m: np.ndarray:
+            first array containing the Bezier Points
+
+        c: np.ndarray:
+            second array containing the Bezier Points
+
+        alpha: float:
+            weight for the first Bezier curve
+
+        beta: float:
+            weight for the first Bezier curve
+
+        t: float:
+            value for which Bezier curves are calculated
+
+        r: int:
+            optional Parameter to calculate only a partial curve if we already have some degree of the bezier points
+
+        interval: Tuple[float,float]:
+            Interval of t used for affine transformation
+
+        Returns
+        -------
+        np.ndarray:
+                point of the weighted combination
+
+        Notes
+        -----
+        The Parameter alpha and beta must hold the following condition: alpha + beta = 1
+        Equation used for computing:
+        math:: \\sum_{j=0}^r (\\alpha \\cdot b_j + \\beta \\cdot c_j)B_j^n(t) =
+        \\alpha \\cdot \\sum_{j=0}^r b_j \\cdot B_j^n(t) + \\beta \\cdot \\sum_{j=0}^r c_j \\cdot B_j^n(t)
+        """
+
+        if alpha + beta != 1:
+            raise ValueError("Alpha and Beta must add up to 1!")
+
+        return m * alpha + c * beta
+
     def __str__(self):
         """
         Returns string represenation as the mathematical bezier curve
@@ -347,6 +358,29 @@ class AbstractBezierCurve(ABC):
         """
         return f"<id {id(self)}, {self.__str__()}>"
 
+    def __call__(self, u):
+        # 3.9
+        a, b = self.interval
+        return self.func((u-a)/(b-a))
+
+    def __mul__(self, other: Union[float, int]):
+        del self.curve  # this is fine and as it should be to reset cached_properties, linters are stupid
+        self._bezier_points *= other
+        return self
+
+    __rmul__ = __mul__
+
+    # TODO write in docstrings that anything other than unit intervals are forbidden
+    def __add__(self, other: AbstractBezierCurve):
+        if not isinstance(other, AbstractBezierCurve):
+            raise TypeError("")
+
+        del self.curve  # this is fine
+
+        self._bezier_points += other._bezier_points
+        self._cnt_ts = max(self._cnt_ts, other._cnt_ts)
+        return self
+
 
 class BezierCurve(AbstractBezierCurve):
     """
@@ -361,32 +395,28 @@ class BezierCurve(AbstractBezierCurve):
     see AbstractBezierCurve
     """
 
-    def init_func(self, m: np.ndarray) -> Callable:
+    def init_func(self) -> Callable[[float], np.ndarray]:
         """
         Method returns the function to calculate all values at once.
 
-        Parameters
-        ----------
-        m: np.ndarray:
-            array containing the Bezier Points
-
         Returns
         -------
         Callable:
             function representing the Bezier Curve
         """
+        m = self._bezier_points
         _, n = m.shape
         m = sy.Matrix(m)
         t = sy.symbols('t')
         for r in range(n):
             m[:, :(n - r - 1)] = (1 - t) * m[:, :(n - r - 1)] + t * m[:, 1:(n - r)]
-        f = sy.lambdify(t, m[:, 0])
-        return np.frompyfunc(f, 1, 1)
+        return sy.lambdify(t, m[:, 0])
 
 
-class BezierCurveThreaded(AbstractBezierCurve):
+# TODO Interval
+class BezierCurveDeCaes(AbstractBezierCurve):
     """
-    Class for creating a 2-dimensional Bezier Curve by using the threaded De Casteljau Algorithm
+    Class for creating a 2-dimensional Bezier Curve by using the De Casteljau Algorithm
 
     Parameters
     ----------
@@ -397,61 +427,152 @@ class BezierCurveThreaded(AbstractBezierCurve):
     see AbstractBezierCurve
     """
 
-    def init_func(self, m: np.ndarray) -> Callable:
+    def init_func(self) -> Callable[[float], np.ndarray]:
         """
         Method returns the function to calculate all values at once.
 
+        Returns
+        -------
+        Callable:
+            function representing the Bezier Curve
+        """
+        return partial(de_caes, self._bezier_points)
+
+
+class BezierCurveBernstein(AbstractBezierCurve):
+
+    def init_func(self) -> Callable[[float], np.ndarray]:
+        return self.bezier_curve_with_bernstein
+
+    def bezier_curve_with_bernstein(self, t: float = 0.5, r: int = 0,
+                                    interval: Tuple[float, float] = (0, 1)) -> np.ndarray:
+        """
+        Method using 5.8 to calculate a point with given bezier points
+
         Parameters
         ----------
-        m: np.ndarray:
-            array containing the Bezier Points
+        t: float:
+            value for which Bezier curve are calculated
+
+        r: int:
+            optional Parameter to calculate only a partial curve if we already have some degree of the bezier points
+
+        interval: Tuple[float,float]:
+            Interval of t used for affine transformation
 
         Returns
         -------
-        Callable:
-            function representing the Bezier Curve
+        np.ndarray:
+                point on the curve
+
+        Notes
+        -----
+        Equation used for computing:
+        math:: b_i^r(t) = \\sum_{j=0}^r b_{i+j} \\cdot B_i^r(t)
         """
-        return self.de_casteljau_threading
+        m = self._bezier_points
+        _, n = m.shape
+        t = (t - interval[0]) / (interval[1] - interval[0])
+        return np.sum([m[:, i] * bernstein_polynomial(n - r - 1, i, t) for i in range(n - r)], axis=0)
+
+
+class HornerBezierCurve(AbstractBezierCurve):
 
-    def de_casteljau_threading(self, ts: np.ndarray = None, cnt_threads: int = 4) -> None:
+    def init_func(self) -> Callable[[float], np.ndarray]:
+        return self.horn_bez
+
+    def horn_bez(self, t: float = 0.5) -> np.ndarray:
         """
-        Method implementing the threading for the De Casteljau algorithm
+        Method using horner like scheme to calculate point on curve with given t
 
         Parameters
         ----------
-        cnt_threads: int
-            number of threads to use
+        t: float:
+            value for which point is calculated
+
+        Returns
+        -------
+        np.ndarray:
+            point calculated with given t
+        """
+        m = self._bezier_points
+        n = m.shape[1] - 1  # need degree of curve (n points means degree = n-1)
+        res = m[:, 0] * (1 - t)
+        for i in range(1, n):
+            res = (res + t ** i * scs.binom(n, i) * m[:, i]) * (1 - t)
+
+        res += t ** n * m[:, n]
+
+        return res
+
 
-        ts: np.ndarray:
-            array containing all ts used to calculate points
+class MonomialBezierCurve(AbstractBezierCurve):
 
+    def init_func(self) -> Callable[[float], np.ndarray]:
         """
-        ts = Tholder(ts)
-        threads = []
-        curve = []
+        TODO Docstring
+        Method calculating monomial representation of given bezier form using 5.27
 
-        for _ in range(cnt_threads):
-            threads.append(CasteljauThread(ts, self._bezier_points))
+        Parameters
+        ----------
+        m: np.ndarray:
+            array containing the Bezier Points
 
-        for t in threads:
-            t.start()
+        Returns
+        -------
+        Callable:
+            bezier function in polynomial form
+
+        Notes
+        -----
+        Equation 5.27 used for computing polynomial form:
+        math:: b^n(t) = \\sum_{j=0}^n \\binom{n}{j} \\Delta^j b_0 t^j
+
+        Initially the method would only compute the polynomial coefficients in an array, and parsing this array with a given
+        t to the horner method we would get a point back. Instead the method uses sympy to calculate a function depending on
+        t. After initial computation, f(t) calculates the value for a given t. Having a function it is simple to map it on
+        an array containing multiple values for t.
+        As a result we do not need to call the horner method for each t individually.
+        """
+        m = self._bezier_points
+        _, n = m.shape
+        diff = self.all_forward_differences()
+        t = sy.symbols('t')
+        res = 0
+        for i in range(n):
+            res += scs.binom(n - 1, i) * diff[:, i] * t ** i
 
-        for t in threads:
-            t.join()
-            tmp = t.res
-            curve = curve + tmp
+        return sy.lambdify(t, res)
 
-        return curve
 
+class BezierCurveApproximation(AbstractBezierCurve):
+    def __init__(self, m: np.ndarray, approx_rounds: int = 5, interval: Tuple[int, int] = (0, 1)) -> None:
+        """
+        We can't actually set the approx_rounds explicitly, because the number of iterations can be
+        changed by using __add__ with any other AbstractBezierCurve
+        """
+        AbstractBezierCurve.__init__(self, m, 2**approx_rounds*m.shape[1], False, interval)
 
-def init() -> None:
-    m = csv_read('test.csv')  # reads csv file with bezier points
-    b1 = BezierCurve(m)
-    print(b1.func(0.5))
-    m = csv_read('test2.csv')  # reads csv file with bezier points
-    b2 = BezierCurve(m)
-    b2.show_funcs([b1])
+    def init_func(self) -> Callable[[float], np.ndarray]:
+        # dummy
+        return partial(de_caes, self._bezier_points)
 
+    serial_execution = None
+    parallel_execution = None
 
-if __name__ == "__main__":
-    init()
+    @cached_property
+    def curve(self) -> Union[Tuple[List[float], List[float]], Tuple[List[float], List[float], List[float]]]:
+        approx_rounds = math.ceil(math.log((self._cnt_ts/self._bezier_points.shape[1]), 2))
+        current = [self._bezier_points]
+        for _ in range(approx_rounds):
+            queue = []
+            for c in current:
+                queue.extend(subdivision(c))
+            current = queue
+
+        ret = np.hstack(current)
+        ret = np.ravel(ret)
+        n = len(ret)
+        if self._dimension == 2:
+            return ret[:n/2], ret[n/2:]
+        return ret[:n/3], ret[n/3:2*n/3], ret[2*n/3:]
diff --git a/curvepy/chapter_5.py b/curvepy/chapter_5.py
deleted file mode 100644
index 991524665e2eacee4fdbf50c07e6c76cd67d27af..0000000000000000000000000000000000000000
--- a/curvepy/chapter_5.py
+++ /dev/null
@@ -1,731 +0,0 @@
-import sys as s
-import numpy as np
-import sympy as sy
-import scipy.special as scs
-from functools import reduce, partial
-from typing import Tuple, Callable, Union, List, Any
-
-
-def bernstein_polynomial_rec(n: int, i: int, t: float = 1) -> float:
-    """
-    Method using 5.8 to calculate a point with given bezier points
-
-    Parameters
-    ----------
-    n: int:
-        degree of the Bernstein Polynomials
-
-    i: int:
-        starting point for calculation
-
-    t: float:
-        value for which Bezier curve are calculated
-
-    Returns
-    -------
-    float:
-        value of Bernstein Polynomial B_i^n(t)
-
-    Notes
-    -----
-    Equation used for computing:
-    Base Case: B_0^0(t) = 1
-    math:: i \\notin \\{0, \\dots, n\\} \\rightarrow B_i^n(t) = 0
-    math:: B_i^n(t) = (1-t) \\cdot B_i^{n-1}(t) + t \\cdot B_{i-1}^{n-1}(t)
-    """
-    if i == n == 0:
-        return 1
-    if not 0 <= i <= n:
-        return 0
-    return (1 - t) * bernstein_polynomial_rec(n - 1, i, t) + t * bernstein_polynomial_rec(n - 1, i - 1, t)
-
-
-def bernstein_polynomial(n: int, i: int, t: float = 1) -> float:
-    """
-    Method using 5.1 to calculate a point with given bezier points
-
-    Parameters
-    ----------
-    n: int:
-        degree of the Bernstein Polynomials
-
-    i: int:
-        starting point for calculation
-
-    t: float:
-        value for which Bezier curve are calculated
-
-    Returns
-    -------
-    float:
-        value of Bernstein Polynomial B_i^n(t)
-
-    Notes
-    -----
-    Equation used for computing:
-    math:: B_i^n(t) = \\binom{n}{i} t^{i} (1-t)^{n-i}
-    """
-    return scs.binom(n, i) * (t ** i) * ((1 - t) ** (n - i))
-
-
-def partial_bernstein_polynomial(n: int, i: int) -> Callable[[float], float]:
-    """
-    Method using 5.1 to calculate a point with given bezier points
-
-    Parameters
-    ----------
-    n: int:
-        degree of the Bernstein Polynomials
-
-    i: int:
-        starting point for calculation
-
-    Returns
-    -------
-    Callable[[float], float]:
-        partial application of 5.1
-
-    Notes
-    -----
-    Equation used for computing:
-    math:: B_i^n(t) = \\binom{n}{i} t^{i} (1-t)^{n-i}
-    """
-    return partial(bernstein_polynomial, n, i)
-
-
-def bezier_curve_with_bernstein(m: np.ndarray, t: float = 0.5, r: int = 0,
-                                interval: Tuple[float, float] = (0, 1)) -> np.ndarray:
-    """
-    Method using 5.8 to calculate a point with given betier points
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        array containing the Bezier Points
-
-    t: float:
-        value for which Bezier curve are calculated
-
-    r: int:
-        optional Parameter to calculate only a partial curve if we already have some degree of the bezier points
-
-    interval: Tuple[float,float]:
-        Interval of t used for affine transformation
-
-    Returns
-    -------
-    np.ndarray:
-            point on the curve
-
-    Notes
-    -----
-    Equation used for computing:
-    math:: b_i^r(t) = \\sum_{j=0}^r b_{i+j} \\cdot B_i^r(t)
-    """
-    _, n = m.shape
-    t = (t - interval[0]) / (interval[1] - interval[0])
-    return np.sum([m[:, i] * bernstein_polynomial(n - r - 1, i, t) for i in range(n - r)], axis=0)
-
-
-def intermediate_bezier_points(m: np.ndarray, r: int, i: int, t: float = 1,
-                               interval: Tuple[float, float] = (0, 1)) -> np.ndarray:
-    """
-    Method using 5.7 an intermediate point of the bezier curve
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        array containing the Bezier Points
-
-    i: int:
-        which intermediate points should be calculated
-
-    t: float:
-        value for which Bezier curve are calculated
-
-    r: int:
-        optional Parameter to calculate only a partial curve
-
-    interval: Tuple[float,float]:
-        Interval of t used for affine transformation
-
-    Returns
-    -------
-    np.ndarray:
-            intermediate point
-
-    Notes
-    -----
-    Equation used for computing:
-    math:: b_i^r(t) = \\sum_{j=0}^r b_{i+j} \\cdot B_i^r(t)
-    """
-    _, n = m.shape
-    t = (t - interval[0]) / (interval[1] - interval[0])
-    return np.sum([m[:, i + j] * bernstein_polynomial(n - 1, j, t) for j in range(r)], axis=0)
-
-
-def barycentric_combination_bezier(m: np.ndarray, c: np.ndarray, alpha: float = 0, beta: float = 1, t: float = 1,
-                                   r: int = 0, interval: Tuple[float, float] = (0, 1)) -> np.ndarray:
-    """
-    Method using 5.13 to calculate the barycentric combination of two given bezier curves
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        first array containing the Bezier Points
-
-    c: np.ndarray:
-        second array containing the Bezier Points
-
-    alpha: float:
-        weight for the first Bezier curve
-
-    beta: float:
-        weight for the first Bezier curve
-
-    t: float:
-        value for which Bezier curves are calculated
-
-    r: int:
-        optional Parameter to calculate only a partial curve if we already have some degree of the bezier points
-
-    interval: Tuple[float,float]:
-        Interval of t used for affine transformation
-
-    Returns
-    -------
-    np.ndarray:
-            point of the weighted combination
-
-    Notes
-    -----
-    The Parameter alpha and beta must hold the following condition: alpha + beta = 1
-    Equation used for computing:
-    math:: \\sum_{j=0}^r (\\alpha \\cdot b_j + \\beta \\cdot c_j)B_j^n(t) =
-    \\alpha \\cdot \\sum_{j=0}^r b_j \\cdot B_j^n(t) + \\beta \\cdot \\sum_{j=0}^r c_j \\cdot B_j^n(t)
-    """
-
-    if alpha + beta != 1:
-        raise Exception("Alpha and Beta must add up to 1!")
-
-    return alpha * bezier_curve_with_bernstein(m, t, r, interval) + beta * bezier_curve_with_bernstein(c, t, r,
-                                                                                                       interval)
-
-
-def horn_bez(m: np.ndarray, t: float = 0.5) -> np.ndarray:
-    """
-    Method using horner like scheme to calculate point on curve with given t
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        array containing the Bezier Points
-
-    t: float:
-        value for which point is calculated
-
-    Returns
-    -------
-    np.ndarray:
-        point calculated with given t
-    """
-    n = m.shape[1] - 1  # need degree of curve (n points means degree = n-1)
-    res = m[:, 0] * (1 - t)
-    for i in range(1, n):
-        res = (res + t ** i * scs.binom(n, i) * m[:, i]) * (1 - t)
-
-    res += t ** n * m[:, n]
-
-    return res
-
-
-def bezier_to_power(m: np.ndarray) -> Callable[[float], np.ndarray]:
-    """
-    Method calculating monomial representation of given bezier form using 5.27
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        array containing the Bezier Points
-
-    Returns
-    -------
-    Callable:
-        bezier function in polynomial form
-
-    Notes
-    -----
-    Equation 5.27 used for computing polynomial form:
-    math:: b^n(t) = \\sum_{j=0}^n \\binom{n}{j} \\Delta^j b_0 t^j
-
-    Initially the method would only compute the polynomial coefficients in an array, and parsing this array with a given
-    t to the horner method we would get a point back. Instead the method uses sympy to calculate a function depending on
-    t. After initial computation, f(t) calculates the value for a given t. Having a function it is simple to map it on
-    an array containing multiple values for t.
-    As a result we do not need to call the horner method for each t individually.
-    """
-    _, n = m.shape
-    diff = all_forward_differences(m)
-    t = sy.symbols('t')
-    res = 0
-    for i in range(n):
-        res += scs.binom(n - 1, i) * diff[:, i] * t ** i
-
-    return sy.lambdify(t, res)
-
-
-def single_forward_difference(m: np.ndarray, i: int = 0, r: int = 0) -> np.ndarray:
-    """
-    Method using equation 5.23 to calculate forward difference of degree r for specific point i
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        array containing the Bezier Points
-
-    i: int:
-        point i for which forward all_forward_differences are calculated
-
-    r: int:
-        degree of forward difference
-
-    Returns
-    -------
-    np.ndarray:
-            forward difference of degree r for point i
-
-    Notes
-    -----
-    Equation used for computing all_forward_differences:
-    math:: \\Delta^r b_i = \\sum_{j=0}^r \\binom{r}{j} (-1)^{r-j} b_{i+j}
-    """
-    _, n = m.shape
-    return np.sum([scs.binom(r, j) * (-1) ** (r - j) * m[:, i + j] for j in range(0, r + 1)], axis=0)
-
-
-def all_forward_differences(m: np.ndarray, i: int = 0) -> np.ndarray:
-    """
-    Method using equation 5.23 to calculate all forward all_forward_differences for a given point i.
-    First entry is first difference, second entry is second difference and so on.
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        array containing the Bezier Points
-
-    i: int:
-        point i for which forward all_forward_differences are calculated
-
-    Returns
-    -------
-    np.ndarray:
-         array holds all forward all_forward_differences for given point i
-
-    Notes
-    -----
-    Equation used for computing all_forward_differences:
-    math:: \\Delta^r b_i = \\sum_{j=0}^r \\binom{r}{j} (-1)^{r-j} b_{i+j}
-    """
-    _, n = m.shape
-    diff = [single_forward_difference(m, i, r) for r in range(0, n)]
-    return np.array(diff).T
-
-
-def derivative_bezier_curve(m: np.ndarray, t: float = 1, r: int = 1) -> np.ndarray:
-    """
-    Method using equation 5.24 to calculate rth derivative of bezier curve at value t
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        array containing the Bezier Points
-
-    t: float:
-        value for which Bezier curves are calculated
-
-    r: int:
-        rth derivative
-
-    Returns
-    -------
-    np.ndarray:
-         point of the rth derivative at value t
-
-    Notes
-    -----
-    Equation used for computing:
-    math:: \\frac{d^r}{dt^r} b^n(t) = \\frac{n!}{(n-r)!} \\cdot \\sum_{j=0}^{n-r} \\Delta^r b_j \\cdot B_j^{n-r}(t)
-    """
-    _, n = m.shape
-    factor = scs.factorial(n) / scs.factorial(n - r)
-    tmp = [factor * single_forward_difference(m, j, r) * bernstein_polynomial(n - r, j, t) for j in range(n - r)]
-    return np.sum(tmp, axis=0)
-
-
-def horner(m: np.ndarray, t: float = 0.5) -> Tuple[Union[float, Any], ...]:
-    """
-    Method using horner's method to calculate point with given t
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        array containing coefficients
-
-    t: float:
-        value for which point is calculated
-
-    Returns
-    -------
-    tuple:
-        point calculated with given t
-    """
-    return tuple(reduce(lambda x, y: t * x + y, m[i, ::-1]) for i in [0, 1])
-
-
-def de_caes_one_step(m: np.ndarray, t: float = 0.5, interval: Tuple[int, int] = (0, 1)) -> np.ndarray:
-    """
-    Method computing one round of de Casteljau
-
-    Parameters
-    ----------
-     m: np.ndarray:
-        array containing coefficients
-
-    t: float:
-        value for which point is calculated
-
-    interval: Tuple[int, int]:
-        if interval != (0,1) we need to transform t
-
-    Returns
-    -------
-    np.ndarray:
-        array containing calculated points with given t
-    """
-    if m.shape[1] < 2:
-        raise Exception("At least two points are needed")
-
-    t1 = (interval[1] - t) / (interval[1] - interval[0]) if interval != (0, 1) else (1 - t)
-    t2 = (t - interval[0]) / (interval[1] - interval[0]) if interval != (0, 1) else t
-
-    m[:, :-1] = t1 * m[:, :-1] + t2 * m[:, 1:]
-    return m
-
-
-def de_caes_n_steps(m: np.ndarray, t: float = 0.5, r: int = 1, interval: Tuple[int, int] = (0, 1)) -> np.ndarray:
-    """
-    Method computing r round of de Casteljau
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        array containing coefficients
-
-    t: float:
-        value for which point is calculated
-
-    r: int:
-        how many rounds of de Casteljau algorithm should be performed
-
-    interval: Tuple[int, int]:
-        if interval != (0,1) we need to transform t
-
-    Returns
-    -------
-    np.ndarray:
-        array containing calculated points with given t
-    """
-
-    if r >= m.shape[1]:
-        raise Exception("Can't perform r rounds!")
-
-    if m.shape[1] < 2:
-        raise Exception("At least two points are needed")
-
-    for i in range(r):
-        m = de_caes_one_step(m, t, interval)
-        if i != r - 1:
-            m = m[:, :-1]
-    return m
-
-
-def de_caes(m: np.ndarray, t: float = 0.5, make_copy: bool = False, interval: Tuple[int, int] = (0, 1)) -> np.ndarray:
-    """
-    Method computing de Casteljau
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        array containing coefficients
-
-    t: float:
-        value for which point is calculated
-
-    make_copy: bool:
-        optional parameter if computation should not be in place
-
-    interval: Tuple[int, int]:
-        if interval != (0,1) we need to transform t
-
-    Returns
-    -------
-    np.ndarray:
-        array containing calculated points with given t
-    """
-
-    if m.shape[1] < 2:
-        raise Exception("At least two points are needed")
-
-    if m.shape[1] < 2:
-        raise Exception("At least two points are needed")
-
-    _, n = m.shape
-    return de_caes_n_steps(m.copy(), t, n, interval) if make_copy else de_caes_n_steps(m, t, n, interval)
-
-
-def de_caes_blossom(m: np.ndarray, ts: List[float], make_copy: bool = False,
-                    interval: Tuple[int, int] = (0, 1)) -> np.ndarray:
-    """
-    Method computing de Casteljau with different values of t in each step
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        array containing coefficients
-
-    ts: List[float]:
-        List containing all ts that are used in calculation
-
-    make_copy: bool:
-        optional parameter if computation should not be in place
-
-    interval: Tuple[int, int]:
-        if interval != (0,1) we need to transform t
-
-    Returns
-    -------
-    np.ndarray:
-        array containing calculated points with given t
-    """
-
-    if m.shape[1] < 2:
-        raise Exception("At least two points are needed")
-
-    if len(ts) >= m.shape[1]:
-        raise Exception("Too many values to use!")
-
-    if not ts:
-        raise Exception("At least one element is needed!")
-
-    c = m.copy() if make_copy else m
-    for i, t in enumerate(ts):
-        c = de_caes_one_step(c, t, interval)
-        if i != len(ts):
-            c = c[:, :-1]
-    return c
-
-
-def subdiv(m: np.ndarray, t: float = 0.5) -> Tuple[np.ndarray, np.ndarray]:
-    """
-    Method using subdivison to calculate right and left polygon with given t using de Casteljau
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        array containing coefficients
-
-    t: float:
-        value for which point is calculated
-
-    Returns
-    -------
-    np.ndarray:
-        right polygon
-    np.ndarray:
-        left polygon
-    """
-    return de_caes(m, t, True), de_caes(m[:, ::-1], 1.0 - t, True)
-
-
-def distance_to_line(p1: np.ndarray, p2: np.ndarray, p_to_check: np.ndarray) -> float:
-    """
-    Method calculating distance of point to line through p1 and p2
-
-    Parameters
-    ----------
-    p1: np.ndarray:
-        beginning point of line
-
-    p2: np.ndarray:
-        end point of line
-
-    p_to_check: np.ndarray:
-        point for which distance is calculated
-
-    Returns
-    -------
-    float:
-        distance from point to line
-
-    Notes
-    -----
-    Given p1 and p2 we can check the distance p3 has to the line going through p1 and p2 as follows:
-    math:: distance(p1,p2,p3) = \\frac{|(x_1-x_1)(y_1-y_3) - (x_1-x_3)(y_2-y_1)|}{//sqrt{(x_2-x_1)^2 + (y_2-y_1)^2}}
-    more information on "https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line"
-    """
-
-    if any((p.shape[0] < 2 for p in [p1, p2, p_to_check])):
-        raise Exception("At least 2 dimensions are needed")
-
-    if p1.shape != p2.shape != p_to_check.shape:
-        raise Exception("points need to be of the same dimension!")
-
-    numerator = abs((p2[0] - p1[0]) * (p1[1] - p_to_check[1]) - (p1[0] - p_to_check[0]) * (p2[1] - p1[1]))
-    denominator = ((p2[0] - p1[0]) ** 2 + (p2[1] - p1[1]) ** 2) ** 0.5
-    return numerator / denominator
-
-
-def check_flat(m: np.ndarray, tol: float = s.float_info.epsilon) -> bool:
-    """
-    Method checking if all points between the first and the last point
-    are less than tol away from line through beginning and end point of bezier curve
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        points of curve
-
-    tol: float:
-        tolerance for distance check
-
-    Returns
-    -------
-    bool:
-        True if all point are less than tol away from line otherwise false
-    """
-    return all(distance_to_line(m[:, 0], m[:, len(m[0]) - 1], m[:, i]) <= tol for i in range(1, len(m[0]) - 1))
-
-
-def min_max_box(m: np.ndarray) -> List[float]:
-    """
-    Method creating the minmaxbox of a given curve
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        points of curve
-
-    Returns
-    -------
-    list:
-        contains the points that describe the minmaxbox
-    """
-    return [m[0, :].min(), m[0, :].max(), m[1, :].min(), m[1, :].max()]
-
-
-def intersect_lines(p1: np.ndarray, p2: np.ndarray, p3: np.ndarray, p4: np.ndarray) -> Union[np.ndarray, None]:
-    """
-    Method checking if line through p1, p2 intersects with line through p3, p4
-
-
-    Parameters
-    ----------
-    p1: np.ndarray:
-        first point of first line
-
-    p2: np.ndarray:
-        second point of first line
-
-    p3: np.ndarray:
-        first point of second line
-
-    p4: np.ndarray:
-        second point of second line
-
-    Returns
-    -------
-    bool:
-        True if all point are less than tol away from line otherwise false
-    """
-
-    if p1.shape != p2.shape != p3.shape != p4.shape:
-        raise Exception("Points need to be of the same dimension!")
-
-    # First we vertical stack the points in an array
-    vertical_stack = np.vstack([p1, p2, p3, p4])
-    # Then we transform them to homogeneous coordinates, to perform a little trick
-    homogeneous = np.hstack((vertical_stack, np.ones((4, 1))))
-    # having our points in this form we can get the lines through the cross product
-    line_1, line_2 = np.cross(homogeneous[0], homogeneous[1]), np.cross(homogeneous[2], homogeneous[3])
-    # when we calculate the cross product of the lines we get intersect point
-    x, y, z = np.cross(line_1, line_2)
-    if z == 0:
-        return None
-    # we divide with z to turn back to 2D space
-    return np.array([x / z, y / z])
-
-
-def intersect(m: np.ndarray, tol: float = s.float_info.epsilon) -> np.ndarray:
-    """
-    Method checks if curve intersects with x-axis
-
-    Parameters
-    ----------
-    m: np.ndarray:
-        points of curve
-
-    tol: float:
-        tolerance for check_flat
-
-    Returns
-    -------
-    np.ndarray:
-        Points where curve and x-axis intersect
-    """
-    box = min_max_box(m)
-    res = np.array([])
-
-    if box[2] * box[3] > 0:
-        # Both y values are positive, ergo curve lies above x_axis
-        return np.array([])
-
-    if check_flat(m, tol):
-        # poly is flat enough, so we can perform intersect of straight lines
-        # since we are assuming poly is a straight line we define a line through first and las point of poly
-        # additionally we create a line which demonstrates the x axis
-        # having these two lines we can check them for intersection
-        p = intersect_lines(m[:, 0], m[:, -1], np.array([0, 0]), np.array([1, 0]))
-        if p is not None:
-            res = np.append(res, p.reshape((2, 1)), axis=1)
-    else:
-        # if poly not flat enough we subdivide and check the resulting polygons for intersection
-        p1, p2 = subdiv(m, 0.5)
-        res = np.append(res, intersect(p1, tol).reshape((2, 1)), axis=1)
-        res = np.append(res, intersect(p2, tol).reshape((2, 1)), axis=1)
-
-    return res
-
-
-def init() -> None:
-    x = [0, 0, 8, 4]
-    y = [0, 2, 2, 0]
-
-    # x_1 = [0]
-    # y_1 = [1]
-    test = np.array([x, y], dtype=float)
-    print(test.shape)
-    # print(bernstein_polynomial(4, 2, 0.5))
-    # print(bezier_curve_with_bernstein(test))
-    _, n = test.shape
-    print(de_caes_blossom(test, [0.5, 0.5, 0.25]))
-    # print(de_caes(test, 0.5))
-    # print(bernstein_polynomial_rec.__doc__)
-    # print(min_max_box(test))
-    # print(np.ndarray([]).size)
-    # print(check_flat(test))
-    # print(horn_bez(test))
-    # print(all_forward_differences(test))
-
-
-if __name__ == "__main__":
-    init()
diff --git a/curvepy/de_caes.py b/curvepy/de_caes.py
new file mode 100644
index 0000000000000000000000000000000000000000..88189afef5294106b1a0a7feb306afd25e37db2d
--- /dev/null
+++ b/curvepy/de_caes.py
@@ -0,0 +1,167 @@
+import concurrent.futures
+from typing import List, Tuple
+import numpy as np
+from multiprocessing import cpu_count
+
+
+def de_caes_one_step(m: np.ndarray, t: float = 0.5, interval: Tuple[int, int] = (0, 1),
+                     make_copy: bool = False) -> np.ndarray:
+    """
+    Method computing one round of de Casteljau
+
+    Parameters
+    ----------
+     m: np.ndarray:
+        array containing coefficients
+
+    t: float:
+        value for which point is calculated
+
+    interval: Tuple[int, int]:
+        if interval != (0,1) we need to transform t
+
+    Returns
+    -------
+    np.ndarray:
+        array containing calculated points with given t
+    """
+    if make_copy:
+        m = m.copy()
+    l, r = interval
+    t2 = (t - l) / (r - l) if interval != (0, 1) else t
+    t1 = 1 - t2
+
+    m[:, :-1] = t1 * m[:, :-1] + t2 * m[:, 1:]
+
+    return m[:, :-1] if m.shape != (2, 1) else m
+
+
+def de_caes_n_steps(m: np.ndarray, t: float = 0.5, r: int = 1, interval: Tuple[int, int] = (0, 1)) -> np.ndarray:
+    """
+    Method computing r round of de Casteljau
+
+    Parameters
+    ----------
+    m: np.ndarray:
+        array containing coefficients
+
+    t: float:
+        value for which point is calculated
+
+    r: int:
+        how many rounds of de Casteljau algorithm should be performed
+
+    interval: Tuple[int, int]:
+        if interval != (0,1) we need to transform t
+
+    Returns
+    -------
+    np.ndarray:
+        array containing calculated points with given t
+    """
+
+    for i in range(r):
+        m = de_caes_one_step(m, t, interval)
+    return m
+
+
+def de_caes(m: np.ndarray, t: float = 0.5, make_copy: bool = False, interval: Tuple[int, int] = (0, 1)) -> np.ndarray:
+    """
+    Method computing de Casteljau
+
+    Parameters
+    ----------
+    m: np.ndarray:
+        array containing coefficients
+
+    t: float:
+        value for which point is calculated
+
+    make_copy: bool:
+        optional parameter if computation should not be in place
+
+    interval: Tuple[int, int]:
+        if interval != (0,1) we need to transform t
+
+    Returns
+    -------
+    np.ndarray:
+        array containing calculated points with given t
+    """
+
+    _, n = m.shape
+    return de_caes_n_steps(m.copy(), t, n, interval) if make_copy else de_caes_n_steps(m, t, n, interval)
+
+
+def de_caes_blossom(m: np.ndarray, ts: List[float], make_copy: bool = False,
+                    interval: Tuple[int, int] = (0, 1)) -> np.ndarray:
+    """
+    Method computing de Casteljau with different values of t in each step
+
+    Parameters
+    ----------
+    m: np.ndarray:
+        array containing coefficients
+
+    ts: List[float]:
+        List containing all ts that are used in calculation
+
+    make_copy: bool:
+        optional parameter if computation should not be in place
+
+    interval: Tuple[int, int]:
+        if interval != (0,1) we need to transform t
+
+    Returns
+    -------
+    np.ndarray:
+        array containing calculated points with given t
+    """
+
+    if m.shape[1] < 2:
+        raise Exception("At least two points are needed")
+
+    if len(ts) >= m.shape[1]:
+        raise Exception("Too many values to use!")
+
+    if not ts:
+        raise Exception("At least one element is needed!")
+
+    c = m.copy() if make_copy else m
+    for t in ts:
+        c = de_caes_one_step(c, t, interval)
+
+    return c
+
+
+def parallel_decaes_unblossomed(m: np.ndarray, ts, interval: Tuple[int, int] = (0, 1)):
+    with concurrent.futures.ThreadPoolExecutor(max_workers=cpu_count() * 2) as executor:
+        return executor.map(lambda t: de_caes(m, t, interval=interval), ts)
+
+
+def subdivision(m: np.ndarray, t: float = 0.5) -> Tuple[np.ndarray, np.ndarray]:
+    left, right = np.zeros(m.shape), np.zeros(m.shape)
+    current = m
+    for i in range(m.shape[1]):
+        left[::, i] = current[::, 0]
+        right[::, i] = current[::, -1]
+        current = de_caes_one_step(current, t, make_copy=True)
+
+    return left, right
+
+
+if __name__ == '__main__':
+    x = [0, 0, 8, 4]
+    y = [0, 2, 2, 0]
+
+    # x_1 = [0]
+    # y_1 = [1]
+    test = np.array([x, y], dtype=float)
+    ptmp = list(parallel_decaes_unblossomed(test, np.linspace(0, 1, 1000)))
+    tmp = [de_caes(test, t, make_copy=True) for t in np.linspace(0, 1, 1000)]
+
+    assert len(ptmp) == len(tmp)
+    print('Länge OK')
+    for a, b in zip(tmp, ptmp):
+        assert all(a == b)
+    print('Alle Elemente gleich')
diff --git a/curvepy/types.py b/curvepy/types.py
index 95232aa9b0acc2c339ec19ad3aa00c7ab834010b..113cf16f05d2e84522f8f1a3fc8031cd3f02c841 100644
--- a/curvepy/types.py
+++ b/curvepy/types.py
@@ -1,9 +1,11 @@
 from enum import Enum
 import numpy as np
 import sys
+import scipy.special as scs
+
 from abc import ABC, abstractmethod
 from typing import Any, Dict, Deque, List, NamedTuple, Tuple, Union, Callable
-from functools import cached_property
+from functools import cached_property, partial
 from collections.abc import Sequence
 from curvepy.utilities import create_straight_line_function
 
@@ -305,3 +307,127 @@ class Triangle:
 
 
 VoronoiRegions2D = Dict[Point2D, Deque[TriangleNode]]
+
+
+def bernstein_polynomial_rec(n: int, i: int, t: float = 1) -> float:
+    """
+    Method using 5.8 to calculate a point with given bezier points
+
+    Parameters
+    ----------
+    n: int:
+        degree of the Bernstein Polynomials
+
+    i: int:
+        starting point for calculation
+
+    t: float:
+        value for which Bezier curve are calculated
+
+    Returns
+    -------
+    float:
+        value of Bernstein Polynomial B_i^n(t)
+
+    Notes
+    -----
+    Equation used for computing:
+    Base Case: B_0^0(t) = 1
+    math:: i \\notin \\{0, \\dots, n\\} \\rightarrow B_i^n(t) = 0
+    math:: B_i^n(t) = (1-t) \\cdot B_i^{n-1}(t) + t \\cdot B_{i-1}^{n-1}(t)
+    """
+    if i == n == 0:
+        return 1
+    if not 0 <= i <= n:
+        return 0
+    return (1 - t) * bernstein_polynomial_rec(n - 1, i, t) + t * bernstein_polynomial_rec(n - 1, i - 1, t)
+
+
+def bernstein_polynomial(n: int, i: int, t: float = 1) -> float:
+    """
+    Method using 5.1 to calculate a point with given bezier points
+
+    Parameters
+    ----------
+    n: int:
+        degree of the Bernstein Polynomials
+
+    i: int:
+        starting point for calculation
+
+    t: float:
+        value for which Bezier curve are calculated
+
+    Returns
+    -------
+    float:
+        value of Bernstein Polynomial B_i^n(t)
+
+    Notes
+    -----
+    Equation used for computing:
+    math:: B_i^n(t) = \\binom{n}{i} t^{i} (1-t)^{n-i}
+    """
+    return scs.binom(n, i) * (t ** i) * ((1 - t) ** (n - i))
+
+
+def partial_bernstein_polynomial(n: int, i: int) -> Callable[[float], float]:
+    """
+    Method using 5.1 to calculate a point with given bezier points
+
+    Parameters
+    ----------
+    n: int:
+        degree of the Bernstein Polynomials
+
+    i: int:
+        starting point for calculation
+
+    Returns
+    -------
+    Callable[[float], float]:
+        partial application of 5.1
+
+    Notes
+    -----
+    Equation used for computing:
+    math:: B_i^n(t) = \\binom{n}{i} t^{i} (1-t)^{n-i}
+    """
+    return partial(bernstein_polynomial, n, i)
+
+
+def intermediate_bezier_points(m: np.ndarray, r: int, i: int, t: float = 0.5,
+                               interval: Tuple[float, float] = (0, 1)) -> np.ndarray:
+    """
+    Method using 5.7 an intermediate point of the bezier curve
+
+    Parameters
+    ----------
+    m: np.ndarray:
+        array containing the Bezier Points
+
+    i: int:
+        which intermediate points should be calculated
+
+    t: float:
+        value for which Bezier curve are calculated
+
+    r: int:
+        optional Parameter to calculate only a partial curve
+
+    interval: Tuple[float,float]:
+        Interval of t used for affine transformation
+
+    Returns
+    -------
+    np.ndarray:
+            intermediate point
+
+    Notes
+    -----
+    Equation used for computing:
+    math:: b_i^r(t) = \\sum_{j=0}^r b_{i+j} \\cdot B_i^r(t)
+    """
+    _, n = m.shape
+    t = (t - interval[0]) / (interval[1] - interval[0])
+    return np.sum([m[:, i + j] * bernstein_polynomial(n - 1, j, t) for j in range(r)], axis=0)
\ No newline at end of file
diff --git a/curvepy/ui_curvepy.py b/curvepy/ui_curvepy.py
index d3df7f5738b41a630a466e48dbc6c7caa6400b24..bd40a8cc12377fde63d176c305547a692b11b379 100644
--- a/curvepy/ui_curvepy.py
+++ b/curvepy/ui_curvepy.py
@@ -66,7 +66,7 @@ class MyOptionWidget(QWidget):
         self.plot_canvas = plot_canvas
 
         self.type = curve_type
-        self.names = ['BezierCurve', 'BezierCurveThreaded', 'BezierCurveBlossoms']
+        self.names = ['BezierCurve', 'BezierCurveDeCaes', 'BezierCurveBlossoms']
 
         self.sld_points_cnt = QSlider(Qt.Horizontal)
         self.sld_thread_cnt = QSlider(Qt.Horizontal)
@@ -149,7 +149,7 @@ class MyTableWidget(QWidget):
         QWidget.__init__(self)
         self.layout = QVBoxLayout(self)
         self.tabs_cnt = 4
-        self.tabs_titles = ['Welcome', 'BezierCurve', 'BezierCurveThreaded', 'BezierCurveBlossoms']
+        self.tabs_titles = ['Welcome', 'BezierCurve', 'BezierCurveDeCaes', 'BezierCurveBlossoms']
 
         # Initialize tab screen
         self.tabs = QTabWidget()
diff --git a/curvepy/utilities.py b/curvepy/utilities.py
index 83cbb3acac59ba27f1ebac09fdd39f8881e57fda..7046f0a145dc275929637a3f5a6fba222778af03 100644
--- a/curvepy/utilities.py
+++ b/curvepy/utilities.py
@@ -1,6 +1,9 @@
-from typing import Any, List, Callable
+from typing import Any, List, Callable, Tuple, Union
 import numpy as np
 import functools
+import sys
+
+from curvepy.de_caes import subdivision
 
 
 def straight_line_point(a: np.ndarray, b: np.ndarray, t: float = 0.5) -> np.ndarray:
@@ -98,10 +101,184 @@ def ratio(left_point: np.ndarray, col_point: np.ndarray, right_point: np.ndarray
     return 0
 
 
+def horner(m: np.ndarray, t: float = 0.5) -> Tuple[Union[float, Any], ...]:
+    """
+    TODO show which problem this is
+    Method using horner's method to calculate point with given t
+
+    Parameters
+    ----------
+    m: np.ndarray:
+        array containing coefficients
+
+    t: float:
+        value for which point is calculated
+
+    Returns
+    -------
+    tuple:
+        point calculated with given t
+    """
+    return tuple(functools.reduce(lambda x, y: t * x + y, m[i, ::-1]) for i in [0, 1])
+
+
+def distance_to_line(p1: np.ndarray, p2: np.ndarray, p_to_check: np.ndarray) -> float:
+    """
+    Method calculating distance of point to line through p1 and p2
+
+    Parameters
+    ----------
+    p1: np.ndarray:
+        beginning point of line
+
+    p2: np.ndarray:
+        end point of line
+
+    p_to_check: np.ndarray:
+        point for which distance is calculated
+
+    Returns
+    -------
+    float:
+        distance from point to line
+
+    Notes
+    -----
+    Given p1 and p2 we can check the distance p3 has to the line going through p1 and p2 as follows:
+    math:: distance(p1,p2,p3) = \\frac{|(x_1-x_1)(y_1-y_3) - (x_1-x_3)(y_2-y_1)|}{//sqrt{(x_2-x_1)^2 + (y_2-y_1)^2}}
+    more information on "https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line"
+    """
+
+    if any((p.shape[0] < 2 for p in [p1, p2, p_to_check])):
+        raise Exception("At least 2 dimensions are needed")
+
+    if p1.shape != p2.shape != p_to_check.shape:
+        raise Exception("points need to be of the same dimension!")
+
+    numerator = abs((p2[0] - p1[0]) * (p1[1] - p_to_check[1]) - (p1[0] - p_to_check[0]) * (p2[1] - p1[1]))
+    denominator = ((p2[0] - p1[0]) ** 2 + (p2[1] - p1[1]) ** 2) ** 0.5
+    return numerator / denominator
+
+
+def check_flat(m: np.ndarray, tol: float = sys.float_info.epsilon) -> bool:
+    """
+    Method checking if all points between the first and the last point
+    are less than tol away from line through beginning and end point of bezier curve
+
+    Parameters
+    ----------
+    m: np.ndarray:
+        points of curve
+
+    tol: float:
+        tolerance for distance check
+
+    Returns
+    -------
+    bool:
+        True if all point are less than tol away from line otherwise false
+    """
+    return all(distance_to_line(m[:, 0], m[:, len(m[0]) - 1], m[:, i]) <= tol for i in range(1, len(m[0]) - 1))
+
+
+def intersect_lines(p1: np.ndarray, p2: np.ndarray, p3: np.ndarray, p4: np.ndarray) -> Union[np.ndarray, None]:
+    """
+    Method checking if line through p1, p2 intersects with line through p3, p4
+
+
+    Parameters
+    ----------
+    p1: np.ndarray:
+        first point of first line
+
+    p2: np.ndarray:
+        second point of first line
+
+    p3: np.ndarray:
+        first point of second line
+
+    p4: np.ndarray:
+        second point of second line
+
+    Returns
+    -------
+    bool:
+        True if all point are less than tol away from line otherwise false
+    """
+
+    if p1.shape != p2.shape != p3.shape != p4.shape:
+        raise Exception("Points need to be of the same dimension!")
+
+    # First we vertical stack the points in an array
+    vertical_stack = np.vstack([p1, p2, p3, p4])
+    # Then we transform them to homogeneous coordinates, to perform a little trick
+    homogeneous = np.hstack((vertical_stack, np.ones((4, 1))))
+    # having our points in this form we can get the lines through the cross product
+    line_1, line_2 = np.cross(homogeneous[0], homogeneous[1]), np.cross(homogeneous[2], homogeneous[3])
+    # when we calculate the cross product of the lines we get intersect point
+    x, y, z = np.cross(line_1, line_2)
+    if z == 0:
+        return None
+    # we divide with z to turn back to 2D space
+    return np.array([x / z, y / z])
+
+
 def flatten_list_of_lists(xss: List[List[Any]]) -> List[Any]:
     return sum(xss, [])
 
 
+def min_max_box(m: np.ndarray) -> List[np.ndarray]:
+    """
+    Method creates minmax box for the corresponding bezier points
+    """
+    box = [m[0, :].min(), m[0, :].max(), m[1, :].min(), m[1, :].max()]
+    if m.shape[0] == 2:
+        return box
+    box.extend([m[2, :].min(), m[2, :].max()])
+    return box
+
+
+def intersect(m: np.ndarray, tol: float = sys.float_info.epsilon) -> np.ndarray:
+    """
+    Method checks if curve intersects with x-axis
+
+    Parameters
+    ----------
+    m: np.ndarray:
+        points of curve
+
+    tol: float:
+        tolerance for check_flat
+
+    Returns
+    -------
+    np.ndarray:
+        Points where curve and x-axis intersect
+    """
+    box = min_max_box(m)
+    res = np.array([])
+
+    if box[2] * box[3] > 0:
+        # Both y values are positive, ergo curve lies above x_axis
+        return np.array([])
+
+    if check_flat(m, tol):
+        # poly is flat enough, so we can perform intersect of straight lines
+        # since we are assuming poly is a straight line we define a line through first and las point of poly
+        # additionally we create a line which demonstrates the x axis
+        # having these two lines we can check them for intersection
+        p = intersect_lines(m[:, 0], m[:, -1], np.array([0, 0]), np.array([1, 0]))
+        if p is not None:
+            res = np.append(res, p.reshape((2, 1)), axis=1)
+    else:
+        # if poly not flat enough we subdivide and check the resulting polygons for intersection
+        p1, p2 = subdivision(m)
+        res = np.append(res, intersect(p1, tol).reshape((2, 1)), axis=1)
+        res = np.append(res, intersect(p2, tol).reshape((2, 1)), axis=1)
+
+    return res
+
+
 def csv_read(file_path: str) -> np.ndarray:
     try:
         with open(file_path, 'r') as csv_file: