basic.cpp

Go to the documentation of this file.

 
/*
 *  GiNaC Copyright (C) 1999-2025 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)
{
}
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(basic, void, {…}
 
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;
}
const basic & basic::operator=(const basic & other) {…}
 
// protected
 
// none (all inlined)
 
// other constructors
 
// none (all inlined)
 
// archiving
 
void basic::read_archive(const archive_node& n, lst& syms)
{ }
void basic::read_archive(const archive_node& n, lst& syms) {…}
 
void basic::archive(archive_node &n) const
{
    n.add_string("class", class_name());
}
void basic::archive(archive_node &n) const {…}
 
// 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(const print_context & c, unsigned level) const {…}
 
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::print_dispatch(const registered_class_info & ri, const print_context & c, unsigned level) const {…}
 
void basic::do_print(const print_context & c, unsigned level) const
{
    c.s << "[" << class_name() << " object]";
}
void basic::do_print(const print_context & c, unsigned level) const {…}
 
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_tree(const print_tree & c, unsigned level) const {…}
 
void basic::do_print_python_repr(const print_python_repr & c, unsigned level) const
{
    c.s << class_name() << "()";
}
void basic::do_print_python_repr(const print_python_repr & c, unsigned level) const {…}
 
void basic::dbgprint() const
{
    this->print(print_dflt(std::cerr));
    std::cerr << std::endl;
}
void basic::dbgprint() const {…}
 
void basic::dbgprinttree() const
{
    this->print(print_tree(std::cerr));
}
void basic::dbgprinttree() const {…}
 
unsigned basic::precedence() const
{
    return 70;
}
unsigned basic::precedence() const {…}
 
bool basic::info(unsigned inf) const
{
    // all possible properties are false for basic objects
    return false;
}
bool basic::info(unsigned inf) const {…}
 
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;
}
size_t basic::nops() const {…}
 
ex basic::op(size_t i) const
{
    throw(std::range_error(std::string("basic::op(): ") + class_name() + std::string(" has no operands")));
}
ex basic::op(size_t i) const {…}
 
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::let_op(size_t i) {…}
 
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[](const ex & index) const {…}
 
ex basic::operator[](size_t i) const
{
    return op(i);
}
ex basic::operator[](size_t i) const {…}
 
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[](const ex & index) {…}
 
ex & basic::operator[](size_t i)
{
    return let_op(i);
}
ex & basic::operator[](size_t 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;
}
bool basic::has(const ex & pattern, unsigned options) const {…}
 
ex basic::map(map_function & f) const
{
    size_t num = nops();
    if (num == 0)
        return *this;
 
    basic *copy = nullptr;
    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 == nullptr)
                copy = duplicate();
            copy->let_op(i) = n;
        }
    }
 
    if (copy) {
        copy->clearflag(status_flags::hash_calculated | status_flags::expanded);
        return *copy;
    } else
        return *this;
}
ex basic::map(map_function & f) const {…}
 
bool basic::is_polynomial(const ex & var) const
{
    return !has(var) || is_equal(ex_to<basic>(var));
}
bool basic::is_polynomial(const ex & var) const {…}
 
int basic::degree(const ex & s) const
{
    return is_equal(ex_to<basic>(s)) ? 1 : 0;
}
int basic::degree(const ex & s) const {…}
 
int basic::ldegree(const ex & s) const
{
    return is_equal(ex_to<basic>(s)) ? 1 : 0;
}
int basic::ldegree(const ex & s) const {…}
 
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::coeff(const ex & s, int n) const {…}
 
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 auto & xi : x) {
                ex key = _ex1;
                ex pre_coeff = xi;
                for (auto & li : l) {
                    int cexp = pre_coeff.degree(li);
                    pre_coeff = pre_coeff.coeff(li, cexp);
                    key *= pow(li, cexp);
                }
                auto ci = cmap.find(key);
                if (ci != cmap.end())
                    ci->second += pre_coeff;
                else
                    cmap.insert(exmap::value_type(key, pre_coeff));
            }
 
            exvector resv;
            for (auto & mi : cmap)
                resv.push_back((mi.first)*(mi.second));
            return dynallocate<add>(resv);
 
        } 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::collect(const ex & s, bool distributed) const {…}
 
ex basic::eval() const
{
    // There is nothing to do for basic objects:
    return hold();
}
ex basic::eval() const {…}
 
struct evalf_map_function : public map_function {
    ex operator()(const ex & e) override { return evalf(e); }
};
struct evalf_map_function : public map_function {…};
 
ex basic::evalf() const
{
    if (nops() == 0)
        return *this;
    else {
        evalf_map_function map_evalf;
        return map(map_evalf);
    }
}
ex basic::evalf() const {…}
 
struct evalm_map_function : public map_function {
    ex operator()(const ex & e) override { return evalm(e); }
} map_evalm;
struct evalm_map_function : public map_function {…};
 
ex basic::evalm() const
{
    if (nops() == 0)
        return *this;
    else
        return map(map_evalm);
}
ex basic::evalm() const {…}
 
struct eval_integ_map_function : public map_function {
    ex operator()(const ex & e) override { return eval_integ(e); }
} map_eval_integ;
struct eval_integ_map_function : public map_function {…};
 
ex basic::eval_integ() const
{
    if (nops() == 0)
        return *this;
    else
        return map(map_eval_integ);
}
ex basic::eval_integ() const {…}
 
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::eval_indexed(const basic & i) const {…}
 
ex basic::add_indexed(const ex & self, const ex & other) const
{
    return self + other;
}
ex basic::add_indexed(const ex & self, const ex & other) const {…}
 
ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const
{
    return self * other;
}
ex basic::scalar_mul_indexed(const ex & self, const numeric & other) const {…}
 
bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
{
    // Do nothing
    return false;
}
bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const {…}
 
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 (auto & it : repl_lst) {
            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;
    }
}
bool basic::match(const ex & pattern, exmap& repl_lst) const {…}
 
ex basic::subs_one_level(const exmap & m, unsigned options) const
{
    if (options & subs_options::no_pattern) {
        ex thisex = *this;  // NB: *this may be deleted here.
        auto it = m.find(thisex);
        if (it != m.end())
            return it->second;
        return thisex;
    } else {
        for (auto & it : m) {
            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_one_level(const exmap & m, unsigned options) const {…}
 
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->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::subs(const exmap & m, unsigned options) const {…}
 
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;
}
ex basic::diff(const symbol & s, unsigned nth) const {…}
 
exvector basic::get_free_indices() const
{
    return exvector(); // return an empty exvector
}
exvector basic::get_free_indices() const {…}
 
ex basic::conjugate() const
{
    return *this;
}
ex basic::conjugate() const {…}
 
ex basic::real_part() const
{
    return real_part_function(*this).hold();
}
ex basic::real_part() const {…}
 
ex basic::imag_part() const
{
    return imag_part_function(*this).hold();
}
ex basic::imag_part() const {…}
 
ex basic::eval_ncmul(const exvector & v) const
{
    return hold_ncmul(v);
}
ex basic::eval_ncmul(const exvector & v) const {…}
 
// protected
 
struct derivative_map_function : public map_function {
    const symbol &s;
    derivative_map_function(const symbol &sym) : s(sym) {}
    ex operator()(const ex & e) override { return diff(e, s); }
};
struct derivative_map_function : public map_function {…};
 
ex basic::derivative(const symbol & s) const
{
    if (nops() == 0)
        return _ex0;
    else {
        derivative_map_function map_derivative(s);
        return map(map_derivative);
    }
}
ex basic::derivative(const symbol & s) const {…}
 
int basic::compare_same_type(const basic & other) const
{
    return compare_pointers(this, &other);
}
int basic::compare_same_type(const basic & other) const {…}
 
bool basic::is_equal_same_type(const basic & other) const
{
    return compare_same_type(other)==0;
}
bool basic::is_equal_same_type(const basic & other) const {…}
 
bool basic::match_same_type(const basic & other) const
{
    // The default is to only consider subexpressions, but not any other
    // attributes
    return true;
}
bool basic::match_same_type(const basic & other) const {…}
 
unsigned basic::return_type() const
{
    return return_types::commutative;
}
unsigned basic::return_type() const {…}
 
return_type_t basic::return_type_tinfo() const
{
    return_type_t rt;
    rt.tinfo = &typeid(*this);
    rt.rl = 0;
    return rt;
}
return_type_t basic::return_type_tinfo() const {…}
 
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;
}
unsigned basic::calchash() const {…}
 
struct expand_map_function : public map_function {
    unsigned options;
    expand_map_function(unsigned o) : options(o) {}
    ex operator()(const ex & e) override { return e.expand(options); }
};
struct expand_map_function : public map_function {…};
 
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);
    }
}
ex basic::expand(unsigned options) const {…}
 
 
// 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);
    }
}
int basic::compare(const basic & other) const {…}
 
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);
}
bool basic::is_equal(const basic & other) const {…}
 
// protected
 
const basic & basic::hold() const
{
    return setflag(status_flags::evaluated);
}
const basic & basic::hold() const {…}
 
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);
}
void basic::ensure_if_modifiable() const {…}
 
// global variables
 
#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