diff --git a/src/shapepy/analytic/__init__.py b/src/shapepy/analytic/__init__.py
index a333e43f..170447d7 100644
--- a/src/shapepy/analytic/__init__.py
+++ b/src/shapepy/analytic/__init__.py
@@ -9,5 +9,7 @@
 * Bezier
 """
 
+from ..tools import To
+from .base import IAnalytic
 from .bezier import Bezier
 from .polynomial import Polynomial
diff --git a/src/shapepy/analytic/base.py b/src/shapepy/analytic/base.py
index ad98a92d..0898a58d 100644
--- a/src/shapepy/analytic/base.py
+++ b/src/shapepy/analytic/base.py
@@ -7,11 +7,11 @@
 
 from abc import ABC, abstractmethod
 from functools import lru_cache
-from typing import Iterable, Set, Union
+from typing import Set, Union
 
-from ..rbool import SubSetR1, WholeR1, from_any
+from ..rbool import SubSetR1
 from ..scalar.reals import Real
-from ..tools import Is
+from ..tools import Is, vectorize
 
 
 @lru_cache(maxsize=None)
@@ -54,9 +54,28 @@ def domain(self) -> SubSetR1:
         raise NotImplementedError
 
     @abstractmethod
-    def __call__(self, node: Real, derivate: int = 0) -> Real:
+    def eval(self, node: Real, derivate: int = 0) -> Real:
+        """
+        Evaluates the given analytic function at given node.
+
+        The optional parameter 'derivate' gives if a derivative is required
+
+        Example
+        -------
+        >>> polynomial = Polynomial([1, 2, 3])  # p(x) = 1 + 2 * x + 3 * x^2
+        >>> polynomial.eval(0)  # p(0) = 1 + 2 * 0 + 3 * 0^2
+        1
+        >>> polynomial.eval(1)  # p(1) = 1 + 2 * 1 + 3 * 1^2
+        6
+        >>> polynomial.eval(1, 1)  # p'(1) = 2 + 6 * 1
+        8
+        """
         raise NotImplementedError
 
+    @vectorize(1, 0)
+    def __call__(self, node: Real) -> Real:
+        return self.eval(node, 0)
+
     @abstractmethod
     def __add__(self, other: Union[Real, IAnalytic]) -> IAnalytic:
         raise NotImplementedError
@@ -66,99 +85,57 @@ def __mul__(self, other: Union[Real, IAnalytic]) -> IAnalytic:
         raise NotImplementedError
 
     @abstractmethod
-    def shift(self, amount: Real) -> IAnalytic:
+    def derivate(self, times: int = 1) -> IAnalytic:
         """
-        Transforms the analytic p(t) into p(t-d) by
-        translating the analytic by 'd' to the right.
+        Derivate the polynomial curve, giving a new one
 
         Example
         -------
-        >>> old_poly = Polynomial([0, 0, 0, 1])
-        >>> print(old_poly)
-        t^3
-        >>> new_poly = shift(poly, 1)  # transform to (t-1)^3
-        >>> print(new_poly)
-        - 1 + 3 * t - 3 * t^2 + t^3
+        >>> poly = Polynomial([1, 2, 5])
+        >>> print(poly)
+        1 + 2 * t + 5 * t^2
+        >>> dpoly = derivate(poly)
+        >>> print(dpoly)
+        2 + 10 * t
         """
         raise NotImplementedError
 
     @abstractmethod
-    def scale(self, amount: Real) -> IAnalytic:
+    def integrate(self, domain: SubSetR1) -> Real:
         """
-        Transforms the analytic p(t) into p(A*t) by
-        scaling the argument of the analytic by 'A'.
+        Derivate the polynomial curve, giving a new one
 
         Example
         -------
-        >>> old_poly = Polynomial([0, 2, 0, 1])
-        >>> print(old_poly)
-        2 * t + t^3
-        >>> new_poly = scale(poly, 2)
-        >>> print(new_poly)
-        4 * t + 8 * t^3
-        """
-        raise NotImplementedError
-
-    @abstractmethod
-    def clean(self) -> IAnalytic:
-        """
-        Cleans the curve, removing the unnecessary coefficients
+        >>> poly = Polynomial([1, 2, 5])
+        >>> print(poly)
+        1 + 2 * t + 5 * t^2
+        >>> dpoly = derivate(poly)
+        >>> print(dpoly)
+        2 + 10 * t
         """
         raise NotImplementedError
 
     @abstractmethod
-    def integrate(self, times: int = 1) -> IAnalytic:
+    def compose(self, function: IAnalytic) -> IAnalytic:
         """
-        Integrates the analytic function
-        """
-        raise NotImplementedError
+        Compose the analytic function: h(x) = f(g(x))
 
-    @abstractmethod
-    def derivate(self, times: int = 1) -> IAnalytic:
-        """
-        Derivates the analytic function
+        Example
+        -------
+        >>> f = Polynomial([0, 0, 1])  # f(x) = x^2
+        >>> g = Polynomial([-1, 1])  # f(x) = x - 1
+        >>> f.compose(g)  # h(x) = (x - 1)^2
+        1 - 2 * x + x^2
         """
         raise NotImplementedError
 
-
-class BaseAnalytic(IAnalytic):
-    """
-    Base class parent of Analytic classes
-    """
-
-    def __init__(self, coefs: Iterable[Real], domain: SubSetR1 = WholeR1()):
-        if not Is.iterable(coefs):
-            raise TypeError("Expected an iterable of coefficients")
-        coefs = tuple(coefs)
-        if len(coefs) == 0:
-            raise ValueError("Cannot receive an empty tuple")
-        self.__coefs = coefs
-        self.domain = domain
-
-    @property
-    def domain(self) -> SubSetR1:
-        return self.__domain
-
-    @domain.setter
-    def domain(self, subset: SubSetR1):
-        self.__domain = from_any(subset)
-
-    @property
-    def ncoefs(self) -> int:
-        """
-        Returns the number of coefficients that determines the analytic
-        """
-        return len(self.__coefs)
-
-    def __iter__(self):
-        yield from self.__coefs
-
-    def __getitem__(self, index):
-        return self.__coefs[index]
-
     def __neg__(self) -> IAnalytic:
         return -1 * self
 
+    def __sub__(self, other: Union[Real, IAnalytic]) -> IAnalytic:
+        return self.__add__(-other)
+
     def __pow__(self, exponent: int) -> IAnalytic:
         if not Is.integer(exponent) or exponent < 0:
             raise ValueError
@@ -171,9 +148,6 @@ def __pow__(self, exponent: int) -> IAnalytic:
             cache[n] = cache[n // 2] * cache[n - n // 2]
         return cache[exponent]
 
-    def __sub__(self, other: Union[Real, IAnalytic]) -> IAnalytic:
-        return self.__add__(-other)
-
     def __rsub__(self, other: Real) -> IAnalytic:
         return (-self).__add__(other)
 
@@ -182,26 +156,3 @@ def __radd__(self, other: Real) -> IAnalytic:
 
     def __rmul__(self, other: Real) -> IAnalytic:
         return self.__mul__(other)
-
-    def __repr__(self) -> str:
-        if self.domain is WholeR1():
-            return str(self)
-        return f"{self.domain}: {self}"
-
-
-def is_analytic(obj: object) -> bool:
-    """
-    Tells if given object is an analytic function
-    """
-    return Is.instance(obj, IAnalytic)
-
-
-def derivate_analytic(ana: IAnalytic) -> IAnalytic:
-    """
-    Computes the derivative of an Analytic function
-    """
-    assert is_analytic(ana)
-    return ana.derivate()
-
-
-Is.analytic = is_analytic
diff --git a/src/shapepy/analytic/bezier.py b/src/shapepy/analytic/bezier.py
index 310e2267..2adb9fa6 100644
--- a/src/shapepy/analytic/bezier.py
+++ b/src/shapepy/analytic/bezier.py
@@ -7,12 +7,10 @@
 from functools import lru_cache
 from typing import Iterable, Tuple, Union
 
-from ..loggers import debug
-from ..rbool import SubSetR1, WholeR1
+from ..rbool import SubSetR1, WholeR1, from_any
 from ..scalar.quadrature import inner
 from ..scalar.reals import Math, Rational, Real
-from ..tools import Is, NotExpectedError, To
-from .base import BaseAnalytic, IAnalytic
+from ..tools import Is, To
 from .polynomial import Polynomial
 
 
@@ -53,141 +51,33 @@ def inverse_caract_matrix(degree: int) -> Tuple[Tuple[Rational, ...], ...]:
     return tuple(map(tuple, matrix))
 
 
-def bezier2polynomial(bezier: Bezier) -> Polynomial:
+def bezier2polynomial(coefs: Iterable[Real]) -> Iterable[Real]:
     """
     Converts a Bezier instance to Polynomial
     """
-    coefs = tuple(bezier)
+    coefs = tuple(coefs)
     matrix = bezier_caract_matrix(len(coefs) - 1)
-    poly_coefs = (inner(weights, bezier) for weights in matrix)
-    return Polynomial(poly_coefs, bezier.domain)
+    return (inner(weights, coefs) for weights in matrix)
 
 
-def polynomial2bezier(polynomial: Polynomial) -> Bezier:
+def polynomial2bezier(coefs: Iterable[Real]) -> Iterable[Real]:
     """
     Converts a Polynomial instance to a Bezier
     """
-    coefs = tuple(polynomial)
+    coefs = tuple(coefs)
     matrix = inverse_caract_matrix(len(coefs) - 1)
-    ctrlpoints = (inner(weights, coefs) for weights in matrix)
-    return Bezier(ctrlpoints, polynomial.domain)
+    return (inner(weights, coefs) for weights in matrix)
 
 
-class Bezier(BaseAnalytic):
+class Bezier(Polynomial):
     """
     Defines the Bezier class, that allows evaluating and operating
     such as adding, subtracting, multiplying, etc
     """
 
-    def __init__(self, coefs: Iterable[Real], domain: SubSetR1 = WholeR1()):
-        super().__init__(coefs, domain)
-        self.__polynomial = bezier2polynomial(self)
-
-    @property
-    def degree(self) -> int:
-        """
-        Returns the degree of the polynomial, which is the
-        highest power of t with a non-zero coefficient.
-        If the polynomial is constant, returns 0.
-        """
-        return self.__polynomial.degree
-
-    def __eq__(self, value: object) -> bool:
-        if not Is.instance(value, IAnalytic):
-            if Is.real(value):
-                return all(ctrlpoint == value for ctrlpoint in self)
-            return NotImplemented
-        if self.domain != value.domain:
-            return False
-        if isinstance(value, Bezier):
-            return self.__polynomial == bezier2polynomial(value)
-        if isinstance(value, Polynomial):
-            return self.__polynomial == value
-        raise NotExpectedError
-
-    def __add__(self, other: Union[Real, Polynomial, Bezier]) -> Bezier:
-        if Is.instance(other, Bezier):
-            other = bezier2polynomial(other)
-        sumpoly = self.__polynomial + other
-        return polynomial2bezier(sumpoly)
-
-    def __mul__(self, other: Union[Real, Polynomial, Bezier]) -> Bezier:
-        if Is.instance(other, Bezier):
-            other = bezier2polynomial(other)
-        mulpoly = self.__polynomial * other
-        return polynomial2bezier(mulpoly)
-
-    def __call__(self, node: Real, derivate: int = 0) -> Real:
-        return self.__polynomial(node, derivate)
-
-    def __str__(self):
-        return str(self.__polynomial)
-
-    @debug("shapepy.analytic.bezier")
-    def clean(self) -> Bezier:
-        """
-        Decreases the degree of the bezier curve if possible
-        """
-        return polynomial2bezier(bezier2polynomial(self).clean())
-
-    @debug("shapepy.analytic.bezier")
-    def scale(self, amount: Real) -> Bezier:
-        """
-        Transforms the polynomial p(t) into p(A*t) by
-        scaling the argument of the polynomial by 'A'.
-
-        p(t) = a0 + a1 * t + ... + ap * t^p
-        p(A * t) = a0 + a1 * (A*t) + ... + ap * (A * t)^p
-                = b0 + b1 * t + ... + bp * t^p
-
-        Example
-        -------
-        >>> old_poly = Polynomial([0, 0, 0, 1])
-        >>> print(old_poly)
-        t^3
-        >>> new_poly = scale(poly, 1)  # transform to (t-1)^3
-        >>> print(new_poly)
-        - 1 + 3 * t - 3 * t^2 + t^3
-        """
-        return polynomial2bezier(bezier2polynomial(self).scale(amount))
-
-    @debug("shapepy.analytic.bezier")
-    def shift(self, amount: Real) -> Bezier:
-        """
-        Transforms the bezier p(t) into p(t-d) by
-        translating the bezier by 'd' to the right.
-        """
-        return polynomial2bezier(bezier2polynomial(self).shift(amount))
-
-    @debug("shapepy.analytic.bezier")
-    def integrate(self, times: int = 1) -> Bezier:
-        """
-        Integrates the bezier analytic
-
-        Example
-        -------
-        >>> poly = Polynomial([1, 2, 5])
-        >>> print(poly)
-        1 + 2 * t + 5 * t^2
-        >>> ipoly = integrate(poly)
-        >>> print(ipoly)
-        t + t^2 + (5/3) * t^3
-        """
-        return polynomial2bezier(bezier2polynomial(self).integrate(times))
-
-    @debug("shapepy.analytic.bezier")
-    def derivate(self, times: int = 1) -> Bezier:
-        """
-        Derivate the bezier curve, giving a new one
-        """
-        return polynomial2bezier(bezier2polynomial(self).derivate(times))
-
-
-def to_bezier(coefs: Iterable[Real]) -> Bezier:
-    """
-    Creates a Bezier instance
-    """
-    return Bezier(coefs).clean()
-
-
-To.bezier = to_bezier
+    def __init__(
+        self, coefs: Iterable[Real], domain: Union[None, SubSetR1] = None
+    ):
+        domain = WholeR1() if domain is None else from_any(domain)
+        poly_coefs = tuple(bezier2polynomial(coefs))
+        super().__init__(poly_coefs, domain)
diff --git a/src/shapepy/analytic/polynomial.py b/src/shapepy/analytic/polynomial.py
index 261e1f8c..c462592d 100644
--- a/src/shapepy/analytic/polynomial.py
+++ b/src/shapepy/analytic/polynomial.py
@@ -7,14 +7,13 @@
 from numbers import Real
 from typing import Iterable, List, Union
 
-from ..loggers import debug
-from ..rbool import move, scale
+from ..rbool import IntervalR1, SubSetR1, WholeR1, from_any, move, scale
 from ..scalar.reals import Math
-from ..tools import Is, To, vectorize
-from .base import BaseAnalytic, IAnalytic
+from ..tools import Is, To
+from .base import IAnalytic
 
 
-class Polynomial(BaseAnalytic):
+class Polynomial(IAnalytic):
     """
     Defines a polynomial with coefficients
 
@@ -37,6 +36,20 @@ class Polynomial(BaseAnalytic):
     5
     """
 
+    def __init__(self, coefs: Iterable[Real], domain: SubSetR1 = WholeR1()):
+        if not Is.iterable(coefs):
+            raise TypeError("Expected an iterable of coefficients")
+        coefs = tuple(coefs)
+        if len(coefs) == 0:
+            raise ValueError("Cannot receive an empty tuple")
+        degree = max((i for i, v in enumerate(coefs) if v * v > 0), default=0)
+        self.__coefs = coefs[: degree + 1]
+        self.__domain = from_any(domain)
+
+    @property
+    def domain(self) -> SubSetR1:
+        return self.__domain
+
     @property
     def degree(self) -> int:
         """
@@ -44,21 +57,31 @@ def degree(self) -> int:
         highest power of t with a non-zero coefficient.
         If the polynomial is constant, returns 0.
         """
-        return self.ncoefs - 1
+        return len(self.__coefs) - 1
+
+    def __iter__(self):
+        yield from self.__coefs
+
+    def __getitem__(self, index):
+        return self.__coefs[index]
 
     def __eq__(self, other: object) -> bool:
-        if Is.instance(other, Polynomial):
-            return tuple(self.clean()) == tuple(other.clean())
-        if Is.real(other):
-            cself = self.clean()
-            return cself.degree == 0 and cself[0] == other
-        return NotImplemented
+        if not Is.instance(other, IAnalytic):
+            if Is.finite(other):
+                return self.degree == 0 and self[0] == other
+            return NotImplemented
+        return (
+            Is.instance(other, Polynomial)
+            and other.degree == self.degree
+            and self.domain == other.domain
+            and tuple(self) == tuple(other)
+        )
 
     def __add__(self, other: Union[Real, Polynomial]) -> Polynomial:
         if not Is.instance(other, IAnalytic):
             coefs = list(self)
             coefs[0] += other
-            return self.__class__(coefs, self.domain)
+            return Polynomial(coefs, self.domain)
         if not Is.instance(other, Polynomial):
             return NotImplemented
         coefs = [0] * (1 + max(self.degree, other.degree))
@@ -66,36 +89,96 @@ def __add__(self, other: Union[Real, Polynomial]) -> Polynomial:
             coefs[i] += coef
         for i, coef in enumerate(other):
             coefs[i] += coef
-        return self.__class__(coefs, self.domain)
+        return Polynomial(coefs, self.domain)
 
     def __mul__(self, other: Union[Real, Polynomial]) -> Polynomial:
         if not Is.instance(other, IAnalytic):
-            return self.__class__((other * coef for coef in self), self.domain)
+            return Polynomial((other * coef for coef in self), self.domain)
         if not Is.instance(other, Polynomial):
             return NotImplemented
         coefs = [0 * self[0]] * (self.degree + other.degree + 1)
         for i, coefi in enumerate(self):
             for j, coefj in enumerate(other):
                 coefs[i + j] += coefi * coefj
-        return self.__class__(coefs, self.domain & other.domain)
+        return Polynomial(coefs, self.domain & other.domain)
 
-    @vectorize(1, 0)
-    def __call__(self, node: Real, derivate: int = 0) -> Real:
+    def eval(self, node: Real, derivate: int = 0) -> Real:
+        if node not in self.domain:
+            raise ValueError(f"Node {node} not in {self.domain}")
         if derivate == 0:
-            node = To.real(node)
-            if self.degree == 0:
-                return self[0]
-            if Is.infinity(node):
-                return self[self.degree] * (
-                    Math.NEGINF
-                    if self.degree % 2 and node == Math.NEGINF
-                    else Math.POSINF
-                )
-            result: Real = 0 * self[0]
-            for coef in self[::-1]:
-                result = node * result + coef
-            return result
-        return self.derivate(derivate)(node)
+            coefs = self.__coefs
+        elif self.degree < derivate:
+            coefs = (0 * self.__coefs[0],)
+        else:
+            coefs = tuple(
+                Math.factorial(n + derivate) // Math.factorial(n) * coef
+                for n, coef in enumerate(self.__coefs[derivate:])
+            )
+        node = To.real(node)
+        if len(coefs) == 1:
+            return coefs[0]
+        if Is.infinity(node):
+            return coefs[len(coefs) - 1] * (
+                Math.NEGINF
+                if self.degree % 2 and node == Math.NEGINF
+                else Math.POSINF
+            )
+        result: Real = 0 * coefs[0]
+        for coef in coefs[::-1]:
+            result = node * result + coef
+        return result
+
+    def derivate(self, times=1):
+        if self.degree < times:
+            return Polynomial([0 * self[0]], self.domain)
+        coefs = (
+            Math.factorial(n + times) // Math.factorial(n) * coef
+            for n, coef in enumerate(self[times:])
+        )
+        return Polynomial(coefs, self.domain)
+
+    def integrate(self, domain):
+        domain = from_any(domain)
+        if domain not in self.domain:
+            raise ValueError(
+                f"Given domain {domain} is not subset of {self.domain}"
+            )
+        if not Is.instance(domain, IntervalR1):
+            raise ValueError(f"Cannot integrate over {domain}")
+        left, right = domain[0], domain[1]
+        potencias = [1]
+        result = 0
+        for i, coef in enumerate(self):
+            result += coef * sum(potencias) / (i + 1)
+            potencias.append(right * potencias[-1])
+            for j in range(i + 1):
+                potencias[j] *= left
+        return result * (right - left)
+
+    def compose(self, function: IAnalytic) -> Polynomial:
+        if not Is.instance(function, Polynomial):
+            raise TypeError(f"Analytic {type(function)} is not suported")
+        if function.degree != 1:
+            raise ValueError("Only polynomials of degree = 1 are allowed")
+        shift_amount, scale_amount = tuple(function)
+        coefs = list(self)
+        domain = self.domain
+        if scale_amount != 1:
+            inv = 1 / scale_amount
+            coefs = list(coef * inv**i for i, coef in enumerate(self))
+            domain = scale(domain, scale_amount)
+        if shift_amount != 0:
+            for i, coef in enumerate(tuple(coefs)):
+                for j in range(i):
+                    value = Math.binom(i, j) * (shift_amount ** (i - j))
+                    if (i + j) % 2:
+                        value *= -1
+                    coefs[j] += coef * value
+            domain = move(domain, shift_amount)
+        return Polynomial(coefs, domain)
+
+    def __repr__(self):
+        return str(self.domain) + ": " + self.__str__()
 
     def __str__(self):
         if self.degree == 0:
@@ -118,117 +201,3 @@ def __str__(self):
                 msg += f"^{i}"
             msgs.append(msg)
         return " ".join(msgs)
-
-    @debug("shapepy.analytic.polynomial")
-    def clean(self) -> Polynomial:
-        """
-        Decreases the degree of the bezier curve if possible
-        """
-        degree = max((i for i, v in enumerate(self) if v * v > 0), default=0)
-        return Polynomial(self[: degree + 1], self.domain)
-
-    @debug("shapepy.analytic.polynomial")
-    def scale(self, amount: Real) -> Polynomial:
-        """
-        Transforms the polynomial p(t) into p(A*t) by
-        scaling the argument of the polynomial by 'A'.
-
-        p(t) = a0 + a1 * t + ... + ap * t^p
-        p(A * t) = a0 + a1 * (A*t) + ... + ap * (A * t)^p
-                = b0 + b1 * t + ... + bp * t^p
-
-        Example
-        -------
-        >>> old_poly = Polynomial([0, 2, 0, 1])
-        >>> print(old_poly)
-        2 * t + t^3
-        >>> new_poly = scale(poly, 2)
-        >>> print(new_poly)
-        4 * t + 8 * t^3
-        """
-        inv = 1 / amount
-        coefs = tuple(coef * inv**i for i, coef in enumerate(self))
-        return Polynomial(coefs, scale(self.domain, amount))
-
-    @debug("shapepy.analytic.polynomial")
-    def shift(self, amount: Real) -> Polynomial:
-        """
-        Transforms the polynomial p(t) into p(t-d) by
-        translating the polynomial by 'd' to the right.
-
-        p(t) = a0 + a1 * t + ... + ap * t^p
-        p(t-d) = a0 + a1 * (t-d) + ... + ap * (t-d)^p
-                = b0 + b1 * t + ... + bp * t^p
-
-        Example
-        -------
-        >>> old_poly = Polynomial([0, 0, 0, 1])
-        >>> print(old_poly)
-        t^3
-        >>> new_poly = shift(poly, 1)  # transform to (t-1)^3
-        >>> print(new_poly)
-        - 1 + 3 * t - 3 * t^2 + t^3
-        """
-        newcoefs = list(self)
-        for i, coef in enumerate(self):
-            for j in range(i):
-                value = Math.binom(i, j) * (amount ** (i - j))
-                if (i + j) % 2:
-                    value *= -1
-                newcoefs[j] += coef * value
-        return Polynomial(newcoefs, move(self.domain, amount))
-
-    @debug("shapepy.analytic.polynomial")
-    def integrate(self, times: int = 1) -> Polynomial:
-        """
-        Integrates the polynomial curve
-
-        Example
-        -------
-        >>> poly = Polynomial([1, 2, 5])
-        >>> print(poly)
-        1 + 2 * t + 5 * t^2
-        >>> ipoly = integrate(poly)
-        >>> print(ipoly)
-        t + t^2 + (5/3) * t^3
-        """
-        polynomial = self
-        for _ in range(times):
-            newcoefs = [0 * polynomial[0]]
-            newcoefs += list(
-                coef / (n + 1) for n, coef in enumerate(polynomial)
-            )
-            polynomial = Polynomial(newcoefs, self.domain)
-        return polynomial
-
-    @debug("shapepy.analytic.polynomial")
-    def derivate(self, times: int = 1) -> Polynomial:
-        """
-        Derivate the polynomial curve, giving a new one
-
-        Example
-        -------
-        >>> poly = Polynomial([1, 2, 5])
-        >>> print(poly)
-        1 + 2 * t + 5 * t^2
-        >>> dpoly = poly.derivate()
-        >>> print(dpoly)
-        2 + 10 * t
-        """
-        if self.degree < times:
-            return Polynomial([0 * self[0]], self.domain)
-        coefs = (
-            Math.factorial(n + times) // Math.factorial(n) * coef
-            for n, coef in enumerate(self[times:])
-        )
-        return Polynomial(coefs, self.domain)
-
-
-def to_polynomial(coeficients: Iterable[Real]) -> Polynomial:
-    """
-    Creates a polynomial instance from given coefficients
-    """
-    return Polynomial(coeficients).clean()
-
-
-To.polynomial = to_polynomial
diff --git a/src/shapepy/analytic/tools.py b/src/shapepy/analytic/tools.py
index dfcb4cb4..4a658c61 100644
--- a/src/shapepy/analytic/tools.py
+++ b/src/shapepy/analytic/tools.py
@@ -18,90 +18,10 @@
 )
 from ..scalar.reals import Math, Real
 from ..tools import Is, NotExpectedError, To
-from .base import IAnalytic, derivate_analytic
-from .bezier import Bezier, bezier2polynomial
+from .base import IAnalytic
 from .polynomial import Polynomial
 
 
-def find_polynomial_roots(
-    polynomial: Polynomial, domain: SubSetR1 = WholeR1()
-) -> SubSetR1:
-    """
-    Finds all the values of t* such p(t*) = 0 inside given domain
-    """
-    assert Is.instance(polynomial, Polynomial)
-    polynomial = polynomial.clean()
-    domain &= polynomial.domain
-    if polynomial.degree == 0:
-        return domain if polynomial[0] == 0 else EmptyR1()
-    if polynomial.degree == 1:
-        numerator = -To.rational(1, 1) * polynomial[0]
-        return SingleR1(numerator / polynomial[1])
-    if polynomial.degree == 2:
-        c, b, a = polynomial
-        delta = b * b - 4 * a * c
-        if delta < 0:
-            return EmptyR1()
-        sqrtdelta = Math.sqrt(delta)
-        half = To.rational(1, 2)
-        x0 = half * (-b - sqrtdelta) / a
-        x1 = half * (-b + sqrtdelta) / a
-        return from_any({x0, x1})
-    roots = np.roots(tuple(polynomial)[::-1])
-    roots = (To.real(np.real(r)) for r in roots if abs(np.imag(r)) < 1e-15)
-    return from_any(set(roots))
-
-
-def where_minimum_polynomial(
-    polynomial: Polynomial, domain: SubSetR1 = WholeR1()
-) -> SubSetR1:
-    """
-    Finds the value of t* such poly(t*) is minimal
-    """
-    assert Is.instance(polynomial, Polynomial)
-    domain &= polynomial.domain
-    if polynomial.degree == 0:
-        return domain
-    if domain == WholeR1() and polynomial.degree % 2:
-        return EmptyR1()
-    relation = {knot: polynomial(knot) for knot in extract_knots(domain)}
-    critical = find_roots(derivate_analytic(polynomial), domain)
-    for knot in extract_knots(critical):
-        relation[knot] = polynomial(knot)
-    minvalue = min(relation.values(), default=float("inf"))
-    relation = {
-        key: value
-        for key, value in relation.items()
-        if value == minvalue and key in domain
-    }
-    return unite(set(relation.keys()))
-
-
-def find_minimum_polynomial(
-    polynomial: Polynomial, domain: SubSetR1 = WholeR1()
-) -> Union[Real, None]:
-    """
-    Finds the minimal value of p(t) in the given domain
-
-    If the minimal does not exist, returns None
-
-    If the polynomial goes to -inf, returns -inf
-    """
-    assert Is.instance(polynomial, Polynomial)
-    if polynomial.degree == 0:
-        return polynomial[0]
-    if domain == WholeR1() and polynomial.degree % 2:
-        return Math.NEGINF
-    relation = {}
-    relation = {knot: polynomial(knot) for knot in extract_knots(domain)}
-    critical = find_roots(derivate_analytic(polynomial), domain)
-    for knot in extract_knots(critical):
-        relation[knot] = polynomial(knot)
-    return min(
-        (val for key, val in relation.items() if key in domain), default=None
-    )
-
-
 @debug("shapepy.analytic.tools")
 def find_roots(analytic: IAnalytic, domain: SubSetR1 = WholeR1()) -> SubSetR1:
     """
@@ -110,10 +30,8 @@ def find_roots(analytic: IAnalytic, domain: SubSetR1 = WholeR1()) -> SubSetR1:
     assert Is.instance(analytic, IAnalytic)
     domain = from_any(domain)
     if Is.instance(analytic, Polynomial):
-        return find_polynomial_roots(analytic, domain)
-    if Is.instance(analytic, Bezier):
-        return find_roots(bezier2polynomial(analytic), domain)
-    raise NotExpectedError
+        return PolynomialFunctions.find_roots(analytic, domain)
+    raise NotExpectedError(f"Invalid analytic: {type(analytic)}")
 
 
 @debug("shapepy.analytic.tools")
@@ -126,10 +44,8 @@ def where_minimum(
     assert Is.instance(analytic, IAnalytic)
     domain = from_any(domain)
     if Is.instance(analytic, Polynomial):
-        return where_minimum_polynomial(analytic, domain)
-    if Is.instance(analytic, Bezier):
-        return where_minimum(bezier2polynomial(analytic), domain)
-    raise NotExpectedError
+        return PolynomialFunctions.where_minimum(analytic, domain)
+    raise NotExpectedError(f"Invalid analytic: {type(analytic)}")
 
 
 @debug("shapepy.analytic.tools")
@@ -142,7 +58,85 @@ def find_minimum(
     assert Is.instance(analytic, IAnalytic)
     domain = from_any(domain)
     if Is.instance(analytic, Polynomial):
-        return find_minimum_polynomial(analytic, domain)
-    if Is.instance(analytic, Bezier):
-        return find_minimum(bezier2polynomial(analytic), domain)
-    raise NotExpectedError
+        return PolynomialFunctions.find_minimum(analytic, domain)
+    raise NotExpectedError(f"Invalid analytic: {type(analytic)}")
+
+
+class PolynomialFunctions:
+    """Static class that stores static functions used for the generics
+    functions above. This class specifics for Polynomial"""
+
+    @staticmethod
+    def find_roots(polynomial: Polynomial, domain: SubSetR1) -> SubSetR1:
+        """
+        Finds all the values of t* such p(t*) = 0 inside given domain
+        """
+        assert Is.instance(polynomial, Polynomial)
+        domain &= polynomial.domain
+        if polynomial.degree == 0:
+            return domain if polynomial[0] == 0 else EmptyR1()
+        if polynomial.degree == 1:
+            numerator = -To.rational(1, 1) * polynomial[0]
+            return SingleR1(numerator / polynomial[1])
+        if polynomial.degree == 2:
+            c, b, a = polynomial
+            delta = b * b - 4 * a * c
+            if delta < 0:
+                return EmptyR1()
+            sqrtdelta = Math.sqrt(delta)
+            half = To.rational(1, 2)
+            x0 = half * (-b - sqrtdelta) / a
+            x1 = half * (-b + sqrtdelta) / a
+            return from_any({x0, x1})
+        roots = np.roots(tuple(polynomial)[::-1])
+        roots = (To.real(np.real(r)) for r in roots if abs(np.imag(r)) < 1e-15)
+        return from_any(set(roots))
+
+    @staticmethod
+    def where_minimum(polynomial: Polynomial, domain: SubSetR1) -> SubSetR1:
+        """
+        Finds the value of t* such poly(t*) is minimal
+        """
+        assert Is.instance(polynomial, Polynomial)
+        domain &= polynomial.domain
+        if polynomial.degree == 0:
+            return domain
+        if domain == WholeR1() and polynomial.degree % 2:
+            return EmptyR1()
+        relation = {knot: polynomial(knot) for knot in extract_knots(domain)}
+        critical = find_roots(polynomial.derivate(), domain)
+        for knot in extract_knots(critical):
+            relation[knot] = polynomial(knot)
+        minvalue = min(relation.values(), default=float("inf"))
+        relation = {
+            key: value
+            for key, value in relation.items()
+            if value == minvalue and key in domain
+        }
+        return unite(set(relation.keys()))
+
+    @staticmethod
+    def find_minimum(
+        polynomial: Polynomial, domain: SubSetR1
+    ) -> Union[Real, None]:
+        """
+        Finds the minimal value of p(t) in the given domain
+
+        If the minimal does not exist, returns None
+
+        If the polynomial goes to -inf, returns -inf
+        """
+        assert Is.instance(polynomial, Polynomial)
+        if polynomial.degree == 0:
+            return polynomial[0]
+        if domain == WholeR1() and polynomial.degree % 2:
+            return Math.NEGINF
+        relation = {}
+        relation = {knot: polynomial(knot) for knot in extract_knots(domain)}
+        critical = find_roots(polynomial.derivate(), domain)
+        for knot in extract_knots(critical):
+            relation[knot] = polynomial(knot)
+        return min(
+            (val for key, val in relation.items() if key in domain),
+            default=None,
+        )
diff --git a/src/shapepy/geometry/factory.py b/src/shapepy/geometry/factory.py
index e1334e8b..d3e40149 100644
--- a/src/shapepy/geometry/factory.py
+++ b/src/shapepy/geometry/factory.py
@@ -104,8 +104,7 @@ def spline_curve(spline_curve) -> JordanCurve:
         """
         beziers = spline_curve.split(spline_curve.knots)
         segments = (
-            FactorySegment.bezier(bezier.ctrlpoints).clean()
-            for bezier in beziers
+            FactorySegment.bezier(bezier.ctrlpoints) for bezier in beziers
         )
         return JordanCurve(map(USegment, segments))
 
diff --git a/src/shapepy/geometry/integral.py b/src/shapepy/geometry/integral.py
index 20830b25..34817b58 100644
--- a/src/shapepy/geometry/integral.py
+++ b/src/shapepy/geometry/integral.py
@@ -4,7 +4,6 @@
 
 from __future__ import annotations
 
-from copy import copy
 from functools import partial
 
 from ..analytic.base import IAnalytic
@@ -37,11 +36,11 @@ def polynomial(curve: Segment, expx: int, expy: int):
         assert Is.instance(curve, Segment)
         xfunc = curve.xfunc
         yfunc = curve.yfunc
-        pcrossdp = xfunc * yfunc.derivate() - yfunc * xfunc.derivate()
+        pcrossdp = xfunc * yfunc.derivate()
+        pcrossdp -= yfunc * xfunc.derivate()
         function = (xfunc**expx) * (yfunc**expy) * pcrossdp
-        assert Is.analytic(function)
-        ipoly = function.integrate()
-        return (ipoly(1) - ipoly(0)) / (expx + expy + 2)
+        assert Is.instance(function, IAnalytic)
+        return function.integrate([0, 1]) / (expx + expy + 2)
 
     @staticmethod
     def turns(curve: Segment, point: Point2D) -> float:
@@ -61,9 +60,8 @@ def turns(curve: Segment, point: Point2D) -> float:
         radius_square = deltax * deltax + deltay * deltay
         if find_minimum(radius_square, [0, 1]) < 1e-6:
             return To.rational(1, 2)
-        crossf = (
-            deltax * copy(deltay).derivate() - deltay * copy(deltax).derivate()
-        )
+        crossf = deltax * deltay.derivate()
+        crossf -= deltay * deltax.derivate()
         function = partial(
             lambda t, cf, rs: cf(t) / rs(t), cf=crossf, rs=radius_square
         )
diff --git a/src/shapepy/geometry/jordancurve.py b/src/shapepy/geometry/jordancurve.py
index 4c96f4d3..abcf1e02 100644
--- a/src/shapepy/geometry/jordancurve.py
+++ b/src/shapepy/geometry/jordancurve.py
@@ -9,6 +9,7 @@
 from copy import copy
 from typing import Iterable, Iterator
 
+from ..analytic import IAnalytic
 from ..loggers import debug
 from ..scalar.reals import Real
 from ..tools import CyclicContainer, Is, pairs, reverse
@@ -201,10 +202,10 @@ def compute_area(jordan: JordanCurve) -> Real:
         segment = usegment.parametrize()
         xfunc = segment.xfunc
         yfunc = segment.yfunc
-        poly = xfunc * yfunc.derivate() - yfunc * xfunc.derivate()
-        assert Is.analytic(poly)
-        ipoly = poly.integrate()
-        total += ipoly(1) - ipoly(0)
+        poly = xfunc * yfunc.derivate()
+        poly -= yfunc * xfunc.derivate()
+        assert Is.instance(poly, IAnalytic)
+        total += poly.integrate([0, 1])
     return total / 2
 
 
diff --git a/src/shapepy/geometry/segment.py b/src/shapepy/geometry/segment.py
index 9f852452..a5c56fda 100644
--- a/src/shapepy/geometry/segment.py
+++ b/src/shapepy/geometry/segment.py
@@ -16,6 +16,7 @@
 from typing import Iterable, Optional, Tuple
 
 from ..analytic.base import IAnalytic
+from ..analytic.bezier import Bezier
 from ..analytic.tools import find_minimum
 from ..loggers import debug
 from ..rbool import IntervalR1, from_any
@@ -40,8 +41,8 @@ def __init__(self, xfunc: IAnalytic, yfunc: IAnalytic):
             raise TypeError
         self.__length = None
         self.__knots = (To.rational(0, 1), To.rational(1, 1))
-        self.__xfunc = xfunc.clean()
-        self.__yfunc = yfunc.clean()
+        self.__xfunc = xfunc
+        self.__yfunc = yfunc
 
     def __str__(self) -> str:
         return f"BS{list(self.knots)}:({self.xfunc}, {self.yfunc})"
@@ -68,8 +69,8 @@ def __contains__(self, point: Point2D) -> bool:
 
     @vectorize(1, 0)
     def __call__(self, node: Real, derivate: int = 0) -> Point2D:
-        xcoord = self.xfunc(node, derivate)
-        ycoord = self.yfunc(node, derivate)
+        xcoord = self.xfunc.eval(node, derivate)
+        ycoord = self.yfunc.eval(node, derivate)
         return cartesian(xcoord, ycoord)
 
     @property
@@ -105,8 +106,8 @@ def derivate(self, times: Optional[int] = 1) -> Segment:
         """
         if not Is.integer(times) or times <= 0:
             raise ValueError(f"Times must be integer >= 1, not {times}")
-        dxfunc = copy(self.xfunc).derivate(times)
-        dyfunc = copy(self.yfunc).derivate(times)
+        dxfunc = self.xfunc.derivate(times)
+        dyfunc = self.yfunc.derivate(times)
         return Segment(dxfunc, dyfunc)
 
     def box(self) -> Box:
@@ -120,12 +121,6 @@ def box(self) -> Box:
         ymax = -find_minimum(-self.yfunc, [0, 1])
         return Box(cartesian(xmin, ymin), cartesian(xmax, ymax))
 
-    def clean(self) -> Segment:
-        """Cleans the segment"""
-        self.__xfunc = self.__xfunc.clean()
-        self.__yfunc = self.__yfunc.clean()
-        return self
-
     def __copy__(self) -> Segment:
         return self.__deepcopy__(None)
 
@@ -137,9 +132,9 @@ def __invert__(self) -> Segment:
         Inverts the direction of the curve.
         If the curve is clockwise, it becomes counterclockwise
         """
-        half = To.rational(1, 2)
-        xfunc = self.__xfunc.shift(-half).scale(-1).shift(half)
-        yfunc = self.__yfunc.shift(-half).scale(-1).shift(half)
+        composition = Bezier([1, 0])
+        xfunc = self.__xfunc.compose(composition)
+        yfunc = self.__yfunc.compose(composition)
         return Segment(xfunc, yfunc)
 
     def split(self, nodes: Iterable[Real]) -> Tuple[Segment, ...]:
@@ -157,8 +152,9 @@ def extract(self, interval: IntervalR1) -> Segment:
             raise TypeError
         knota, knotb = interval[0], interval[1]
         denom = 1 / (knotb - knota)
-        nxfunc = copy(self.xfunc).shift(-knota).scale(denom)
-        nyfunc = copy(self.yfunc).shift(-knota).scale(denom)
+        composition = Bezier([(-knota) * denom, (1 - knota) * denom])
+        nxfunc = self.xfunc.compose(composition)
+        nyfunc = self.yfunc.compose(composition)
         return Segment(nxfunc, nyfunc)
 
 
@@ -168,10 +164,8 @@ def compute_length(segment: Segment) -> Real:
     Computes the length of the jordan curve
     """
     domain = (0, 1)
-    dpsquare: IAnalytic = (
-        segment.xfunc.derivate() ** 2 + segment.yfunc.derivate() ** 2
-    )
-    assert Is.analytic(dpsquare)
+    dpsquare = segment.xfunc.derivate() ** 2 + segment.yfunc.derivate() ** 2
+    assert Is.instance(dpsquare, IAnalytic)
     if dpsquare == dpsquare(0):  # Check if it's constant
         return (domain[1] - domain[0]) * Math.sqrt(dpsquare(0))
     integrator = IntegratorFactory.clenshaw_curtis(3)
diff --git a/tests/analytic/test_bezier.py b/tests/analytic/test_bezier.py
index c82e679b..2a14f6ba 100644
--- a/tests/analytic/test_bezier.py
+++ b/tests/analytic/test_bezier.py
@@ -36,23 +36,7 @@ def test_degree():
     bezier = Bezier([1, 2])
     assert bezier.degree == 1
     bezier = Bezier([1, 2, 3])
-    assert bezier.degree == 2
-
-
-@pytest.mark.order(4)
-@pytest.mark.dependency(depends=["test_build"])
-def test_coefficients():
-    bezier = Bezier([0])
-    assert bezier[0] == 0
-    bezier = Bezier([1])
-    assert bezier[0] == 1
-    bezier = Bezier([1, 2])
-    assert bezier[0] == 1
-    assert bezier[1] == 2
-    bezier = Bezier([1, 2, 3])
-    assert bezier[0] == 1
-    assert bezier[1] == 2
-    assert bezier[2] == 3
+    assert bezier.degree == 1
 
 
 @pytest.mark.order(4)
@@ -194,26 +178,22 @@ def test_conversions():
     for _ in range(ntests):
         degree = np.random.randint(0, 6)
         ctrlpoints = tuple(np.random.randint(-3, 4, degree + 1))
-        bezier = Bezier(ctrlpoints)
         for _ in range(4):
-            bezier = bezier2polynomial(bezier)
-            bezier = polynomial2bezier(bezier)
-            assert bezier == Bezier(ctrlpoints)
+            poly_coefs = bezier2polynomial(ctrlpoints)
+            bezier_coefs = tuple(polynomial2bezier(poly_coefs))
+            assert bezier_coefs == ctrlpoints
 
 
 @pytest.mark.order(4)
 @pytest.mark.dependency(depends=["test_build", "test_matrices"])
 def test_clean():
-    ctrlpoints = [1, 2, 3, 4]
-    bezier = Bezier(ctrlpoints)
-    assert bezier.clean() == Bezier([1, 4])
+    assert Bezier([1, 2, 3, 4]) == Bezier([1, 4])
 
 
 @pytest.mark.order(4)
 @pytest.mark.dependency(
     depends=[
         "test_build",
-        "test_coefficients",
         "test_degree",
         "test_matrices",
         "test_compare",
diff --git a/tests/analytic/test_derivate.py b/tests/analytic/test_derivate.py
index be2a21a6..3759867a 100644
--- a/tests/analytic/test_derivate.py
+++ b/tests/analytic/test_derivate.py
@@ -1,7 +1,7 @@
 import pytest
 
+from shapepy.analytic.bezier import Bezier
 from shapepy.analytic.polynomial import Polynomial
-from shapepy.tools import To
 
 
 @pytest.mark.order(9)
@@ -19,43 +19,37 @@ def test_begin():
 @pytest.mark.order(9)
 @pytest.mark.dependency(depends=["test_begin"])
 def test_polynomial():
-    poly = To.polynomial([0])
-    assert poly.derivate(1) == 0
-    assert poly.derivate(2) == 0
+    poly = Polynomial([0])
+    assert poly.derivate() == 0
 
-    poly = To.polynomial([3])
-    assert poly.derivate(1) == 0
-    assert poly.derivate(2) == 0
+    poly = Polynomial([3])
+    assert poly.derivate() == 0
 
-    poly = To.polynomial([1, 1, 1, 1, 1])
+    poly = Polynomial([1, 1, 1, 1, 1])
     assert poly.derivate(1) == Polynomial([1, 2, 3, 4])
     assert poly.derivate(2) == Polynomial([2, 6, 12])
     assert poly.derivate(3) == Polynomial([6, 24])
 
-    assert poly(0, 1) == 1
-    assert poly(1, 1) == 10
-    assert poly(0, 2) == 2
-    assert poly(1, 2) == 20
+    assert poly.eval(0, 1) == 1
+    assert poly.eval(1, 1) == 10
+    assert poly.eval(0, 2) == 2
+    assert poly.eval(1, 2) == 20
 
 
 @pytest.mark.order(9)
 @pytest.mark.dependency(depends=["test_begin"])
 def test_bezier():
-    bezier = To.bezier([0])
-    assert bezier.derivate(1) == 0
-    assert bezier.derivate(2) == 0
+    bezier = Bezier([0])
+    assert bezier.derivate() == 0
 
-    bezier = To.bezier([3])
-    assert bezier.derivate(1) == 0
-    assert bezier.derivate(2) == 0
+    bezier = Bezier([3])
+    assert bezier.derivate() == 0
 
-    bezier = To.bezier([1, 1, 1, 1, 1])
-    assert bezier.derivate(1) == 0
-    assert bezier.derivate(2) == 0
-    assert bezier.derivate(3) == 0
+    bezier = Bezier([1, 1, 1, 1, 1])
+    assert bezier.derivate() == 0
 
-    assert bezier(0, 1) == 0
-    assert bezier(1, 1) == 0
+    assert bezier.eval(0, 1) == 0
+    assert bezier.eval(1, 1) == 0
 
 
 @pytest.mark.order(9)
diff --git a/tests/analytic/test_integrate.py b/tests/analytic/test_integrate.py
index 7797e701..16f0b878 100644
--- a/tests/analytic/test_integrate.py
+++ b/tests/analytic/test_integrate.py
@@ -2,7 +2,6 @@
 
 from shapepy.analytic.bezier import Bezier
 from shapepy.analytic.polynomial import Polynomial
-from shapepy.tools import To
 
 
 @pytest.mark.order(9)
@@ -21,34 +20,30 @@ def test_begin():
 @pytest.mark.dependency(depends=["test_begin"])
 def test_polynomial():
     poly = Polynomial([0])
-    assert poly.integrate(1) == 0
-    assert poly.integrate(2) == 0
+    assert poly.integrate([0, 1]) == 0
 
     poly = Polynomial([3])
-    assert poly.integrate(1) == Polynomial([0, 3])
-    assert poly.integrate(2) == Polynomial([0, 0, 3 / 2])
+    assert poly.integrate([0, 1]) == 3
+    assert poly.integrate([0, 2]) == 6
 
     poly = Polynomial([6, 24, 60])
-    assert poly.integrate(1) == Polynomial([0, 6, 12, 20])
-    assert poly.integrate(2) == Polynomial([0, 0, 3, 4, 5])
-    assert poly.integrate(3) == Polynomial([0, 0, 0, 1, 1, 1])
+    assert poly.integrate([0, 1]) == 38
+    assert poly.integrate([0, 2]) == 220
 
 
 @pytest.mark.order(9)
 @pytest.mark.dependency(depends=["test_begin"])
 def test_bezier():
-    bezier = To.bezier([0])
-    assert bezier.integrate(1) == 0
-    assert bezier.integrate(2) == 0
-
-    bezier = To.bezier([3])
-    assert bezier.integrate(1) == Bezier([0, 3])
-    assert bezier.integrate(2) == Bezier([0, 0, 1.5])
-
-    bezier = To.bezier([6, 6, 6, 6, 6])
-    assert bezier.integrate(1) == Bezier([0, 6])
-    assert bezier.integrate(2) == Bezier([0, 0, 3])
-    assert bezier.integrate(3) == Bezier([0, 0, 0, 1])
+    bezier = Bezier([0])
+    assert bezier.integrate([0, 1]) == 0
+
+    bezier = Bezier([3])
+    assert bezier.integrate([0, 0.5]) == 3 / 2
+    assert bezier.integrate([0.5, 1]) == 3 / 2
+
+    bezier = Bezier([6, 12, 6])
+    assert bezier.integrate([0, 0.5]) == 4
+    assert bezier.integrate([0.5, 1]) == 4
 
 
 @pytest.mark.order(9)
diff --git a/tests/analytic/test_random.py b/tests/analytic/test_random.py
index cded3be7..1318e004 100644
--- a/tests/analytic/test_random.py
+++ b/tests/analytic/test_random.py
@@ -176,7 +176,7 @@ def test_shift():
     tsample = np.linspace(-1, 1, 17)
     for funca in generator_analytic(100):
         amount = random.randint(-5, 5)
-        funcb = funca.shift(amount)
+        funcb = funca.compose(Polynomial([amount, 1]))
         valuesa = funca(tsample)
         valuesb = funcb(amount + tsample)
         np.testing.assert_allclose(valuesa, valuesb)
@@ -200,7 +200,7 @@ def test_scale():
     tsample = np.linspace(-1, 1, 17)
     for funca in generator_analytic(100):
         amount = random.randint(2, 5)
-        funcb = funca.scale(amount)
+        funcb = funca.compose(Polynomial([0, amount]))
         valuesa = funca(tsample)
         valuesb = funcb(amount * tsample)
         np.testing.assert_allclose(valuesa, valuesb, atol=1e-12, rtol=1)