support fac in formula generation (#460)

Author: loader3229

CHANGELOG.md

@@ -2,6 +2,10 @@ To install or update LODA, please follow the [installation instructions](https:/
 
 ## [Unreleased]
 
+### Enhancements
+
+* Support `fac` in formula generation
+
 ### v25.9.26
 
 ### Enhancements

src/cmd/test.cpp

@@ -771,12 +771,12 @@ void Test::incEval() {
     i++;
   }
   // OEIS sequence test cases
-  std::vector<size_t> ids = {
-      8,     45,    142,    165,    178,    204,   246,   253,   278,
-      280,   407,   542,    933,    1075,   1091,  1117,  1304,  1353,
-      1360,  1519,  1541,   1542,   1609,   2081,  3411,  7661,  7981,
-      8581,  10362, 11218,  12866,  14979,  22564, 25774, 49349, 57552,
-      79309, 80493, 122593, 130487, 247309, 302643};
+  std::vector<size_t> ids = {8,     45,    142,    178,    204,    246,   253,
+                             278,   280,   407,    542,    933,    1075,  1091,
+                             1117,  1304,  1353,   1360,   1519,   1541,  1542,
+                             1609,  2081,  3411,   7661,   7981,   8581,  10362,
+                             11218, 12866, 14979,  22564,  25774,  49349, 57552,
+                             79309, 80493, 122593, 130487, 247309, 302643};
   for (auto id : ids) {
     checkEvaluator(settings, id, "", EVAL_INCREMENTAL, true);
   }

src/form/expression.cpp

@@ -370,7 +370,9 @@ bool Expression::needsBrackets(bool isRoot, Expression::Type parentType) const {
     return false;
   }
   if (type == Expression::Type::FACTORIAL &&
-      parentType == Expression::Type::FACTORIAL) {
+      (parentType == Expression::Type::FACTORIAL ||
+       parentType == Expression::Type::PRODUCT ||
+       parentType == Expression::Type::SUM)) {
     return false;
   }
   if (parentType == Expression::Type::FACTORIAL && !isRoot) {

src/form/formula_gen.cpp

@@ -113,6 +113,43 @@ bool FormulaGenerator::bitfunc(Operation::Type type, const Expression& a,
   return true;
 }
 
+// Express falling/rising factorial using standard factorial
+bool FormulaGenerator::facToExpression(const Expression& a,
+                                             const Expression& b, Expression& res) const {
+
+  if (ExpressionUtil::canBeNegative(a, offset) ||
+      ExpressionUtil::canBeNegative(b, offset)) {
+    return false;
+  }
+  
+  Expression num(Expression::Type::FACTORIAL);
+  Expression denom(Expression::Type::FACTORIAL);
+  // Falling factorial: a!/(a+b)!
+  // If b <= 0: (a)!/(a+b)!
+  // If b > 0: rising factorial: (a+b-1)!/(a-1)!
+  if (b.type == Expression::Type::CONSTANT && b.value <= 0) {
+    // falling factorial
+    num.children = {a};
+    denom.children = {sum({a, b})};
+  } else {
+    // rising factorial
+    // (a+b-1)!/(a-1)!
+    num.children = {sum({a, sum({b, constant(-1)})})};
+    denom.children = {sum({a, constant(-1)})};
+  }
+  // simplify immediately
+  auto& d = denom.children.front();
+  ExpressionUtil::normalize(d);
+  if (d.type == Expression::Type::CONSTANT &&
+      (d.value == Number::ZERO || d.value == Number::ONE)) {
+    res =  num;
+  } else {
+    res =  divToFraction(num, denom);
+  }
+  
+  return true;
+}
+
 bool FormulaGenerator::update(const Operation& op) {
   auto source = operandToExpression(op.source);
   auto target = operandToExpression(op.target);
@@ -168,6 +205,10 @@ bool FormulaGenerator::update(const Operation& op) {
       res = func("binomial", {prevTarget, source});
       break;
     }
+    case Operation::Type::FAC: {
+      okay = facToExpression(prevTarget, source, res);
+      break;
+    }
     case Operation::Type::LOG: {
       res = func("logint", {prevTarget, source});
       break;

src/form/formula_gen.hpp

@@ -40,6 +40,9 @@ class FormulaGenerator {
   bool bitfunc(Operation::Type type, const Expression& a, const Expression& b,
                Expression& res) const;
 
+  bool facToExpression(const Expression& a, const Expression& b,
+               Expression& res) const;
+
   bool update(const Operation& op);
 
   bool update(const Program& p);

tests/formula/formula.txt

@@ -11,7 +11,7 @@ A000071: a(n) = b(n)-1, b(n) = b(n-1)+b(n-2), b(1) = 1, b(0) = 0
 A000073: a(n) = a(n-1)+a(n-2)+a(n-3), a(2) = 1, a(1) = 0, a(0) = 0
 A000096: a(n) = binomial(n+2,2)-1
 A000142: a(n) = n*a(n-1), a(0) = 1
-A000165: a(n) = b(2*n), b(n) = n*b(n-2), b(1) = 1, b(0) = 1
+A000165: a(n) = 2^n*n!
 A000168: a(n) = truncate((2*A151383(n))/(n+2)), A151383(n) = floor((floor(binomial(2*n,n)/(n+1))*3^(n+1))/3)
 A000212: a(n) = floor((n^2)/3)
 A000247: a(n) = 2^n-n-2
@@ -57,6 +57,7 @@ A001911: a(n) = b(n)-2, b(n) = b(n-1)+b(n-2), b(1) = 3, b(0) = 2
 A001923: a(n) = n^n+a(n-1), a(0) = 0
 A002015: a(n) = (n^2)%100
 A002309: a(n) = (2*n-1)^4+a(n-1), a(2) = 82, a(1) = 1, a(0) = 0
+A002378: a(n) = floor(((n+1)!)/((n-1)!))
 A002489: a(n) = n^(n^2)
 A002535: a(n) = 9*a(n-2)+2*a(n-1), a(2) = 11, a(1) = 1, a(0) = 1
 A004273: a(n) = max(2*n-1,0)

tests/programs/oeis/000/A000165.asm

@@ -1,10 +1,9 @@
 ; A000165: Double factorial of even numbers: (2n)!! = 2^n*n!.
 ; 1,2,8,48,384,3840,46080,645120,10321920,185794560,3715891200,81749606400,1961990553600,51011754393600,1428329123020800,42849873690624000,1371195958099968000,46620662575398912000,1678343852714360832000,63777066403145711616000,2551082656125828464640000,107145471557284795514880000,4714400748520531002654720000,216862434431944426122117120000,10409396852733332453861621760000,520469842636666622693081088000000,27064431817106664380040216576000000,1461479318123759876522171695104000000
 
-mov $1,1
-mul $0,2
-lpb $0
-  mul $1,$0
-  sub $0,2
-lpe
+mov $2,1
+fac $2,$0
+mov $1,2
+pow $1,$0
+mul $1,$2
 mov $0,$1

tests/programs/oeis/002/A002378.asm

@@ -1,6 +1,4 @@
 ; A002378: Oblong (or promic, pronic, or heteromecic) numbers: a(n) = n*(n+1).
 ; 0,2,6,12,20,30,42,56,72,90,110,132,156,182,210,240,272,306,342,380,420,462,506,552,600,650,702,756,812,870,930,992,1056,1122,1190,1260,1332,1406,1482,1560,1640,1722,1806,1892,1980,2070,2162,2256,2352,2450,2550,2652,2756,2862,2970,3080,3192,3306,3422,3540,3660,3782,3906,4032,4160,4290,4422,4556,4692,4830,4970,5112,5256,5402,5550,5700,5852,6006,6162,6320
 
-mov $1,$0
-add $0,1
-mul $0,$1
+fac $0,2

tests/programs/oeis/007/A007531.asm

@@ -0,0 +1,4 @@
+; A007531: a(n) = n*(n-1)*(n-2) (or n!/(n-3)!).
+; 0,0,0,6,24,60,120,210,336,504,720,990,1320,1716,2184,2730,3360,4080,4896,5814,6840,7980,9240,10626,12144,13800,15600,17550,19656,21924,24360,26970,29760,32736,35904,39270,42840,46620,50616,54834,59280,63960,68880,74046,79464,85140,91080,97290,103776,110544,117600,124950,132600,140556,148824,157410,166320,175560,185136,195054,205320,215940,226920,238266,249984,262080,274560,287430,300696,314364,328440,342930,357840,373176,388944,405150,421800,438900,456456,474474
+
+fac $0,-3