basic.cpp

Go to the documentation of this file.

/*
 *  GiNaC Copyright (C) 1999-2011 Johannes Gutenberg University Mainz, Germany
 *
 *  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.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */

#include "basic.h"
#include "ex.h"
#include "numeric.h"
#include "power.h"
#include "add.h"
#include "symbol.h"
#include "lst.h"
#include "ncmul.h"
#include "relational.h"
#include "operators.h"
#include "wildcard.h"
#include "archive.h"
#include "utils.h"
#include "hash_seed.h"
#include "inifcns.h"

#include <iostream>
#include <stdexcept>
#include <typeinfo>

namespace GiNaC {

GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(basic, void,
  print_func<print_context>(&basic::do_print).
  print_func<print_tree>(&basic::do_print_tree).
  print_func<print_python_repr>(&basic::do_print_python_repr))


// default constructor, destructor, copy constructor and assignment operator

// public


basic::basic(const basic & other) : flags(other.flags & ~status_flags::dynallocated), hashvalue(other.hashvalue)
{
}

const basic & basic::operator=(const basic & other)
{
    unsigned fl = other.flags & ~status_flags::dynallocated;
    if (typeid(*this) != typeid(other)) {
        // The other object is of a derived class, so clear the flags as they
        // might no longer apply (especially hash_calculated). Oh, and don't
        // copy the tinfo_key: it is already set correctly for this object.
        fl &= ~(status_flags::evaluated | status_flags::expanded | status_flags::hash_calculated);
    } else {
        // The objects are of the exact same class, so copy the hash value.
        hashvalue = other.hashvalue;
    }
    flags = fl;
    set_refcount(0);
    return *this;
}

// protected

// none (all inlined)

// other constructors

// none (all inlined)

// archiving

void basic::read_archive(const archive_node& n, lst& syms)
{ }

void basic::archive(archive_node &n) const
{
    n.add_string("class", class_name());
}

// new virtual functions which can be overridden by derived classes

// public

void basic::print(const print_context & c, unsigned level) const
{
    print_dispatch(get_class_info(), c, level);
}

void basic::print_dispatch(const registered_class_info & ri, const print_context & c, unsigned level) const
{
    // Double dispatch on object type and print_context type
    const registered_class_info * reg_info = &ri;
    const print_context_class_info * pc_info = &c.get_class_info();

next_class:
    const std::vector<print_functor> & pdt = reg_info->options.get_print_dispatch_table();

next_context:
    unsigned id = pc_info->options.get_id();
    if (id >= pdt.size() || !(pdt[id].is_valid())) {

        // Method not found, try parent print_context class
        const print_context_class_info * parent_pc_info = pc_info->get_parent();
        if (parent_pc_info) {
            pc_info = parent_pc_info;
            goto next_context;
        }

        // Method still not found, try parent class
        const registered_class_info * parent_reg_info = reg_info->get_parent();
        if (parent_reg_info) {
            reg_info = parent_reg_info;
            pc_info = &c.get_class_info();
            goto next_class;
        }

        // Method still not found. This shouldn't happen because basic (the
        // base class of the algebraic hierarchy) registers a method for
        // print_context (the base class of the print context hierarchy),
        // so if we end up here, there's something wrong with the class
        // registry.
        throw (std::runtime_error(std::string("basic::print(): method for ") + class_name() + "/" + c.class_name() + " not found"));

    } else {

        // Call method
        pdt[id](*this, c, level);
    }
}

void basic::do_print(const print_context & c, unsigned level) const
{
    c.s << "[" << class_name() << " object]";
}

void basic::do_print_tree(const print_tree & c, unsigned level) const
{
    c.s << std::string(level, ' ') << class_name() << " @" << this
        << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec;
    if (nops())
        c.s << ", nops=" << nops();
    c.s << std::endl;
    for (size_t i=0; i<nops(); ++i)
        op(i).print(c, level + c.delta_indent);
}

void basic::do_print_python_repr(const print_python_repr & c, unsigned level) const
{
    c.s << class_name() << "()";
}

void basic::dbgprint() const
{
    this->print(print_dflt(std::cerr));
    std::cerr << std::endl;
}

void basic::dbgprinttree() const
{
    this->print(print_tree(std::cerr));
}

unsigned basic::precedence() const
{
    return 70;
}

bool basic::info(unsigned inf) const
{
    // all possible properties are false for basic objects
    return false;
}

size_t basic::nops() const
{
    // iterating from 0 to nops() on atomic objects should be an empty loop,
    // and accessing their elements is a range error.  Container objects should
    // override this.
    return 0;
}

ex basic::op(size_t i) const
{
    throw(std::range_error(std::string("basic::op(): ") + class_name() + std::string(" has no operands")));
}

ex & basic::let_op(size_t i)
{
    ensure_if_modifiable();
    throw(std::range_error(std::string("basic::let_op(): ") + class_name() + std::string(" has no operands")));
}

ex basic::operator[](const ex & index) const
{
    if (is_exactly_a<numeric>(index))
        return op(static_cast<size_t>(ex_to<numeric>(index).to_int()));

    throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
}

ex basic::operator[](size_t i) const
{
    return op(i);
}

ex & basic::operator[](const ex & index)
{
    if (is_exactly_a<numeric>(index))
        return let_op(ex_to<numeric>(index).to_int());

    throw(std::invalid_argument(std::string("non-numeric indices not supported by ") + class_name()));
}

ex & basic::operator[](size_t i)
{
    return let_op(i);
}

bool basic::has(const ex & pattern, unsigned options) const
{
    exmap repl_lst;
    if (match(pattern, repl_lst))
        return true;
    for (size_t i=0; i<nops(); i++)
        if (op(i).has(pattern, options))
            return true;
    
    return false;
}

ex basic::map(map_function & f) const
{
    size_t num = nops();
    if (num == 0)
        return *this;

    basic *copy = NULL;
    for (size_t i=0; i<num; i++) {
        const ex & o = op(i);
        const ex & n = f(o);
        if (!are_ex_trivially_equal(o, n)) {
            if (copy == NULL)
                copy = duplicate();
            copy->let_op(i) = n;
        }
    }

    if (copy) {
        copy->setflag(status_flags::dynallocated);
        copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
        return *copy;
    } else
        return *this;
}

bool basic::is_polynomial(const ex & var) const
{
    return !has(var) || is_equal(ex_to<basic>(var));
}

int basic::degree(const ex & s) const
{
    return is_equal(ex_to<basic>(s)) ? 1 : 0;
}

int basic::ldegree(const ex & s) const
{
    return is_equal(ex_to<basic>(s)) ? 1 : 0;
}

ex basic::coeff(const ex & s, int n) const
{
    if (is_equal(ex_to<basic>(s)))
        return n==1 ? _ex1 : _ex0;
    else
        return n==0 ? *this : _ex0;
}

ex basic::collect(const ex & s, bool distributed) const
{
    ex x;
    if (is_a<lst>(s)) {

        // List of objects specified
        if (s.nops() == 0)
            return *this;
        if (s.nops() == 1)
            return collect(s.op(0));

        else if (distributed) {

            x = this->expand();
            if (! is_a<add>(x))
                return x; 
            const lst& l(ex_to<lst>(s));

            exmap cmap;
            cmap[_ex1] = _ex0;
            for (const_iterator xi=x.begin(); xi!=x.end(); ++xi) {
                ex key = _ex1;
                ex pre_coeff = *xi;
                for (lst::const_iterator li=l.begin(); li!=l.end(); ++li) {
                    int cexp = pre_coeff.degree(*li);
                    pre_coeff = pre_coeff.coeff(*li, cexp);
                    key *= pow(*li, cexp);
                }
                exmap::iterator ci = cmap.find(key);
                if (ci != cmap.end())
                    ci->second += pre_coeff;
                else
                    cmap.insert(exmap::value_type(key, pre_coeff));
            }

            exvector resv;
            for (exmap::const_iterator mi=cmap.begin(); mi != cmap.end(); ++mi)
                resv.push_back((mi->first)*(mi->second));
            return (new add(resv))->setflag(status_flags::dynallocated);

        } else {

            // Recursive form
            x = *this;
            size_t n = s.nops() - 1;
            while (true) {
                x = x.collect(s[n]);
                if (n == 0)
                    break;
                n--;
            }
        }

    } else {

        // Only one object specified
        for (int n=this->ldegree(s); n<=this->degree(s); ++n)
            x += this->coeff(s,n)*power(s,n);
    }
    
    // correct for lost fractional arguments and return
    return x + (*this - x).expand();
}

ex basic::eval(int level) const
{
    // There is nothing to do for basic objects:
    return hold();
}

struct evalf_map_function : public map_function {
    int level;
    evalf_map_function(int l) : level(l) {}
    ex operator()(const ex & e) { return evalf(e, level); }
};

ex basic::evalf(int level) const
{
    if (nops() == 0)
        return *this;
    else {
        if (level == 1)
            return *this;
        else if (level == -max_recursion_level)
            throw(std::runtime_error("max recursion level reached"));
        else {
            evalf_map_function map_evalf(level - 1);
            return map(map_evalf);
        }
    }
}

struct evalm_map_function : public map_function {
    ex operator()(const ex & e) { return evalm(e); }
} map_evalm;

ex basic::evalm() const
{
    if (nops() == 0)
        return *this;
    else
        return map(map_evalm);
}

struct eval_integ_map_function : public map_function {
    ex operator()(const ex & e) { return eval_integ(e); }
} map_eval_integ;

ex basic::eval_integ() const
{
    if (nops() == 0)
        return *this;
    else
        return map(map_eval_integ);
}

ex basic::eval_indexed(const basic & i) const
 // this function can't take a "const ex & i" because that would result
 // in an infinite eval() loop
{
    // There is nothing to do for basic objects
    return i.hold();
}

ex basic::add_indexed(const ex & self, const ex & other) const
{
    return self + other;
}

ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
{
    return self * other;
}

bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
{
    // Do nothing
    return false;
}

bool basic::match(const ex & pattern, exmap& repl_lst) const
{
/*
    Sweet sweet shapes, sweet sweet shapes,
    That's the key thing, right right.
    Feed feed face, feed feed shapes,
    But who is the king tonight?
    Who is the king tonight?
    Pattern is the thing, the key thing-a-ling,
    But who is the king of Pattern?
    But who is the king, the king thing-a-ling,
    Who is the king of Pattern?
    Bog is the king, the king thing-a-ling,
    Bog is the king of Pattern.
    Ba bu-bu-bu-bu bu-bu-bu-bu-bu-bu bu-bu
    Bog is the king of Pattern.
*/

    if (is_exactly_a<wildcard>(pattern)) {

        // Wildcard matches anything, but check whether we already have found
        // a match for that wildcard first (if so, the earlier match must be
        // the same expression)
        for (exmap::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
            if (it->first.is_equal(pattern))
                return is_equal(ex_to<basic>(it->second));
        }
        repl_lst[pattern] = *this;
        return true;

    } else {

        // Expression must be of the same type as the pattern
        if (typeid(*this) != typeid(ex_to<basic>(pattern)))
            return false;

        // Number of subexpressions must match
        if (nops() != pattern.nops())
            return false;

        // No subexpressions? Then just compare the objects (there can't be
        // wildcards in the pattern)
        if (nops() == 0)
            return is_equal_same_type(ex_to<basic>(pattern));

        // Check whether attributes that are not subexpressions match
        if (!match_same_type(ex_to<basic>(pattern)))
            return false;

        // Even if the expression does not match the pattern, some of
        // its subexpressions could match it. For example, x^5*y^(-1)
        // does not match the pattern $0^5, but its subexpression x^5
        // does. So, save repl_lst in order to not add bogus entries.
        exmap tmp_repl = repl_lst;
        // Otherwise the subexpressions must match one-to-one
        for (size_t i=0; i<nops(); i++)
            if (!op(i).match(pattern.op(i), tmp_repl))
                return false;

        // Looks similar enough, match found
        repl_lst = tmp_repl;
        return true;
    }
}

ex basic::subs_one_level(const exmap & m, unsigned options) const
{
    exmap::const_iterator it;

    if (options & subs_options::no_pattern) {
        ex thisex = *this;
        it = m.find(thisex);
        if (it != m.end())
            return it->second;
        return thisex;
    } else {
        for (it = m.begin(); it != m.end(); ++it) {
            exmap repl_lst;
            if (match(ex_to<basic>(it->first), repl_lst))
                return it->second.subs(repl_lst, options | subs_options::no_pattern);
            // avoid infinite recursion when re-substituting the wildcards
        }
    }

    return *this;
}

ex basic::subs(const exmap & m, unsigned options) const
{
    size_t num = nops();
    if (num) {

        // Substitute in subexpressions
        for (size_t i=0; i<num; i++) {
            const ex & orig_op = op(i);
            const ex & subsed_op = orig_op.subs(m, options);
            if (!are_ex_trivially_equal(orig_op, subsed_op)) {

                // Something changed, clone the object
                basic *copy = duplicate();
                copy->setflag(status_flags::dynallocated);
                copy->clearflag(status_flags::hash_calculated | status_flags::expanded);

                // Substitute the changed operand
                copy->let_op(i++) = subsed_op;

                // Substitute the other operands
                for (; i<num; i++)
                    copy->let_op(i) = op(i).subs(m, options);

                // Perform substitutions on the new object as a whole
                return copy->subs_one_level(m, options);
            }
        }
    }

    // Nothing changed or no subexpressions
    return subs_one_level(m, options);
}

ex basic::diff(const symbol & s, unsigned nth) const
{
    // trivial: zeroth derivative
    if (nth==0)
        return ex(*this);
    
    // evaluate unevaluated *this before differentiating
    if (!(flags & status_flags::evaluated))
        return ex(*this).diff(s, nth);
    
    ex ndiff = this->derivative(s);
    while (!ndiff.is_zero() &&    // stop differentiating zeros
           nth>1) {
        ndiff = ndiff.diff(s);
        --nth;
    }
    return ndiff;
}

exvector basic::get_free_indices() const
{
    return exvector(); // return an empty exvector
}

ex basic::conjugate() const
{
    return *this;
}

ex basic::real_part() const
{
    return real_part_function(*this).hold();
}

ex basic::imag_part() const
{
    return imag_part_function(*this).hold();
}

ex basic::eval_ncmul(const exvector & v) const
{
    return hold_ncmul(v);
}

// protected

struct derivative_map_function : public map_function {
    const symbol &s;
    derivative_map_function(const symbol &sym) : s(sym) {}
    ex operator()(const ex & e) { return diff(e, s); }
};

ex basic::derivative(const symbol & s) const
{
    if (nops() == 0)
        return _ex0;
    else {
        derivative_map_function map_derivative(s);
        return map(map_derivative);
    }
}

int basic::compare_same_type(const basic & other) const
{
    return compare_pointers(this, &other);
}

bool basic::is_equal_same_type(const basic & other) const
{
    return compare_same_type(other)==0;
}

bool basic::match_same_type(const basic & other) const
{
    // The default is to only consider subexpressions, but not any other
    // attributes
    return true;
}

unsigned basic::return_type() const
{
    return return_types::commutative;
}

return_type_t basic::return_type_tinfo() const
{
    return_type_t rt;
    rt.tinfo = &typeid(*this);
    rt.rl = 0;
    return rt;
}

unsigned basic::calchash() const
{
    unsigned v = make_hash_seed(typeid(*this));
    for (size_t i=0; i<nops(); i++) {
        v = rotate_left(v);
        v ^= this->op(i).gethash();
    }

    // store calculated hash value only if object is already evaluated
    if (flags & status_flags::evaluated) {
        setflag(status_flags::hash_calculated);
        hashvalue = v;
    }

    return v;
}

struct expand_map_function : public map_function {
    unsigned options;
    expand_map_function(unsigned o) : options(o) {}
    ex operator()(const ex & e) { return e.expand(options); }
};

ex basic::expand(unsigned options) const
{
    if (nops() == 0)
        return (options == 0) ? setflag(status_flags::expanded) : *this;
    else {
        expand_map_function map_expand(options);
        return ex_to<basic>(map(map_expand)).setflag(options == 0 ? status_flags::expanded : 0);
    }
}


// non-virtual functions in this class

// public

int basic::compare(const basic & other) const
{
#ifdef GINAC_COMPARE_STATISTICS
    compare_statistics.total_basic_compares++;
#endif
    const unsigned hash_this = gethash();
    const unsigned hash_other = other.gethash();
    if (hash_this<hash_other) return -1;
    if (hash_this>hash_other) return 1;
#ifdef GINAC_COMPARE_STATISTICS
    compare_statistics.compare_same_hashvalue++;
#endif

    const std::type_info& typeid_this = typeid(*this);
    const std::type_info& typeid_other = typeid(other);
    if (typeid_this == typeid_other) {
//      int cmpval = compare_same_type(other);
//      if (cmpval!=0) {
//          std::cout << "hash collision, same type: " 
//                    << *this << " and " << other << std::endl;
//          this->print(print_tree(std::cout));
//          std::cout << " and ";
//          other.print(print_tree(std::cout));
//          std::cout << std::endl;
//      }
//      return cmpval;
#ifdef GINAC_COMPARE_STATISTICS
        compare_statistics.compare_same_type++;
#endif
        return compare_same_type(other);
    } else {
//      std::cout << "hash collision, different types: " 
//                << *this << " and " << other << std::endl;
//      this->print(print_tree(std::cout));
//      std::cout << " and ";
//      other.print(print_tree(std::cout));
//      std::cout << std::endl;
        return (typeid_this.before(typeid_other) ? -1 : 1);
    }
}

bool basic::is_equal(const basic & other) const
{
#ifdef GINAC_COMPARE_STATISTICS
    compare_statistics.total_basic_is_equals++;
#endif
    if (this->gethash()!=other.gethash())
        return false;
#ifdef GINAC_COMPARE_STATISTICS
    compare_statistics.is_equal_same_hashvalue++;
#endif
    if (typeid(*this) != typeid(other))
        return false;
    
#ifdef GINAC_COMPARE_STATISTICS
    compare_statistics.is_equal_same_type++;
#endif
    return is_equal_same_type(other);
}

// protected

const basic & basic::hold() const
{
    return setflag(status_flags::evaluated);
}

void basic::ensure_if_modifiable() const
{
    if (get_refcount() > 1)
        throw(std::runtime_error("cannot modify multiply referenced object"));
    clearflag(status_flags::hash_calculated | status_flags::evaluated);
}

// global variables

int max_recursion_level = 1024;


#ifdef GINAC_COMPARE_STATISTICS
compare_statistics_t::~compare_statistics_t()
{
    std::clog << "ex::compare() called " << total_compares << " times" << std::endl;
    std::clog << "nontrivial compares: " << nontrivial_compares << " times" << std::endl;
    std::clog << "basic::compare() called " << total_basic_compares << " times" << std::endl;
    std::clog << "same hashvalue in compare(): " << compare_same_hashvalue << " times" << std::endl;
    std::clog << "compare_same_type() called " << compare_same_type << " times" << std::endl;
    std::clog << std::endl;
    std::clog << "ex::is_equal() called " << total_is_equals << " times" << std::endl;
    std::clog << "nontrivial is_equals: " << nontrivial_is_equals << " times" << std::endl;
    std::clog << "basic::is_equal() called " << total_basic_is_equals << " times" << std::endl;
    std::clog << "same hashvalue in is_equal(): " << is_equal_same_hashvalue << " times" << std::endl;
    std::clog << "is_equal_same_type() called " << is_equal_same_type << " times" << std::endl;
    std::clog << std::endl;
    std::clog << "basic::gethash() called " << total_gethash << " times" << std::endl;
    std::clog << "used cached hashvalue " << gethash_cached << " times" << std::endl;
}

compare_statistics_t compare_statistics;
#endif

} // namespace GiNaC