Dear GiNaC team, during debugging a complicated Jacobian generated with my FEM framework pyoomph via GiNaC, I was stuck, since the expressions were quite lenghty. I tracked down an issue with the C code generation of inverse powers ( x**^-1) when combining both an int (-1) and a float (-1.0) in the same expressions. It is now manually fixed in pyoomph y pre-parsing the expression, but of course, it should be patched in GiNaC itself. Below please find a code to reproduce it. Best regards, Christian PS: Thanks for GiNaC and CLN! It was the best decision to use it in my code! // GiNaC: print_csrc silently drops the reciprocal of a factor whose exponent is an inexact -1, i.e. // numeric(-1.0) rather than numeric(-1). // Seen on 1.8.7 (Ubuntu libginac-dev, shared, against libcln.so.6) and on 1.8.10 built from source. // // GiNaC compares numbers by value, so numeric(-1) and numeric(-1.0) are is_equal and hash alike; // the two expressions below are one and the same expression as far as every computation is // concerned. mul::do_print_csrc, however, asks two different questions about the exponent and gets // inconsistent answers (ginac/mul.cpp, do_print_csrc): // // // If the first argument is a negative integer power, it gets printed as "1.0/<expr>" // if (it == seq.begin() && it->coeff.info(info_flags::negint)) { ... c.s << "1.0/"; } // // // If the exponent is 1 or -1, it is left out // if (it->coeff.is_equal(_ex1) || it->coeff.is_equal(_ex_1)) // it->rest.print(c, precedence()); // ... // // Separator is "/" for negative integer powers, "*" otherwise // if (it->coeff.info(info_flags::negint)) c.s << "/"; else c.s << "*"; // // numeric::info(info_flags::negint) requires a CLN integer and is therefore false for -1.0, while // is_equal(_ex_1) compares by value and is therefore true. So neither the "1.0/" prefix nor the "/" // separator is emitted, and the exponent is left out all the same: x^(-1.0) inside a product is // printed as a plain multiplication by x. #include <ginac/ginac.h> #include <iostream> #include <sstream> #include <string> using namespace GiNaC; template <typename Context> static std::string csrc(const ex &e) { std::ostringstream os; e.print(Context(os)); return os.str(); } // `as_printed` is what the emitted C source literally says, written out as an expression again, so // that the two can be evaluated side by side. It is 0 where the output is correct. static void show(const std::string &what, const ex &e, const symbol &x, const ex &as_printed = 0) { std::cout << " " << what << "\n" << " expression : " << e << "\n" << " C source : " << csrc<print_csrc_double>(e) << "\n" << " expression @x=2 : " << e.subs(x == 2).evalf() << "\n"; if (!as_printed.is_zero()) std::cout << " that C @x=2 : " << as_printed.subs(x == 2).evalf() << " <-- WRONG\n"; std::cout << "\n"; } int main() { symbol x("x"); const ex exact = 3 * pow(1 + x, -1); // exponent is a CLN integer const ex inexact = 3 * pow(1 + x, numeric(-1.0)); // exponent is a CLN float const ex minus_one = numeric(-1), minus_one_point_zero = numeric(-1.0); std::cout << "The two exponents, and hence the two expressions, are equal to GiNaC:\n" << " -1 is_equal -1.0 = " << minus_one.is_equal(minus_one_point_zero) << "\n" << " equal hashes = " << (minus_one.gethash() == minus_one_point_zero.gethash()) << "\n" << " (exact - inexact).is_zero() = " << (exact - inexact).is_zero() << "\n" << " info(negint) of -1 and -1.0 = " << minus_one.info(info_flags::negint) << " and " << minus_one_point_zero.info(info_flags::negint) << " <-- the inconsistency\n\n"; std::cout << "They are not printed as the same C code:\n\n"; show("3*(1+x)^(-1) correct", exact, x); show("3*(1+x)^(-1.0) the reciprocal is gone", inexact, x, 3 * (1 + x)); std::cout << "Only an exponent of exactly -1 is affected, and only inside a product:\n\n"; show("(1+x)^(-1.0) correct - power::do_print_csrc handles it", pow(1 + x, numeric(-1.0)), x); show("3*(1+x)^(-2.0) correct - falls through to pow()", 3 * pow(1 + x, numeric(-2.0)), x); show("3*(1+x)^(-1/2) correct - not an integer either way", 3 * pow(1 + x, numeric(-1, 2)), x); show("x*(1+x)^(-1.0) the reciprocal is gone", x * pow(1 + x, numeric(-1.0)), x, x * (1 + x)); std::cout << "The other C source contexts print it the same way:\n" << " print_csrc_double : " << csrc<print_csrc_double>(inexact) << "\n" << " print_csrc_float : " << csrc<print_csrc_float>(inexact) << "\n" << " print_csrc_cl_N : " << csrc<print_csrc_cl_N>(inexact) << "\n" << "for comparison, the same expression with an exact exponent:\n" << " print_csrc_double : " << csrc<print_csrc_double>(exact) << "\n" << " print_csrc_float : " << csrc<print_csrc_float>(exact) << "\n" << " print_csrc_cl_N : " << csrc<print_csrc_cl_N>(exact) << "\n\n"; return 0; }