diff --git a/CHANGELOG.md b/CHANGELOG.md index 971ad1aa..db7ca6bd 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/form/formula_gen.cpp b/src/form/formula_gen.cpp index dfd3295d..6c897625 100644 --- a/src/form/formula_gen.cpp +++ b/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)})}; } diff --git a/tests/formula/formula.txt b/tests/formula/formula.txt index a3507da7..eea92521 100644 --- a/tests/formula/formula.txt +++ b/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 diff --git a/tests/formula/pari-function.txt b/tests/formula/pari-function.txt index b0f337f4..b2852522 100644 --- a/tests/formula/pari-function.txt +++ b/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) diff --git a/tests/programs/oeis/006/A006430.asm b/tests/programs/oeis/006/A006430.asm new file mode 100644 index 00000000..29ba57fc --- /dev/null +++ b/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