diff --git a/include/fintamath/functions/IFunction.hpp b/include/fintamath/functions/IFunction.hpp
index b70be33a0..e409f99ab 100644
--- a/include/fintamath/functions/IFunction.hpp
+++ b/include/fintamath/functions/IFunction.hpp
@@ -65,8 +65,7 @@ class IFunction : public IMathObject {
 
   static std::unique_ptr<IMathObject> makeExprChecked(const IFunction &function, const ArgumentsRefVector &args);
 
-  // TODO: uncomment
-  // static std::unique_ptr<IMathObject> makeExpr(const IFunction &function, const ArgumentsRefVector &args);
+  static std::unique_ptr<IMathObject> makeExpr(const IFunction &function, const ArgumentsPtrVector &args);
 
   static bool isExpression(const ArgumentRef &arg);
 
diff --git a/include/fintamath/functions/powers/Root.hpp b/include/fintamath/functions/powers/Root.hpp
new file mode 100644
index 000000000..c55a3b5f1
--- /dev/null
+++ b/include/fintamath/functions/powers/Root.hpp
@@ -0,0 +1,20 @@
+#pragma once
+
+#include "fintamath/functions/IFunction.hpp"
+#include "fintamath/numbers/INumber.hpp"
+
+namespace fintamath {
+
+class Root : public IFunctionCRTP<INumber, Root, INumber, INumber> {
+public:
+  Root() = default;
+
+  std::string toString() const override {
+    return "root";
+  }
+
+protected:
+  std::unique_ptr<IMathObject> call(const ArgumentsRefVector &argsVect) const override;
+};
+
+}
diff --git a/include/fintamath/numbers/IntegerFunctions.hpp b/include/fintamath/numbers/IntegerFunctions.hpp
index 48d6c7b72..c84f81192 100644
--- a/include/fintamath/numbers/IntegerFunctions.hpp
+++ b/include/fintamath/numbers/IntegerFunctions.hpp
@@ -48,6 +48,8 @@ Integer factorial(const Integer &rhs);
 
 Integer factorial(const Integer &rhs, size_t order);
 
+std::map<Integer, Integer> factors(Integer rhs, Integer limit = -1);
+
 Integer combinations(const Integer &totalNumber, const Integer &choosedNumber);
 
 Integer multinomialCoefficient(const Integer &totalNumber, const std::vector<Integer> &groupNumbers);
diff --git a/src/fintamath/config/ExpressionConfig.cpp b/src/fintamath/config/ExpressionConfig.cpp
index 9a2294ca0..b0cb8f28a 100644
--- a/src/fintamath/config/ExpressionConfig.cpp
+++ b/src/fintamath/config/ExpressionConfig.cpp
@@ -57,6 +57,7 @@
 #include "fintamath/functions/other/Rad.hpp"
 #include "fintamath/functions/powers/Exp.hpp"
 #include "fintamath/functions/powers/Pow.hpp"
+#include "fintamath/functions/powers/Root.hpp"
 #include "fintamath/functions/powers/Sqrt.hpp"
 #include "fintamath/functions/trigonometry/Acos.hpp"
 #include "fintamath/functions/trigonometry/Acot.hpp"
@@ -308,6 +309,11 @@ struct ExpressionConfig {
       return makeExpr(Pow(), args.front(), std::make_shared<Rational>(1, 2));
     });
 
+    // TODO! uncomment
+    // Expression::registerFunctionExpressionMaker<Root>([](const ArgumentsPtrVector &args) {
+    //   return makeExpr(Pow(), args.front(), makeExpr(Div(), std::make_shared<Integer>(1), args.back()));
+    // });
+
     Expression::registerFunctionExpressionMaker<Sin>([](const ArgumentsPtrVector &args) {
       return std::make_shared<TrigonometryExpression>(Sin(), args.front());
     });
diff --git a/src/fintamath/config/FunctionConfig.cpp b/src/fintamath/config/FunctionConfig.cpp
index d4cf500cf..4c23181ab 100644
--- a/src/fintamath/config/FunctionConfig.cpp
+++ b/src/fintamath/config/FunctionConfig.cpp
@@ -9,10 +9,9 @@ std::unique_ptr<IMathObject> IFunction::makeExprChecked(const IFunction &functio
   return fintamath::makeExprChecked(function, args);
 }
 
-// TODO: uncomment
-// std::unique_ptr<IMathObject> IFunction::makeExpr(const IFunction &function, const ArgumentsRefVector &args) {
-//   return fintamath::makeExpr(function, args);
-// }
+std::unique_ptr<IMathObject> IFunction::makeExpr(const IFunction &function, const ArgumentsPtrVector &args) {
+  return fintamath::makeExpr(function, args)->clone();
+}
 
 bool IFunction::isExpression(const ArgumentRef &arg) {
   return is<IExpression>(arg);
diff --git a/src/fintamath/config/ParserConfig.cpp b/src/fintamath/config/ParserConfig.cpp
index 5ca5862ab..be27ce4c6 100644
--- a/src/fintamath/config/ParserConfig.cpp
+++ b/src/fintamath/config/ParserConfig.cpp
@@ -49,6 +49,7 @@
 #include "fintamath/functions/other/Rad.hpp"
 #include "fintamath/functions/powers/Exp.hpp"
 #include "fintamath/functions/powers/Pow.hpp"
+#include "fintamath/functions/powers/Root.hpp"
 #include "fintamath/functions/powers/Sqrt.hpp"
 #include "fintamath/functions/trigonometry/Acos.hpp"
 #include "fintamath/functions/trigonometry/Acot.hpp"
@@ -172,6 +173,7 @@ struct ParserConfig {
     IFunction::registerType<Lg>();
     IFunction::registerType<Exp>();
     IFunction::registerType<Sqrt>();
+    IFunction::registerType<Root>();
     IFunction::registerType<Sin>();
     IFunction::registerType<Cos>();
     IFunction::registerType<Tan>();
diff --git a/src/fintamath/expressions/binary/PowExpression.cpp b/src/fintamath/expressions/binary/PowExpression.cpp
index cfa45223c..a983bbe25 100644
--- a/src/fintamath/expressions/binary/PowExpression.cpp
+++ b/src/fintamath/expressions/binary/PowExpression.cpp
@@ -7,6 +7,7 @@
 #include "fintamath/functions/arithmetic/Mul.hpp"
 #include "fintamath/functions/arithmetic/Neg.hpp"
 #include "fintamath/functions/powers/Pow.hpp"
+#include "fintamath/functions/powers/Root.hpp"
 #include "fintamath/functions/powers/Sqrt.hpp"
 #include "fintamath/numbers/Integer.hpp"
 #include "fintamath/numbers/IntegerFunctions.hpp"
@@ -23,218 +24,216 @@ std::string PowExpression::toString() const {
     const Integer &numerator = val->numerator();
     const Integer &denominator = val->denominator();
 
-    if (denominator == 2) {
-      if (numerator == 1) {
+    if (numerator == 1) {
+      if (denominator == 2) {
         return functionToString(Sqrt(), {lhsChild});
       }
 
-      PowExpression res(std::make_shared<PowExpression>(lhsChild, Rational(1, 2).clone()), numerator.clone());
-      return res.IBinaryExpression::toString();
+      return functionToString(Root(), {lhsChild, denominator.clone()});
     }
+
+    return IBinaryExpression::toString();
   }
 
-  return IBinaryExpression::toString();
-}
+  std::shared_ptr<IFunction> PowExpression::getOutputFunction() const {
+    if (*rhsChild == Rational(1, 2)) {
+      return std::make_shared<Sqrt>();
+    }
 
-std::shared_ptr<IFunction> PowExpression::getOutputFunction() const {
-  if (*rhsChild == Rational(1, 2)) {
-    return std::make_shared<Sqrt>();
+    return IBinaryExpression::getFunction();
   }
 
-  return IBinaryExpression::getFunction();
-}
+  ArgumentPtr PowExpression::preciseSimplify() const {
+    if (*rhsChild == Rational(1, 2)) {
+      auto preciseExpr = cast<PowExpression>(clone());
+      preciseSimplifyChild(preciseExpr->lhsChild);
+      return preciseExpr;
+    }
 
-ArgumentPtr PowExpression::preciseSimplify() const {
-  if (*rhsChild == Rational(1, 2)) {
-    auto preciseExpr = cast<PowExpression>(clone());
-    preciseSimplifyChild(preciseExpr->lhsChild);
-    return preciseExpr;
+    return IBinaryExpression::preciseSimplify();
   }
 
-  return IBinaryExpression::preciseSimplify();
-}
+  PowExpression::SimplifyFunctionsVector PowExpression::getFunctionsForPreSimplify() const {
+    static const PowExpression::SimplifyFunctionsVector simplifyFunctions = {
+        &PowExpression::negSimplify, //
+        &PowExpression::powSimplify, //
+    };
+    return simplifyFunctions;
+  }
 
-PowExpression::SimplifyFunctionsVector PowExpression::getFunctionsForPreSimplify() const {
-  static const PowExpression::SimplifyFunctionsVector simplifyFunctions = {
-      &PowExpression::negSimplify, //
-      &PowExpression::powSimplify, //
-  };
-  return simplifyFunctions;
-}
+  PowExpression::SimplifyFunctionsVector PowExpression::getFunctionsForPostSimplify() const {
+    static const PowExpression::SimplifyFunctionsVector simplifyFunctions = {
+        &PowExpression::numSimplify,     //
+        &PowExpression::polynomSimplify, //
+    };
+    return simplifyFunctions;
+  }
 
-PowExpression::SimplifyFunctionsVector PowExpression::getFunctionsForPostSimplify() const {
-  static const PowExpression::SimplifyFunctionsVector simplifyFunctions = {
-      &PowExpression::numSimplify,     //
-      &PowExpression::polynomSimplify, //
-  };
-  return simplifyFunctions;
-}
+  // Use bites representation for generate all partitions of numbers, using stars and bars method
+  // https://en.wikipedia.org/wiki/Stars_and_bars_(combinatorics)
+  // https://en.wikipedia.org/wiki/Partition_(number_theory)
+  // This method generate first bit number in sequence of numbers for each partition
+  Integer PowExpression::generateFirstNum(const Integer &countOfOne) {
+    Integer n = 0;
 
-// Use bites representation for generate all partitions of numbers, using stars and bars method
-// https://en.wikipedia.org/wiki/Stars_and_bars_(combinatorics)
-// https://en.wikipedia.org/wiki/Partition_(number_theory)
-// This method generate first bit number in sequence of numbers for each partition
-Integer PowExpression::generateFirstNum(const Integer &countOfOne) {
-  Integer n = 0;
+    for (int i = 0; i < countOfOne; i++) {
+      n = n << 1 | 1;
+    }
 
-  for (int i = 0; i < countOfOne; i++) {
-    n = n << 1 | 1;
+    return n;
   }
 
-  return n;
-}
+  // Use alghorithm for generating next number in sequence of partitions
+  // https://en.wikipedia.org/wiki/Partition_(number_theory)
+  // https://en.wikipedia.org/wiki/Combinatorial_number_system#Applications
+  Integer PowExpression::generateNextNumber(Integer n) {
+    Integer u = n & -n;
+    Integer v = u + n;
+    n = v + (((v ^ n) / u) >> 2);
+    return n;
+  }
 
-// Use alghorithm for generating next number in sequence of partitions
-// https://en.wikipedia.org/wiki/Partition_(number_theory)
-// https://en.wikipedia.org/wiki/Combinatorial_number_system#Applications
-Integer PowExpression::generateNextNumber(Integer n) {
-  Integer u = n & -n;
-  Integer v = u + n;
-  n = v + (((v ^ n) / u) >> 2);
-  return n;
-}
+  std::vector<Integer> PowExpression::getPartition(Integer bitNumber, const Integer &variableCount) {
+    std::vector<Integer> result;
+    Integer counter = 0;
 
-std::vector<Integer> PowExpression::getPartition(Integer bitNumber, const Integer &variableCount) {
-  std::vector<Integer> result;
-  Integer counter = 0;
+    while (int64_t(result.size()) < variableCount) {
+      if (bitNumber % 2 == 1) {
+        counter++;
+      }
 
-  while (int64_t(result.size()) < variableCount) {
-    if (bitNumber % 2 == 1) {
-      counter++;
-    }
+      if (bitNumber % 2 == 0) {
+        result.emplace_back(counter);
+        counter = 0;
+      }
 
-    if (bitNumber % 2 == 0) {
-      result.emplace_back(counter);
-      counter = 0;
+      bitNumber >>= 1;
     }
 
-    bitNumber >>= 1;
+    return result;
   }
 
-  return result;
-}
+  // Use multinomial theorem for exponentiation of sum
+  //  https://en.wikipedia.org/wiki/Multinomial_theorem
+  ArgumentPtr PowExpression::sumPolynomSimplify(const ArgumentPtr &expr, const Integer &powValue) {
+    auto sumExpr = cast<IExpression>(expr);
+    ArgumentsPtrVector polynom;
 
-// Use multinomial theorem for exponentiation of sum
-//  https://en.wikipedia.org/wiki/Multinomial_theorem
-ArgumentPtr PowExpression::sumPolynomSimplify(const ArgumentPtr &expr, const Integer &powValue) {
-  auto sumExpr = cast<IExpression>(expr);
-  ArgumentsPtrVector polynom;
+    if (sumExpr && is<Add>(sumExpr->getFunction())) {
+      polynom = sumExpr->getChildren();
+    }
+    else {
+      return {};
+    }
 
-  if (sumExpr && is<Add>(sumExpr->getFunction())) {
-    polynom = sumExpr->getChildren();
-  }
-  else {
-    return {};
-  }
+    ArgumentsPtrVector newPolynom;
+    Integer variableCount = int64_t(polynom.size());
 
-  ArgumentsPtrVector newPolynom;
-  Integer variableCount = int64_t(polynom.size());
+    Integer bitNumber = generateFirstNum(powValue);
 
-  Integer bitNumber = generateFirstNum(powValue);
+    for (int i = 0; i < combinations(powValue + variableCount - 1, powValue); i++) {
+      std::vector<Integer> vectOfPows = getPartition(bitNumber, variableCount);
+      bitNumber = generateNextNumber(bitNumber);
 
-  for (int i = 0; i < combinations(powValue + variableCount - 1, powValue); i++) {
-    std::vector<Integer> vectOfPows = getPartition(bitNumber, variableCount);
-    bitNumber = generateNextNumber(bitNumber);
+      ArgumentsPtrVector mulExprPolynom;
+      mulExprPolynom.emplace_back(std::make_shared<Integer>(multinomialCoefficient(powValue, vectOfPows)));
 
-    ArgumentsPtrVector mulExprPolynom;
-    mulExprPolynom.emplace_back(std::make_shared<Integer>(multinomialCoefficient(powValue, vectOfPows)));
+      for (size_t j = 0; j < size_t(variableCount); j++) {
+        auto powExpr = makeExpr(Pow(), polynom[j], std::make_shared<Integer>(std::move(vectOfPows[j])));
+        mulExprPolynom.emplace_back(powExpr);
+      }
 
-    for (size_t j = 0; j < size_t(variableCount); j++) {
-      auto powExpr = makeExpr(Pow(), polynom[j], std::make_shared<Integer>(std::move(vectOfPows[j])));
-      mulExprPolynom.emplace_back(powExpr);
+      ArgumentPtr mulExpr = makeExpr(Mul(), mulExprPolynom);
+      newPolynom.emplace_back(mulExpr);
     }
 
-    ArgumentPtr mulExpr = makeExpr(Mul(), mulExprPolynom);
-    newPolynom.emplace_back(mulExpr);
+    return makeExpr(Add(), newPolynom)->toMinimalObject();
   }
 
-  return makeExpr(Add(), newPolynom)->toMinimalObject();
-}
+  ArgumentPtr PowExpression::negSimplify(const IFunction & /*func*/, const ArgumentPtr &lhs, const ArgumentPtr &rhs) {
+    if (auto lhsExpr = cast<IExpression>(lhs); lhsExpr && is<Neg>(lhsExpr->getFunction())) {
+      auto lhsMul = makeExpr(Pow(), std::make_shared<Integer>(-1), rhs);
+      auto rhsMul = makeExpr(Pow(), lhsExpr->getChildren().front(), rhs);
+      return makeExpr(Mul(), lhsMul, rhsMul)->toMinimalObject();
+    }
 
-ArgumentPtr PowExpression::negSimplify(const IFunction & /*func*/, const ArgumentPtr &lhs, const ArgumentPtr &rhs) {
-  if (auto lhsExpr = cast<IExpression>(lhs); lhsExpr && is<Neg>(lhsExpr->getFunction())) {
-    auto lhsMul = makeExpr(Pow(), std::make_shared<Integer>(-1), rhs);
-    auto rhsMul = makeExpr(Pow(), lhsExpr->getChildren().front(), rhs);
-    return makeExpr(Mul(), lhsMul, rhsMul)->toMinimalObject();
+    return {};
   }
 
-  return {};
-}
+  ArgumentPtr PowExpression::powSimplify(const IFunction & /*func*/, const ArgumentPtr &lhs, const ArgumentPtr &rhs) {
+    if (auto lhsExpr = cast<IExpression>(lhs); lhsExpr && is<Pow>(lhsExpr->getFunction())) {
+      auto lhsPow = lhsExpr->getChildren().front();
+      auto rhsPow = makeExpr(Mul(), lhsExpr->getChildren().back(), rhs);
+      return makeExpr(Pow(), lhsPow, rhsPow)->toMinimalObject();
+    }
 
-ArgumentPtr PowExpression::powSimplify(const IFunction & /*func*/, const ArgumentPtr &lhs, const ArgumentPtr &rhs) {
-  if (auto lhsExpr = cast<IExpression>(lhs); lhsExpr && is<Pow>(lhsExpr->getFunction())) {
-    auto lhsPow = lhsExpr->getChildren().front();
-    auto rhsPow = makeExpr(Mul(), lhsExpr->getChildren().back(), rhs);
-    return makeExpr(Pow(), lhsPow, rhsPow)->toMinimalObject();
+    return {};
   }
 
-  return {};
-}
+  ArgumentPtr PowExpression::numSimplify(const IFunction & /*func*/, const ArgumentPtr &lhs, const ArgumentPtr &rhs) {
+    auto lhsInt = cast<Integer>(lhs);
+    auto rhsInt = cast<Integer>(rhs);
 
-ArgumentPtr PowExpression::numSimplify(const IFunction & /*func*/, const ArgumentPtr &lhs, const ArgumentPtr &rhs) {
-  auto lhsInt = cast<Integer>(lhs);
-  auto rhsInt = cast<Integer>(rhs);
+    if (rhsInt) {
+      if (*rhsInt == 0) {
+        return std::make_shared<Integer>(1);
+      }
 
-  if (rhsInt) {
-    if (*rhsInt == 0) {
-      return std::make_shared<Integer>(1);
-    }
+      if (*rhsInt == 1 || (lhsInt && *lhsInt == 1)) {
+        return lhs;
+      }
 
-    if (*rhsInt == 1 || (lhsInt && *lhsInt == 1)) {
-      return lhs;
-    }
+      if (*rhsInt == -1) {
+        return makeExpr(Div(), std::make_shared<Integer>(1), lhs)->toMinimalObject();
+      }
 
-    if (*rhsInt == -1) {
-      return makeExpr(Div(), std::make_shared<Integer>(1), lhs)->toMinimalObject();
+      if (*rhsInt < 0) {
+        return makeExpr(Div(), std::make_shared<Integer>(1), makeExpr(Pow(), lhs, makeExpr(Neg(), rhsInt)))
+            ->toMinimalObject();
+      }
     }
 
-    if (*rhsInt < 0) {
-      return makeExpr(Div(), std::make_shared<Integer>(1), makeExpr(Pow(), lhs, makeExpr(Neg(), rhsInt)))
-          ->toMinimalObject();
+    if (lhsInt && *lhsInt == 0) {
+      return lhs;
     }
-  }
 
-  if (lhsInt && *lhsInt == 0) {
-    return lhs;
+    return {};
   }
 
-  return {};
-}
+  ArgumentPtr PowExpression::polynomSimplify(const IFunction & /*func*/, const ArgumentPtr &lhs,
+                                             const ArgumentPtr &rhs) {
+    if (auto res = mulSimplify(lhs, rhs)) {
+      return res;
+    }
 
-ArgumentPtr PowExpression::polynomSimplify(const IFunction & /*func*/, const ArgumentPtr &lhs, const ArgumentPtr &rhs) {
-  if (auto res = mulSimplify(lhs, rhs)) {
-    return res;
-  }
+    if (auto res = sumSimplify(lhs, rhs)) {
+      return res;
+    }
 
-  if (auto res = sumSimplify(lhs, rhs)) {
-    return res;
+    return {};
   }
 
-  return {};
-}
+  ArgumentPtr PowExpression::mulSimplify(const ArgumentPtr &lhs, const ArgumentPtr &rhs) {
+    if (auto mulExpr = cast<IExpression>(lhs); mulExpr && is<Mul>(mulExpr->getFunction())) {
+      ArgumentsPtrVector args = mulExpr->getChildren();
 
-ArgumentPtr PowExpression::mulSimplify(const ArgumentPtr &lhs, const ArgumentPtr &rhs) {
-  if (auto mulExpr = cast<IExpression>(lhs); mulExpr && is<Mul>(mulExpr->getFunction())) {
-    ArgumentsPtrVector args = mulExpr->getChildren();
+      for (auto &arg : args) {
+        arg = makeExpr(Pow(), arg, rhs->clone());
+      }
 
-    for (auto &arg : args) {
-      arg = makeExpr(Pow(), arg, rhs->clone());
+      return makeExpr(Mul(), args)->toMinimalObject();
     }
 
-    return makeExpr(Mul(), args)->toMinimalObject();
+    return {};
   }
 
-  return {};
-}
-
-ArgumentPtr PowExpression::sumSimplify(const ArgumentPtr &lhs, const ArgumentPtr &rhs) {
-  if (const auto rhsInt = cast<Integer>(rhs)) {
-    if (auto sumExpr = sumPolynomSimplify(lhs, *rhsInt)) {
-      return sumExpr;
+  ArgumentPtr PowExpression::sumSimplify(const ArgumentPtr &lhs, const ArgumentPtr &rhs) {
+    if (const auto rhsInt = cast<Integer>(rhs)) {
+      if (auto sumExpr = sumPolynomSimplify(lhs, *rhsInt)) {
+        return sumExpr;
+      }
     }
-  }
-
-  return {};
-}
 
+    return {};
+  }
 }
diff --git a/src/fintamath/functions/powers/Pow.cpp b/src/fintamath/functions/powers/Pow.cpp
index da92b22c7..08c4082a6 100644
--- a/src/fintamath/functions/powers/Pow.cpp
+++ b/src/fintamath/functions/powers/Pow.cpp
@@ -1,5 +1,6 @@
 #include "fintamath/functions/powers/Pow.hpp"
 
+#include "fintamath/functions/powers/Root.hpp"
 #include "fintamath/functions/powers/Sqrt.hpp"
 #include "fintamath/numbers/Integer.hpp"
 #include "fintamath/numbers/IntegerFunctions.hpp"
@@ -29,10 +30,13 @@ std::unique_ptr<IMathObject> Pow::call(const ArgumentsRefVector &argsVect) const
         return pow(lhs, numerator).toMinimalObject();
       }
 
-      if (denominator == 2) { // TODO: implement nth root
+      // TODO! uncomment
+      // return Pow()(*Root()(lhs, denominator), numerator);
+
+      // TODO! remove
+      if (denominator == 2) {
         return Pow()(*Sqrt()(lhs), numerator);
       }
-
       return pow(convert<Real>(lhs), convert<Real>(rhs)).toMinimalObject();
     });
 
diff --git a/src/fintamath/functions/powers/Root.cpp b/src/fintamath/functions/powers/Root.cpp
new file mode 100644
index 000000000..07a8158a5
--- /dev/null
+++ b/src/fintamath/functions/powers/Root.cpp
@@ -0,0 +1,64 @@
+#include "fintamath/functions/powers/Root.hpp"
+
+#include "fintamath/functions/arithmetic/Add.hpp"
+#include "fintamath/numbers/Integer.hpp"
+#include "fintamath/numbers/IntegerFunctions.hpp"
+#include "fintamath/numbers/Rational.hpp"
+#include "fintamath/numbers/Real.hpp"
+#include "fintamath/numbers/RealFunctions.hpp"
+
+namespace fintamath {
+
+std::map<Integer, Rational> rootSimpl(const Integer &lhs, const Integer &rhs) {
+  std::map<Integer, Rational> factorPowers;
+  std::map<Integer, Integer> factorRates = factors(lhs);
+
+  for (const auto &factorRate : factorRates) {
+    factorPowers
+  }
+
+  return factorPowers;
+}
+
+std::unique_ptr<IMathObject> Root::call(const ArgumentsRefVector &argsVect) const {
+  static const auto multiRoot = [] {
+    static MultiMethod<std::unique_ptr<IMathObject>(const INumber &, const INumber &)> outMultiRoot;
+
+    outMultiRoot.add<Integer, Integer>([](const Integer &lhs, const Integer &rhs) {
+      if (rhs == 2) { // TODO: implement perfect nth-roots minimization
+        Integer remainder;
+        Integer lhsSqrt = sqrt(lhs, remainder);
+
+        if (remainder == 0) {
+          return lhsSqrt.clone();
+        }
+      }
+
+      auto rootVals = rootSimpl(lhs, rhs);
+
+      return makeExpr(Add(), {outsideRoot.clone(), makeExpr(Root(), {insideRoot.clone(), rhs.clone()})});
+    });
+
+    outMultiRoot.add<Rational, Integer>([](const Rational &lhs, const Integer &rhs) {
+      return nullptr;
+    });
+
+    outMultiRoot.add<Rational, Integer>([](const Rational &lhs, const Integer &rhs) {
+      return nullptr;
+    });
+
+    return outMultiRoot;
+  }();
+
+  const auto &lhs = cast<INumber>(argsVect.front().get());
+
+  if (const auto rhs = cast<Integer>(argsVect.back().get()))
+
+    if (rhs < Integer(0)) {
+      return multiRoot(lhs, rhs);
+    }
+
+  return Pow()();
+}
+
+}
diff --git a/src/fintamath/functions/powers/Sqrt.cpp b/src/fintamath/functions/powers/Sqrt.cpp
index 83635ffe8..66e93e375 100644
--- a/src/fintamath/functions/powers/Sqrt.cpp
+++ b/src/fintamath/functions/powers/Sqrt.cpp
@@ -9,6 +9,8 @@
 namespace fintamath {
 
 std::unique_ptr<IMathObject> Sqrt::call(const ArgumentsRefVector &argsVect) const {
+  // TODO! move this logic to Root
+
   static const auto multiSqrt = [] {
     static MultiMethod<std::unique_ptr<IMathObject>(const INumber &)> outMultiSqrt;
 
diff --git a/src/fintamath/numbers/IntegerFunctions.cpp b/src/fintamath/numbers/IntegerFunctions.cpp
index 35822c763..43bf08f0b 100644
--- a/src/fintamath/numbers/IntegerFunctions.cpp
+++ b/src/fintamath/numbers/IntegerFunctions.cpp
@@ -86,6 +86,40 @@ Integer factorial(const Integer &rhs, size_t order) {
   return res;
 }
 
+std::map<Integer, Integer> factors(Integer rhs, Integer limit) {
+  if (rhs <= 0) {
+    throw UndefinedFunctionException("factors", {rhs.toString()});
+  }
+
+  std::map<Integer, Integer> factorRates;
+
+  while (rhs % 2 == 0) {
+    factorRates[2]++;
+    rhs /= 2;
+  }
+
+  {
+    Integer rhsSqrt = sqrt(rhs);
+
+    if (limit < 0 || limit > rhsSqrt) {
+      limit = rhsSqrt;
+    }
+  }
+
+  for (Integer i = 3; i <= limit; i += 2) {
+    while (rhs % i == 0) {
+      factorRates[i]++;
+      rhs = rhs / i;
+    }
+  }
+
+  if (rhs > 1) {
+    factorRates[rhs]++;
+  }
+
+  return factorRates;
+}
+
 // Use number of combinations formula.
 // https://en.wikipedia.org/wiki/Combination#Number_of_k-combinations.
 Integer combinations(const Integer &totalNumber, const Integer &choosedNumber) {
diff --git a/tests/src/functions/powers/RootTests.cpp b/tests/src/functions/powers/RootTests.cpp
new file mode 100644
index 000000000..723d3615c
--- /dev/null
+++ b/tests/src/functions/powers/RootTests.cpp
@@ -0,0 +1,73 @@
+#include "gtest/gtest.h"
+
+#include "fintamath/exceptions/UndefinedException.hpp"
+#include "fintamath/functions/powers/Root.hpp"
+
+#include "fintamath/literals/Variable.hpp"
+#include "fintamath/numbers/Rational.hpp"
+#include "fintamath/numbers/Real.hpp"
+
+using namespace fintamath;
+
+const Root f;
+
+TEST(RootTests, toStringTest) {
+  EXPECT_EQ(f.toString(), "root");
+}
+
+TEST(RootTests, getFunctionTypeTest) {
+  EXPECT_EQ(f.getFunctionType(), IFunction::Type::Binary);
+}
+
+TEST(RootTests, callTest) {
+  EXPECT_EQ(f(Integer(100), Integer(2))->toString(), "10");
+  EXPECT_EQ(f(Integer(25), Integer(2))->toString(), "5");
+  EXPECT_EQ(f(Integer(100), Integer(2))->toString(), "10");
+  EXPECT_EQ(f(Integer(144), Integer(2))->toString(), "12");
+  EXPECT_EQ(f(Integer("10000000000000000000000000000000000000000000000000000"), Integer(2))->toString(),
+            "100000000000000000000000000");
+  EXPECT_EQ(f(Integer(10), Integer(2))->toString(),
+            "3.1622776601683793319988935444327185337195551393252168268575048527925944386392382");
+  EXPECT_EQ(f(Integer(35), Integer(2))->toString(),
+            "5.916079783099616042567328291561617048415501230794340322879719669142822459105653");
+  EXPECT_EQ(f(Integer(4212), Integer(2))->toString(),
+            "64.899922958351807276145982814468927032523338329214431828788155012089005069274188");
+  EXPECT_EQ(f(Integer("992188888888"), Integer(2))->toString(),
+            "996086.7878292533601855411744463733737437739712744225743190752774953323543572221");
+  // TODO! pass this test
+  // EXPECT_EQ(f(Integer("68732648273642987365932706179432649827364"), Integer(2))->toString(),
+  //           "262169121510606178721.58566765529332797469887558257026225414531129497927294552571");
+
+  // EXPECT_EQ(f(Rational(25))->toString(), "5");
+  // EXPECT_EQ(f(Rational(25, 169))->toString(), "5/13");
+  // EXPECT_EQ(f(Rational(144, 49))->toString(), "12/7");
+  // EXPECT_EQ(f(Rational(144, 49))->toString(), "12/7");
+  // EXPECT_EQ(f(Rational("1.44"))->toString(), "6/5");
+  // EXPECT_EQ(f(Rational(1, 10))->toString(),
+  //           "0.31622776601683793319988935444327185337195551393252168268575048527925944386392382");
+  // EXPECT_EQ(f(Rational(25, 48))->toString(),
+  //           "0.72168783648703220563643597562744681955950218908765859502325290810497209037866668");
+  // EXPECT_EQ(f(Rational(26, 49))->toString(),
+  //           "0.72843135908468354714688915843182599850911013515708520108356725777513337600888917");
+  // EXPECT_EQ(f(Rational(45, 67))->toString(),
+  //           "0.81953754706222955297929294102425245239518786681123328404647027909565181499392632");
+  // EXPECT_EQ(f(Rational("68732648273642987365932706179432649827364.144"))->toString(),
+  //           "262169121510606178721.5856676552933279746991502144551506779910090326304113931067");
+
+  // EXPECT_EQ(f(Real(144))->toString(), "12");
+  // EXPECT_EQ(f(Real(2))->toString(),
+  // "1.414213562373095048801688724209698078569671875376948073176679737990732478462107");
+
+  EXPECT_EQ(f(Variable("a"), Integer(2))->toString(), "root(a, 2)");
+
+  EXPECT_THROW(f(Integer(-10)), UndefinedException);
+  EXPECT_THROW(f(Rational(-9289, 10)), UndefinedException);
+  EXPECT_THROW(f(Real(-9289)), UndefinedException);
+
+  EXPECT_THROW(f(), InvalidInputFunctionException);
+  EXPECT_THROW(f(Integer(1), Integer(1), Integer(1)), InvalidInputFunctionException);
+}
+
+TEST(RootTests, toStringTest1) {
+  EXPECT_EQ(f.toString(), "root");
+}