Fix PARI formula generation for rising factorial when base can be zero

CHANGELOG.md

@@ -10,6 +10,7 @@
 ### Bugfixes
 
 * Fix parsing error during stats generation
+* Fix PARI formula generation for rising factorial when the base can be zero (e.g., A006430)
 
 ### Enhancements
 

src/form/formula_gen.cpp

@@ -124,15 +124,19 @@ Expression FormulaGenerator::cellFunc(int64_t index) const {
   return ExpressionUtil::newFunction(getCellName(index));
 }
 
+// Maximum size of product expansion for factorial
+constexpr int64_t MAX_PRODUCT_EXPANSION = 10;
+
 // Express falling/rising factorial using standard factorial
 bool FormulaGenerator::facToExpression(const Expression& a, const Expression& b,
                                        const Operand& aOp, const Operand& bOp,
                                        const RangeMap* ranges,
                                        Expression& res) const {
-
-  // Optimize for small constants: convert to PRODUCT instead of FACTORIAL
-  if (b.type == Expression::Type::CONSTANT && b.value >= Number(-3) &&
-      b.value <= Number(3)) {
+  // Check if b is a constant within the product expansion threshold.
+  // This range check also ensures safe conversion to int64_t.
+  if (b.type == Expression::Type::CONSTANT &&
+      b.value >= Number(-MAX_PRODUCT_EXPANSION) &&
+      b.value <= Number(MAX_PRODUCT_EXPANSION)) {
     auto k = b.value.asInt();
 
     // Trivial case: k = 0 -> result is 1
@@ -147,7 +151,8 @@ bool FormulaGenerator::facToExpression(const Expression& a, const Expression& b,
       return true;
     }
 
-    // Build product expression
+    // Use product expansion for small constants.
+    // This avoids factorial division issues when a can be zero.
     res = Expression(Expression::Type::PRODUCT);
     if (k > 0) {
       // Rising factorial: a * (a+1) * (a+2) * ... * (a+k-1)
@@ -171,13 +176,14 @@ bool FormulaGenerator::facToExpression(const Expression& a, const Expression& b,
     }
     return true;
   }
+
   // TODO: can we relax the negativity check for b?
   if (canBeNegativeWithRanges(a, aOp, ranges) ||
       canBeNegativeWithRanges(b, bOp, ranges)) {
     return false;
   }
 
-  // General case
+  // General case: use factorial division
   Expression num(Expression::Type::FACTORIAL);
   Expression denom(Expression::Type::FACTORIAL);
 
@@ -191,6 +197,7 @@ bool FormulaGenerator::facToExpression(const Expression& a, const Expression& b,
   } else {
     // rising factorial
     // (a+b-1)!/(a-1)!
+    // Note: we already checked that a-1 >= 0 above for small constants
     num.children = {sum({a, sum({b, constant(-1)})})};
     denom.children = {sum({a, constant(-1)})};
   }

tests/formula/formula.txt

@@ -77,6 +77,7 @@ A005408: a(n) = 2*n+1
 A005600: a(n) = (3*min(n-3,6)-((n%4)==0)+1)%10
 A005843: a(n) = 2*n
 A006221: a(n) = n*(n*(34*n+51)+27)+5
+A006430: a(n) = truncate((floor(((-(1==n)+n)*(-(1==n)+n+1)*(-(1==n)+n+2)*(-(1==n)+n+3))/(-(1==n)+n+1))*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*(23*n-23*(1==n)+963)+17544)+147952)+481675)-1052153)-7850914)-2900162)+60869272)+37067400)-179920800))/79833600)
 A007395: a(n) = 2
 A007531: a(n) = n*(n-2)*(n-1)
 A007583: a(n) = 2*floor((4^n)/3)+1

tests/formula/pari-function.txt

@@ -5,3 +5,4 @@ A000058: (a(n) = b(n)+1); (b(n) = if(n==0,1,local(l1=b(n-1)); l1*(l1+1)))
 A000247: a(n) = 2^n-n-2
 A001286: a(n) = floor(((n-1)*n!)/2)
 A001715: (a(n) = truncate(b(n)/6)); (b(n) = if(n==0,1,n*b(n-1)))
+A006430: a(n) = truncate((floor(((-(1==n)+n)*(-(1==n)+n+1)*(-(1==n)+n+2)*(-(1==n)+n+3))/(-(1==n)+n+1))*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*((-(1==n)+n)*(23*n-23*(1==n)+963)+17544)+147952)+481675)-1052153)-7850914)-2900162)+60869272)+37067400)-179920800))/79833600)

tests/programs/oeis/006/A006430.asm

@@ -0,0 +1,36 @@
+; A006430: Number of loopless tree-rooted planar maps with 5 vertices and n faces and no isthmuses.
+; 0,5,360,7350,73700,474588,2292790,9046807,30676440,92393015,252872984,639382605,1512137536,3377126024,7176513960,14599539314,28575632350,54036739617,99069119952,176618150000,306965183268,521265871700,866527603370,1412513294049,2261193633300,3559534488211,5516605536480,8426235973615,12696729614240,18889491891600,27768822919308,40365602495540,58058144150922,82674136440925,116618332018840,163031500997274,225987148221380,310733619012236,423990500605854,574309684127035,772513103084624
+
+#offset 1
+
+mov $3,1
+equ $3,$0
+sub $0,$3
+mov $1,$0
+mov $2,$0
+add $2,1
+fac $0,4
+div $0,$2
+sub $2,1
+mul $1,23
+add $1,963
+mul $1,$2
+add $1,17544
+mul $1,$2
+add $1,147952
+mul $1,$2
+add $1,481675
+mul $1,$2
+sub $1,1052153
+mul $1,$2
+sub $1,7850914
+mul $1,$2
+sub $1,2900162
+mul $1,$2
+add $1,60869272
+mul $1,$2
+add $1,37067400
+mul $1,$2
+sub $1,179920800
+mul $0,$1
+div $0,79833600