Implement Factorial of Real numbers

fintarin · Aug 1, 2023 · 8df1e11 · 8df1e11
1 parent a7bff04
commit 8df1e11
Show file tree

Hide file tree

Showing 11 changed files with 198 additions and 36 deletions.
diff --git a/include/fintamath/functions/other/Factorial.hpp b/include/fintamath/functions/other/Factorial.hpp
@@ -7,11 +7,19 @@
 
 namespace fintamath {
 
+class Integer;
+class Rational;
+class Real;
+
 class Factorial : public IOperatorCRTP<INumber, Factorial, INumber> {
 public:
   Factorial() : IOperatorCRTP(IOperator::Priority::PostfixUnary) {
   }
 
+  Factorial(size_t inOrder) : Factorial() {
+    setOrder(inOrder);
+  }
+
   std::string toString() const override {
     return std::string(order, '!');
   }
@@ -28,6 +36,15 @@ class Factorial : public IOperatorCRTP<INumber, Factorial, INumber> {
 protected:
   std::unique_ptr<IMathObject> call(const ArgumentsRefVector &argsVect) const override;
 
+private:
+  static std::unique_ptr<IMathObject> multiFactorialSimpl(const INumber &lhs, size_t order);
+
+  static std::unique_ptr<IMathObject> factorialSimpl(const Integer &rhs, size_t order);
+
+  static std::unique_ptr<IMathObject> factorialSimpl(const Rational &rhs, size_t order);
+
+  static std::unique_ptr<IMathObject> factorialSimpl(const Real &rhs, size_t order);
+
 private:
   size_t order = 1;
 };

diff --git a/include/fintamath/numbers/Integer.hpp b/include/fintamath/numbers/Integer.hpp
@@ -26,6 +26,10 @@ class Integer : public IIntegerCRTP<Integer> {
 
   explicit Integer(std::string str);
 
+  template <typename T, typename = std::enable_if_t<std::is_arithmetic_v<T>>>
+  explicit Integer(T val) : backend(val) {
+  }
+
   Integer(int64_t val);
 
   std::string toString() const override;

diff --git a/include/fintamath/numbers/RealFunctions.hpp b/include/fintamath/numbers/RealFunctions.hpp
@@ -56,4 +56,6 @@ Real getE();
 
 Real getPi();
 
+Real tgamma(const Real &rhs);
+
 }
diff --git a/src/fintamath/expressions/ExpressionUtils.cpp b/src/fintamath/expressions/ExpressionUtils.cpp
@@ -34,13 +34,16 @@ std::string functionToString(const IFunction &func, const ArgumentsPtrVector &ar
   return result;
 }
 
-std::string binaryOperatorChildToString(const IOperator &oper, const ArgumentPtr &child) {
+std::string operatorChildToString(const IOperator &oper, const ArgumentPtr &child) {
   std::string childStr = child->toString();
   std::shared_ptr<IOperator> childOper;
 
   if (const auto childExpr = cast<IExpression>(child)) {
     childOper = cast<IOperator>(childExpr->getOutputFunction());
   }
+  else if (oper.getFunctionType() == IFunction::Type::Unary && childStr.front() == Neg().toString().front()) {
+    childOper = std::make_shared<Neg>();
+  }
   else if (is<Rational>(child)) {
     childOper = std::make_shared<Div>();
   }
@@ -58,7 +61,7 @@ std::string binaryOperatorChildToString(const IOperator &oper, const ArgumentPtr
       return putInBrackets(childStr);
     }
 
-    if (lhsOperPriority > operPriority || (lhsOperPriority == operPriority && !oper.isAssociative())) {
+    if (lhsOperPriority >= operPriority || (lhsOperPriority == operPriority && !oper.isAssociative())) {
       return putInBrackets(childStr);
     }
   }
@@ -76,34 +79,18 @@ std::string binaryOperatorToString(const IOperator &oper, const ArgumentPtr &lhs
     operStr = ' ' + operStr + ' ';
   }
 
-  std::string lhsStr = binaryOperatorChildToString(oper, lhs);
-  std::string rhsStr = binaryOperatorChildToString(oper, rhs);
+  std::string lhsStr = operatorChildToString(oper, lhs);
+  std::string rhsStr = operatorChildToString(oper, rhs);
 
   return lhsStr + operStr + rhsStr;
 }
 
 std::string prefixUnaryOperatorToString(const IOperator &oper, const ArgumentPtr &rhs) {
-  std::string result = oper.toString();
-
-  if (const auto child = cast<IExpression>(rhs)) {
-    if (is<IOperator>(child->getOutputFunction())) {
-      return result + putInBrackets(rhs->toString());
-    }
-  }
-
-  return result + rhs->toString();
+  return oper.toString() + operatorChildToString(oper, rhs);
 }
 
 std::string postfixUnaryOperatorToString(const IOperator &oper, const ArgumentPtr &rhs) {
-  std::string result = rhs->toString();
-
-  if (const auto child = cast<IExpression>(rhs)) {
-    if (is<IOperator>(child->getOutputFunction())) {
-      return putInBrackets(result) + oper.toString();
-    }
-  }
-
-  return result + oper.toString();
+  return operatorChildToString(oper, rhs) + oper.toString();
 }
 
 bool hasVariables(const std::shared_ptr<const IExpression> &expr) {

diff --git a/src/fintamath/functions/other/Factorial.cpp b/src/fintamath/functions/other/Factorial.cpp
@@ -1,18 +1,75 @@
 #include "fintamath/functions/other/Factorial.hpp"
 
-#include "fintamath/exceptions/UndefinedException.hpp"
+#include "fintamath/literals/constants/ComplexInf.hpp"
 #include "fintamath/numbers/IntegerFunctions.hpp"
+#include "fintamath/numbers/Rational.hpp"
+#include "fintamath/numbers/Real.hpp"
+#include "fintamath/numbers/RealFunctions.hpp"
 
 namespace fintamath {
 
 std::unique_ptr<IMathObject> Factorial::call(const ArgumentsRefVector &argsVect) const {
-  const auto &rhs = argsVect.front().get();
+  const auto &rhs = cast<INumber>(argsVect.front().get());
+  return multiFactorialSimpl(rhs, order);
+}
+
+std::unique_ptr<IMathObject> Factorial::multiFactorialSimpl(const INumber &lhs, size_t order) {
+  static const auto multiFactorial = [] {
+    static MultiMethod<std::unique_ptr<IMathObject>(const INumber &, const Integer &)> outMultiAbs;
+
+    outMultiAbs.add<Integer, Integer>([](const Integer &inRhs, const Integer &inOrder) {
+      return factorialSimpl(inRhs, size_t(inOrder));
+    });
+
+    outMultiAbs.add<Rational, Integer>([](const Rational &inRhs, const Integer &inOrder) {
+      return factorialSimpl(inRhs, size_t(inOrder));
+    });
+
+    outMultiAbs.add<Real, Integer>([](const Real &inRhs, const Integer &inOrder) {
+      return factorialSimpl(inRhs, size_t(inOrder));
+    });
+
+    return outMultiAbs;
+  }();
 
-  if (!is<Integer>(rhs)) {
-    throw UndefinedUnaryOperatorException(toString(), rhs.toString(), UndefinedUnaryOperatorException::Type::Postfix);
+  return multiFactorial(lhs, Integer(order));
+}
+
+std::unique_ptr<IMathObject> Factorial::factorialSimpl(const Integer &rhs, size_t order) {
+  if (rhs < 0) {
+    if (order != 1) {
+      return makeExpr(Factorial(order), rhs);
+    }
+
+    return ComplexInf().clone();
   }
 
-  return factorial(cast<Integer>(rhs), order).toMinimalObject();
+  return factorial(rhs, order).toMinimalObject();
+}
+
+std::unique_ptr<IMathObject> Factorial::factorialSimpl(const Rational &rhs, size_t order) {
+  if (rhs.denominator() == 1) {
+    return factorialSimpl(rhs.numerator(), order);
+  }
+
+  if (order != 1) {
+    return makeExpr(Factorial(order), rhs);
+  }
+
+  return factorialSimpl(Real(rhs), order);
+}
+
+std::unique_ptr<IMathObject> Factorial::factorialSimpl(const Real &rhs, size_t order) {
+  if (order != 1) {
+    return makeExpr(Factorial(order), rhs);
+  }
+
+  try {
+    return tgamma(rhs + 1).toMinimalObject();
+  }
+  catch (const UndefinedException &) {
+    return ComplexInf().clone();
+  }
 }
 
 }
diff --git a/src/fintamath/numbers/RealFunctions.cpp b/src/fintamath/numbers/RealFunctions.cpp
@@ -208,4 +208,16 @@ Real getPi() {
   return {cpp_dec_float_100(get_constant_pi<cpp_dec_float_100::backend_type>())};
 }
 
+Real tgamma(const Real &rhs) {
+  try {
+    return boost::math::tgamma(rhs.getBackend());
+  }
+  catch (const std::domain_error &) {
+    throw UndefinedFunctionException("tgamma", {rhs.toString()});
+  }
+  catch (const std::overflow_error &) {
+    throw UndefinedFunctionException("tgamma", {rhs.toString()});
+  }
+}
+
 }
diff --git a/tests/src/FintamathTests.cpp b/tests/src/FintamathTests.cpp
@@ -43,7 +43,7 @@ TEST(FintamathTests, fintamathTests) {
   expr = Expression("(4x^4 + 1 + 3x^3 + 2x) / (x^2 + x + 2)");
   EXPECT_EQ(expr.toString(), "4 x^2 - x - 7 + (11 x + 15)/(x^2 + x + 2)");
 
-  expr = log(2, 256) + ln(pow(E(), 2));
+  expr = log(2, 256) + ln(pow(e(), 2));
   EXPECT_EQ(expr.toString(), "10");
 
   expr = sqrt(Expression(8));

diff --git a/tests/src/expressions/ExpressionTests.cpp b/tests/src/expressions/ExpressionTests.cpp
@@ -112,6 +112,11 @@ TEST(ExpressionTests, stringConstructorTest) {
   EXPECT_EQ(Expression("30!!!!!!").toString(), "933120");
   EXPECT_EQ(Expression("30!!!!!!!").toString(), "198720");
 
+  EXPECT_EQ(Expression("(2/3)!").toString(), "(2/3)!");
+  EXPECT_EQ(Expression("E!").toString(), "E!");
+  EXPECT_EQ(Expression("(2/3)!!").toString(), "(2/3)!!");
+  EXPECT_EQ(Expression("(-1)!!").toString(), "(-1)!!");
+
   EXPECT_EQ(Expression("sqrt144").toString(), "12");
   EXPECT_EQ(Expression("sqrt0").toString(), "0");
   EXPECT_EQ(Expression("sqrt(25)").toString(), "5");
@@ -738,6 +743,7 @@ TEST(ExpressionTests, stringConstructorTest) {
   EXPECT_EQ(Expression("log(1, 10)").toString(), "ComplexInf");
   EXPECT_EQ(Expression("log(10, 0)").toString(), "-Inf");
   EXPECT_EQ(Expression("log(1/10, 0)").toString(), "Inf");
+  EXPECT_EQ(Expression("(-1)!").toString(), "ComplexInf");
 
   EXPECT_EQ(Expression("0*Inf").toString(), "Undefined");
   EXPECT_EQ(Expression("0*-Inf").toString(), "Undefined");
@@ -1064,10 +1070,6 @@ TEST(ExpressionTests, stringConstructorNegativeTest) {
   EXPECT_THROW(Expression("max()"), InvalidInputException);
   EXPECT_THROW(Expression("max(True, False)"), InvalidInputException);
 
-  EXPECT_THROW(Expression("(-1)!"), UndefinedException);
-  EXPECT_THROW(Expression("(2/3)!"), UndefinedException);
-  EXPECT_THROW(Expression("(-1)!!"), UndefinedException);
-  EXPECT_THROW(Expression("(2/3)!!"), UndefinedException);
   EXPECT_THROW(Expression("sqrt(-1)"), UndefinedException);
   EXPECT_THROW(Expression("sqrt(-1)"), UndefinedException);
   EXPECT_THROW(Expression("(-1)^(2/3)"), UndefinedException);
@@ -1078,7 +1080,6 @@ TEST(ExpressionTests, stringConstructorNegativeTest) {
   EXPECT_THROW(Expression("lg(-1)"), UndefinedException);
   EXPECT_THROW(Expression("log(-1, 1)"), UndefinedException);
   // TODO constants
-  // EXPECT_THROW(Expression("E!"), UndefinedException);
   // EXPECT_THROW(Expression("tan(Pi/2)"), UndefinedException);
   // EXPECT_THROW(Expression("cot(0)"), UndefinedException);
   // EXPECT_THROW(Expression("asin(2)"), UndefinedException);
@@ -1186,6 +1187,15 @@ TEST(ExpressionTests, preciseTest) {
   EXPECT_EQ(Expression("((x - z)^2 / 8) * (x / y)").precise().toString(),
             "(0.125 x^3)/y + (-0.25 x^2 z)/y + (0.125 x z^2)/y");
 
+  EXPECT_EQ(Expression("(2/3)!").precise().toString(),
+            "0.90274529295093361129685868543634252367955151070452913226268164530918864360116169");
+  EXPECT_EQ(Expression("E!").precise().toString(),
+            "4.2608204763570033817001212246457024649334243739593219749116048935993443487275001");
+  EXPECT_EQ(Expression("(2/3)!!").precise().toString(),
+            "0.66666666666666666666666666666666666666666666666666666666666666666666666666666667!!");
+  EXPECT_EQ(Expression("(1/1000000000000000000000000000000000000000)!!").precise().toString(), "(1*10^-39)!!");
+  EXPECT_EQ(Expression("(-1)!!").precise().toString(), "(-1)!!");
+
   EXPECT_EQ(Expression("ln(x)").precise().toString(), "ln(x)");
   EXPECT_EQ(Expression("sqrt(x)").precise().toString(), "sqrt(x)");
   EXPECT_EQ(Expression("root(x, 3)").precise().toString(), "root(x, 3)");

diff --git a/tests/src/functions/other/FactorialTests.cpp b/tests/src/functions/other/FactorialTests.cpp
@@ -7,6 +7,7 @@
 #include "fintamath/functions/arithmetic/UnaryPlus.hpp"
 #include "fintamath/literals/Variable.hpp"
 #include "fintamath/numbers/Rational.hpp"
+#include "fintamath/numbers/Real.hpp"
 
 using namespace fintamath;
 
@@ -27,12 +28,40 @@ TEST(FactorialTests, getOperatorPriorityTest) {
 TEST(FactorialTests, callTest) {
   EXPECT_EQ(f(Integer(0))->toString(), "1");
   EXPECT_EQ(f(Integer(1))->toString(), "1");
+  EXPECT_EQ(f(Integer(2))->toString(), "2");
+  EXPECT_EQ(f(Integer(3))->toString(), "6");
+  EXPECT_EQ(f(Integer(4))->toString(), "24");
+  EXPECT_EQ(f(Integer(5))->toString(), "120");
   EXPECT_EQ(f(Integer(10))->toString(), "3628800");
+  EXPECT_EQ(f(Integer(500))->toString(),
+            "1220136825991110068701238785423046926253574342803192842192413588385845373153881997605496447502203281863013"
+            "6164771482035841633787220781772004807852051593292854779075719393306037729608590862704291745478824249127263"
+            "4430567017327076946106280231045264421887878946575477714986349436778103764427403382736539747138647787849543"
+            "8489595537537990423241061271326984327745715546309977202781014561081188373709531016356324432987029563896628"
+            "9116589747695720879269288712817800702651745077684107196243903943225364226052349458501299185715012487069615"
+            "6814162535905669342381300885624924689156412677565448188650659384795177536089400574523894033579847636394490"
+            "5313062323749066445048824665075946735862074637925184200459369692981022263971952597190945217823331756934581"
+            "5085523328207628200234026269078983424517120062077146409794561161276291459512372299133401695523638509428855"
+            "9201872743379517301458635757082835578015873543276888868012039988238470215146760544540766353598417443048012"
+            "8938313896881639487469658817504506926365338175055478128640000000000000000000000000000000000000000000000000"
+            "000000000000000000000000000000000000000000000000000000000000000000000000000");
 
-  EXPECT_EQ(f(Variable("a"))->toString(), "a!");
+  EXPECT_EQ(f(Rational(10, 1))->toString(), "3628800");
+  EXPECT_EQ(f(Rational(1, 10))->toString(),
+            "0.9513507698668731836292487177265402192550578626088377343050000770434265383322821");
+  EXPECT_EQ(f(Rational(-1, 10))->toString(),
+            "1.0686287021193193548973053356944807781698387850609731790493706839815721770254476");
+
+  EXPECT_EQ(f(Real("0.1"))->toString(),
+            "0.9513507698668731836292487177265402192550578626088377343050000770434265383322821");
+
+  // TODO!!! add tests with double, triple, etc. factorials
 
-  EXPECT_THROW(f(Integer(-10)), UndefinedUnaryOperatorException);
-  EXPECT_THROW(f(Rational(1, 10)), UndefinedUnaryOperatorException);
+  EXPECT_EQ(f(Integer(-10))->toString(), "ComplexInf");
+  EXPECT_EQ(f(Rational(-10))->toString(), "ComplexInf");
+  EXPECT_EQ(f(Real(-10))->toString(), "ComplexInf");
+
+  EXPECT_EQ(f(Variable("a"))->toString(), "a!");
 
   EXPECT_THROW(f(), InvalidInputFunctionException);
   EXPECT_THROW(f(Integer(1), Integer(1), Integer(1)), InvalidInputFunctionException);

diff --git a/tests/src/numbers/IntegerTests.cpp b/tests/src/numbers/IntegerTests.cpp
@@ -24,6 +24,13 @@ TEST(IntegerTests, stringConstructorTest) {
   EXPECT_THROW(Integer("+"), InvalidInputException);
 }
 
+TEST(IntegerTests, templateConstructorTest) {
+  EXPECT_EQ(Integer(std::numeric_limits<uint64_t>::max()).toString(), "18446744073709551615");
+  EXPECT_EQ(Integer(std::numeric_limits<int64_t>::min()).toString(), "-9223372036854775808");
+  EXPECT_EQ(Integer(std::numeric_limits<uint8_t>::max()).toString(), "255");
+  EXPECT_EQ(Integer(std::numeric_limits<int8_t>::min()).toString(), "-128");
+}
+
 TEST(IntegerTests, intConstructorTest) {
   EXPECT_EQ(Integer(10), 10);
 

diff --git a/tests/src/numbers/RealFunctionsTests.cpp b/tests/src/numbers/RealFunctionsTests.cpp
@@ -248,3 +248,40 @@ TEST(RealFunctionsTests, getETest) {
 TEST(RealFunctionsTests, getEPiTest) {
   EXPECT_EQ(getPi().toString(), "3.141592653589793238462643383279502884197169399375105820974944592307816406286209");
 }
+
+TEST(RealFunctionsTests, tgammaTest) {
+  EXPECT_EQ(tgamma(Real("1")).toString(), "1");
+  EXPECT_EQ(tgamma(Real("2")).toString(), "1");
+  EXPECT_EQ(tgamma(Real("3")).toString(), "2");
+  EXPECT_EQ(tgamma(Real("4")).toString(), "6");
+  EXPECT_EQ(tgamma(Real("5")).toString(), "24");
+  EXPECT_EQ(tgamma(Real("10")).toString(), "362880");
+  EXPECT_EQ(tgamma(Real("1000000")).toString(),
+            "8.2639316883312400623766461031726662911353479789638730451677758855633796110356451*10^5565702");
+
+  EXPECT_EQ(tgamma(Real("0.5")).toString(),
+            "1.7724538509055160272981674833411451827975494561223871282138077898529112845910322");
+  EXPECT_EQ(tgamma(Real("1.5")).toString(),
+            "0.88622692545275801364908374167057259139877472806119356410690389492645564229551609");
+  EXPECT_EQ(tgamma(Real("1.23345")).toString(),
+            "0.90997105032312305741807535506003700168381988861711231924267338164564458593112273");
+  EXPECT_EQ(tgamma(Real("10.888")).toString(),
+            "2790176.1744474603546712580801633763601209060179121811338563253860175273379208134");
+  EXPECT_EQ(tgamma(Real("-0.5")).toString(),
+            "-3.5449077018110320545963349666822903655950989122447742564276155797058225691820644");
+  EXPECT_EQ(tgamma(Real("-1.5")).toString(),
+            "2.3632718012073547030642233111215269103967326081631828376184103864705483794547096");
+  EXPECT_EQ(tgamma(Real("-1.23345")).toString(),
+            "4.1813173753338157686838484929306398465161945343856820450249490692138002199035393");
+  EXPECT_EQ(tgamma(Real("-10.888")).toString(),
+            "-3.0005493180448293721869822512345553711440749496581172835311328514510130843563454*10^-7");
+
+  EXPECT_THROW(tgamma(Real("0")), UndefinedFunctionException);
+  EXPECT_THROW(tgamma(Real("-1")), UndefinedFunctionException);
+  EXPECT_THROW(tgamma(Real("-2")), UndefinedFunctionException);
+  EXPECT_THROW(tgamma(Real("-3")), UndefinedFunctionException);
+  EXPECT_THROW(tgamma(Real("-32352")), UndefinedFunctionException);
+
+  EXPECT_THROW(tgamma(Real("1000000000")), UndefinedFunctionException);
+  EXPECT_THROW(tgamma(Real("-1000000000")), UndefinedFunctionException);
+}