Latex Equation Compiler Code

/*************************************************************************
    func.cpp  -  class func, extending the GiNaC function clas
    *** added in 0.5
                             -------------------
    begin                : Fri Feb 13 2004
    copyright            : (C) 2004 by Jan Rheinlaender
    email                : jrheinlaender@users.sourceforge.net
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include "../config/config.h" // *** added in 1.0
// Note: this include must come before "func.h" for cross-compiling with mingw!
#include <cmath> // *** added in 1.2
#include "func.h"
#include "printing.h"
#include "msgdriver.h" // *** added 1.3.1
#include "equation.h" // *** added in 1.0
#include "eqc.h" // *** added in 1.3.0
//#include "../config/config.h" // *** added in 1.3.1
#include <sstream>
#ifndef _MSC_VER
  #include <stdint.h> // *** added in 1.3.1 for hash function
#endif
//#include "optstack.h" // *** added in 1.2

// *** changed basic to exprseq in 1.4.1
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(func, exprseq, print_func<print_context>(&func::do_print).print_func<imathprint>(&func::do_print_imath));

// Define static members------------------------------------------------------------------
// Note that the map functions must be defined BEFORE func_initializer, because func::init()
// accesses this map!
  std::map<const std::string, funcrec*> func::functions; // *** changed to funcmap in 1.0, to std::map in 1.4.2
  std::map<const std::string, std::string> func::hard_names;
  std::map<const std::string, std::string> func::hard_names_rev; // *** added in 1.0
  std::map<std::string, unsigned> func::hint_map = std::map<std::string, unsigned>();
  std::map<std::string, std::string> func::func_inv = std::map<std::string, std::string>();
  bool func_init::called = false;
  func_init func_init::func_initializer;
//  diff_t func::diff_type = dfdt; // *** added in 0.6, *** replaced by option o_difftype in 1.4.1

func_init::func_init() {
  if (called) {
    // Removed in 1.3.1, msg_error might not be initialized yet 
    //    msg_error.prio(0) << "Attempt to call func::init() more than one time!" << endline;
    return;
  } else
    called = true;
  func::init();
}

// Required constructors and destructors and other GiNaC-specific methods-----------------
// Default constructor etc.
// *** changed in 1.2.1
//func::func() : inherited(TINFO_func), name("") {} // func::func()
func::func() : name("") {
  //tinfo_key = &func::tinfo_static; *** removed in 1.3.1
} // func::func

void func::checkargs() { // *** made separate function in 0.6
  if (!seq.empty()) {
    if ((functions[name] != NULL) && (seq.size() != functions[name]->vars.size())) {
      msg::warn(0) << "Warning: Number of arguments does not match for " << name
       << ". Using default arguments." << endline;
      seq = functions[name]->vars;
    }
    numargs = seq.size();
  } else {
    if (functions[name] == NULL) {
      msg::warn(0) << "Warning: Number of arguments not known for " << name
        << ". Setting to 1" << endline;
      numargs = 1;
    } else {
      numargs = functions[name]->vars.size();
    }
  }

  MSG_INFO(3) << "Constructor of function " << name << endline;
} // func::checkargs()

func::func(const std::string &n, const exvector &args, bool discardable)
    : exprseq(args, discardable), name(n) {
//  tinfo_key = &func::tinfo_static; // *** changed in 1.2.1, *** removed in 1.3.1
  checkargs();
} // func::func()

func::func(const std::string &n, const expression &e)
    : exprseq(make_exvector(e)), name(n) {
//  tinfo_key = &func::tinfo_static; // *** changed in 1.2.1, *** removed in 1.3.1
  checkargs();
} // func::func()

func::func(const std::string &n, const exprseq &es)
    : exprseq(es), name(n) {
//  tinfo_key = &func::tinfo_static; // *** changed in 1.2.1, *** removed in 1.3.1
  checkargs();
} // func::func()

func::func(const std::string &n, exvector &v, bool discardable)
    : exprseq(v, discardable), name(n) { // *** changed for GiNaC 1.2.1 in 0.7
//  tinfo_key = &func::tinfo_static; // *** changed in 1.2.1, *** removed in 1.3.1
  checkargs();
} // func::func()

func::func(const std::string &n, std::auto_ptr<exvector> vp)
  : exprseq(vp), name(n) { // *** added in 0.7
//  tinfo_key = &func::tinfo_static; // *** changed TINFO_function to TINFO_func in 0.9, changed in 1.2.1, *** removed in 1.3.1
  checkargs();
} // func::func()

func::func(const std::string &n, const function &f)
  : name(n) {
  // TODO: Is it a problem that we can't do ": exprseq(...)" here?
  // *** changed to const_iterator in 0.7
  for (func::const_iterator i = f.begin(); i != f.end(); i++) seq.push_back(*i);
//  tinfo_key = &func::tinfo_static; // *** changed in 1.2.1, *** removed in 1.3.1
  checkargs();
} // func::func()

// *** removed const qualifier for sym_lst in 0.7, removed in 1.3.1
/*func::func(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) {
  ex temp;
  n.find_string("func.name", name);
  //n.find_ex("func.seq", temp, sym_lst); seq = ex_to<lst>(temp);
  n.find_ex("func.definition", temp, sym_lst); functions[name]->definition = temp;
}
*/

// *** removed const qualifier for sym_lst in 0.7, changed to read_archive() format in 1.3.1
void func::read_archive(const archive_node &n, lst &sym_lst) {
  inherited::read_archive(n, sym_lst);
  ex temp;
  n.find_string("func.name", name);
  //n.find_ex("func.seq", temp, sym_lst); seq = ex_to<lst>(temp);
  n.find_ex("func.definition", temp, sym_lst); functions[name]->definition = temp;
//  return (new func(n, sym_lst))->setflag(status_flags::dynallocated);
}

void func::archive(archive_node &n) const {
  inherited::archive(n);
  n.add_string("func.name", name);
  //n.add_ex("func.seq", lst(seq));
  n.add_ex("func.definition", functions[name]->definition); // *** changed in 0.9
}
GINAC_BIND_UNARCHIVER(func);

void func::printseq(const print_context & c, const std::string &openbracket, char delim,
    const std::string &closebracket, unsigned this_precedence, unsigned upper_precedence) const {
  // modified from exprseq::printseq to use print_ltx()
  if (this_precedence <= upper_precedence)
    c.s << openbracket;

  if (!seq.empty()) {
    exvector::const_iterator it = seq.begin();
    while (it != seq.end()) {
      if (is_a<print_latex>(c)) {
        print_ltx(*it, c.s);
      } else {
        it->print(c, this_precedence);
      }
      if (it < seq.end() - 1) c.s << delim;
      ++it;
    }
  }

  if (this_precedence <= upper_precedence)
    c.s << closebracket;
} // func::printseq()

void func::do_print(const print_context &c, unsigned level) const {
  // *** changed name to do_print in 0.8
  print(c, 1, level);
} // func::do_print()

void func::do_print_imath(const imathprint &c, unsigned level) const {
  // *** added in 1.4.1 Is this necessary or could do_print differentiate the methods?
  print_imath(c, 1, level);
} // func::do_print_imath()

// Added for workaround *** 1.4.1
static void mindex_print_imath(const ex &e, const ex &row, const ex &col, const imathprint &c);
static void partialdiff_print_imath(const ex &e, const ex &s, const ex &n, const imathprint &c);
static void sum_print_imath(const ex &lower, const ex &higher, const ex &e, const imathprint &c);

// *** added in 1.4.1
void func::print_imath(const imathprint &c, const ex &p, unsigned level) const {
    if ((functions[name] != NULL) && functions[name]->hard &&
        (functions[name]->hints & hint_map["print"])) {
      // Fall through to the hard-coded print function
      // TODO: The hardcoded print function for \ln prints \log!!
      // Work-around because REGISTER_FUNCTION does not accept mindex_print_imath() etc
      if (name == "mindex") {
        mindex_print_imath(seq[0], seq[1], seq[2], c);
      } else if ((name == "diff") || (name == "\\diff")) {
        partialdiff_print_imath(seq[0], seq[1], seq[2], c);
      } else if ((name == "sum") || (name == "\\sum")) {
        sum_print_imath(seq[0], seq[1], seq[2], c);
      } else {
        unsigned s = functions[name]->serial;
        function(s, seq).print(c);
      }
      return;
    }
    if (p != 1) { // Print the exponent of a trigonometric function directly after the function name
      if (c.poptions->get(o_eqparse).boolean) { // Option eqparse disables this feature
        c.s << "{" << name;
        if (seq.empty()) { c.s << "}"; return; }
      } else {
        c.s << name << "^";
        print_ltx(p, c.s);
        if (seq.empty()) return;
      }
    } else {
      c.s << name;
      if (seq.empty()) return;
    }
    c.s << "{";

    if (functions[name] == NULL) {
      printseq(c, " left(", ',', " right)");
    } else { // Some functions get special treatment, depending on their hints
      if (functions[name]->hints & hint_map["nobracket"]) { // do not print brackets
        printseq(c, "", ',', "");
      } else if ((functions[name]->hints & hint_map["trig"]) &&
                 (nops() == 1) && (is_a<numeric>(op(0)) || is_a<symbol>(op(0)) || is_a<power>(op(0)) ||
                              is_a<Unit>(op(0)) || is_a<func>(op(0)))) { // do not print brackets
          print_ltx(seq[0], c.s);
      } else {
        printseq(c, " left(", ',', " right)");
      }
    }
    c.s << "}";
    if ((p != 1) && c.poptions->get(o_eqparse).boolean) {
      // Print the exponent of a trigonometric function. Option eqparse enables this feature
      c.s << "}^";
      print_ltx(p, c.s);
    }
} // func::print(imathprint)

void func::print(const print_context &c, const ex &p, unsigned level) const {
  // Note that the only time when p!=1 will be when this function is called directly
  // from print_ltx() and the function is trigonometric
  if (is_a<print_tree>(c)) {
    c.s << std::string(level, ' ') << class_name() << " " << name
        << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
        << ", nops=" << nops() << std::endl;
    unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
    for (unsigned i=0; i < seq.size(); ++i)
            seq[i].print(c, level + delta_indent);
    c.s << std::string(level + delta_indent, ' ') << "=====" << std::endl;
  } else if (is_a<print_csrc>(c)) {
    // Print function name in lowercase
    std::string lname = name;
    for (unsigned i = 0; i < lname.size(); i++) lname[i] = tolower(lname[i]);
    c.s << lname << "(";
    // Print arguments, separated by commas
    // TODO: use func::printseq() here?
    exvector::const_iterator it = seq.begin();
    while (it != seq.end()) {
      it->print(c);
      ++it;
      if (it != seq.end()) c.s << ",";
    }
    c.s << ")";
  } else if (is_a<print_latex>(c)) {
    if ((functions[name] != NULL) && functions[name]->hard &&
        (functions[name]->hints & hint_map["print"])) { // *** added in 1.0
      // Fall through to the hard-coded print function
      // TODO: The hardcoded print function for \ln prints \log!!
      unsigned s = functions[name]->serial;
      function(s, seq).print(c);
      return;
    }
    if (p != 1) { // Print the exponent of a trigonometric function directly after the function name
      if (optstack::options->get(o_eqparse).boolean) { // Option eqparse disables this feature *** added in 1.0
        c.s << "{" << name;
        if (seq.empty()) { c.s << "}"; return; }
      } else {
        c.s << name << "^";
        print_ltx(p, c.s);
        if (seq.empty()) return;
      }
    } else {
      c.s << name;
      if (seq.empty()) return;
    }
    c.s << "{";

    if (functions[name] == NULL) {
      printseq(c, "\\left(", ',', "\\right)");
    } else { // Some functions get special treatment, depending on their hints
      if (functions[name]->hints & hint_map["nobracket"]) { // do not print brackets
        printseq(c, "", ',', "");
      } else if ((functions[name]->hints & hint_map["trig"]) &&
                 (nops() == 1) && (is_a<numeric>(op(0)) || is_a<symbol>(op(0)) || is_a<power>(op(0)) ||
                              is_a<Unit>(op(0)) || is_a<func>(op(0)))) { // do not print brackets
          print_ltx(seq[0], c.s);
      } else {// *** changed brackets to \left( and \right) in 0.6
        printseq(c, "\\left(", ',', "\\right)");
      }
    }
    c.s << "}";
    if ((p != 1) && optstack::options->get(o_eqparse).boolean) {
      // Print the exponent of a trigonometric function
      // Option eqparse enables this feature *** added in 1.0
      c.s << "}^";
      print_ltx(p, c.s);
    }
  } else {
    c.s << name;
    printseq(c, "(", ',', ")", exprseq::precedence(), func::precedence());
  }
} // func::print(print_context)

// *** added in 1.4.1
void func::print_diff_line(const unsigned ndiff, const print_context& c) const {
  // TODO: All this is mostly copied from func::print()
  if ((functions[name] != NULL) && functions[name]->hard && (functions[name]->hints & hint_map["print"]))
    msg::warn(0) << "Warning: Cannot use hardcoded print function for difftype line" << endline;

  c.s << name;
  for (unsigned i = 0; i < ndiff; i++) c.s << "'";
  if (seq.empty()) return;
   
  c.s << "{";

  if (functions[name] == NULL) {
    is_a<print_latex>(c) ? printseq(c, "\\left(", ',', "\\right)") : printseq(c, " left(", ',', " right)");
  } else { // Some functions get special treatment, depending on their hints
    if (functions[name]->hints & hint_map["nobracket"]) { // do not print brackets
      printseq(c, "", ',', "");
    } else if ((functions[name]->hints & hint_map["trig"]) &&
               (nops() == 1) && (is_a<numeric>(op(0)) || is_a<symbol>(op(0)) || is_a<power>(op(0)) ||
                            is_a<Unit>(op(0)) || is_a<func>(op(0)))) { // do not print brackets
        print_ltx(seq[0], c.s);
    } else {// *** changed brackets to \left( and \right) in 0.6
      is_a<print_latex>(c) ? printseq(c, "\\left(", ',', "\\right)") : printseq(c, " left(", ',', " right)");
    }
  }
  c.s << "}";
} // func::print_diff_line()

ex func::expand_definition() const throw (std::runtime_error, std::logic_error) {
  // *** removed argument in 0.6
  if (functions[name] == NULL)
    throw (std::runtime_error("Warning: Function " + name + " is unknown! Please register it first!"));

  if (functions[name]->definition.is_empty() && !functions[name]->hard) {
    throw(std::logic_error("Warning: Function " + name + " has no definition! Not expanded."));
  }

  // *** since 0.6, size mismatch in arguments can not happen any more
  if (seq.empty()) { // There is nothing to expand *** changed in 0.6
    if (functions[name]->hard) // *** changed in 0.9
      return this->setflag(status_flags::expanded);
    else
      return functions[name]->definition;
  }

  // *** 0.6: This call used 'seq', but is that correct?? Removed expand() and call_func()
  if (functions[name]->hard)
    return(func(name, seq).setflag(status_flags::expanded));
    // TODO: Why not this->setflag(...)?

  // Substitute the arguments in the definition
  // *** changed to use exmap and subs_options::algebraic in 0.8
  exmap subst;
  for (unsigned i = 0; i < seq.size(); i++)
    subst[functions[name]->vars[i]] = seq[i];
  ex result = functions[name]->definition.subs(subst, subs_options::algebraic);
  if (msg::info().checkprio(1)) {
    msg::info() << "Expanding user function " << name << "(" << subst << ") = "
                     << functions[name]->definition << endline;
    msg::info()  << "Expansion result: " << result << endline;
  }
  return(result);
} // func::expand_definition()

ex func::expand(unsigned options) const {
  MSG_INFO(2)  << "Expanding " << name << endline;
  func result;
  
  // *** changed for mindex in 1.4.1
  if (options & expand_options::expand_function_args) { // Only expand arguments when asked to do so
    if (name == "mindex") { // Special treatment, because we don't want the indixes to became floats...
      std::cout << "Expanding mindex: " << *this << std::endl;
      result = *this;
    } else {
      result = ex_to<func>(inherited::expand(options));
    }
  } else {
    result = *this;
  }
  
  if ((functions[name] != NULL) && !functions[name]->definition.is_empty()) { // is there a function definition that can be expanded?
    try {
      return (result.expand_definition());//TODO: should we do another expand() on this result,
        // since the definition itself has not been expanded yet?
    } catch (std::runtime_error &e) {
      msg::error(0) << e.what() << endline;
    } catch (std::logic_error &e) {
      msg::error(0)  << e.what() << endline;
    }
    return result;
  } else {
    return (options == 0) ? result.setflag(status_flags::expanded) : ex(result);
  }
} // func::expand()

ex func::eval(int level) const {
  // If the function is hardcoded, drop through to the GiNaC::function::eval() method
  MSG_INFO(3) << "Doing eval of " << name << endline;
  if ((functions[name] != NULL) && functions[name]->hard) {
    if (seq.size() == 1) {
      if (is_exactly_a<func>(seq[0]) && functions[ex_to<func>(seq[0]).name]->hard) { // *** added in 0.9, added 'hard' in 1.0
        // Take advantage of the hard-coded GiNac eval rules, e.g. tan(atan(x)) = x
        ex func_arg = function(functions[ex_to<func>(seq[0]).name]->serial, ex_to<func>(seq[0]).seq);
        ex result = function(functions[name]->serial, func_arg).eval(level);
        if (is_a<function>(result)) return this->hold(); // Nothing happened in eval...
        return result; // TODO: Could the result possibly be a GiNaC function???
      } else if ((name == "\\ln") && is_a<constant>(seq[0]) && (ex_to<constant>(seq[0]) == Euler_number)) { // handle ln e = 1 *** added in 1.3.0
        return numeric(1);
      }
    } 
    if (!seq.empty()) { // *** added this distinction in 0.9
      // Take advantage of the hard-coded GiNaC eval rules, e.g. sin(-2) = -sin(2)
      ex result = function(functions[name]->serial, seq).eval(level);
      return (is_a<function>(result) ? func(name, ex_to<function>(result)).hold() : result);
      // Don't introduce any GiNaC::function into the system!
    } else {
      return (this->hold());
    }
  }

  func result = *this;
  if (level > 1) result.seq = evalchildren(level);

  //const function_options &opt = functions[name]->opts;
  // Canonicalize argument order according to the symmetry properties
  // Not possible because symtree is declared protected in class function_options!
  //  if ((result.seq.size() > 1) && !(opt.symtree.is_zero())) {
  //    exvector v = result.seq;
  //    GINAC_ASSERT(is_a<symmetry>(opt.symtree));
  //    int sig = canonicalize(v.begin(), ex_to<symmetry>(opt.symtree()));
  //    if (sig != INT_MAX) {
  //      // Something has changed while sorting arguments, more evaluations later
  //      if (sig == 0) return(ex(0));
  //      return ex(sig) * thisexprseq(v);
  //    }
  //  }

  if (functions[name]->hints & hint_map["expand"])  // Expand function immediately *** added in 0.6
    if (!functions[name]->definition.is_empty()) {
      return result.expand_definition(); // no full expansion! Just the function, not the arguments
    }

  // Nothing else happens here, since user-defined eval functions are not implemented!
  return result.hold(); // This sets statusflags::evaluated!
} // func::eval()

ex func::reduce_double_funcs::operator()(const ex &e) { // *** added in 1.0
  MSG_INFO(2) << "Reducing pairs of function/inverse function in " << e << endline;
    
  if (is_a<function>(e)) {
    function f = ex_to<function>(e);
    std::cout << "Found G-function " << f.get_name() << ". This should not happen!!!!" << std::endl;
  }
  if (is_a<func>(e)) {
    func f = ex_to<func>(e);
        
    if (f.seq.size() == 0) return e; // Nothing to reduce *** changed to seq.size() in 1.0
    if (f.numargs > 1) return e.map(*this); // map children

      
    if (is_a<func>(f.seq[0])) { // The argument is a function
      func argf = ex_to<func>(ex_to<func>(f.seq[0]).map(*this));
      std::string fnamebare = f.name; // *** added in 1.4.1 for iMath
      std::string argfnamebare = argf.name;
      if (fnamebare[0] == '\\') fnamebare.erase(0,1);
      if (argfnamebare[0] == '\\') argfnamebare.erase(0,1);

      if (func_inv[fnamebare] == argfnamebare) {
        if (argf.seq.size() > 1) { // The function has several arguments
          return lst(argf.seq.begin(), argf.seq.end());
        } else {
          return argf.seq[0];
        }
      }
    }
  }

  return e.map(*this);
} // reduce_double_funcs::operator()

ex func::replace_function_by_func::operator()(const ex &e) { // *** added in 1.0
  MSG_INFO(2) << "Replacing all GiNaC functions by EQC funcs in " << e << endline;

  if (is_a<function>(e)) {
    std::string fname = hard_names_rev[ex_to<function>(e).get_name()]; // *** added in 1.4.1 for iMath
    if (!is_a_func(fname)) fname = "\\" + fname;
    ex result = func(fname, ex_to<function>(e));
    return result.map(*this);
  }
  return e.map(*this);
} // replace_functions_by_func::operator()

ex func::evalf (int level) const throw(std::runtime_error) { // *** removed argument args in 0.6
  MSG_INFO(3)  << "Evaluating " << *this << endline;
  // If the function is hardcoded, drop through to the GiNaC::function::evalf() method
  if ((functions[name] != NULL) && functions[name]->hard) {
    ex result = function(functions[name]->serial, seq).evalf(level);
    return (is_a<function>(result) ? func(name, ex_to<function>(result)) :  result);
  }

  // Evaluate children first
  exvector eseq;
  if (level == 1)
    eseq = seq;
  else if (level == -max_recursion_level)
    throw(std::runtime_error("max recursion level reached"));
  else
    eseq.reserve(seq.size());

  --level;
  exvector::const_iterator it = seq.begin();
  while (it != seq.end()) {
    eseq.push_back(it->evalf(level));
    ++it;
  }

  if (functions[name] != NULL) {
    if (!functions[name]->definition.is_empty()) {
      // Note: The code in GiNaC::functions.cpp::evalf() seems to drop eseq here!
      //if (level == 1) // evaluate only one level. Is this correct?
      //  return (expand_definition(args));
      try {
        return (func(name, eseq).expand_definition().evalf(level));
      } catch (std::exception &e) { // The function cannot be expanded
        msg::error(0)  << e.what() << endline;
      }
    }
  }

  return func(name, eseq).hold();
} // evalf()

// The hash functions are copied from GiNaCs utils.h and hash_seed.h and crc32.h *** in 1.0 and again in 1.3.1
#ifdef HAVE_STDINT_H
typedef uintptr_t p_int;
#else
typedef unsigned long p_int;
#endif
/**
 * CRC32 hash function. Shamelessly stolen from GNU coreutils (cksum.c)
 */
static unsigned const crctab[256] =
{
  0x00000000,
  0x04c11db7, 0x09823b6e, 0x0d4326d9, 0x130476dc, 0x17c56b6b,
  0x1a864db2, 0x1e475005, 0x2608edb8, 0x22c9f00f, 0x2f8ad6d6,
  0x2b4bcb61, 0x350c9b64, 0x31cd86d3, 0x3c8ea00a, 0x384fbdbd,
  0x4c11db70, 0x48d0c6c7, 0x4593e01e, 0x4152fda9, 0x5f15adac,
  0x5bd4b01b, 0x569796c2, 0x52568b75, 0x6a1936c8, 0x6ed82b7f,
  0x639b0da6, 0x675a1011, 0x791d4014, 0x7ddc5da3, 0x709f7b7a,
  0x745e66cd, 0x9823b6e0, 0x9ce2ab57, 0x91a18d8e, 0x95609039,
  0x8b27c03c, 0x8fe6dd8b, 0x82a5fb52, 0x8664e6e5, 0xbe2b5b58,
  0xbaea46ef, 0xb7a96036, 0xb3687d81, 0xad2f2d84, 0xa9ee3033,
  0xa4ad16ea, 0xa06c0b5d, 0xd4326d90, 0xd0f37027, 0xddb056fe,
  0xd9714b49, 0xc7361b4c, 0xc3f706fb, 0xceb42022, 0xca753d95,
  0xf23a8028, 0xf6fb9d9f, 0xfbb8bb46, 0xff79a6f1, 0xe13ef6f4,
  0xe5ffeb43, 0xe8bccd9a, 0xec7dd02d, 0x34867077, 0x30476dc0,
  0x3d044b19, 0x39c556ae, 0x278206ab, 0x23431b1c, 0x2e003dc5,
  0x2ac12072, 0x128e9dcf, 0x164f8078, 0x1b0ca6a1, 0x1fcdbb16,
  0x018aeb13, 0x054bf6a4, 0x0808d07d, 0x0cc9cdca, 0x7897ab07,
  0x7c56b6b0, 0x71159069, 0x75d48dde, 0x6b93dddb, 0x6f52c06c,
  0x6211e6b5, 0x66d0fb02, 0x5e9f46bf, 0x5a5e5b08, 0x571d7dd1,
  0x53dc6066, 0x4d9b3063, 0x495a2dd4, 0x44190b0d, 0x40d816ba,
  0xaca5c697, 0xa864db20, 0xa527fdf9, 0xa1e6e04e, 0xbfa1b04b,
  0xbb60adfc, 0xb6238b25, 0xb2e29692, 0x8aad2b2f, 0x8e6c3698,
  0x832f1041, 0x87ee0df6, 0x99a95df3, 0x9d684044, 0x902b669d,
  0x94ea7b2a, 0xe0b41de7, 0xe4750050, 0xe9362689, 0xedf73b3e,
  0xf3b06b3b, 0xf771768c, 0xfa325055, 0xfef34de2, 0xc6bcf05f,
  0xc27dede8, 0xcf3ecb31, 0xcbffd686, 0xd5b88683, 0xd1799b34,
  0xdc3abded, 0xd8fba05a, 0x690ce0ee, 0x6dcdfd59, 0x608edb80,
  0x644fc637, 0x7a089632, 0x7ec98b85, 0x738aad5c, 0x774bb0eb,
  0x4f040d56, 0x4bc510e1, 0x46863638, 0x42472b8f, 0x5c007b8a,
  0x58c1663d, 0x558240e4, 0x51435d53, 0x251d3b9e, 0x21dc2629,
  0x2c9f00f0, 0x285e1d47, 0x36194d42, 0x32d850f5, 0x3f9b762c,
  0x3b5a6b9b, 0x0315d626, 0x07d4cb91, 0x0a97ed48, 0x0e56f0ff,
  0x1011a0fa, 0x14d0bd4d, 0x19939b94, 0x1d528623, 0xf12f560e,
  0xf5ee4bb9, 0xf8ad6d60, 0xfc6c70d7, 0xe22b20d2, 0xe6ea3d65,
  0xeba91bbc, 0xef68060b, 0xd727bbb6, 0xd3e6a601, 0xdea580d8,
  0xda649d6f, 0xc423cd6a, 0xc0e2d0dd, 0xcda1f604, 0xc960ebb3,
  0xbd3e8d7e, 0xb9ff90c9, 0xb4bcb610, 0xb07daba7, 0xae3afba2,
  0xaafbe615, 0xa7b8c0cc, 0xa379dd7b, 0x9b3660c6, 0x9ff77d71,
  0x92b45ba8, 0x9675461f, 0x8832161a, 0x8cf30bad, 0x81b02d74,
  0x857130c3, 0x5d8a9099, 0x594b8d2e, 0x5408abf7, 0x50c9b640,
  0x4e8ee645, 0x4a4ffbf2, 0x470cdd2b, 0x43cdc09c, 0x7b827d21,
  0x7f436096, 0x7200464f, 0x76c15bf8, 0x68860bfd, 0x6c47164a,
  0x61043093, 0x65c52d24, 0x119b4be9, 0x155a565e, 0x18197087,
  0x1cd86d30, 0x029f3d35, 0x065e2082, 0x0b1d065b, 0x0fdc1bec,
  0x3793a651, 0x3352bbe6, 0x3e119d3f, 0x3ad08088, 0x2497d08d,
  0x2056cd3a, 0x2d15ebe3, 0x29d4f654, 0xc5a92679, 0xc1683bce,
  0xcc2b1d17, 0xc8ea00a0, 0xd6ad50a5, 0xd26c4d12, 0xdf2f6bcb,
  0xdbee767c, 0xe3a1cbc1, 0xe760d676, 0xea23f0af, 0xeee2ed18,
  0xf0a5bd1d, 0xf464a0aa, 0xf9278673, 0xfde69bc4, 0x89b8fd09,
  0x8d79e0be, 0x803ac667, 0x84fbdbd0, 0x9abc8bd5, 0x9e7d9662,
  0x933eb0bb, 0x97ffad0c, 0xafb010b1, 0xab710d06, 0xa6322bdf,
  0xa2f33668, 0xbcb4666d, 0xb8757bda, 0xb5365d03, 0xb1f740b4
};

static unsigned crc32(const char* c, const unsigned len, const unsigned crcinit)
{
	unsigned crc = crcinit ^ 0xFFFFFFFFU;
	for (unsigned i = 0; i < len; ++i)
		crc = (crc << 8) ^ crctab[((crc >> 24) ^ c[i]) & 0xFF];
	crc = ~crc & 0xFFFFFFFFU;
	return crc;
}

/**
 * We need a hash function which gives different values for objects of
 * different types. Hence we need some unique integer for each type.
 * Fortunately, standard C++ RTTI class `type_info' stores a pointer to
 * mangled type name. Normally this pointer is the same for all objects
 * of the same type (although it changes from run to run), so it can be
 * used for computing hashes. However, on some platforms (such as woe32)
 * the pointer returned by type_info::name() might be different even for
 * objects of the same type! Hence we need to resort to comparing string
 * representation of the (mangled) type names. This is quite expensive,
 * so we compare crc32 hashes of those strings. We might got more hash
 * collisions (and slower evaluation as a result), but being a bit slower
 * is much better than being wrong.
 */
#ifdef _WIN32
#define GINAC_HASH_USE_MANGLED_NAME 1
#endif
#ifndef GINAC_HASH_USE_MANGLED_NAME
unsigned golden_ratio_hash(p_int n);
unsigned make_hash_seed(const std::type_info& tinfo)
{
	// this pointer is the same for all objects of the same type.
	// Hence we can use that pointer 
	const void* mangled_name_ptr = (const void*)tinfo.name();
	unsigned v = golden_ratio_hash((p_int)mangled_name_ptr);
	return v;
}
#else
unsigned make_hash_seed(const std::type_info& tinfo)
{
	const char* mangled_name = tinfo.name();
	return crc32(mangled_name, std::strlen(mangled_name), 0);
}
#endif

/** Rotate bits of unsigned value by one bit to the left.
  * This can be necesary if the user wants to define its own hashes. */
unsigned rotate_left(unsigned n)
{
        return (n & 0x80000000U) ? (n << 1 | 0x00000001U) : (n << 1);
}
 
unsigned golden_ratio_hash(p_int n)
{
        // This function works much better when fast arithmetic with at
        // least 64 significant bits is available.
        if (sizeof(long) >= 8) {
        // So 'long' has 64 bits.  Excellent!  We prefer it because it might be
        // more efficient than 'long long'.
        unsigned long l = n * 0x4f1bbcddUL;
        return (unsigned)l;
        }
#ifdef HAVE_LONG_LONG
        else if (sizeof(long long) >= 8) {
        // This requires 'long long' (or an equivalent 64 bit type)---which is,
        // unfortunately, not ANSI-C++-compliant.
        // (Yet C99 demands it, which is reason for hope.)
        unsigned long long l = n * 0x4f1bbcddULL;
        return (unsigned)l;
        }
#endif   
        // Without a type with 64 significant bits do the multiplication manually
        // by splitting n up into the lower and upper two bytes.
        const unsigned n0 = (n & 0x0000ffffU);
        const unsigned n1 = (n & 0xffff0000U) >> 16;
        return (n0 * 0x0000bcddU) + ((n1 * 0x0000bcddU + n0 * 0x00004f1bU) << 16);
}


unsigned func::calchash(void) const {
  // *** copied from GiNac::function in 1.0 and adapted in 1.3.1
  unsigned ser;
  if (functions[name] == NULL) {
    //msg::error(0) << "Error: Function is not registered" << endline;
    ser = 0;
  } else
    ser = functions[name]->serial;

  unsigned v = golden_ratio_hash(make_hash_seed(typeid(*this)) ^ ser); // *** changed in 1.2.1 to user (unsigned int) = (p_int)
  for (size_t i=0; i<nops(); i++) {
    v = rotate_left(v);
    v ^= this->op(i).gethash();
  }

  if (flags & status_flags::evaluated) {
    setflag(status_flags::hash_calculated);
    hashvalue = v;
  }
  MSG_INFO(3) << "Hash value of " << name << " (serial " << ser << "): " << v << endline;
  return v;
} // func::calchash()

ex func::thiscontainer(const exvector &v) const {
  return func(name, v);
} // func::thisexprseq()

ex func::thiscontainer(std::auto_ptr<exvector> vp) const {
  return func(name, vp);
} // func::thisexprseq()

ex func::series(const relational &r, int order, unsigned options) const {
  // Not implemented yet. Returns basic::series
  return basic::series(r, order, options); // *** added options in 1.0
//if (registered_functions()[serial].series_f==0) return basic::series(r, order);

//  ex res;
//  current_serial = serial;
//  if (registered_functions()[serial].series_use_exvector_args) {
//    try {
//      res = ((series_funcp_exvector)(registered_functions()[serial].series_f))(seq, r, order
//    } catch (do_taylor) {
//      res = basic::series(r, order, options);
//    }
//    return res;
//  }
} // func::series()

ex func::derivative(const symbol & s) const throw (std::runtime_error) { // *** added in 0.6
  if (functions[name] == NULL) {
    throw std::runtime_error("Function " + name + " is not registered yet!");
  }
  MSG_INFO(2) << "Calculating derivative of " << *this << " to " << s << endline;

  if (name == "\\diff") { // *** added in 1.0
    if (ex_to<symbol>(seq[1]) == s)
      return func("\\diff", partialdiff(seq[0], seq[1], seq[2] + 1));
  } else if (name == "diff") { // *** added in 1.4.1 for iMath
    if (ex_to<symbol>(seq[1]) == s)
      return func("diff", partialdiff(seq[0], seq[1], seq[2] + 1));
  }

  if (functions[name]->hard) { // Fall through to GiNaC::function::diff(). derivative() is protected!
    ex result = function(functions[name]->serial, seq).diff(s);
    replace_function_by_func replace_functions;
    return replace_functions(result); // *** changed in 1.0
  }

  ex result;
  ex arg_diff;

  if (seq.size() == 0) { // *** added in 1.0, removed ex_to<symbol> in 1.0
    if (name == s.get_name()) return 1; // Handle case dr() / dr()
    for (unsigned i = 0; i < numargs; i++) {
      result += pderivative(i) * functions[name]->vars[i].diff(s);
      MSG_INFO(2) << "Derivative of " << *this << " stage " << i << ": " << result << endline;
    }
  } else {
    for (unsigned i = 0; i < seq.size(); i++) {
      arg_diff = seq[i].diff(s);
      if (!arg_diff.is_zero()) result += pderivative(i)*arg_diff;
    }
  }
  return result;
} // func::derivative()

int func::compare_same_type(const basic &other) const {
  // This function is important because it is used to simplify expressions!
  // If (*this == other) then the expression "*this / other" will be 1.
  // TODO: What is the point of comparing > and < ?
  MSG_INFO(3)  << "Comparing " << *this << " and " << other << endline;
  const func &o = static_cast<const func &>(other);
  if (name == o.name) { // *** added diff_num in 0.6
    return exprseq::compare_same_type(o);
  }
  else if (name < o.name)
    return -1;
  else
    return 1;
}

bool func::is_equal_same_type(const basic & other) const { // *** added in 0.6
  // This function seems to be used for substitution: If a.is_equal_same_type(b),
  // then b==c will be substituted
  MSG_INFO(3)  << "Checking equality of " << *this << " and " << other << endline;
  const func & o = static_cast<const func &>(other);

  if (name != o.name)
    return false;
  else
    return (exprseq::is_equal_same_type(o));
} // func::is_equal_same_type()

bool func::match_same_type(const basic &other) const { // *** added in 0.6
  MSG_INFO(3) << "Checking match of " << *this << " and " << other << endline;
  const func & o = static_cast<const func &>(other);
  return name == o.name;
}

unsigned func::return_type(void) const { // *** added in 0.6
  MSG_INFO(3) << "Return type of " << *this << " requested." << endline;
  if (seq.empty())
    return return_types::commutative;
  else
    return seq.begin()->return_type();
}

/*tinfo_t func::return_type_tinfo(void) const { // *** added in 0.6, changed to tinfo_t in 1.2.1, removed in 1.3.1
  MSG_INFO(3)  << "tinfo of " << *this << " requested." << endline;
  if (seq.empty())
    return tinfo_key;
  else
    return seq.begin()->return_type_tinfo();
} // func::return_type_tinfo()
*/

ex func::pderivative(unsigned diff_param) const throw(std::logic_error, std::runtime_error) { // *** added in 0.6
  if (functions[name] == NULL) {
    throw std::runtime_error("Function " + name + " is not registered yet!");
  }
  MSG_INFO(2) << "Calculating partial derivative of " << *this <<
                        " to " << functions[name]->vars[diff_param] << endline;

  ex result;
  // We assume that no hardcoded functions are called with this method!
  if (functions[name]->definition.is_empty() || !(functions[name]->hints & hint_map["defdiff"])) { // *** added defdiff in 1.4.1
    if (is_a_func("\\diff")) { // ***changed for iMath
      result = func("\\diff", partialdiff(*this, functions[name]->vars[diff_param], 1));
    } else {
      result = func("diff", partialdiff(*this, functions[name]->vars[diff_param], 1));
    }
  } else {
    if (diff_param >= numargs)
      throw std::logic_error("The requested dependant variable does not exist in " + name + "()");

    ex diffvar = functions[name]->vars[diff_param]; // *** added diffvar in 1.0
    if (is_a<symbol>(diffvar)) {
      result = functions[name]->definition.diff(ex_to<symbol>(diffvar));
    } else { // The diffvar is a (pure) function
      result = diff_to_func(functions[name]->definition, ex_to<func>(diffvar));
    }

    if (!seq.empty()) { // Substitute the function arguments in the result *** added in 1.0
      exmap subs_map;
      for (unsigned i = 0; i < numargs; i++) subs_map[functions[name]->vars[i]] = seq[i];
      result = result.subs(subs_map);
    }
  }

  MSG_INFO(2) << "Partial derivative #" << diff_param << " of " << *this << ": "
                     << result << endline;
  return result;
} // func::pderivative()

//--------------------------------------------------------------------------------------

bool func::initialized = false;


void func::remove(const std::string& fname) { // *** added in 1.4.1
  functions[fname] = NULL;
  functions.erase(fname);
} // func::delete()

void func::clear() {
  // *** changed to funcmap in 1.0, added fdel in 1.0, to std::map in 1.4.2
  std::map<const std::string, funcrec*>::iterator it_func, fdel;
  MSG_INFO(2) << "Clearing functions..." << endline;
  if (msg::info().checkprio(3)) { // *** added in 1.0
    msg::info() << "List of functions" << endline;
    for (it_func = functions.begin(); it_func != functions.end(); it_func++)
      if (it_func->second != NULL) {
        msg::info() << it_func->first << " = ";
        if (it_func->second->hard)
          msg::info() << "[hard-coded]";
        else
          msg::info() << it_func->second->definition;
        msg::info() << endline;
      }
  }

  it_func = functions.begin();
  while (it_func != functions.end()) {
    if ((it_func->second != NULL) && (!it_func->second->hints & hint_map["lib"])) { // *** changed in 1.0
      MSG_INFO(3) << "Deleting " << it_func->first << endline;
      //delete(it_func->second); // TODO: ERASE_HACK *** added in 1.0
      it_func->second = NULL;
      fdel = it_func; // *** added in 1.0
      it_func++;
      functions.erase(fdel);
    } else
      it_func++;
  }
} // func::clear()

void func::clearall() {
  functions.clear(); // *** changed to clear() in 0.8
} // func::clearall()

void func::print_functions() {
  std::map<const std::string, funcrec*>::iterator it_func = functions.begin(); // *** changed to funcmap in 1.0, to std::map in 1.4.2
  while (it_func != functions.end()) {
    MSG_INFO(0) << it_func->first << "(" << it_func->second->vars << ")";
    if (!it_func->second->hard || (it_func->second->definition.is_empty())) // *** changed in 0.6
      MSG_INFO(0) << " = " << it_func->second->definition;
    MSG_INFO(0) << " serial: " << it_func->second->serial << endline; // *** added in 0.6
    it_func++;
  }
} // print_functions()

const unsigned func::hint(const std::string &s) { // *** added in 0.6
    std::string str = s;
    
    if (s == "no_bracket")  {
      msg::warn(0) << "Warning: Replace deprecated function hint 'no_bracket' with 'nobracket'" << endline;
      str = "nobracket"; // *** added in 1.3.1, changed syntax, but handles legacy files
    }
    
    if (hint_map[str] == 0) {
      msg::warn(0) << "Warning: The hint " << str << " is not defined." << endline;
      return(0);
    } else {
      return hint_map[str];
    }
  } // func::hint()

  const unsigned max_unsigned = 1 << (sizeof(unsigned) * 8 - 1); // *** TODO: 1 << ...*8 gives warning

void func::registr(const std::string &n, const lst &args, const symbol &defarg, const unsigned h)
    throw(std::invalid_argument) {
  // *** added hints in 0.6
  if (functions[n] != NULL)
    throw (std::invalid_argument("Function " + n + " already exists"));

  functions[n] = new funcrec;
  functions[n]->hints = h;

  try { // Is the function hard-coded into GiNaC?
    std::string nbare = n; // *** added in 1.4.1 for iMath
    if (nbare[0] == '\\') nbare.erase(0,1);
    functions[n]->serial = function::find_function(
      ((hard_names[nbare] == "") ? n : hard_names[nbare]), args.nops());
    functions[n]->hard = true;
    MSG_INFO(1) << "Function " <<
      ((hard_names[nbare] == "") ? n : hard_names[nbare]) << " is hard-coded into GiNaC." << endline;
  } catch (std::exception &e) {
    functions[n]->serial = max_unsigned - functions.size(); // *** changed in 1.0
    functions[n]->hard = false;
  }

  unsigned j = 0;
  for (lst::const_iterator i = args.begin(); i != args.end(); i++) {
    // *** changed to const_iterator in 0.7, arguments may now be (pure) functions since 1.0
    if (is_a<symbol>(*i) || (is_a<func>(*i) && ex_to<func>(*i).is_pure())) {
      functions[n]->vars.push_back(*i);
    } else {
      msg::warn(0) << "Warning: Argument " << j+1 << " of " << n << " is no symbol! "
           << "Using argument 'x' instead." << endline;
      functions[n]->vars.push_back(defarg);
    }
    j++;
  }
  MSG_INFO(1) << "Registered function " << n <<
    "(" << functions[n]->vars << ((h & hint_map["lib"]) ? "), which is a library function." : ").")  << endline;
} // func::registr()

void func::define(const std::string &n, const expression &def)
    throw(std::invalid_argument) {
  if (functions[n] == NULL)
   throw (std::invalid_argument("Function " + n + " does not exist! Please register it first."));

  // TODO: handle the case that the function has been differentiated!
  functions[n]->definition = def;
  MSG_INFO(1) << "Defined function " << n <<
    "(" << functions[n]->vars << ")  = " << def << endline;
} // define()

const bool func::is_a_func(const std::string &fname) {
  return (functions[fname] != NULL);
}

const bool func::is_trig() const throw(std::runtime_error) {
  if (functions[name] == NULL)
    throw std::runtime_error("Function is not registered");
  return (functions[name]->hints & hint_map["trig"]);
} // is_trig()

const bool func::is_lib() const throw(std::runtime_error) {
  if (functions[name] == NULL)
    throw std::runtime_error("Function is not registered");
  return (functions[name]->hints & hint_map["lib"]);
} // is_lib()

// ------------------------------------------------------------------------------
// hard-coded functions
// *** round, ceil and floor added in 1.2

static ex round_eval(const ex &e, const ex &n) throw(std::runtime_error) {
  if (!is_a<numeric>(e.evalf()))
    return func("\\round", round(e, n)).hold();
  numeric num = ex_to<numeric>(e.evalf());
  if (!num.info(info_flags::real))
    throw std::runtime_error("Can only round real numbers");
  
	numeric digits;
  if (!is_a<numeric>(n)) {
		ex ndig = n.evalf(); // Maybe n becomes a numeric now? Note: Doing evalf() right away converts integers to floats
		
		if (!is_a<numeric>(ndig)) {
    	throw std::runtime_error("Number of digits to round to must be an integer");
		} else {
			digits = ex_to<numeric>(ndig);
		}
	} else {
	  digits = ex_to<numeric>(n);
	}
  // *** changed nonnegint to nonnegative in 1.2, decimal exponents don't hurt pow()...
	if (!digits.info(info_flags::nonnegative))
    throw std::runtime_error("Number of digits to round to must be a positive integer");
  cln::cl_F m = std::pow(double(10), digits.to_double());
  return numeric(truncate1(cln::the<cln::cl_R>(num.to_cl_N()) * m + csgn(num) * cln::cl_F(0.5)) / m);
} // round_eval()

static ex round_deriv(const ex &e, const ex &n, unsigned deriv_param) {
  msg::error(0) << "Error: \\round cannot be differentiated" << endline;
  return func("\\round", round(e, n)).hold();
} // round_deriv()

static void round_print_ltx(const ex &e, const ex &n, const print_context &c) {
  c.s << "\\round{";
  print_ltx(e, c.s);
  c.s << ", ";
  print_ltx(n, c.s);
  c.s << "}";
} // round_print_ltx()

REGISTER_FUNCTION(round, eval_func(round_eval).
                       print_func<print_latex>(round_print_ltx).
                       derivative_func(round_deriv));

static ex floor_eval(const ex &e, const ex &n) throw(std::runtime_error) {
  if (!is_a<numeric>(e.evalf()))
    return func("\\floor", floor(e, n)).hold();
  numeric num = ex_to<numeric>(e.evalf());
  if (!num.info(info_flags::real))
    throw std::runtime_error("Can only floor real numbers");
		
	numeric digits;
  if (!is_a<numeric>(n)) {
		ex ndig = n.evalf(); // Maybe n becomes a numeric now? Note: Doing evalf() right away converts integers to floats
		
		if (!is_a<numeric>(ndig)) {
    	throw std::runtime_error("Number of digits to floor to must be an integer");
		} else {
			digits = ex_to<numeric>(ndig);
		}
	} else {
	  digits = ex_to<numeric>(n);
	}
 // *** changed nonnegint to nonnegative in 1.2, decimal exponents don't hurt pow()... 
	if (!digits.info(info_flags::nonnegative))
    throw std::runtime_error("Number of digits to floor to must be a positive integer");
  cln::cl_F m = std::pow(double(10), digits.to_double());
  return numeric(floor1(cln::the<cln::cl_R>(num.to_cl_N()) * m) / m);
} // floor_eval()

static ex floor_deriv(const ex &e, const ex &n, unsigned deriv_param) {
  msg::error(0) << "Error: \\floor cannot be differentiated" << endline;
  return func("\\floor", floor(e, n)).hold();
} // floor_deriv()

static void floor_print_ltx(const ex &e, const ex &n, const print_context &c) {
  c.s << "\\floor{";
  print_ltx(e, c.s);
  c.s << ", ";
  print_ltx(n, c.s);
  c.s << "}";
} // floor_print_ltx()

REGISTER_FUNCTION(floor, eval_func(floor_eval).
                       print_func<print_latex>(floor_print_ltx).
                       derivative_func(floor_deriv));

static ex ceil_eval(const ex &e, const ex &n) throw(std::runtime_error) {
  if (!is_a<numeric>(e.evalf()))
    return func("\\ceil", ceil(e, n)).hold();
  numeric num = ex_to<numeric>(e.evalf());
  if (!num.info(info_flags::real))
    throw std::runtime_error("Can only ceil real numbers");
		
	numeric digits;
  if (!is_a<numeric>(n)) {
		ex ndig = n.evalf(); // Maybe n becomes a numeric now? Note: Doing evalf() right away converts integers to floats
		
		if (!is_a<numeric>(ndig)) {
    	throw std::runtime_error("Number of digits to ceil to must be an integer");
		} else {
			digits = ex_to<numeric>(ndig);
		}
	} else {
	  digits = ex_to<numeric>(n);
	}
	
// *** changed nonnegint to nonnegative in 1.2, decimal exponents don't hurt pow()...	
    if (!digits.info(info_flags::nonnegative))
          throw std::runtime_error("Number of digits to ceil to must be a positive integer");
  cln::cl_F m = std::pow(double(10), digits.to_double());
  return numeric(ceiling1(cln::the<cln::cl_R>(num.to_cl_N()) * m) / m);
} // ceil_eval()

static ex ceil_deriv(const ex &e, const ex &n, unsigned deriv_param) {
  msg::error(0) << "Error: \\ceil cannot be differentiated" << endline;
  return func("\\ceil", ceil(e, n)).hold();
} // ceil_deriv()

static void ceil_print_ltx(const ex &e, const ex &n, const print_context &c) {
  c.s << "\\ceil{";
  print_ltx(e, c.s);
  c.s << ", ";
  print_ltx(n, c.s);
  c.s << "}";
} // ceil_print_ltx()

REGISTER_FUNCTION(ceil, eval_func(ceil_eval).
                       print_func<print_latex>(ceil_print_ltx).
                       derivative_func(ceil_deriv));

// Sum function *** added in 1.0
static ex sum_eval(const ex &lower, const ex &higher, const ex &e) {
  if (!is_a<relational>(lower))
    throw std::runtime_error("Lower bound must be an equation");
  if (!is_a<symbol>((ex_to<relational>(lower)).lhs()))
    throw std::runtime_error("Lower bound must assign a value to a symbol");
  if (lower.is_equal(higher)) return e;
  // TODO: check if higher < lower
  if (e.is_zero()) return 0;
  if (func::is_a_func("\\sum")) { // added in 1.4.1 for imath
    return func("\\sum", sum(lower, higher, e)).hold();
  } else {
    return func("sum", sum(lower, higher, e)).hold();
  }
} // sum_eval()

static ex sum_deriv(const ex &lower, const ex &higher, const ex &e, unsigned deriv_param) {
  msg::error(0) << "Error: sum operator cannot be differentiated" << endline;
  if (func::is_a_func("\\sum")) { // added in 1.4.1 for imath
    return func("\\sum", sum(lower, higher, e)).hold();
  } else {
    return func("sum", sum(lower, higher, e)).hold();  
  }
} // sum_deriv()

static void sum_print_ltx(const ex &lower, const ex &higher, const ex &e, const print_context &c) {
  c.s << (optstack::options->get(o_eqparse).boolean ? "\\sum{" : "\\sum_{");
  optstack::options->save(o_eqalign);
  option_value o;
  o.align = none;
  optstack::options->set(o_eqalign, o);
  print_ltx(lower, c.s); // Don't print any alignment characters in the subscript of the \sum!
  optstack::options->restore(o_eqalign);
  c.s << (optstack::options->get(o_eqparse).boolean ? ";" : "}^{");
  print_ltx(higher, c.s);
  c.s << (optstack::options->get(o_eqparse).boolean ? ";" : "} ");
  if (is_a<add>(e) || is_a<mul>(e)) c.s << "\\left(";
  print_ltx(e, c.s);
  if (is_a<add>(e) || is_a<mul>(e)) c.s << "\\right)";
  if (optstack::options->get(o_eqparse).boolean) c.s << "}";
} // sum_print_ltx()

static void sum_print_imath(const ex &lower, const ex &higher, const ex &e, const imathprint &c) {
  c.s << "sum from {";
/*  optstack::options->save(o_eqalign); // Don't print any alignment characters in the subscript of the \sum!
  option_value o;
  o.align = none;
  optstack::options->set(o_eqalign, o);*/
  lower.print(c);
//  optstack::options->restore(o_eqalign);
  c.s << "} to {";
  higher.print(c);
  c.s << "} {";
  if (is_a<add>(e) || is_a<mul>(e)) c.s << "left(";
  e.print(c);
  if (is_a<add>(e) || is_a<mul>(e)) c.s << "right)";
  c.s << "}";
} // sum_print_imath()

REGISTER_FUNCTION(sum, eval_func(sum_eval).
                       print_func<print_latex>(sum_print_ltx).
                       print_func<imathprint>(sum_print_ltx). // TODO: This should be sum_print_imath
                       derivative_func(sum_deriv));

// Differentiation functions *** added in 1.0
// TODO: Perhaps this should be a new class, not a function?
//static ex partialdiff_evalf(const ex &e, const ex &s, const ex &n) throw (std::logic_error) {
//  int ndiff = (is_a<numeric>(n) ? numeric_to_int(ex_to<numeric>(n)) : -1);
//
//  MSG_INFO(2) << "Calculating derivative #" << n << " of " << e << " to " << s << endline;
//  if (ndiff < 0)
//    throw std::logic_error("Can only find positive integer derivatives");
//  if (is_a<symbol>(s))
//    return e.diff(ex_to<symbol>(s), ndiff);
//  else if (is_a<func>(s))
//    return diff_to_func(e, ex_to<func>(s), ndiff);
//  else
//    throw std::logic_error("Can only find derivatives with respect to a symbol or a pure function");
//} // partialdiff_evalf()

static ex partialdiff_eval(const ex &e, const ex &s, const ex &n) {
  if (n.is_equal(0)) return e;
  if (e.is_equal(s)) return 1;
  if (is_a<symbol>(e) || is_a<numeric>(e)) return 0;
  if (func::is_a_func("\\diff")) { // *** changed in 1.4.1 for iMath
    return func("\\diff", partialdiff(e, s, n)).hold();
  } else {
    return func("diff", partialdiff(e, s, n)).hold();
  }
} // partialdiff_eval()

static ex partialdiff_deriv(const ex &e, const ex &s, const ex &n, unsigned deriv_param) {
  msg::error(0) << "Error: \\diff cannot be differentiated" << endline;
  if (func::is_a_func("\\diff")) { // *** changed in 1.4.1 for iMath
    return func("\\diff", partialdiff(e, s, n)).hold(); // *** added hold() in 1.2
  } else {
    return func("diff", partialdiff(e, s, n)).hold();  
  }
} // partialdiff_deriv()

// *** added in 1.4.1
static void partialdiff_print_imath(const ex &e, const ex &s, const ex &n, const imathprint &c) {
  int ndiff;
  ndiff = (is_a<numeric>(n) ? ex_to<numeric>(n).to_int() : -1);

  std::string difftype = *imathprint(c).poptions->get(o_difftype).str;
  if (difftype == "dot") {
    if (!is_a<func>(e) || ex_to<func>(e).get_numargs() != 1) {
      msg::warn(0) << "Warning: Diff type 'dot' makes no sense with " << e << endline;
    } else if ((ndiff < 0) || (ndiff > 2)) {
      msg::warn(0) << "Warning: Diff type 'dot' is not implemented for diff level " << n << endline;
    } else {
      c.s << ((ndiff == 1) ? " dot " : " ddot ");
      if (!(is_a<symbol>(e) || is_a<func>(e))) c.s << " left(";
      e.print(c);
      if (!(is_a<symbol>(e) || is_a<func>(e))) c.s << " right)";
    }
  } else if (difftype == "line") {
      if (!is_a<func>(e) || ex_to<func>(e).get_numargs() != 1) {
        msg::warn(0) << "Warning: Diff type 'line' makes no sense with " << e << endline;
      } else if (ndiff < 0) {
        msg::warn(0) << "Warning: Diff type 'line' is not implemented for diff level " << n << endline;
      } else {
        if (!(is_a<symbol>(e) || is_a<func>(e))) c.s << " left(";
        if (is_a<func>(e)) {
          ex_to<func>(e).print_diff_line(ndiff, c);
        } else {
          e.print(c);
        }
        if (!(is_a<symbol>(e) || is_a<func>(e))) c.s << " right)";
        if (!is_a<func>(e)) 
          for (int i = 0; i < ndiff; i++) c.s << "'";
      }
  } else {
      // Everything else is printed in 'dfdt' style
      std::string d_sign;
      if (is_a<func>(e) && ex_to<func>(e).get_numargs() != 1)
        d_sign = "partial";
      else
        d_sign = "d"; // TODO: Using "d" skews the formula formatting in OOo

      c.s << "{{" << d_sign;
      if (ndiff != 1) c.s << "^{" << n << "}";
      if (!(is_a<symbol>(e) || is_a<func>(e))) c.s << " left("; else c.s << " ";
      e.print(c);
      if (!(is_a<symbol>(e) || is_a<func>(e))) c.s << " right)";
      c.s << "} over {" << d_sign << " ";
      s.print(c);
      if (ndiff != 1) c.s << "^{" << n << "}";
      c.s << "}}";
  }
} // partialdiff_print_imath()


static void partialdiff_print_ltx(const ex &e, const ex &s, const ex &n, const print_context &c) {
  int ndiff;
  ndiff = (is_a<numeric>(n) ? ex_to<numeric>(n).to_int() : -1);

  // *** replaced by option o_difftype in 1.4.1
  std::string difftype = *optstack::options->get(o_difftype).str;
  if (difftype == "dot") {
      if (!is_a<func>(e) || ex_to<func>(e).get_numargs() != 1) {
        msg::warn(0) << "Warning: Diff type 'dot' makes no sense with " << e << endline;
      } else if ((ndiff < 0) || (ndiff > 2)) {
        msg::warn(0) << "Warning: Diff type 'dot' is not implemented for diff level " << n << endline;
      } else {
        c.s << ((ndiff == 1) ? "\\dot " : "\\ddot ");
        if (!(is_a<symbol>(e) || is_a<func>(e))) c.s << "\\left(";
        print_ltx(e, c.s);
        if (!(is_a<symbol>(e) || is_a<func>(e))) c.s << "\\right)";
      }
  } else if (difftype == "line") {
      if (!is_a<func>(e) || ex_to<func>(e).get_numargs() != 1) {
        msg::warn(0) << "Warning: Diff type 'line' makes no sense with " << e << endline;
      } else if (ndiff < 0) {
        msg::warn(0) << "Warning: Diff type 'line' is not implemented for diff level " << n << endline;
      } else {
        if (!(is_a<symbol>(e) || is_a<func>(e))) c.s << "\\left(";
        if (is_a<func>(e)) {
          ex_to<func>(e).print_diff_line(ndiff, c);
        } else {
          print_ltx(e, c.s);
        }
        if (!(is_a<symbol>(e) || is_a<func>(e))) c.s << "\\right)";
        if (!is_a<func>(e))
          for (int i = 0; i < ndiff; i++) c.s << "'";
      }
  } else {
      // Everything else is printed in 'dfdt' style
      std::string d_sign;
      if (is_a<func>(e) && ex_to<func>(e).get_numargs() != 1)
        d_sign = "\\delta";
      else
        d_sign = "\\text{d}";

      c.s << "\\frac{" << d_sign;
      if (ndiff != 1) c.s << "^{" << n << "}";
      if (!(is_a<symbol>(e) || is_a<func>(e))) c.s << "\\left("; else c.s << "{}";
      print_ltx(e, c.s);
      if (!(is_a<symbol>(e) || is_a<func>(e))) c.s << "\\right)";
      c.s << "}{" << d_sign << "{}";
      print_ltx(s, c.s);
      if (ndiff != 1) c.s << "^{" << n << "}";
      c.s << "}";
  }
} // partialdiff_print_ltx()

// TODO: implement series expansion?

REGISTER_FUNCTION(partialdiff, eval_func(partialdiff_eval).
                         //evalf_func(partialdiff_evalf).
                         print_func<print_latex>(partialdiff_print_ltx).
                         print_func<imathprint>(partialdiff_print_ltx). // TODO: partialdiff_print_imath should work here!!
                         derivative_func(partialdiff_deriv));

// Helper function, find derivative with respect to a (pure) function
// Warning: No check is done whether s really is a function!
ex diff_to_func(const ex &e, const func &s, const unsigned n) { // *** added in 1.0
  symbol diffsym(s.get_name()); // Create a temporary symbol
  MSG_INFO(2) << "Differentiating " << e << " to function " << s << endline;
  ex result = e.subs(s == diffsym).diff(diffsym, n);
  MSG_INFO(2) << "Intermediate result " << result << endline;
  func::replace_function_by_func replace_functions; // get rid of GiNaC functions that might have been introduced *** added in 1.3.0
  return replace_functions(result.subs(diffsym == s));
} // diff_to_func()

 ex func::expand_sum::operator()(const ex &e) throw (std::logic_error) { // *** added in 1.0
  MSG_INFO(2) << "Expanding sum in " << e << endline;

  if (is_a<func>(e) && ((ex_to<func>(e).name == "\\sum") || (ex_to<func>(e).name == "sum"))) { // *** added sum in 1.4.1 for imath
    func f = ex_to<func>(e);

    // Handle children first
    exvector eseq;
    expand_sum expand_s;
    eseq.reserve(f.seq.size());
    for (exvector::const_iterator i = f.seq.begin(); i != f.seq.end(); i++)
      eseq.push_back(expand_s(*i));
    f.seq = eseq;

    MSG_INFO(0) << "Lower bound: " << f.seq[0] << endline;
    symbol var = ex_to<symbol>(ex_to<equation>(f.seq[0]).lhs());
    expression lbound = ex_to<equation>(f.seq[0]).rhs();
    expression hbound = f.seq[1];
    MSG_INFO(0) << "Summing up " << f.seq[2] << " from " << var << " = " << lbound
                       << " to " << hbound << endline;
    
    int l, h;
    if (!is_a<numeric>(lbound)) { // *** Changed in 1.4.1 because VAL(x) makes all numbers become floats
      throw std::logic_error("Lower bound of sum must be an integer");
    } else {
      l = numeric_to_int(ex_to<numeric>(lbound));
    }
    if (!is_a<numeric>(hbound)) {// *** changed in 1.4.1
      throw std::logic_error("Higher bound of sum must be an integer");
    } else {
      h = numeric_to_int(ex_to<numeric>(hbound));
    }
    
    expression result;
    while (l <= h) {
      MSG_INFO(0) << "Summing up: current value: " << result << endline;
      result = result + expression(f.seq[2].subs(var == l));
      l++;
    }
    return result;
  }

  return e.map(*this);
} // expand_sum::operator()

ex func::expand_partialdiff::operator()(const ex &e) throw (std::logic_error) { // *** added in 1.0
  MSG_INFO(2) << "Expanding \\diff in " << e << endline;

  if (is_a<func>(e) && ((ex_to<func>(e).name == "\\diff") || (ex_to<func>(e).name == "diff"))) { // *** added "diff" in 1.4.1 for iMath
    func f = ex_to<func>(e);

    // Handle children first
    exvector eseq;
    expand_partialdiff expand_diff;
    eseq.reserve(f.seq.size());
    for (exvector::const_iterator i = f.seq.begin(); i != f.seq.end(); i++)
      eseq.push_back(expand_diff(*i));
    f.seq = eseq;

    int ndiff = (is_a<numeric>(f.seq[2]) ? numeric_to_int(ex_to<numeric>(f.seq[2])) : -1);
    MSG_INFO(2) << "Calculating derivative #" << f.seq[2] << " of " << f.seq[0]
                       << " to " << f.seq[1] << endline;
    if (ndiff < 0)
      throw std::logic_error("Can only find positive integer derivatives");
    if (is_a<symbol>(f.seq[1]))
      return f.seq[0].diff(ex_to<symbol>(f.seq[1]), ndiff);
    else if (is_a<func>(f.seq[1]))
      return diff_to_func(f.seq[0], ex_to<func>(f.seq[1]), ndiff);
    else
      throw std::logic_error("Can only find derivatives with respect to a symbol or a pure function");
  }

  return e.map(*this);
} // expand_partialdiff::operator()

static ex mindex_eval(const ex &e, const ex &r, const ex &c) {
  // Immediately evaluate if e is a matrix and r and c are integers
  MSG_INFO(3) << "mindex() eval: " << e << ", " << r << ", " << c << endline;
  
  if (is_a<matrix>(e)) {
    matrix m = ex_to<matrix>(e);
    int row = 0;
    int col = 0;
    
    if (is_a<numeric>(r)) {
      numeric rnum = ex_to<numeric>(r);

      if  (rnum.info(info_flags::posint)) {
        row = rnum.to_int();
      } else { // try harder, since non-integer indixes make no sense
        row = numeric_to_int(rnum);
      }
      row--; // Adjust index to count from 0
      if (row < 0) {
        MSG_WARN(0) << "Warning: Row index out of bounds for " << m << endline;
        return func("mindex", mindex(e, r, c)).hold();
      }
    } else if (is_a<wildcard>(r)) {
      row = -1;
    } else {
      return func("mindex", mindex(e, r, c)).hold();
    }
    
    MSG_INFO(2) << "mindex() row: " << row << ", rows(): " << m.rows() << endline;
    
    if (is_a<numeric>(c)) {
      numeric cnum = ex_to<numeric>(c);
      
      if (cnum.info(info_flags::posint)) {
        col = cnum.to_int();
      } else {
        col = numeric_to_int(cnum);
      }
      col--; // Adjust index to count from 0
      if (col < 0) {
        MSG_WARN(0) << "Warning: Column index out of bounds for " << m << endline;
        func("mindex", mindex(e, r, c)).hold();
      }
    } else if (is_a<wildcard>(c)) {
      col = -1;
    } else {
      return func("mindex", mindex(e, r, c)).hold();
    }
    
    MSG_INFO(2) << "mindex() col: " << col << ", cols(): " << m.cols() << endline;
   
    if (row == -1) {
      if (col == -1) {
        return e; // m[wild, wild] returns complete matrix
      } else {
        matrix result(m.rows(), 1);
        if (col > (int)m.cols()) {
          MSG_WARN(0) << "Warning: Column index out of bounds for " << m << endline;
          return func("mindex", mindex(e, r, c)).hold();
        }
        for (unsigned i = 0; i < m.rows(); i++) result(i, 0) = m(i, col);
        return result;
      }
    } else {
      if (col == -1) {
        if (m.rows() == 1) {
          if (row > (int)m.cols()) {
            MSG_WARN(0) << "Warning: Index out of bounds for " << m << endline;
            return func("mindex", mindex(e, r, c)).hold();
          }
          return m(0, row); // Special treatment for vector: Return one element only
        } else {
          matrix result(1, m.cols());
          if (row > (int)m.rows()) {
            MSG_WARN(0) << "Warning: Row index out of bounds for " << m << endline;
            return func("mindex", mindex(e, r, c)).hold();
          }
          for (unsigned i = 0; i < m.cols(); i++) result(0, i) = m(row, i);
          return result;
        }
      } else {
        return m(row, col);
      }
    }
  } else {
    return func("mindex", mindex(e, r, c)).hold();
  }
} // mindex_eval()

static ex mindex_deriv(const ex &e, const ex &r, const ex &c, unsigned deriv_param) {
  msg::error(0) << "Error: \\mindex cannot be differentiated" << endline;
  return func("mindex", mindex(e, r, c)).hold();  
} // mindex_deriv()

static void mindex_print_ltx(const ex &e, const ex &row, const ex &col, const print_context &c) {
  std::cout << "mindex_print_ltx" << std::endl;
  print_ltx(e, c.s);
  c.s << "[";
  if (is_a<func>(row) && (ex_to<func>(row).get_name() == "wild")) {
    c.s << "*,";
  } else {
    print_ltx(row, c.s);
    c.s << ",";
  }
  if (is_a<func>(col) && (ex_to<func>(col).get_name() == "wild")) {
    // Just leave away column index so that it looks nicer
  } else {
    print_ltx(col, c.s);
  }
  c.s << "]";
} // mindex_print_ltx()

static void mindex_print_imath(const ex &e, const ex &row, const ex &col, const imathprint &c) {
  e.print(c);
  c.s << "[";
  if (is_a<wildcard>(row)) {
    c.s << "*,";
  } else {
    row.print(c);
  }
  if (is_a<wildcard>(col)) {
    // Just leave away column index so that it looks nicer
  } else {
    c.s << ",";
    col.print(c);
  }
  c.s << "]";
} // mindex_print_imath()

REGISTER_FUNCTION(mindex, eval_func(mindex_eval).
                         //evalf_func(mindex_evalf).
                         print_func<print_latex>(mindex_print_ltx).
                         print_func<imathprint>(mindex_print_ltx). // TODO: This should be mindex_print_imath, but that is not possible
                         derivative_func(mindex_deriv));

// new input/output operators
/*message &operator<<(message &ms, const func &f) {
  std::ostringstream os;
  os << f;
  ms << os.str();
  return (ms);
} // operator<<
*/

// Moved to end because partialdiff needs to have been declared

void func::init() { // *** moved here from ltxeqlex.ll and adapted in 0.5
  if (initialized) return;
#ifdef CPPU_ENV
  functions = std::map<const std::string, funcrec*>(); 
  hard_names =  std::map<const std::string, std::string>();
  hard_names_rev = std::map<const std::string, std::string>();
  hint_map = std::map<std::string, unsigned>();
  func_inv = std::map<std::string, std::string>();  
//                diff_t func::diff_type = dfdt;
#endif

  // Initialize map of functions known by Latex and/or GiNaC.
  // TODO: complete this list
  // *** removed the \\ to make the map useable by iMath, too. In 1.4.1
  MSG_INFO(2) << "Initializing function names..." << endline;

  //hard_names[""] = "abs";
  hard_names["arccos"] = "acos";
  hard_names["arcosh"] = "acosh";
  hard_names["arcsin"] = "asin";
  hard_names["arsinh"] = "asinh";
  hard_names["arctan"] = "atan";
  //hard_names["", aa, "atan2", 2);
  hard_names["artanh"] = "atanh";
  //hard_names["arg"]    = "";
  hard_names["ceil"] = "ceil"; // *** added in 1.2
  hard_names["cos"] = "cos";
  hard_names["cosh"] = "cosh";
  hard_names["diff"] = "partialdiff"; // *** added in 1.4.1 for iMath
  //hard_names[""] = "csgn";
  //hard_names["deg"] = "";
  //hard_names["det"]    = "";
  //hard_names["dim"]    = "";
  hard_names["exp"] = "exp";
  hard_names["floor"] = "floor"; // *** added in 1.2
  //hard_names["gcd"]    = "";
  //hard_names["hom"]    = "";
  //hard_names["inf"]    = "";
  //hard_names["ker"]    = "";
  //hard_names["lg"]     = "";
  //hard_names["lim"]    = "";
  //hard_names["liminf"] = "";
  //hard_names["limsup"] = "";
  hard_names["ln"] = "log";
  //hard_names["max"]    = "";
  //hard_names["min"]    = "";
  //hard_names["Pr"]     = "";
  hard_names["round"] = "round"; // *** added in 1.2
  hard_names["sin"] = "sin";
  hard_names["sinh"] = "sinh";
  hard_names["sum"] = "sum";
  //hard_names["sup"]    = "";
  hard_names["tan"] = "tan";
  hard_names["tanh"] = "tanh";
  hard_names["wild"] = "wild";
  hard_names["mindex"] = "mindex";
  // !!!!!!!!!!!! Remember to add these functions to mathconstants.tex and init.imath!!!!!!!!

  // Create reverse map
  // *** made const_iterator iterator, std::string const std::string in 1.4.1 for MSVC compile
  for (std::map<const std::string, std::string>::iterator i = hard_names.begin();
       i != hard_names.end(); i++) hard_names_rev[i->second] = i->first;

  // Initialize hint_map *** added in 0.6, 'print' added in 1.0
  hint_map["lib"] = 1;     // A library function, will not be deleted by \clearequations
  hint_map["trig"] = 2;    // For special printing of trigonometric functions
  hint_map["expand"] = 4;  // Function is to be expanded immediately (e.g. \square)
  hint_map["nobracket"] = 8;   // Function does not need brackets when printed (e.g. \sqrt{...})
  hint_map["print"] = 16;  // Function has a hard-coded printing function for latex
  hint_map["none"] = 32; // No effect, added for iMath *** added in 1.4.1
  hint_map["defdiff"] = 64; // Differentiate definition, not function itself *** added in 1.4.1

  // Initialize list of inverse functions
  // *** Removed the \\ to make the list compatible with iMath in 1.4.1
  func_inv["sin"]  = "arcsin";
  func_inv["sinh"] = "arsinh";
  func_inv["cos"]  = "arccos";
  func_inv["cosh"] = "arcosh";
  func_inv["tan"]  = "arctan";
  func_inv["tanh"] = "artanh";
  func_inv["cot"]  = "arccot";
  func_inv["ln"]   = "exp";
  for (std::map<std::string, std::string>::iterator i = func_inv.begin(); // *** made const_iterator iterator in 1.4.1
       i != func_inv.end(); i++)
    func_inv[i->second] = i->first;

  initialized = true;
} // func::init()