[GiNaC-list] "Internal error..." report -- probably misused map.

10 Aug 2006

      Hi,

I am working on a program that generates code for other libraries
or programs.  One of the target libraries is CVODE, a solver for
differential equations (part of the SUNDIALS suite).

In CVODE, numbers are represented with a typedef 'realtype', which
can be chosen to be float, double, or one or two other options.
To ensure consistent use of realtype, constants are wrapped in the
macro RCONST. I am trying to automate this in ginac.  I am using
ginac version 1.3.4 (compiled with gcc 4.0.2 on Ubuntu Linux).

In some of my first attempts to do this, I occasionally got an
"internal error" message.  The following code is a simplified
version that causes the error:

--- ginacwrap1.cpp --------------------------------------------
#include <iostream>
#include <ginac/ginac.h>

using namespace std;
using namespace GiNaC;

DECLARE_FUNCTION_1P(RCONST)
REGISTER_FUNCTION(RCONST,dummy())

ex wrap1(const ex & e)
    {
    if (is_a<numeric>(e))
        return RCONST(e);
    else
        return e.map(wrap1);
    }

int main(void)
    {
    symbol x("x");

    ex p = cos(1+x);   // Internal error: ... line 93 has been reached!
    // ex p = 2*x;        // Internal error: ... line 30 has been reached!
    // ex p = 2+sin(2);   // Works: wrap1(p) = RCONST(2)+sin(RCONST(2))
    ex w = wrap1(p);
    cout << "wrap1(p) = " << w << endl;
    }
---------------------------------------------------------------
This results in an "internal error":

  Internal error: statement in file ./real/elem/cl_R_mul.cc, line 93 has been reached!!
  Please send the authors of the program a description how you produced this error!

If the expression p is changed to 2*x (see the comments), a similar
internal error is reported.

After more work, I ended up with the code shown below, which seems to
do what I want (but any suggestions for a better method would be
appreciated).  I had to define my own print function for the RCONST
function, because the default print function converted the name to
lower case when csrc was used.  One remaining nuisance:  In csrc mode,
ginac prints the expression -x as "-1.0*x".  Is there a way to have it
output "-x" in csrc mode?

Best regards,

Warren Weckesser

--- ginacwrap2.cpp --------------------------------------------
#include <iostream>
#include <ginac/ginac.h>

using namespace std;
using namespace GiNaC;

static void RCONST_print_csrc_double(const ex & arg, const print_context & c)
     {
     c.s << "RCONST(";
     arg.print(c);
     c.s << ")";
     }

DECLARE_FUNCTION_1P(RCONST)
REGISTER_FUNCTION(RCONST,print_func<print_csrc_double>(RCONST_print_csrc_double))

ex wrap2(const ex& e)
    {
    if (is_a<numeric>(e))
        return RCONST(e);
    if (is_a<symbol>(e))
        return e;
    if (is_a<add>(e))
        {
        ex s = 0;
        for (size_t k = 0; k < e.nops(); ++k)
            s = s + wrap2(e.op(k));
        return s;
        }
    if (is_a<mul>(e))
        {
        ex p = 1;
        for (size_t k = 0; k < e.nops(); ++k)
            p = p * wrap2(e.op(k));
        return p;
        }
    return e.map(wrap2);
    }

int main(void)
    {
    symbol x("x");

    cout << csrc;
    ex p = -exp(-x)+x*x;
    cout << "p = " << p << endl;
    ex w = wrap2(p);
    cout << "wrap2(p) = " << w << endl;
    }
----------------------------------------------------------------