Latex Equation Compiler Code

/***************************************************************************
                          utils.cpp  -  Utility functions for EQC
                             -------------------
    begin                : Mon Oct 22 2001
    copyright            : (C) 2001 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 <cmath> // *** added in 1.0
#include <algorithm>
#include <stdexcept> // *** changed to stdexcept in 0.7
#include "../config/config.h" // *** added in 1.3.1
#ifdef HAVE_SSTREAM // *** added in 1.3.1
  #include <sstream>
#else
  #include <strstream>
#endif
#include "utils.h"
#include "operands.h"
#include "msgdriver.h" // *** added in 1.3.1
#include "unit.h" // *** added in 1.0

void set_print(std::set<const char*, ltstr> &s, std::ostream &os) {
  copy(s.begin(), s.end(), std::ostream_iterator<const char*>(os, " "));
}

/*message &operator<<(message &ms, std::set<const char*, ltstr> &s) {
  std::ostringstream os;
  copy(s.begin(), s.end(), std::ostream_iterator<const char*>(os, " "));
  ms << os.str();
  return ms;
} // message::operator<<()
*/
bool set_has(std::set<const char*, ltstr> &s, const std::string &element) {
  const char *temp[1] = {element.c_str()};
  std::set<const char*, ltstr> e(temp, temp+1);
  std::set<const char*, ltstr> result;

  set_intersection(s.begin(), s.end(), e.begin(), e.end(), inserter(result, result.begin()), ltstr());
  return !result.empty();
}

std::set<const char *, ltstr> set_add(const std::set<const char*, ltstr> &s, const std::string &element) {
  // element has to be copied in this complicated way because otherwise it will disappear later on - why?
  std::string *elem = new std::string (element.c_str());
  const char *temp[1] = {elem->c_str()};
  std::set<const char*, ltstr> e(temp, temp+1);
  std::set<const char*, ltstr> result;
  set_union(s.begin(), s.end(), e.begin(), e.end(), inserter(result, result.begin()), ltstr());
  return result;
}

bool is_negex(const ex &e) {
  // *** changed in 1.0 to use info_flags::negative
  //0) << "Checking for less than zero: " << e << endline;
  return ((is_a<numeric>(e) && (e.info(info_flags::negative))) ||
          (is_a<expairseq>(e) && (e.op(e.nops()).info(info_flags::negative))));
} // is_negex()

bool is_negpower(const ex &p) {
  if (!is_a<power>(p)) return false;
  return (is_negex(ex_to<power>(p).op(1)));
} // is_negpower()

bool is_equal_one(const numeric &n) { // *** added in 1.0
  //0) << "Testing " << (n-1).to_double() << " for equality to zero" << endline;
  return (abs(n - 1).to_double() < std::pow(double(10), double(-Digits)));
} // is_equal_one()

bool is_equal_zero(const numeric &n) { // *** added in 1.0
  return (abs(n).to_double() < std::pow(double(10), double(-Digits)));
} // is_equal_zero()

// Test whether the argument is a quantity (i.e., a numeric multiplied
// with any number of units or powers of units)
const bool is_quantity (const expression &quantity) {
  // *** changed ex to expression in 0.4, rewrite in 1.0 for efficiency reasons
  MSG_INFO(2) << "Testing expression " << quantity << endline;

  if (is_a<mul>(quantity)) {
    for (const_iterator i = quantity.begin(); i != quantity.end(); i++) {
      if (is_a<power>(*i)) {
        if (!is_a<Unit>(get_basis(ex_to<power>(*i)))) return false;
        if (!get_exp(ex_to<power>(*i)).info(info_flags::real)) return false;
      } else if (!(is_a<numeric>(*i) || is_a<Unit>(*i))) {
        return false;
      }
    }
    return true;
  } else if (is_a<matrix>(quantity)) {
    for (unsigned i = 0; i < ex_to<matrix>(quantity).rows(); i++)
      for (unsigned j = 0; j < ex_to<matrix>(quantity).cols(); j++)
        if (!is_quantity(ex_to<matrix>(quantity)(i,j))) return false;
    return true;
  }

  return(is_a<numeric>(quantity) || is_a<Unit>(quantity));
} // is_quantity()

lst make_lst_from_matrix(const expression &e) { // *** added in 1.3.0
  lst result;
  if (is_a<matrix>(e)) {
    matrix m = ex_to<matrix>(e);
    if (m.rows() > 1) throw std::invalid_argument("Matrix may have only one row");
    
    for (unsigned i = 0; i < m.cols(); i++)
      result.append(m(0,i));
    MSG_INFO(1) << "Create lst from matrix " << m << ", result: " << result << endline;
  } else if (is_a<lst>(e)) {
    result = ex_to<lst>(e);
 } else {
   result.append(e);
 }
 return result;
}

exvector make_exvector(const expression &e) { // *** added in 0.6
  exvector result;
  if (is_a<lst>(e)) {
    lst l = ex_to<lst>(e);
    for (lst::const_iterator i = l.begin(); i != l.end(); i++) {
      // *** changed to const_iterator in 0.7
      result.push_back(*i);
    }
  } else {
    result.push_back(e);
  }
  return result;
} // make_exvector()

matrix make_matrix(const expression &e) { // *** added in 1.3.0
  if (is_a<matrix>(e)) {
    return ex_to<matrix>(e);
  } else if (is_a<lst>(e)) {
    return matrix (1, ex_to<lst>(e).nops(), ex_to<lst>(e));
  } else {
    exvector elements;
    elements.push_back(e);
    return matrix (1, elements.size(), elements);
  }
} // make_matrix() 
  
// Make an equation label printable in Latex without needing
// a verbatim environment (which the align environment will not allow in a tag!)
const std::string ltx_label(const std::string &l) { // *** added in 0.8
  std::string result = l;
  unsigned pos = 0;

  while ((pos = result.find('_', pos)) < result.size()) {
    result.insert(pos, 1, '\\');
    pos = pos + 2;
  }
  pos = 0;
  while ((pos = result.find('^', pos)) < result.size()) {
    result.insert(pos, 1, '\\');
    pos = pos + 2;
  }

  return result;
} // ltx_label()

const unsigned numeric_to_uint(const numeric &n) { // *** added in 1.0
  return cln::cl_I_to_uint(cln::floor1(cln::cl_F(n.to_double())));
}

const int numeric_to_int(const numeric &n) { // *** added in 1.0
  return cln::cl_I_to_int(cln::floor1(cln::cl_F(n.to_double())));
}

// Map functions ================
// *** added in 1.0

// Helper function for expand_real_powers
// Reduce double powers in units. This can be done safely since a unit is
// always positive (otherwise sqrt(x^2) could be +x or -x)
const power check_double_power(const ex &e, const ex &exp) {
  if (is_a<power>(e) && is_a<Unit>(get_basis(ex_to<power>(e))))
    return power(get_basis(ex_to<power>(e)), get_exp(ex_to<power>(e)) * exp);
  else
    return power(e, exp);
} // check_double_power()

ex expand_real_powers::operator()(const ex &e) { // *** added in 1.0
  MSG_INFO(3) << "Special expanding powers in " << e << endline;

  if (is_a<power>(e)) {
    expand_real_powers expand_power;
    ex the_basis = expand_power(get_basis(ex_to<power>(e)));
    ex the_exp   = expand_power(get_exp(ex_to<power>(e)));
    //0) << "Basis: " << the_basis << ", exponent " << the_exp << endline;

    if (is_a<mul>(the_basis) && !the_exp.info(info_flags::integer)) {
      expression result = numeric(1);
      for (const_iterator i = the_basis.begin(); i != the_basis.end(); i++) {
        result = result * check_double_power(*i, the_exp);
      }
      return result;
    } else if (is_a<power>(the_basis)) {
      return check_double_power(the_basis, the_exp);
    }
  }

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

ex reduce_double_powers::operator()(const ex &e) {
  MSG_INFO(2) << "Reducing double powers in " << e << endline;

  if (is_a<power>(e)) {
    reduce_double_powers reduce_powers;
    ex the_basis = reduce_powers(get_basis(ex_to<power>(e)));
    ex the_exp   = reduce_powers(get_exp(ex_to<power>(e)));
    if (is_a<power>(the_basis)) {
      power ebasis = ex_to<power>(the_basis);
      return expression(power(get_basis(ebasis), get_exp(ebasis) * the_exp));
    } else {
      return power(the_basis, the_exp);
    }
  }
  return e.map(*this);
} // reduce_double_powers::operator()

ex gather_sqrt::operator()(const ex &e) { // *** added in 1.2
	MSG_INFO(0) << "Gathering square roots in " << e << endline;
		
	if (is_a<mul>(e)) {
		ex sqrts = numeric(1);
		ex non_sqrt = numeric(1);
		gather_sqrt gather_sqrts;
			
		for (const_iterator i = e.begin(); i != e.end(); i++) {
			if ((is_a<power>(*i)) && (get_exp(ex_to<power>(*i)) == numeric(1, 2))) {
				ex the_basis = gather_sqrts(get_basis(ex_to<power>(*i)));
				sqrts *= the_basis;
			} else {
				non_sqrt *= gather_sqrts(*i);
			}
		}
			
		return sqrt(sqrts) * non_sqrt;
	} 
			
	return e.map(*this);
} // gather_sqrt::operator()

// Stream output operators ==============

std::ostream &operator<<(std::ostream &os, const exvector &v) {
  // *** moved here from func.cpp in 0.7
  // *** added check for empty() in 0.8
  if (v.empty()) {
    os << "[empty]";
  } else {
    if (v.size() > 1)
      copy(v.begin(), v.end() - 1, std::ostream_iterator<const ex>(os, "; "));
    os << *(v.end() - 1);
  }
  return (os);
} // operator<<

/*message &operator<<(message &ms, const GiNaC::basic &o) {
  std::ostringstream os;
  os << o; 
  ms << os.str();
  return ms; 
}

message &operator<<(message &ms, const GiNaC::lst &l) { 
  std::ostringstream os;
  os << l; 
  ms << os.str();
  return ms; 
}

message &operator<<(message &ms, const GiNaC::ex &e) { 
  std::ostringstream os;
  os << e; 
  ms << os.str();
  return ms;
}

message &operator<<(message &ms, const cln::cl_R &x) { 
  std::ostringstream os;
  os << x; 
  ms << os.str();
  return ms; 
}

message &operator<<(message &ms, const GiNaC::exvector &v) {
  // *** added in 0.7, added check for empty(), changed to use copy in 0.8
  if (v.empty()) {
    ms << "[empty]";
  } else {
    std::ostringstream os;
    if (v.size() > 1)
      copy(v.begin(), v.end() - 1, std::ostream_iterator<const ex>(os, "; "));
    os << *(v.end() - 1);
    ms << os.str();
  }
  return ms;
} // operator<<

message &operator<<(message &ms, const GiNaC::exmap &m) { // *** added in 0.8
  if (m.empty()) {
    ms << "[empty]";
  } else if (m.size() > 1) {
    for (exmap::const_iterator i = m.begin(); i != --m.end(); i++)
      ms << i->first << " = " << i->second << "; ";
    ms  << (--m.end())->first << " = " << (--m.end())->second;
  } else {
    ms << m.begin()->first << " = " << m.begin()->second;
  }
  return ms;
} // operator<<

message &operator<<(message &ms, const std::vector<expression> &v) { // *** added in 1.0
  if (v.empty()) {
    ms << "[empty]";
  } else if (v.size() > 1) {
    for (std::vector<expression>::const_iterator i = v.begin(); i != --v.end(); i++)
      ms << *i << "; ";
    ms  << *(--v.end());
  } else {
    ms << *(v.begin());
  }
  return ms;
} // operator<<
*/

// Arithmetic with lists *** added in 1.3.0
/*matrix operator*(const lst &l, const ex &e) {
  matrix op1(1, l.nops(), l);
  if (is_a<numeric>(e)) {
    return op1.mul(ex_to<numeric>(e));
  } else {
    return op1.mul(make_matrix(e));
  }
} // operator*(const lst&, const ex&)

ex operator*(const ex &e, const lst &l) {
  matrix op2(1, l.nops(), l);
  return (e * op2);
} // operator*(const ex&, const lst&)
*/

/* This function can be used to calculated the coefficients of 
	(x - x0) (x - x1) .. (x - xn)
	which again is required for parabolic interpolation with Lagrange's formula
	(see Bronstein p. 517)
	The coefficients are (example for order 4)
	c0 = x0 x1 x2 x3
	c1 = x0 x1 x2 + x0 x1 x3 + x1 x2 x3
	c2 = x0 x1 + x0 x2 + x0 x3 + x1 x2 + x1 x3 + x2 x3
	c3 = x0 + x1 + x2 + x3
	Since not all coefficients are always required, the number of coefficients that should be 
	calculated is passed as a parameter.
	
	The number of these so-called combinations is
	(n over m) = n! / (m! (n -m)!
	where n is the order and m the number of elements to be combined (Bronstein 139).
	
	The combinations are coded as a binary, for example, 1101 means x0 * x1 * x3. They
	follow the following algorithm (example for order 5)
	11111 -> combination of all four variables
	11110, 11101, 11011, 10111, 01111 -> (5 over 4) = 5 combinations of four variables
		created by cycling a zero through the previous number
	11100, 11010, 10110, 01110; 11001, 10101, 01101; 10011, 01011; 00111 -> (5 over 3) = 10 permutations of three variables
		created by cycling a zero through the previous numbers, beginning at position 4 with the first, 3 with the second...
	10000, 01000; 00100; 00010; 00001 -> permutations of one variable
	11000, 10100, 10010, 10001; 01100, 01010, 01001; 00110, 00101; 00011 -> (5 over 2) = 10 permutations of two variables
		created by cycling a one through the previous numbers
*/
// const std::vector< std::vector<long long unsigned> > permutations(const unsigned order. const unsigned coefficients)
// 	throw(std::invalid_argument) {
// 	if (order > sizeof(long long unsigned))
// 		throw std::invalid_argument("Maximum order exceeded");
// 	
// 	// The initial permutation
// 	std::vector< std::vector<long long unsigned> > result;
// 	result[0][0] = -1; // TODO: this sets all bits to 1 ??
// 	
// 	for (unsigned o = 1; o < order; o++) {
// 		// work on the permutations of the previous order
// 		for (std::vector<long long unsigned>::const_iterator p = result[o].begin(); p != result[o].end(); p++) {
// 			
// 
// }