/Users/mourrain/Devel/mmx/basix/src/syntactic.cpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : syntactic.cpp
* DESCRIPTION: Syntactic expressions with basic simplification
* COPYRIGHT  : (C) 2008  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#include <basix/string.hpp>
#include <basix/identifiers.hpp>
#include <basix/literal.hpp>
#include <basix/compound.hpp>
#include <basix/vector_sort.hpp>
namespace mmx {

/******************************************************************************
* Basic syntactic operations
******************************************************************************/

syntactic::syntactic (const int& i): rep (as_string (i)) {}
bool is_atom (const syntactic& g) { return is<literal> (*g); }
string as_string (const syntactic& g) { return literal_to_string (*g); }

syntactic
apply (const vector<syntactic>& v) {
  VERIFY (N(v) != 0, "non-empty vector expected");
  if (is<literal> (*v[0])) {
    string s= as_string (as<literal> (*v[0]));
    if (N(s) == 1) {
      if (s[0] == '-' && N(v) == 2) return -v[1];
      if (s[0] == '+' && N(v) == 3) return v[1] + v[2];
      if (s[0] == '-' && N(v) == 3) return v[1] - v[2];
      if (s[0] == '*' && N(v) == 3) return v[1] * v[2];
      if (s[0] == '/' && N(v) == 3) return v[1] / v[2];
      if (s[0] == '^' && N(v) == 3) return pow (v[1], v[2]);
    }
  }
  return as_syntactic (vector_to_compound (*((vector<generic>*) (void*) &v)));
}

syntactic apply (const syntactic& f) {
  return apply (vec (f)); }
syntactic apply (const syntactic& f, const syntactic& a1) {
  return apply (vec (f, a1)); }
syntactic apply (const syntactic& f, const syntactic& a1,
                 const syntactic& a2) {
  return apply (vec (f, a1, a2)); }
syntactic apply (const syntactic& f, const syntactic& a1,
                 const syntactic& a2, const syntactic& a3) {
  return apply (vec (f, a1, a2, a3)); }
syntactic apply (const syntactic& f,
                 const syntactic& a1, const syntactic& a2,
                 const syntactic& a3, const syntactic& a4) {
  return apply (vec (f, a1, a2, a3, a4)); }
syntactic apply (const syntactic& f,
                 const syntactic& a1, const syntactic& a2,
                 const syntactic& a3, const syntactic& a4,
                 const syntactic& a5) {
  return apply (vec (f, a1, a2, a3, a4, a5)); }
syntactic apply (const syntactic& f,
                 const syntactic& a1, const syntactic& a2,
                 const syntactic& a3, const syntactic& a4,
                 const syntactic& a5, const syntactic& a6) {
  return apply (vec (f, a1, a2, a3, a4, a5, a6)); }
syntactic apply (const syntactic& f, const vector<syntactic>& a) {
  return apply (cons (f, a)); }

bool is_func (const syntactic& g, const char* f) {
  return is_func (*g, f); }
bool is_func (const syntactic& g, const char* f, nat n) {
  return is_func (*g, f, n); }
bool is_func (const syntactic& g, const syntactic& f) {
  return is_func (*g, *f); }
bool is_func (const syntactic& g, const syntactic& f, nat n) {
  return is_func (*g, *f, n); }
bool is_func (const syntactic& g, const generic& f) {
  return is_func (*g, f); }
bool is_func (const syntactic& g, const generic& f, nat n) {
  return is_func (*g, f, n); }

vector<syntactic>
components (const syntactic& g) {
  vector<syntactic> v= fill<syntactic> (N(g));
  for (nat i=0; i<N(g); i++) v[i]= g[i];
  return v;
}

vector<syntactic>
arguments (const syntactic& g) {
  return cdr (components (g));
}

/******************************************************************************
* Constants
******************************************************************************/

void set_smallest (syntactic& x) { x= as_syntactic (GEN_SMALLEST); }
void set_largest (syntactic& x) { x= as_syntactic (GEN_LARGEST); }
void set_accuracy (syntactic& x) { x= as_syntactic (GEN_ACCURACY); }
void set_pi (syntactic& x) { x= as_syntactic (GEN_PI); }
void set_log2 (syntactic& x) { x= log (as_syntactic (2)); }
void set_euler (syntactic& x) { x= as_syntactic (GEN_EULER); }
void set_catalan (syntactic& x) { x= as_syntactic (GEN_CATALAN); }
void set_imaginary (syntactic& x) { x= as_syntactic (GEN_I); }
void set_nan (syntactic& x) { x= as_syntactic (GEN_NAN); }
void set_fuzz (syntactic& x) { x= as_syntactic (GEN_FUZZ); }
void set_infinity (syntactic& x) { x= as_syntactic (GEN_INFINITY); }
void set_maximal (syntactic& x) { x= as_syntactic (GEN_INFINITY); }
void set_minimal (syntactic& x) { x= -as_syntactic (GEN_INFINITY); }

/******************************************************************************
* Additive arithmetic with simplifications
******************************************************************************/

static void
signed_decompose (const syntactic& g, syntactic& abs_g, int& sgn_g) {
  if (is_atom (g)) {
    const string s= as_string (g);
    if (N(s) > 1 && s[0] == '-' &&
        (s[1] == '.' || (s[1] >= '0' && s[1] <= '9'))) {
      abs_g= syntactic (s (1, N(s)));
      sgn_g= -1;
    }
    else {
      abs_g= g;
      sgn_g= 1;
    }
  }
  else if (is_func (g, GEN_MINUS, 1)) {
    signed_decompose (g[1], abs_g, sgn_g);
    sgn_g= -sgn_g;
  }
  else if (is_func (g, GEN_TIMES, 2) || is_func (g, GEN_OVER, 2)) {
    syntactic a1, a2;
    int     s1, s2;
    signed_decompose (g[1], a1, s1);
    signed_decompose (g[2], a2, s2);
    if (is_func (g, GEN_TIMES, 2)) abs_g= a1 * a2;
    else abs_g= a1 / a2;
    sgn_g= s1 * s2;
  }
  else {
    abs_g= g;
    sgn_g= 1;
  }
}

syntactic
operator - (const syntactic& g) {
  if (exact_eq (g, 0)) return g;
  if (is_func (g, GEN_MINUS, 1))
    return g[1];
  if (is_func (g, GEN_PLUS, 2))
    return (-g[1]) - g[2];
  if (is_func (g, GEN_MINUS, 2))
    return (-g[1]) + g[2];
  syntactic abs_g; int sgn_g;
  signed_decompose (g, abs_g, sgn_g);
  if (sgn_g > 0) return syn (GEN_MINUS, abs_g);
  else return abs_g;
}

syntactic
operator + (const syntactic& g1, const syntactic& g2) {
  if (exact_eq (g1, 0)) return g2;
  if (exact_eq (g2, 0)) return g1;
  if (is_func (g2, GEN_MINUS, 1))
    return g1 - g2[1];
  if (is_func (g2, GEN_PLUS, 2))
    return (g1 + g2[1]) + g2[2];
  if (is_func (g2, GEN_MINUS, 2))
    return (g1 + g2[1]) - g2[2];
  syntactic abs_g2; int sgn_g2;
  signed_decompose (g2, abs_g2, sgn_g2);
  if (sgn_g2 < 0) return syn (GEN_MINUS, g1, abs_g2);
  else return syn (GEN_PLUS, g1, abs_g2);
}

syntactic
operator - (const syntactic& g1, const syntactic& g2) {
  if (exact_eq (g2, 0)) return g1;
  if (is_func (g2, GEN_MINUS, 1))
    return g1 + g2[1];
  if (is_func (g2, GEN_PLUS, 2))
    return (g1 - g2[1]) - g2[2];
  if (is_func (g2, GEN_MINUS, 2))
    return (g1 - g2[1]) + g2[2];
  syntactic abs_g2; int sgn_g2;
  signed_decompose (g2, abs_g2, sgn_g2);
  if (exact_eq (g1, 0)) {
    if (sgn_g2 > 0) return syn (GEN_MINUS, abs_g2);
    else return abs_g2;    
  }
  else {
    if (sgn_g2 > 0) return syn (GEN_MINUS, g1, abs_g2);
    else return syn (GEN_PLUS, g1, abs_g2);
  }
}

static syntactic
migrate_sub (const syntactic& g, const syntactic& h) {
  if (is_func (g, GEN_MINUS, 2))
    return syn (GEN_MINUS, migrate_sub (g[1], h), g[2]);
  syntactic abs_g; int sgn_g;
  signed_decompose (g, abs_g, sgn_g);
  return syn (GEN_MINUS, h, abs_g);
}

static syntactic
migrate_sub (const syntactic& g, nat found) {
  if (found == 0) return migrate_sub (g[1], g[2]);
  return syn (g[0], migrate_sub (g[1], found-1), g[2]);
}

syntactic
migrate_negate (const syntactic& g, nat maximal) {
  int i= 0, found= -1;
  syntactic h= g;
  while (is_func (h, GEN_PLUS, 2) || is_func (h, GEN_MINUS, 2)) {
    if (is_func (h, GEN_PLUS, 2)) found= i;
    if (maximal == 0) return g;
    h= h[1];
    maximal--;
    i++;
  }
  if (found == -1) return g;
  syntactic abs_h; int sgn_h;
  signed_decompose (h, abs_h, sgn_h);
  if (sgn_h >= 0) return g;
  return migrate_sub (g, found);
}

/******************************************************************************
* Multiplicative arithmetic with simplifications
******************************************************************************/

static bool frac_flag= false;

bool
set_frac_flag (bool new_val) {
  bool ret= frac_flag;
  frac_flag= new_val;
  return ret;
}

syntactic
operator * (const syntactic& g1, const syntactic& g2) {
  if (frac_flag && is_func (g1, GEN_OVER, 2)) {
    if (is_func (g2, GEN_OVER, 2))
      return syn (GEN_OVER, g1[1] * g2[1], g1[2] * g2[2]);
    else return syn (GEN_OVER, g1[1] * g2, g1[2]);
  }
  if (frac_flag && is_func (g2, GEN_OVER, 2))
    return syn (GEN_OVER, g1 * g2[1], g2[2]);
  if (exact_eq (g1, 0)) return 0;
  if (exact_eq (g2, 0)) return 0;
  if (exact_eq (g1, 1)) return g2;
  if (exact_eq (g2, 1)) return g1;
  if (exact_eq (g1, -1)) return -g2;
  if (exact_eq (g2, -1)) return -g1;
  return syn (GEN_TIMES, g1, g2);
}

syntactic
operator / (const syntactic& g1, const syntactic& g2) {
  if (frac_flag && is_func (g2, GEN_OVER, 2))
    return g1 * syn (GEN_OVER, g2[2], g2[1]);
  if (exact_eq (g2, 1)) return g1;
  if (exact_eq (g2, -1)) return -g1;
  if (frac_flag) return g1 * syn (GEN_OVER, 1, g2);
  else return syn (GEN_OVER, g1, g2);
}

syntactic
square (const syntactic& g) {
  return syn (GEN_POWER, g, syntactic (2));
}

syntactic
invert (const syntactic& g) {
  return syntactic (1) / g;
}

syntactic
pow (const syntactic& g1, const syntactic& g2) {
  if (exact_eq (g1, 0)) return 0;
  if (exact_eq (g1, 1)) return 1;
  if (exact_eq (g2, 0)) return 1;
  if (exact_eq (g2, 1)) return g1;
  if (exact_eq (g2, syn (GEN_OVER, syntactic (1), syntactic (2))))
    return syn (GEN_SQRT, g1);
  if (frac_flag) {
    syntactic abs; int sgn;
    signed_decompose (g2, abs, sgn);
    if (sgn < 0) return syn (GEN_OVER, 1, pow (g1, abs));
  }
  return syn (GEN_POWER, g1, g2);
}

syntactic
pow (const syntactic& g1, const int& g2) {
  return pow (g1, flatten (g2));
}

syntactic
pow (const int& g1, const syntactic& g2) {
  return pow (flatten (g1), g2);
}

/******************************************************************************
* Sorting big sums and big products
******************************************************************************/

bool is_numeric (const string& s);

struct sum_less_op {
  static bool
  op (const syntactic& s1, const syntactic& s2) {
    generic g1= as_generic (s1);
    generic g2= as_generic (s2);
    if ((is_func (g1, GEN_TIMES, 2) || is_func (g1, GEN_OVER, 2)) &&
        is<literal> (g1[1]) && is_numeric (as_string (as<literal> (g1[1]))))
      g1= (is_func (g1, GEN_TIMES, 2)? g1[2]: gen (GEN_OVER, 1, g1[2]));
    if ((is_func (g2, GEN_TIMES, 2) || is_func (g2, GEN_OVER, 2)) &&
        is<literal> (g2[1]) && is_numeric (as_string (as<literal> (g2[1]))))
      g2= (is_func (g2, GEN_TIMES, 2)? g2[2]: gen (GEN_OVER, 1, g2[2]));
    return big_small_compare (g1, g2) <= 0;
  }
};

syntactic
ordered_sum (const vector<syntactic>& v) {
  vector<syntactic> w= copy (v);
  sort_leq<sum_less_op> (w);
  syntactic r= 0;
  for (nat i=0; i<N(w); i++) r= r + w[i];
  return migrate_negate (r, 5);
}

struct product_less_op {
  static bool
  op (const syntactic& s1, const syntactic& s2) {
    generic g1= as_generic (s1);
    generic g2= as_generic (s2);
    if (is_func (g1, GEN_OVER, 2) && exact_eq (g1[1], 1)) g1= g1[2];
    if (is_func (g2, GEN_OVER, 2) && exact_eq (g2[1], 1)) g2= g2[2];
    if (is_func (g1, GEN_POWER, 2) &&
        is<literal> (g1[2]) && is_numeric (as_string (as<literal> (g1[2]))))
      g1= g1[1];
    if (is_func (g2, GEN_POWER, 2) &&
        is<literal> (g2[2]) && is_numeric (as_string (as<literal> (g2[2]))))
      g2= g2[1];
    return small_big_compare (g1, g2) <= 0;
  }
};

syntactic
ordered_product (const vector<syntactic>& v) {
  vector<syntactic> w= copy (v);
  sort_leq<product_less_op> (w);
  bool old= set_frac_flag (true);
  syntactic r= 1;
  for (nat i=0; i<N(w); i++) r= r * w[i];
  (void) set_frac_flag (old);
  return r;
}

/******************************************************************************
* Gcd
******************************************************************************/

syntactic gcd (const syntactic& g1, const syntactic& g2) {
  if (exact_eq (g1, 0)) return g2;
  if (exact_eq (g2, 0)) return g1;
  return syn (GEN_GCD, g1, g2);
}

syntactic lcm (const syntactic& g1, const syntactic& g2) {
  if (exact_eq (g1, 0) || exact_eq (g2, 0)) return 0;
  return syn (GEN_LCM, g1, g2);
}

syntactic xgcd (const syntactic& g1, const syntactic& g2) {
  if (exact_eq (g1, 0))
    return exact_eq (g2, 0) ?
      syn (GEN_SQTUPLE, syn (GEN_ROW, 1, 0     ), syn (GEN_ROW, 0, 1)) :
      syn (GEN_SQTUPLE, syn (GEN_ROW, 0, 1 / g2), syn (GEN_ROW, 1, 0));
  return exact_eq (g2, 0) ?
    syn (GEN_SQTUPLE, syn (GEN_ROW, 1 / g1, 0), syn (GEN_ROW, 0  , 1 )) :
    syn (GEN_SQTUPLE, syn (GEN_ROW, 1 / g1, 0), syn (GEN_ROW, -g2, g1));
}

/******************************************************************************
* Other functions
******************************************************************************/

syntactic sqrt (const syntactic& g) { return syn (GEN_SQRT, g); }
syntactic exp (const syntactic& g) { return syn (GEN_EXP, g); }
syntactic log (const syntactic& g) { return syn (GEN_LOG, g); }
syntactic cos (const syntactic& g) { return syn (GEN_COS, g); }
syntactic sin (const syntactic& g) { return syn (GEN_SIN, g); }
syntactic tan (const syntactic& g) { return syn (GEN_TAN, g); }
syntactic acos (const syntactic& g) { return syn (GEN_ARCCOS, g); }
syntactic asin (const syntactic& g) { return syn (GEN_ARCSIN, g); }
syntactic atan (const syntactic& g) { return syn (GEN_ARCTAN, g); }
syntactic cosh (const syntactic& g) { return syn (GEN_CH, g); }
syntactic sinh (const syntactic& g) { return syn (GEN_SH, g); }
syntactic tanh (const syntactic& g) { return syn (GEN_TH, g); }
syntactic acosh (const syntactic& g) { return syn (GEN_ARGCH, g); }
syntactic asinh (const syntactic& g) { return syn (GEN_ARGSH, g); }
syntactic atanh (const syntactic& g) { return syn (GEN_ARGTH, g); }

syntactic trig (const syntactic& g) {
  return syn (GEN_SQTUPLE, cos (g), sin (g)); }
syntactic hypot (const syntactic& g1, const syntactic& g2) {
  return syn (GEN_HYPOT, g1, g2); }
syntactic atan2 (const syntactic& g1, const syntactic& g2) {
  return syn (GEN_ARCTAN2, g1, g2); }

syntactic Re (const syntactic& g1) { return syn (GEN_RE, g1); }
syntactic Im (const syntactic& g1) { return syn (GEN_IM, g1); }
syntactic abs (const syntactic& g1) { return syn (GEN_ABS, g1); }
syntactic arg (const syntactic& g1) { return syn (GEN_ARG, g1); }
syntactic conj (const syntactic& g1) { return syn (GEN_CONJ, g1); }
syntactic gaussian (const syntactic& g1, const syntactic& g2) {
  return g1 + g2 * Imaginary (syntactic); }
syntactic polar (const syntactic& g1, const syntactic& g2) {
  return g1 * exp (g2 * Imaginary (syntactic)); }

syntactic center (const syntactic& g1) { return syn (GEN_CENTER, g1); }
syntactic radius (const syntactic& g1) { return syn (GEN_RADIUS, g1); }
syntactic numerator (const syntactic& g1) {
  return syn (GEN_NUMERATOR, g1); }
syntactic denominator (const syntactic& g1) {
  return syn (GEN_DENOMINATOR, g1); }

syntactic operator << (const syntactic& x1, const syntactic& x2) {
  return syn (GEN_LESSLESS, x1, x2); }
syntactic operator >> (const syntactic& x1, const syntactic& x2) {
  return syn (GEN_GTRGTR, x1, x2); }
syntactic lshiftz (const syntactic& g, const syntactic& sh) {
  return syn ("lshiftz", g, sh); }
syntactic rshiftz (const syntactic& g, const syntactic& sh) {
  return syn ("rshiftz", g, sh); }

syntactic derive (const syntactic& g) {
  return syn (GEN_DERIVE, g); }
syntactic derive (const syntactic& g, const syntactic& v) {
  return syn (GEN_DERIVE, g, v); }
syntactic xderive (const syntactic& g) {
  return syn (GEN_XDERIVE, g); }
syntactic xderive (const syntactic& g, const syntactic& v) {
  return syn (GEN_XDERIVE, g, v); }
syntactic integrate (const syntactic& g) {
  return syn (GEN_INTEGRATE, g); }
syntactic integrate (const syntactic& g, const syntactic& v) {
  return syn (GEN_INTEGRATE, g, v); }

syntactic sqrt_init (const syntactic& x, const syntactic& c) {
  return syn ("sqrt_init", x, c); }
syntactic log_init (const syntactic& x, const syntactic& c) {
  return syn ("log_init", x, c); }
syntactic acos_init (const syntactic& x, const syntactic& c) {
  return syn ("acos_init", x, c); }
syntactic asin_init (const syntactic& x, const syntactic& c) {
  return syn ("asin_init", x, c); }
syntactic atan_init (const syntactic& x, const syntactic& c) {
  return syn ("atan_init", x, c); }
syntactic integrate_init (const syntactic& g, const syntactic& c) {
  return syn ("integrate_init", g, c); }
syntactic integrate_init (const syntactic& g, const syntactic& v,
                          const syntactic& c) {
  return syn ("integrate_init", g, v, c); }
syntactic solve_lde_init (const syntactic& x, const syntactic& c) {
  return syn ("solve_lde_init", x, c); }

syntactic access (const syntactic& g, const syntactic& i) {
  return syn (GEN_ACCESS, g, i); }
syntactic access (const syntactic& g, const syntactic& i, const syntactic& j) {
  return syn (GEN_ACCESS, g, i, j); }
syntactic inject (const syntactic& x, const syntactic& y, const syntactic& z) {
  return syn ("inject", x, y, z); }
syntactic project (const syntactic& x, const syntactic& y) {
  return syn ("project", x, y); }

/******************************************************************************
* Flattening related routines
******************************************************************************/

nat      the_series_prec = 10;
nat      the_series_zt   = 25;
bool     the_series_expr = false;

syntactic&
the_series_var () {
  static syntactic v ("z");
  return v;
}

syntactic&
the_polynomial_var () {
  static syntactic v ("x");
  return v;
}

syntactic flatten (const bool& b) {
  return b? syntactic ("true"): syntactic ("false"); }
syntactic flatten (const char& i) {
  return syntactic (as_string (i)); }
syntactic flatten (const signed char& i) {
  return syntactic (as_string (i)); }
syntactic flatten (const unsigned char& i) {
  return syntactic (as_string (i)); }
syntactic flatten (const short int& i) {
  return syntactic (as_string (i)); }
syntactic flatten (const short unsigned int& i) {
  return syntactic (as_string (i)); }
syntactic flatten (const int& i) {
  return syntactic (as_string (i)); }
syntactic flatten (const unsigned int& i) {
  return syntactic (as_string (i)); }
syntactic flatten (const long int& i) {
  return syntactic (as_string (i)); }
syntactic flatten (const long unsigned int& i) {
  return syntactic (as_string (i)); }
syntactic flatten (const long long int& i) {
  return syntactic (as_string (i)); }
syntactic flatten (const long long unsigned int& i) {
  return syntactic (as_string (i)); }
syntactic flatten (const float& x) {
  return syntactic (as_string (x)); }
syntactic flatten (const double& x) {
  return syntactic (as_string (x)); }
syntactic flatten (const long double& x) {
  return syntactic (as_string (x)); }
syntactic flatten (char* const & s) {
  return syn ("$text", syntactic (quote (s))); }
syntactic flatten (const string& s) {
  return syn ("$text", syntactic (quote (s))); }

} // namespace mmx