use PARI wrappers in plain formulas (#313)

Author: ckrause

CHANGELOG.md

@@ -8,6 +8,7 @@ To install or update LODA, please follow the [installation instructions](https:/
 
 ### Enhancements
 
+* Use floor/truncate functions in plain formulas
 * Simplify formulas using Gaussian elimination
 * Simplify sums and products in formulas
 * Extended internal maintenance command

src/form/expression_util.cpp

@@ -298,8 +298,10 @@ bool ExpressionUtil::canBeNegative(const Expression& e) {
           std::any_of(e.children.begin(), e.children.end(),
                       [](const Expression& c) { return !canBeNegative(c); })) {
         return false;
+      } else if (e.name == "binomial") {
+        break;  // infer from children
       } else {
-        return true;
+        return true;  // unknown function
       }
     case Expression::Type::POWER:
       if (e.children.size() == 2 &&

src/form/formula_gen.cpp

@@ -57,6 +57,18 @@ Expression FormulaGenerator::operandToExpression(Operand op) const {
   throw std::runtime_error("internal error");  // unreachable
 }
 
+Expression divToFraction(const Expression& numerator,
+                         const Expression& denominator) {
+  Expression frac(Expression::Type::FRACTION, "", {numerator, denominator});
+  std::string func = "floor";
+  if (ExpressionUtil::canBeNegative(numerator) ||
+      ExpressionUtil::canBeNegative(denominator)) {
+    func = "truncate";
+  }
+  Expression wrapper(Expression::Type::FUNCTION, func, {frac});
+  return wrapper;
+}
+
 bool FormulaGenerator::update(const Operation& op) {
   auto source = operandToExpression(op.source);
   auto target = operandToExpression(op.target);
@@ -89,15 +101,34 @@ bool FormulaGenerator::update(const Operation& op) {
       break;
     }
     case Operation::Type::DIV: {
-      res = Expression(Expression::Type::FRACTION, "", {prevTarget, source});
+      res = divToFraction(prevTarget, source);
       break;
     }
     case Operation::Type::POW: {
       res = Expression(Expression::Type::POWER, "", {prevTarget, source});
+      if (ExpressionUtil::canBeNegative(source)) {
+        Expression wrapper(Expression::Type::FUNCTION, "truncate", {res});
+        res = wrapper;
+      }
       break;
     }
     case Operation::Type::MOD: {
-      res = Expression(Expression::Type::MODULUS, "", {prevTarget, source});
+      auto c1 = prevTarget;
+      auto c2 = source;
+      if (ExpressionUtil::canBeNegative(c1) ||
+          ExpressionUtil::canBeNegative(c2)) {
+        Expression wrapper(Expression::Type::SUM);
+        wrapper.newChild(c1);
+        wrapper.newChild(Expression::Type::PRODUCT);
+        auto frac = divToFraction(c1, c2);
+        wrapper.children[1].newChild(
+            Expression(Expression::Type::CONSTANT, "", Number(-1)));
+        wrapper.children[1].newChild(c2);
+        wrapper.children[1].newChild(frac);
+        res = wrapper;
+      } else {
+        res = Expression(Expression::Type::MODULUS, "", {c1, c2});
+      }
       break;
     }
     case Operation::Type::BIN: {

src/form/pari.cpp

@@ -43,31 +43,7 @@ bool convertExprToPari(Expression& expr) {
       return false;
     }
   }
-  if (expr.type == Expression::Type::FRACTION) {
-    convertFracToPari(expr);
-  } else if (expr.type == Expression::Type::POWER) {
-    if (ExpressionUtil::canBeNegative(expr.children.at(1))) {
-      Expression wrapper(Expression::Type::FUNCTION, "truncate", {expr});
-      expr = wrapper;
-    }
-  } else if (expr.type == Expression::Type::MODULUS) {
-    auto c1 = expr.children.at(0);
-    auto c2 = expr.children.at(1);
-    if (ExpressionUtil::canBeNegative(c1) ||
-        ExpressionUtil::canBeNegative(c2)) {
-      Expression wrapper(Expression::Type::SUM);
-      wrapper.newChild(c1);
-      wrapper.newChild(Expression::Type::PRODUCT);
-      Expression frac(Expression::Type::FRACTION, "", {c1, c2});
-      convertFracToPari(frac);
-      wrapper.children[1].newChild(
-          Expression(Expression::Type::CONSTANT, "", Number(-1)));
-      wrapper.children[1].newChild(c2);
-      wrapper.children[1].newChild(frac);
-      expr = wrapper;
-    }
-  } else if (expr.type == Expression::Type::FUNCTION &&
-             expr.name == "binomial") {
+  if (expr.type == Expression::Type::FUNCTION && expr.name == "binomial") {
     // TODO: check feedback from PARI team to avoid this limitation
     if (ExpressionUtil::canBeNegative(expr.children.at(1))) {
       return false;

tests/formula/program-formula.txt

@@ -1,5 +1,5 @@
 A000004: a(n) = 0
-A000008: a(n) = b(n+1), b(n) = b(n-1)+binomial((n+1)/2-((2*n+2)/5)+1,2), b(0) = 0
+A000008: a(n) = b(n+1), b(n) = b(n-1)+binomial(-floor((2*n+2)/5)+floor((n+1)/2)+1,2), b(0) = 0
 A000012: a(n) = 1
 A000023: a(n) = (n-1)*(2*a(n-2)+a(n-1))-a(n-1), a(3) = -2, a(2) = 2, a(1) = -1, a(0) = 1
 A000027: a(n) = n+1
@@ -12,9 +12,9 @@ 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
-A000168: a(n) = (2*A151383(n))/(n+2), A151383(n) = ((binomial(2*n,n)/(n+1))*3^(n+1))/3
-A000212: a(n) = (n^2)/3
-A000276: a(n) = (2*c(n+1)*(n+3))/2, b(n) = b(n-1)*(n+1), b(2) = 6, b(1) = 2, b(0) = 1, c(n) = c(n-1)*(n+1)+b(n-1), c(2) = 5, c(1) = 1, c(0) = 0
+A000168: a(n) = truncate((2*A151383(n))/(n+2)), A151383(n) = truncate((floor(binomial(2*n,n)/(n+1))*3^(n+1))/3)
+A000212: a(n) = floor((n^2)/3)
+A000276: a(n) = truncate((2*c(n+1)*(n+3))/2), b(n) = b(n-1)*(n+1), b(2) = 6, b(1) = 2, b(0) = 1, c(n) = c(n-1)*(n+1)+b(n-1), c(2) = 5, c(1) = 1, c(0) = 0
 A000278: a(n) = (-a(n-2))^2+a(n-1), a(1) = 1, a(0) = 0
 A000280: a(n) = a(n-2)^3+a(n-1), a(1) = 1, a(0) = 0
 A000284: a(n) = a(n-1)^3+a(n-2), a(3) = 2, a(2) = 1, a(1) = 1, a(0) = 0
@@ -23,17 +23,17 @@ A000290: a(n) = n^2
 A000295: a(n) = 2^n-n-1
 A000321: a(n) = -a(n-2)*(2*n-2)-a(n-1), a(2) = -1, a(1) = -1, a(0) = 1
 A000407: a(n) = a(n-1)*(4*n+2), a(1) = 6, a(0) = 1
-A000463: a(n) = gcd(n+1,n/2+1)*(n/2+1)
+A000463: a(n) = gcd(n+1,floor(n/2)+1)*(floor(n/2)+1)
 A000466: a(n) = 4*n^2-1
 A000538: a(n) = n^4+a(n-1), a(0) = 0
 A000561: a(n) = (n+3)*(binomial(-3*n-1,2)+1)
-A000601: a(n) = b(n+4), b(n) = ((n-2)^2)/4+b(n-3), b(2) = 0, b(1) = 0, b(0) = 0
+A000601: a(n) = b(n+4), b(n) = b(n-3)+floor(((n-2)^2)/4), b(2) = 0, b(1) = 0, b(0) = 0
 A000748: a(n) = -3*a(n-1)-3*a(n-2), a(2) = 6, a(1) = -3, a(0) = 1
 A000792: a(n) = min(n+1,(n+1)%3)*c(n+1)+b(n+1), b(n) = 2*b(n-3)+2*c(n-3), b(5) = 2, b(4) = 2, b(3) = 2, b(2) = 1, b(1) = 1, b(0) = 1, c(n) = b(n-3)+c(n-3), c(5) = 1, c(4) = 1, c(3) = 1, c(2) = 0, c(1) = 0, c(0) = 0
 A000918: a(n) = 2^n-2
 A000931: a(n) = a(n-2)+a(n-3), a(5) = 1, a(4) = 0, a(3) = 1, a(2) = 0, a(1) = 0, a(0) = 1
 A000932: a(n) = n*a(n-2)+a(n-1), a(2) = 3, a(1) = 1, a(0) = 1
-A000933: a(n) = (b(max(n-2,0))+5)/6, b(n) = b(n-1)+n-1, b(0) = 0
+A000933: a(n) = truncate((b(max(n-2,0))+5)/6), b(n) = b(n-1)+n-1, b(0) = 0
 A001044: a(n) = A000142(n)^2, A000142(n) = n*A000142(n-1), A000142(0) = 1
 A001054: a(n) = a(n-2)*a(n-1)-1, a(3) = -2, a(2) = -1, a(1) = 1, a(0) = 0
 A001107: a(n) = n*(4*n-3)
@@ -48,7 +48,7 @@ A001542: a(n) = 3*a(n-1)+2*b(n-1), a(1) = 2, a(0) = 0, b(n) = 4*a(n-1)+3*b(n-1),
 A001611: a(n) = A000045(n)+1, A000045(n) = A000045(n-1)+A000045(n-2), A000045(1) = 1, A000045(0) = 0
 A001630: a(n) = 2*a(n-1)-a(n-5), a(6) = 12, a(5) = 6, a(4) = 3, a(3) = 2, a(2) = 1, a(1) = 0, a(0) = 0
 A001687: a(n) = a(n-2)+a(n-5), a(5) = 1, a(4) = 0, a(3) = 1, a(2) = 0, a(1) = 1, a(0) = 0
-A001715: a(n) = b(n+3)/6, b(n) = n*b(n-1), b(0) = 1
+A001715: a(n) = truncate(b(n+3)/6), b(n) = n*b(n-1), b(0) = 1
 A001906: a(n) = 3*a(n-1)-a(n-2), a(3) = 8, a(2) = 3, a(1) = 1, a(0) = 0
 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
@@ -59,14 +59,14 @@ A004273: a(n) = max(2*n-1,0)
 A004277: a(n) = max(2*n,1)
 A005408: a(n) = 2*n+1
 A005843: a(n) = 2*n
-A007583: a(n) = 2*((4^n)/3)+1
+A007583: a(n) = 2*floor((4^n)/3)+1
 A008785: a(n) = (n+4)^n
 A008999: a(n) = 2*a(n-1)+a(n-4), a(4) = 17, a(3) = 8, a(2) = 4, a(1) = 2, a(0) = 1
 A014731: a(n) = 4*b(n)^2, b(n) = 4*b(n-1)+b(n-2), b(1) = -2, b(0) = -1
 A021019: a(n) = 6*min(n,1)
 A022322: a(n) = b(n+2), b(n) = c(n-1), b(2) = 1, b(1) = 6, b(0) = 0, c(n) = c(n-1)+c(n-2)+1, c(3) = 10, c(2) = 8, c(1) = 1, c(0) = 6
 A022958: a(n) = -n+2
-A037536: a(n) = b(n+1), b(n) = (3*c(n-1))/13, b(1) = 1, b(0) = 0, c(n) = 3*c(n-1), c(1) = 24, c(0) = 8
+A037536: a(n) = b(n+1), b(n) = truncate((3*c(n-1))/13), b(1) = 1, b(0) = 0, c(n) = 3*c(n-1), c(1) = 24, c(0) = 8
 A048745: a(n) = 2*a(n-1)+a(n-2)+3, a(1) = 5, a(0) = 1
 A062815: a(n) = b(n+1), b(n) = n*n^n+b(n-1), b(0) = 0
 A078013: a(n) = b(n-2), a(2) = 0, a(1) = 0, a(0) = 1, b(n) = -b(n-3)+b(n-1), b(2) = -1, b(1) = -1, b(0) = 0
@@ -76,7 +76,7 @@ A123111: a(n) = ((n+1)^5+(n+1)^3+n+2)*(n+1)^2+1
 A134271: a(n) = b(max(n-2,0))+min(n-2,0), b(n) = 2*b(n-3)+b(n-2), b(2) = 4, b(1) = 0, b(0) = 2
 A154118: a(n) = 2*a(n-1)+5, a(1) = 2, a(0) = 1
 A157823: a(n) = (min(n-1,0)-1)*(b(max(n-1,0))+3)+3, b(n) = 2*b(n-1), b(0) = 1
-A173673: a(n) = (b(n)+min(n,n%2))/2, b(n) = b(n-2)+b(n-4), b(4) = 1, b(3) = 1, b(2) = 1, b(1) = 0, b(0) = 0
+A173673: a(n) = truncate((b(n)+min(n,n%2))/2), b(n) = b(n-2)+b(n-4), b(4) = 1, b(3) = 1, b(2) = 1, b(1) = 0, b(0) = 0
 A275698: a(n) = a(n-1)*(a(n-1)-3)+3, a(2) = 13, a(1) = 5, a(0) = 2
 A293006: a(n) = b(max(n-1,0)), b(n) = 2*b(n-1)+2*b(n-2)+4, b(1) = 2, b(0) = 0
 A318791: a(n) = (3*n-39)*(3*n-38)+41