Merge branch 'master' of git://ffmssmsc.jinr.ru:443/varg/ginac
authorJens Vollinga <jensv@nikhef.nl>
Tue, 30 Sep 2008 09:47:25 +0000 (11:47 +0200)
committerJens Vollinga <jensv@nikhef.nl>
Tue, 30 Sep 2008 09:47:25 +0000 (11:47 +0200)
51 files changed:
check/Makefile.am
check/error_report.hpp [new file with mode: 0644]
check/exam_cra.cpp [new file with mode: 0644]
check/exam_mod_gcd.cpp [new file with mode: 0644]
check/match_bug.cpp [new file with mode: 0644]
check/parser_bugs.cpp [new file with mode: 0644]
check/time_parser.cpp
doc/examples/Makefile.am
doc/examples/derivative.cpp [new file with mode: 0644]
doc/examples/ginac-examples.texi
doc/tutorial/ginac.texi
ginac/Makefile.am
ginac/basic.cpp
ginac/basic.h
ginac/ex.cpp
ginac/ex.h
ginac/expairseq.cpp
ginac/expairseq.h
ginac/ginac.h
ginac/indexed.cpp
ginac/inifcns_nstdsums.cpp
ginac/input_lexer.h [deleted file]
ginac/input_lexer.ll [deleted file]
ginac/input_parser.yy [deleted file]
ginac/matrix.cpp
ginac/mul.cpp
ginac/ncmul.cpp
ginac/parser/lexer.cpp
ginac/parser/parse_context.cpp
ginac/parser/parse_context.hpp
ginac/parser/parser.cpp
ginac/parser/parser.hpp
ginac/parser/parser_compat.cpp [new file with mode: 0644]
ginac/polynomial/cra_garner.cpp [new file with mode: 0644]
ginac/polynomial/cra_garner.hpp [new file with mode: 0644]
ginac/polynomial/debug.hpp [new file with mode: 0644]
ginac/polynomial/gcd_euclid.tcc [new file with mode: 0644]
ginac/polynomial/mod_gcd.cpp [new file with mode: 0644]
ginac/polynomial/mod_gcd.hpp [new file with mode: 0644]
ginac/polynomial/normalize.tcc [new file with mode: 0644]
ginac/polynomial/remainder.tcc [new file with mode: 0644]
ginac/polynomial/ring_traits.hpp [new file with mode: 0644]
ginac/polynomial/upoly.hpp [new file with mode: 0644]
ginac/polynomial/upoly_io.cpp [new file with mode: 0644]
ginac/polynomial/upoly_io.hpp [new file with mode: 0644]
ginac/power.cpp
ginac/structure.h
ginac/wildcard.cpp
ginac/wildcard.h
ginsh/ginsh_parser.yy
tools/Makefile.am

index b656c631bc5628f1b44865964f87994f22c70a18..5c24aee640be7408798afe1949800e815fc89f3e 100644 (file)
@@ -7,6 +7,8 @@ CHECKS = check_numeric \
 
 EXAMS = exam_paranoia \
        exam_heur_gcd \
+       match_bug \
+       parser_bugs \
        exam_numeric_archive \
        exam_numeric \
        exam_powerlaws  \
@@ -25,7 +27,9 @@ EXAMS = exam_paranoia \
        exam_archive  \
        exam_structure  \
        exam_hashmap  \
-       exam_misc
+       exam_misc \
+       exam_mod_gcd \
+       exam_cra
 
 TIMES = time_dennyfliegner \
        time_gammaseries \
@@ -74,6 +78,12 @@ exam_paranoia_LDADD = ../ginac/libginac.la
 exam_heur_gcd_SOURCES = heur_gcd_bug.cpp 
 exam_heur_gcd_LDADD = ../ginac/libginac.la
 
+match_bug_SOURCES = match_bug.cpp error_report.hpp
+match_bug_LDADD = ../ginac/libginac.la
+
+parser_bugs_SOURCES = parser_bugs.cpp
+parser_bugs_LDADD = ../ginac/libginac.la
+
 exam_numeric_archive_SOURCES = numeric_archive.cpp
 exam_numeric_archive_LDADD = ../ginac/libginac.la
 
@@ -132,6 +142,11 @@ exam_hashmap_LDADD = ../ginac/libginac.la
 exam_misc_SOURCES = exam_misc.cpp
 exam_misc_LDADD = ../ginac/libginac.la
 
+exam_mod_gcd_SOURCES = exam_mod_gcd.cpp
+exam_mod_gcd_LDADD = ../ginac/libginac.la
+
+exam_cra_SOURCES = exam_cra.cpp
+exam_cra_LDADD = ../ginac/libginac.la
 
 time_dennyfliegner_SOURCES = time_dennyfliegner.cpp \
                             randomize_serials.cpp timer.cpp timer.h
@@ -233,6 +248,6 @@ time_parser_SOURCES = time_parser.cpp \
 time_parser_LDADD = ../ginac/libginac.la
 
 
-AM_CPPFLAGS = -I$(srcdir)/../ginac -I../ginac
+AM_CPPFLAGS = -I$(srcdir)/../ginac -I../ginac -DIN_GINAC
 
 CLEANFILES = exam.gar
diff --git a/check/error_report.hpp b/check/error_report.hpp
new file mode 100644 (file)
index 0000000..6b2dfe6
--- /dev/null
@@ -0,0 +1,17 @@
+#ifndef GINAC_CHECK_ERROR_REPORT_HPP
+#define GINAC_CHECK_ERROR_REPORT_HPP
+#include <sstream>
+#include <stdexcept>
+
+#define bug_on(cond, what)                             \
+do {                                                   \
+if (cond) {                                            \
+       std::ostringstream err_stream;                  \
+       err_stream << __FILE__ << ':' << __LINE__       \
+                  << what;                             \
+       throw std::logic_error(err_stream.str());       \
+}                                                      \
+} while (0)
+
+#endif // GINAC_CHECK_ERROR_REPORT_HPP
+
diff --git a/check/exam_cra.cpp b/check/exam_cra.cpp
new file mode 100644 (file)
index 0000000..11e1ab7
--- /dev/null
@@ -0,0 +1,120 @@
+#include <iostream>
+#include <cln/integer.h>
+#include <cln/integer_io.h>
+#include <cln/random.h>
+#include <cln/numtheory.h>
+#include "polynomial/cra_garner.hpp"
+#include <map>
+#include <vector>
+#include <stdexcept>
+using namespace cln;
+using namespace std;
+
+/// Generate a sequences of primes p_i such that \prod_i p_i < limit
+static std::vector<cln::cl_I>
+make_random_moduli(const cln::cl_I& limit);
+
+static std::vector<cln::cl_I>
+calc_residues(const cln::cl_I& x, const std::vector<cln::cl_I>& moduli);
+
+static void dump(const std::vector<cln::cl_I>& v);
+
+/// Make @a n random relatively prime moduli, each < limit, make a 
+/// random number x < \prod_{i=0}{n-1}, calculate residues, and 
+/// compute x' by chinese remainder algorithm. Check if the result
+/// of computation matches the original value x.
+static void run_test_once(const cln::cl_I& lim)
+{
+       std::vector<cln::cl_I> moduli = make_random_moduli(lim);
+       cln::cl_I x = random_I(lim) + 1;
+
+       if (x > (lim >> 1))
+               x = x - lim;
+
+       std::vector<cln::cl_I> residues = calc_residues(x, moduli);
+       cln::cl_I x_test;
+
+       bool error = false;
+       try {
+               x_test = integer_cra(residues, moduli);
+       } catch (std::exception& oops) {
+               std::cerr << "Oops: " << oops.what() << std::endl;
+               error = true;
+       }
+
+       if (x != x_test)
+               error = true;
+
+       if (error) {
+               std::cerr << "Expected x = " << x << ", got " <<
+                       x_test << " instead" << std::endl;
+               std::cerr << "moduli = ";
+               dump(moduli);
+               std::cerr << std::endl;
+               std::cerr << "residues = ";
+               dump(residues);
+               std::cerr << std::endl;
+               throw std::logic_error("bug in integer_cra?");
+       }
+}
+
+static void run_test(const cln::cl_I& limit, const std::size_t ntimes)
+{
+       for (std::size_t i = 0; i < ntimes; ++i)
+               run_test_once(limit);
+}
+
+int main(int argc, char** argv)
+{
+       typedef std::map<cln::cl_I, std::size_t> map_t;
+       map_t the_map;
+       // Run 1024 tests with native 32-bit numbers
+       the_map[cln::cl_I(std::numeric_limits<int>::max())] = 1024;
+
+       // Run 512 tests with native 64-bit integers
+       if (sizeof(long) > sizeof(int)) 
+               the_map[cln::cl_I(std::numeric_limits<long>::max())] = 512;
+
+       // Run 32 tests with a bit bigger numbers
+       the_map[cln::cl_I("987654321098765432109876543210")] = 32;
+
+       std::cout << "examining Garner's integer chinese remainder algorithm " << std::flush;
+
+       for (map_t::const_iterator i = the_map.begin(); i != the_map.end(); ++i)
+               run_test(i->first, i->second);
+
+       return 0;
+}
+
+static std::vector<cln::cl_I>
+calc_residues(const cln::cl_I& x, const std::vector<cln::cl_I>& moduli)
+{
+       std::vector<cln::cl_I> residues(moduli.size());
+       for (std::size_t i = moduli.size(); i-- != 0; )
+               residues[i] = mod(x, moduli[i]);
+       return residues;
+}
+
+static std::vector<cln::cl_I>
+make_random_moduli(const cln::cl_I& limit)
+{
+       std::vector<cln::cl_I> moduli;
+       cln::cl_I prod(1);
+       cln::cl_I next = random_I(std::min(limit >> 1, cln::cl_I(128)));
+       do {
+               cln::cl_I tmp = nextprobprime(next);
+               next = tmp + random_I(cln::cl_I(10)) + 1;
+               prod = prod*tmp;
+               moduli.push_back(tmp);
+       } while (prod < limit);
+       return moduli;
+}
+
+static void dump(const std::vector<cln::cl_I>& v)
+{
+       std::cerr << "[ ";
+       for (std::size_t i = 0; i < v.size(); ++i)
+               std::cerr << v[i] << " ";
+       std::cerr << "]";
+}
+
diff --git a/check/exam_mod_gcd.cpp b/check/exam_mod_gcd.cpp
new file mode 100644 (file)
index 0000000..4407d90
--- /dev/null
@@ -0,0 +1,85 @@
+#include "polynomial/upoly.hpp"
+#include "polynomial/upoly_io.hpp"
+#include "polynomial/mod_gcd.hpp"
+#include "ginac.h"
+#include <cln/random.h>
+#include <string>
+#include <iostream>
+#include <map>
+using namespace GiNaC;
+
+static upoly ex_to_upoly(const ex& e, const symbol& x);
+static ex upoly_to_ex(const upoly& p, const symbol& x);
+
+// make a univariate polynomial \in Z[x] of degree deg
+static upoly make_random_upoly(const std::size_t deg);
+
+static void run_test_once(const std::size_t deg)
+{
+       static const symbol xsym("x");
+
+       const upoly a = make_random_upoly(deg);
+       const upoly b = make_random_upoly(deg);
+
+       upoly g;
+       mod_gcd(g, a, b);
+
+       ex ea = upoly_to_ex(a, xsym);
+       ex eb = upoly_to_ex(b, xsym);
+       ex eg = gcd(ea, eb);
+
+       const upoly g_check = ex_to_upoly(eg, xsym);
+       if (g != g_check) {
+               std::cerr << "a = " << a << std::endl;
+               std::cerr << "b = " << b << std::endl;
+               std::cerr << "mod_gcd(a, b) = " << g << std::endl;
+               std::cerr << "sr_gcd(a, b) = " << g_check << std::endl;
+               throw std::logic_error("bug in mod_gcd");
+       }
+}
+
+int main(int argc, char** argv)
+{
+       std::cout << "examining modular gcd. ";
+       std::map<std::size_t, std::size_t> n_map;
+       // run 256 tests with polynomials of degree 10
+       n_map[10] = 256;
+       // run 32 tests with polynomials of degree 100
+       n_map[100] = 32;
+       std::map<std::size_t, std::size_t>::const_iterator i = n_map.begin();
+       for (; i != n_map.end(); ++i) {
+               for (std::size_t k = 0; k < i->second; ++k)
+                       run_test_once(i->first);
+       }
+       return 0;
+}
+
+static upoly ex_to_upoly(const ex& e, const symbol& x)
+{
+       upoly p(e.degree(x) + 1);
+       for (int i = 0; i <= e.degree(x); ++i)
+               p[i] = cln::the<cln::cl_I>(ex_to<numeric>(e.coeff(x, i)).to_cl_N());
+       return p;
+}
+
+static ex upoly_to_ex(const upoly& p, const symbol& x)
+{
+       exvector tv(p.size());
+       for (std::size_t i = 0; i < p.size(); ++i)
+               tv[i] = pow(x, i)*numeric(p[i]);
+       return (new add(tv))->setflag(status_flags::dynallocated);
+}
+
+static upoly make_random_upoly(const std::size_t deg)
+{
+       static const cln::cl_I biggish("98765432109876543210");
+       upoly p(deg + 1);
+       for (std::size_t i = 0; i <= deg; ++i)
+               p[i] = cln::random_I(biggish);
+
+       // Make sure the leading coefficient is non-zero
+       while (zerop(p[deg])) 
+               p[deg] = cln::random_I(biggish);
+       return p;
+}
+
diff --git a/check/match_bug.cpp b/check/match_bug.cpp
new file mode 100644 (file)
index 0000000..df67629
--- /dev/null
@@ -0,0 +1,66 @@
+/**
+ * @file match_bug.cpp
+ *
+ * Check for bug in GiNaC::ex::match() described here:
+ * http://www.ginac.de/pipermail/ginac-devel/2006-April/000942.html
+ */
+#include "ginac.h"
+#include "error_report.hpp"
+#include <iostream>
+using namespace GiNaC;
+
+/*
+ * basic::match(lst&) used to have an obscure side effect: repl_lst
+ * could be modified even if the match failed! Although this "feature"
+ * was documented it happend to be very confusing *even for GiNaC
+ * developers*, see 
+ * http://www.ginac.de/pipermail/ginac-devel/2006-April/000942.html
+ *
+ * It was fixed in 192ed7390b7b2b705ad100e3db0a92eedd2b20ad. Let's make
+ * sure it will be never re-added:
+ */
+static void failed_match_have_side_effects()
+{
+       symbol x("x");
+       ex e = pow(x, 5);
+       ex pattern = pow(wild(0), -1);
+       // obviously e does NOT match the pattern
+       exmap repls;
+       bool match_p = e.match(pattern, repls);
+       bug_on(match_p, "match(" << e << ", " << pattern << ") says \"Yes\"");
+       bug_on(repls.size() != 0,
+               "failed match have side effects: repls = " << repls);
+}
+
+/*
+ * As a consequence of the bug described above pattern matching can wrongly
+ * fail. In particular, x^5*y^(-1) fails to match ($0)^(-1)*x^($2).
+ *
+ * The first thing that is attempted to match is x^5 with $0^(-1). This match
+ * will fail. However repl_lst will contain $0 == x as a side effect. This
+ * repl_lst will prevent the match of y^(-1) to ($0)^(-1) to succeed.
+ *
+ * This issue was worked around by 73f0ce4cf8d91f073f35a45443f5fbe886921c5c.
+ * Now we have a real fix (192ed7390b7b2b705ad100e3db0a92eedd2b20ad), but
+ * let's add a check.
+ */
+static void match_false_negative()
+{
+       symbol x("x"), y("y");
+       ex e = pow(x, 5)*pow(y, -1);
+       ex pattern = pow(wild(0), -1)*pow(x, wild(2));
+       exmap repls;
+       bool match_p = e.match(pattern, repls);
+       bug_on(!match_p, "false negative: " << e << " did not match "
+                       << pattern);
+}
+
+int main(int argc, char** argv)
+{
+       std::cout << "checking for historical bugs in match()... " << std::flush;
+       failed_match_have_side_effects();
+       match_false_negative();
+       std::cout << "not found. ";
+       return 0;
+}
+
diff --git a/check/parser_bugs.cpp b/check/parser_bugs.cpp
new file mode 100644 (file)
index 0000000..af32c5e
--- /dev/null
@@ -0,0 +1,94 @@
+/// @file parser_a_b.cpp Check for some silly bugs in the parser.
+#include "ginac.h"
+#include <string>
+#include <iostream>
+#include <stdexcept>
+#include <sstream>
+using namespace GiNaC;
+
+// - a - b was misparsed as -a + b due to a bug in parser::parse_unary_expr()
+static int check1(std::ostream& err_str)
+{
+       const std::string srep("-a-b");
+       parser reader;
+       ex e = reader(srep);
+       ex a = reader.get_syms()["a"];
+       ex b = reader.get_syms()["b"];
+       ex g = - a - b;
+       ex d = (e - g).expand();
+       if (!d.is_zero()) {
+               err_str << "\"" << srep << "\" was misparsed as \""
+                       << e << "\"" << std::endl;
+               return 1;
+       }
+       return 0;
+}
+
+/// Parser was rejecting the valid expression '5 - (3*x)/10'.
+static int check2(std::ostream& err_str)
+{
+       const std::string srep("5-(3*x)/10");
+       parser reader;
+       ex e = reader(srep);
+       ex x = reader.get_syms()["x"];
+       ex g = 5 - (3*x)/10;
+       ex d = (e - g).expand();
+       if (!d.is_zero()) {
+               err_str << "\"" << srep << "\" was misparsed as \""
+                       << e << "\"" << std::endl;
+               return 1;
+       }
+       return 0;
+}
+
+/// parse_literal_expr forget to consume the token, so parser get
+/// very confused.
+static int check3(std::ostream& err_str)
+{
+       const std::string srep("5-(2*I)/3");
+       parser reader;
+       ex e = reader(srep);
+       ex g = numeric(5) - (numeric(2)*I)/3;
+       ex d = (e - g).expand();
+       if (!d.is_zero()) {
+               err_str << "\"" << srep << "\" was misparsed as \""
+                       << e << "\"" << std::endl;
+               return 1;
+       }
+       return 0;
+}
+
+/// parser happily accepted various junk like 'x^2()+1'
+static int check4(std::ostream& err_str)
+{
+       const std::string junk("x^2()+1");
+       parser reader;
+       ex e;
+       try {
+               e = reader(junk);
+               err_str << "parser accepts junk: \"" << junk << "\"" << std::endl;
+               return 1;
+       } catch (parse_error& err) {
+               // Ok, parser rejects the nonsense.
+               return 0;
+       }
+}
+
+int main(int argc, char** argv)
+{
+       std::cout << "checking for parser bugs. " << std::flush;
+       std::ostringstream err_str;
+       int errors = 0;
+       errors += check1(err_str);
+       errors += check2(err_str);
+       errors += check3(err_str);
+       errors += check4(err_str);
+       if (errors) {
+               std::cout << "Yes, unfortunately:" << std::endl;
+               std::cout << err_str.str();
+       } else {
+               std::cout << "Not found. ";
+       }
+       return errors;
+}
+
index ca35fb629fcb7ee98e4ec81ddef770b30cd84fe5..4e6081b194f923f2e87c57184e40f158344151dc 100644 (file)
@@ -6,7 +6,6 @@
 #include <vector>
 #include <stdexcept>
 #include "ginac.h"
-#include "parser/parser.hpp"
 #include "timer.h"
 extern void randomify_symbol_serials();
 using namespace std;
@@ -22,33 +21,14 @@ static string prepare_str(const unsigned n, const char x = 'x')
        return s.str();
 }
 
-void benchmark_and_cmp(const string& srep, double& t_new, double& t_old)
+static double benchmark_and_cmp(const string& srep)
 {
        parser the_parser;
        timer RSD10;
        RSD10.start();
        ex e = the_parser(srep);
-       t_new = RSD10.read();
-       RSD10.stop();
-       if (t_new > 2.0)
-               cout << '.' << flush;
-
-       symtab syms = the_parser.get_syms();
-       const symbol x = find_or_insert_symbol("x", syms, true);
-       lst sl;
-       sl = x;
-       RSD10.start();
-       ex e2(srep, sl);
-       t_old = RSD10.read();
-       
-       if (t_old > 2.0)
-               cout << '.' << flush;
-
-       ex dif = (e - e2).expand();
-       if (!dif.is_zero()) {
-               cerr << "Got a difference: " << dif << endl;
-               throw std::logic_error("New and old parser give different results");
-       }
+       const double t = RSD10.read();
+       return t;
 }
 
 int main(int argc, char** argv)
@@ -60,20 +40,18 @@ int main(int argc, char** argv)
        if (argc > 1)
                n_max = atoi(argv[1]);
 
-       vector<double> times_new, times_old;
+       vector<double> times;
        vector<unsigned> ns;
        for (unsigned n = n_min; n <= n_max; n = n << 1) {
-               double t_new, t_old;
                string srep = prepare_str(n);
-               benchmark_and_cmp(srep, t_new, t_old);
-               times_new.push_back(t_new);
-               times_old.push_back(t_old);
+               const double t = benchmark_and_cmp(srep);
+               times.push_back(t);
                ns.push_back(n);
        }
 
        cout << "OK" << endl;
-       cout << "# terms  new parser, s  old parser, s" << endl;
-       for (size_t i = 0; i < times_new.size(); i++)
-               cout << " " << ns[i] << '\t' << times_new[i] << '\t' << times_old[i] << endl;
+       cout << "# terms  time, s" << endl;
+       for (size_t i = 0; i < times.size(); i++)
+               cout << " " << ns[i] << '\t' << times[i] << endl;
        return 0;
 }
index ef8d4d1421ad57b38b27b84691df47ccf4f3eee0..2a5a10888d18e5cc97374c409f34fe16c60d8674 100644 (file)
@@ -1,7 +1,7 @@
 ## Process this file with automake to produce Makefile.in
 
 EXFILES = archive1.cpp compile1.cpp compile2.cpp compile3.cpp lanczos.cpp \
- mystring.cpp
+ mystring.cpp derivative.cpp
 
 TEXI = ginac-examples.texi
 
diff --git a/doc/examples/derivative.cpp b/doc/examples/derivative.cpp
new file mode 100644 (file)
index 0000000..c1a09b9
--- /dev/null
@@ -0,0 +1,30 @@
+// Input expression containing variable 'x' and compute its derivative
+// with respect to 'x'.
+// Example from the tutorial (chapter Input/Output, section `Expression
+// input').
+#include <iostream>
+#include <string>
+#include <stdexcept>
+#include <ginac/ginac.h>
+using namespace std;
+using namespace GiNaC;
+
+int main()
+{
+       cout << "Enter an expression containing 'x': " << flush;
+       parser reader;
+
+       try {
+               ex e = reader(cin);
+               symtab table = reader.get_syms();
+
+               symbol x = table.find("x") != table.end() ? 
+                          ex_to<symbol>(table["x"]) : symbol("x");
+
+               cout << "The derivative of " << e << " with respect to x is ";
+               cout << e.diff(x) << "." << endl;
+       } catch (exception &p) {
+               cerr << p.what() << endl;
+       }
+}
+
index 8870dafb84a8fde391a5937067176b602774f84a..2a109daa4c25496d559fb4d61913de88e1462b9f 100644 (file)
@@ -30,6 +30,13 @@ This is a collection of code examples using GiNaC.
 
 Two expression are stored in an archive on the disk and are restored again.
 
+@section Input and output of expressions
+
+@subsection Expression input @uref{derivative.cpp, (source)}
+
+Interactively input expression and compute its derivative with respect
+to the ``x'' variable.
+
 @chapter Monte Carlo Integration
 
 @section Example showing compile_ex @uref{compile1.cpp, (source)}
index 2f160c7bfb16459bd1d18fd7c46fb9c43273ad27..e17d4dfa8f4f59e974bcb0645bfff81f19e00b94 100644 (file)
@@ -4413,14 +4413,14 @@ matches a given pattern. This is done by the function
 
 @example
 bool ex::match(const ex & pattern);
-bool ex::match(const ex & pattern, lst & repls);
+bool ex::match(const ex & pattern, exmap& repls);
 @end example
 
 This function returns @code{true} when the expression matches the pattern
 and @code{false} if it doesn't. If used in the second form, the actual
-subexpressions matched by the wildcards get returned in the @code{repls}
-object as a list of relations of the form @samp{wildcard == expression}.
-If @code{match()} returns false,  @code{repls} remains unmodified.
+subexpressions matched by the wildcards get returned in the associative
+array @code{repls} with @samp{wildcard} as a key. If @code{match()}
+returns false,  @code{repls} remains unmodified.
 
 The matching algorithm works as follows:
 
@@ -4559,7 +4559,7 @@ Again some examples in @command{ginsh} for illustration (in @command{ginsh},
 The method
 
 @example
-bool ex::find(const ex & pattern, lst & found);
+bool ex::find(const ex & pattern, exset& found);
 @end example
 
 works a bit like @code{has()} but it doesn't stop upon finding the first
@@ -6380,22 +6380,72 @@ and have the @samp{x} and @samp{y} correspond to the symbols @code{x} and
 @code{y} you defined in your program and there is no way to specify the
 desired symbols to the @code{>>} stream input operator.
 
-Instead, GiNaC lets you construct an expression from a string, specifying the
-list of symbols to be used:
+Instead, GiNaC lets you read an expression from a stream or a string,
+specifying the mapping between the input strings and symbols to be used:
 
 @example
 @{
-    symbol x("x"), y("y");
-    ex e("2*x+sin(y)", lst(x, y));
+    symbol x, y;
+    symtab table;
+    table["x"] = x;
+    table["y"] = y;
+    parser reader(table);
+    ex e = reader("2*x+sin(y)");
 @}
 @end example
 
 The input syntax is the same as that used by @command{ginsh} and the stream
-output operator @code{<<}. The symbols in the string are matched by name to
-the symbols in the list and if GiNaC encounters a symbol not specified in
-the list it will throw an exception.
+output operator @code{<<}. Matching between the input strings and expressions
+is given by @samp{table}. The @samp{table} in this example instructs GiNaC
+to substitute any input substring ``x'' with symbol @code{x}. Likewise,
+the substring ``y'' will be replaced with symbol @code{y}. It's also possible
+to map input (sub)strings to arbitrary expressions:
 
-With this constructor, it's also easy to implement interactive GiNaC programs:
+@example
+@{
+    symbol x, y;
+    symtab table;
+    table["x"] = x+log(y)+1;
+    parser reader(table);
+    ex e = reader("5*x^3 - x^2");
+    // e = 5*(x+log(y)+1)^3 + (x+log(y)+1)^2
+@}
+@end example
+
+If no mapping is specified for a particular string GiNaC will create a symbol
+with corresponding name. Later on you can obtain all parser generated symbols
+with @code{get_syms()} method:
+
+@example
+@{
+    parser reader;
+    ex e = reader("2*x+sin(y)");
+    symtab table = reader.get_syms();
+    symbol x = reader["x"];
+    symbol y = reader["y"];
+@}
+@end example
+
+Sometimes you might want to prevent GiNaC from inserting these extra symbols
+(for example, you want treat an unexpected string in the input as an error).
+
+@example
+@{
+       symtab table;
+       table["x"] = symbol();
+       parser reader(table);
+       parser.strict = true;
+       ex e;
+       try @{
+               e = reader("2*x+sin(y)");
+       @} catch (parse_error& err) @{
+               cerr << err.what() << endl;
+               // prints "unknown symbol "y" in the input"
+       @}
+@}
+@end example
+
+With this parser, it's also easy to implement interactive GiNaC programs:
 
 @example
 #include <iostream>
@@ -6407,19 +6457,19 @@ using namespace GiNaC;
 
 int main()
 @{
-    symbol x("x");
-    string s;
-
-    cout << "Enter an expression containing 'x': ";
-    getline(cin, s);
-
-    try @{
-        ex e(s, lst(x));
-        cout << "The derivative of " << e << " with respect to x is ";
-        cout << e.diff(x) << ".\n";
-    @} catch (exception &p) @{
-        cerr << p.what() << endl;
-    @}
+       cout << "Enter an expression containing 'x': " << flush;
+       parser reader;
+
+       try @{
+               ex e = reader(cin);
+               symtab table = reader.get_syms();
+               symbol x = table.find("x") != table.end() ? 
+                          ex_to<symbol>(table["x"]) : symbol("x");
+               cout << "The derivative of " << e << " with respect to x is ";
+               cout << e.diff(x) << "." << endl;
+       @} catch (exception &p) @{
+               cerr << p.what() << endl;
+       @}
 @}
 @end example
 
index 17b272f230f142d9dba9597b81ea3f19ab2fe719..b683d9b410c752e69a1266fda9d1f9d32342d1f5 100644 (file)
@@ -8,16 +8,29 @@ libginac_la_SOURCES = add.cpp archive.cpp basic.cpp clifford.cpp color.cpp \
   integral.cpp lst.cpp matrix.cpp mul.cpp ncmul.cpp normal.cpp numeric.cpp \
   operators.cpp power.cpp registrar.cpp relational.cpp remember.cpp \
   pseries.cpp print.cpp symbol.cpp symmetry.cpp tensor.cpp \
-  utils.cpp wildcard.cpp input_parser.yy input_lexer.ll \
-  input_lexer.h remember.h tostring.h utils.h compiler.h \
+  utils.cpp wildcard.cpp \
+  remember.h tostring.h utils.h compiler.h \
   parser/parse_binop_rhs.cpp \
   parser/parser.cpp \
   parser/parse_context.cpp \
   parser/builtin_fcns.cpp \
   parser/lexer.cpp \
   parser/lexer.hpp \
-  parser/debug.hpp
-  
+  parser/parser_compat.cpp \
+  parser/debug.hpp \
+polynomial/mod_gcd.cpp \
+polynomial/cra_garner.cpp \
+polynomial/gcd_euclid.tcc \
+polynomial/remainder.tcc \
+polynomial/normalize.tcc \
+polynomial/upoly.hpp \
+polynomial/ring_traits.hpp \
+polynomial/mod_gcd.hpp \
+polynomial/cra_garner.hpp \
+polynomial/upoly_io.hpp \
+polynomial/upoly_io.cpp \
+polynomial/debug.hpp
+
 libginac_la_LDFLAGS = -version-info $(LT_VERSION_INFO) -release $(LT_RELEASE)
 libginac_la_LIBADD = $(DL_LIBS)
 ginacincludedir = $(includedir)/ginac
@@ -30,9 +43,7 @@ ginacinclude_HEADERS = ginac.h add.h archive.h assertion.h basic.h class_info.h
   parser/parser.hpp \
   parser/parse_context.hpp
 
-AM_LFLAGS = -Pginac_yy -olex.yy.c
-AM_YFLAGS = -p ginac_yy -d
-EXTRA_DIST = function.pl input_parser.h version.h.in \
+EXTRA_DIST = function.pl version.h.in \
 parser/default_reader.tpl parser/builtin_fcns.def
 
 # Files produced by autogen(1) from templates
index 7a0633de19a18221f241292b21742416cdb29dce..e8c893ba26dd59a963e00fd1056daf34f4f47bf8 100644 (file)
@@ -290,7 +290,7 @@ ex & basic::operator[](size_t i)
  *  but e.has(x+y) is false. */
 bool basic::has(const ex & pattern, unsigned options) const
 {
-       lst repl_lst;
+       exmap repl_lst;
        if (match(pattern, repl_lst))
                return true;
        for (size_t i=0; i<nops(); i++)
@@ -534,9 +534,9 @@ bool basic::contract_with(exvector::iterator self, exvector::iterator other, exv
 }
 
 /** Check whether the expression matches a given pattern. For every wildcard
- *  object in the pattern, an expression of the form "wildcard == matching_expression"
- *  is added to repl_lst. */
-bool basic::match(const ex & pattern, lst & repl_lst) const
+ *  object in the pattern, a pair with the wildcard as a key and matching 
+ *  expression as a value is added to repl_lst. */
+bool basic::match(const ex & pattern, exmap& repl_lst) const
 {
 /*
        Sweet sweet shapes, sweet sweet shapes,
@@ -559,11 +559,11 @@ bool basic::match(const ex & pattern, lst & repl_lst) const
                // 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 (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
-                       if (it->op(0).is_equal(pattern))
-                               return is_equal(ex_to<basic>(it->op(1)));
+               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.append(pattern == *this);
+               repl_lst[pattern] = *this;
                return true;
 
        } else {
@@ -589,7 +589,7 @@ bool basic::match(const ex & pattern, lst & repl_lst) const
                // 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.
-               lst tmp_repl = repl_lst;
+               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))
@@ -614,9 +614,10 @@ ex basic::subs_one_level(const exmap & m, unsigned options) const
                return thisex;
        } else {
                for (it = m.begin(); it != m.end(); ++it) {
-                       lst repl_lst;
+                       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 it->second.subs(repl_lst, options | subs_options::no_pattern);
+                       // avoid infinite recursion when re-substituting the wildcards
                }
        }
 
index 348523791a19afa50abd1c84006c49fe003758c8..ffc508dbee0dd48ae4b4d6e18b18b96a4f68b1da 100644 (file)
@@ -163,7 +163,7 @@ public:
 
        // pattern matching
        virtual bool has(const ex & other, unsigned options = 0) const;
-       virtual bool match(const ex & pattern, lst & repl_lst) const;
+       virtual bool match(const ex & pattern, exmap & repls) const;
 protected:
        virtual bool match_same_type(const basic & other) const;
 public:
index 23a30833b5a436375918e34358a4208376cf00d6..25d5607fdf2d6ba92ce2ebe1d0ae524c025f2318 100644 (file)
@@ -32,7 +32,6 @@
 #include "power.h"
 #include "lst.h"
 #include "relational.h"
-#include "input_lexer.h"
 #include "utils.h"
 
 namespace GiNaC {
@@ -95,7 +94,7 @@ ex ex::diff(const symbol & s, unsigned nth) const
 /** Check whether expression matches a specified pattern. */
 bool ex::match(const ex & pattern) const
 {
-       lst repl_lst;
+       exmap repl_lst;
        return bp->match(pattern, repl_lst);
 }
 
@@ -103,12 +102,10 @@ bool ex::match(const ex & pattern) const
  *  the "found" list. If the expression itself matches the pattern, the
  *  children are not further examined. This function returns true when any
  *  matches were found. */
-bool ex::find(const ex & pattern, lst & found) const
+bool ex::find(const ex & pattern, exset& found) const
 {
        if (match(pattern)) {
-               found.append(*this);
-               found.sort();
-               found.unique();
+               found.insert(*this);
                return true;
        }
        bool any_found = false;
@@ -555,17 +552,6 @@ basic & ex::construct_from_double(double d)
        return *bp;
 }
 
-ptr<basic> ex::construct_from_string_and_lst(const std::string &s, const ex &l)
-{
-       set_lexer_string(s);
-       set_lexer_symbols(l);
-       ginac_yyrestart(NULL);
-       if (ginac_yyparse())
-               throw (std::runtime_error(get_parser_error()));
-       else
-               return parsed_ex.bp;
-}
-       
 //////////
 // static member variables
 //////////
index 783c89ee2fc204cdaed11f95d1477e1009180db0..4a603ce03778d441a9bfb1209cfdb0b3995bce2a 100644 (file)
@@ -141,9 +141,9 @@ public:
 
        // pattern matching
        bool has(const ex & pattern, unsigned options = 0) const { return bp->has(pattern, options); }
-       bool find(const ex & pattern, lst & found) const;
+       bool find(const ex & pattern, exset& found) const;
        bool match(const ex & pattern) const;
-       bool match(const ex & pattern, lst & repl_lst) const { return bp->match(pattern, repl_lst); }
+       bool match(const ex & pattern, exmap & repls) const { return bp->match(pattern, repls); }
 
        // substitutions
        ex subs(const exmap & m, unsigned options = 0) const;
@@ -702,7 +702,7 @@ inline ex imag_part(const ex & thisex)
 inline bool has(const ex & thisex, const ex & pattern, unsigned options = 0)
 { return thisex.has(pattern, options); }
 
-inline bool find(const ex & thisex, const ex & pattern, lst & found)
+inline bool find(const ex & thisex, const ex & pattern, exset& found)
 { return thisex.find(pattern, found); }
 
 inline bool is_polynomial(const ex & thisex, const ex & vars)
@@ -762,7 +762,7 @@ inline ex diff(const ex & thisex, const symbol & s, unsigned nth = 1)
 inline ex series(const ex & thisex, const ex & r, int order, unsigned options = 0)
 { return thisex.series(r, order, options); }
 
-inline bool match(const ex & thisex, const ex & pattern, lst & repl_lst)
+inline bool match(const ex & thisex, const ex & pattern, exmap& repl_lst)
 { return thisex.match(pattern, repl_lst); }
 
 inline ex simplify_indexed(const ex & thisex, unsigned options = 0)
index 3ed9e2892a23040c8bbb8e2ce4bf3ebfed36705f..67aa4f8ef681bf3453eb83fecc0181076edba70e 100644 (file)
@@ -391,7 +391,7 @@ bool expairseq::is_polynomial(const ex & var) const
        return true;
 }
 
-bool expairseq::match(const ex & pattern, lst & repl_lst) const
+bool expairseq::match(const ex & pattern, exmap & repl_lst) const
 {
        // This differs from basic::match() because we want "a+b+c+d" to
        // match "d+*+b" with "*" being "a+c", and we want to honor commutativity
@@ -427,20 +427,10 @@ bool expairseq::match(const ex & pattern, lst & repl_lst) const
                                continue;
                        exvector::iterator it = ops.begin(), itend = ops.end();
                        while (it != itend) {
-                               lst::const_iterator last_el = repl_lst.end();
-                               --last_el;
                                if (it->match(p, repl_lst)) {
                                        ops.erase(it);
                                        goto found;
                                }
-                               while(true) {
-                                       lst::const_iterator next_el = last_el;
-                                       ++next_el;
-                                       if(next_el == repl_lst.end())
-                                               break;
-                                       else
-                                               repl_lst.remove_last();
-                               }
                                ++it;
                        }
                        return false; // no match found
@@ -458,11 +448,11 @@ found:            ;
                        for (size_t i=0; i<num; i++)
                                vp->push_back(split_ex_to_pair(ops[i]));
                        ex rest = thisexpairseq(vp, default_overall_coeff());
-                       for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
-                               if (it->op(0).is_equal(global_wildcard))
-                                       return rest.is_equal(it->op(1));
+                       for (exmap::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it) {
+                               if (it->first.is_equal(global_wildcard))
+                                       return rest.is_equal(it->second);
                        }
-                       repl_lst.append(global_wildcard == rest);
+                       repl_lst[global_wildcard] = rest;
                        return true;
 
                } else {
index 6abd10074eb2036bf55c76ab8724416423d4284f..5f707f7679e19e9780da6586d3e0613deca1dd02 100644 (file)
@@ -84,7 +84,7 @@ public:
        ex eval(int level=0) const;
        ex to_rational(exmap & repl) const;
        ex to_polynomial(exmap & repl) const;
-       bool match(const ex & pattern, lst & repl_lst) const;
+       bool match(const ex & pattern, exmap& repl_lst) const;
        ex subs(const exmap & m, unsigned options = 0) const;
        ex conjugate() const;
        bool is_polynomial(const ex & var) const;
index ad96e6403c057f6ee01d4f5375711a6ed343111f..2d2a9e9d5280916acaff2f55ecc99410f172d58b 100644 (file)
 
 #include "excompiler.h"
 
+#ifndef IN_GINAC
+#include "parser.hpp"
+#else
+#include "parser/parser.hpp"
+#endif
+
 #ifdef __MAKECINT__
 #pragma link C++ nestedclass;
 #pragma link C++ nestedtypedef;
index c8b710ae7f1e63690cbb0356e2e1ccc10159f1e0..aa8500399e270facb7fa72a59a6c9277c26214fa 100644 (file)
@@ -1409,7 +1409,7 @@ exvector get_all_dummy_indices_safely(const ex & e)
        else if (is_a<mul>(e) || is_a<ncmul>(e)) {
                exvector dummies;
                exvector free_indices;
-               for (int i=0; i<e.nops(); ++i) {
+               for (std::size_t i = 0; i < e.nops(); ++i) {
                        exvector dummies_of_factor = get_all_dummy_indices_safely(e.op(i));
                        dummies.insert(dummies.end(), dummies_of_factor.begin(),
                                dummies_of_factor.end());
@@ -1425,7 +1425,7 @@ exvector get_all_dummy_indices_safely(const ex & e)
        }
        else if(is_a<add>(e)) {
                exvector result;
-               for(int i=0; i<e.nops(); ++i) {
+               for(std::size_t i = 0; i < e.nops(); ++i) {
                        exvector dummies_of_term = get_all_dummy_indices_safely(e.op(i));
                        sort(dummies_of_term.begin(), dummies_of_term.end());
                        exvector new_vec;
index d85044d04bd1f311b2d26ad5ab8ccf0572d1ebe7..84a78c42c8e45d6c7b2ee76a9475fd42922f59d5 100644 (file)
@@ -501,19 +501,6 @@ cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl
 }
 
 
-// converts parameter types and calls multipleLi_do_sum (convenience function for G_numeric)
-cln::cl_N mLi_do_summation(const lst& m, const lst& x)
-{
-       std::vector<int> m_int;
-       std::vector<cln::cl_N> x_cln;
-       for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
-               m_int.push_back(ex_to<numeric>(*itm).to_int());
-               x_cln.push_back(ex_to<numeric>(*itx).to_cl_N());
-       }
-       return multipleLi_do_sum(m_int, x_cln);
-}
-
-
 // forward declaration for Li_eval()
 lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf);
 
@@ -990,199 +977,228 @@ ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2
             + shuffle_G(a02, a1, a2_removed, pendint, a_old, scale, gsyms);
 }
 
-
 // handles the transformations and the numerical evaluation of G
 // the parameter x, s and y must only contain numerics
-ex G_numeric(const lst& x, const lst& s, const ex& y)
+static cln::cl_N
+G_numeric(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
+         const cln::cl_N& y);
+
+// do acceleration transformation (hoelder convolution [BBB])
+// the parameter x, s and y must only contain numerics
+static cln::cl_N
+G_do_hoelder(std::vector<cln::cl_N> x, /* yes, it's passed by value */
+            const std::vector<int>& s, const cln::cl_N& y)
 {
-       // check for convergence and necessary accelerations
-       bool need_trafo = false;
-       bool need_hoelder = false;
-       int depth = 0;
-       for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
-               if (!(*it).is_zero()) {
-                       ++depth;
-                       if (abs(*it) - y < -pow(10,-Digits+1)) {
-                               need_trafo = true;
+       cln::cl_N result;
+       const std::size_t size = x.size();
+       for (std::size_t i = 0; i < size; ++i)
+               x[i] = x[i]/y;
+
+       for (std::size_t r = 0; r <= size; ++r) {
+               cln::cl_N buffer(1 & r ? -1 : 1);
+               cln::cl_RA p(2);
+               bool adjustp;
+               do {
+                       adjustp = false;
+                       for (std::size_t i = 0; i < size; ++i) {
+                               if (x[i] == cln::cl_RA(1)/p) {
+                                       p = p/2 + cln::cl_RA(3)/2;
+                                       adjustp = true;
+                                       continue;
+                               }
                        }
-                       if (abs((abs(*it) - y)/y) < 0.01) {
-                               need_hoelder = true;
+               } while (adjustp);
+               cln::cl_RA q = p/(p-1);
+               std::vector<cln::cl_N> qlstx;
+               std::vector<int> qlsts;
+               for (std::size_t j = r; j >= 1; --j) {
+                       qlstx.push_back(cln::cl_N(1) - x[j-1]);
+                       if (instanceof(x[j-1], cln::cl_R_ring) &&
+                           realpart(x[j-1]) > 1 && realpart(x[j-1]) <= 2) {
+                               qlsts.push_back(s[j-1]);
+                       } else {
+                               qlsts.push_back(-s[j-1]);
                        }
                }
-       }
-       if (x.op(x.nops()-1).is_zero()) {
-               need_trafo = true;
-       }
-       if (depth == 1 && x.nops() == 2 && !need_trafo) {
-               return -Li(x.nops(), y / x.op(x.nops()-1)).evalf();
-       }
-       
-       // do acceleration transformation (hoelder convolution [BBB])
-       if (need_hoelder) {
-               
-               ex result;
-               const int size = x.nops();
-               lst newx;
-               for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
-                       newx.append(*it / y);
+               if (qlstx.size() > 0) {
+                       buffer = buffer*G_numeric(qlstx, qlsts, 1/q);
                }
-               
-               for (int r=0; r<=size; ++r) {
-                       ex buffer = pow(-1, r);
-                       ex p = 2;
-                       bool adjustp;
-                       do {
-                               adjustp = false;
-                               for (lst::const_iterator it = newx.begin(); it != newx.end(); ++it) {
-                                       if (*it == 1/p) {
-                                               p += (3-p)/2; 
-                                               adjustp = true;
-                                               continue;
-                                       }
-                               }
-                       } while (adjustp);
-                       ex q = p / (p-1);
-                       lst qlstx;
-                       lst qlsts;
-                       for (int j=r; j>=1; --j) {
-                               qlstx.append(1-newx.op(j-1));
-                               if (newx.op(j-1).info(info_flags::real) && newx.op(j-1) > 1 && newx.op(j-1) <= 2) {
-                                       qlsts.append( s.op(j-1));
-                               } else {
-                                       qlsts.append( -s.op(j-1));
-                               }
-                       }
-                       if (qlstx.nops() > 0) {
-                               buffer *= G_numeric(qlstx, qlsts, 1/q);
-                       }
-                       lst plstx;
-                       lst plsts;
-                       for (int j=r+1; j<=size; ++j) {
-                               plstx.append(newx.op(j-1));
-                               plsts.append(s.op(j-1));
-                       }
-                       if (plstx.nops() > 0) {
-                               buffer *= G_numeric(plstx, plsts, 1/p);
-                       }
-                       result += buffer;
+               std::vector<cln::cl_N> plstx;
+               std::vector<int> plsts;
+               for (std::size_t j = r+1; j <= size; ++j) {
+                       plstx.push_back(x[j-1]);
+                       plsts.push_back(s[j-1]);
                }
-               return result;
+               if (plstx.size() > 0) {
+                       buffer = buffer*G_numeric(plstx, plsts, 1/p);
+               }
+               result = result + buffer;
        }
-       
-       // convergence transformation
-       if (need_trafo) {
-
-               // sort (|x|<->position) to determine indices
-               std::multimap<ex,int> sortmap;
-               int size = 0;
-               for (int i=0; i<x.nops(); ++i) {
-                       if (!x[i].is_zero()) {
-                               sortmap.insert(std::pair<ex,int>(abs(x[i]), i));
-                               ++size;
-                       }
-               }
-               // include upper limit (scale)
-               sortmap.insert(std::pair<ex,int>(abs(y), x.nops()));
-
-               // generate missing dummy-symbols
-               int i = 1;
-               // holding dummy-symbols for the G/Li transformations
-               exvector gsyms;
-               gsyms.push_back(symbol("GSYMS_ERROR"));
-               ex lastentry;
-               for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
-                       if (it != sortmap.begin()) {
-                               if (it->second < x.nops()) {
-                                       if (x[it->second] == lastentry) {
-                                               gsyms.push_back(gsyms.back());
-                                               continue;
-                                       }
-                               } else {
-                                       if (y == lastentry) {
-                                               gsyms.push_back(gsyms.back());
-                                               continue;
-                                       }
+       return result;
+}
+
+// convergence transformation, used for numerical evaluation of G function.
+// the parameter x, s and y must only contain numerics
+static cln::cl_N
+G_do_trafo(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
+          const cln::cl_N& y)
+{
+       // sort (|x|<->position) to determine indices
+       typedef std::multimap<cln::cl_R, std::size_t> sortmap_t;
+       sortmap_t sortmap;
+       std::size_t size = 0;
+       for (std::size_t i = 0; i < x.size(); ++i) {
+               if (!zerop(x[i])) {
+                       sortmap.insert(std::make_pair(abs(x[i]), i));
+                       ++size;
+               }
+       }
+       // include upper limit (scale)
+       sortmap.insert(std::make_pair(abs(y), x.size()));
+
+       // generate missing dummy-symbols
+       int i = 1;
+       // holding dummy-symbols for the G/Li transformations
+       exvector gsyms;
+       gsyms.push_back(symbol("GSYMS_ERROR"));
+       cln::cl_N lastentry(0);
+       for (sortmap_t::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
+               if (it != sortmap.begin()) {
+                       if (it->second < x.size()) {
+                               if (x[it->second] == lastentry) {
+                                       gsyms.push_back(gsyms.back());
+                                       continue;
                                }
-                       }
-                       std::ostringstream os;
-                       os << "a" << i;
-                       gsyms.push_back(symbol(os.str()));
-                       ++i;
-                       if (it->second < x.nops()) {
-                               lastentry = x[it->second];
                        } else {
-                               lastentry = y;
+                               if (y == lastentry) {
+                                       gsyms.push_back(gsyms.back());
+                                       continue;
+                               }
                        }
                }
+               std::ostringstream os;
+               os << "a" << i;
+               gsyms.push_back(symbol(os.str()));
+               ++i;
+               if (it->second < x.size()) {
+                       lastentry = x[it->second];
+               } else {
+                       lastentry = y;
+               }
+       }
 
-               // fill position data according to sorted indices and prepare substitution list
-               Gparameter a(x.nops());
-               lst subslst;
-               int pos = 1;
-               int scale;
-               for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
-                       if (it->second < x.nops()) {
-                               if (s[it->second] > 0) {
-                                       a[it->second] = pos;
-                               } else {
-                                       a[it->second] = -pos;
-                               }
-                               subslst.append(gsyms[pos] == x[it->second]);
+       // fill position data according to sorted indices and prepare substitution list
+       Gparameter a(x.size());
+       exmap subslst;
+       std::size_t pos = 1;
+       int scale;
+       for (sortmap_t::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
+               if (it->second < x.size()) {
+                       if (s[it->second] > 0) {
+                               a[it->second] = pos;
                        } else {
-                               scale = pos;
-                               subslst.append(gsyms[pos] == y);
+                               a[it->second] = -int(pos);
                        }
-                       ++pos;
+                       subslst[gsyms[pos]] = numeric(x[it->second]);
+               } else {
+                       scale = pos;
+                       subslst[gsyms[pos]] = numeric(y);
                }
+               ++pos;
+       }
 
-               // do transformation
-               Gparameter pendint;
-               ex result = G_transform(pendint, a, scale, gsyms);
-               // replace dummy symbols with their values
-               result = result.eval().expand();
-               result = result.subs(subslst).evalf();
-               
-               return result;
+       // do transformation
+       Gparameter pendint;
+       ex result = G_transform(pendint, a, scale, gsyms);
+       // replace dummy symbols with their values
+       result = result.eval().expand();
+       result = result.subs(subslst).evalf();
+       if (!is_a<numeric>(result))
+               throw std::logic_error("G_do_trafo: G_transform returned non-numeric result");
+       
+       cln::cl_N ret = ex_to<numeric>(result).to_cl_N();
+       return ret;
+}
+
+// handles the transformations and the numerical evaluation of G
+// the parameter x, s and y must only contain numerics
+static cln::cl_N
+G_numeric(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
+         const cln::cl_N& y)
+{
+       // check for convergence and necessary accelerations
+       bool need_trafo = false;
+       bool need_hoelder = false;
+       std::size_t depth = 0;
+       for (std::size_t i = 0; i < x.size(); ++i) {
+               if (!zerop(x[i])) {
+                       ++depth;
+                       const cln::cl_N x_y = abs(x[i]) - y;
+                       if (instanceof(x_y, cln::cl_R_ring) &&
+                           realpart(x_y) < cln::least_negative_float(cln::float_format(Digits - 2)))
+                               need_trafo = true;
+
+                       if (abs(abs(x[i]/y) - 1) < 0.01)
+                               need_hoelder = true;
+               }
        }
+       if (zerop(x[x.size() - 1]))
+               need_trafo = true;
+
+       if (depth == 1 && x.size() == 2 && !need_trafo)
+               return - Li_projection(2, y/x[1], cln::float_format(Digits));
+       
+       // do acceleration transformation (hoelder convolution [BBB])
+       if (need_hoelder)
+               return G_do_hoelder(x, s, y);
+       
+       // convergence transformation
+       if (need_trafo)
+               return G_do_trafo(x, s, y);
 
        // do summation
-       lst newx;
-       lst m;
+       std::vector<cln::cl_N> newx;
+       newx.reserve(x.size());
+       std::vector<int> m;
+       m.reserve(x.size());
        int mcount = 1;
-       ex sign = 1;
-       ex factor = y;
-       for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
-               if ((*it).is_zero()) {
+       int sign = 1;
+       cln::cl_N factor = y;
+       for (std::size_t i = 0; i < x.size(); ++i) {
+               if (zerop(x[i])) {
                        ++mcount;
                } else {
-                       newx.append(factor / (*it));
-                       factor = *it;
-                       m.append(mcount);
+                       newx.push_back(factor/x[i]);
+                       factor = x[i];
+                       m.push_back(mcount);
                        mcount = 1;
                        sign = -sign;
                }
        }
 
-       return sign * numeric(mLi_do_summation(m, newx));
+       return sign*multipleLi_do_sum(m, newx);
 }
 
 
 ex mLi_numeric(const lst& m, const lst& x)
 {
        // let G_numeric do the transformation
-       lst newx;
-       lst s;
-       ex factor = 1;
+       std::vector<cln::cl_N> newx;
+       newx.reserve(x.nops());
+       std::vector<int> s;
+       s.reserve(x.nops());
+       cln::cl_N factor(1);
        for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
                for (int i = 1; i < *itm; ++i) {
-                       newx.append(0);
-                       s.append(1);
+                       newx.push_back(cln::cl_N(0));
+                       s.push_back(1);
                }
-               newx.append(factor / *itx);
-               factor /= *itx;
-               s.append(1);
+               const cln::cl_N xi = ex_to<numeric>(*itx).to_cl_N();
+               newx.push_back(factor/xi);
+               factor = factor/xi;
+               s.push_back(1);
        }
-       return pow(-1, m.nops()) * G_numeric(newx, s, _ex1);
+       return numeric(cln::cl_N(1 & m.nops() ? - 1 : 1)*G_numeric(newx, s, cln::cl_N(1)));
 }
 
 
@@ -1210,7 +1226,8 @@ static ex G2_evalf(const ex& x_, const ex& y)
        if (x.op(0) == y) {
                return G(x_, y).hold();
        }
-       lst s;
+       std::vector<int> s;
+       s.reserve(x.nops());
        bool all_zero = true;
        for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
                if (!(*it).info(info_flags::numeric)) {
@@ -1220,16 +1237,21 @@ static ex G2_evalf(const ex& x_, const ex& y)
                        all_zero = false;
                }
                if ( !ex_to<numeric>(*it).is_real() && ex_to<numeric>(*it).imag() < 0 ) {
-                       s.append(-1);
+                       s.push_back(-1);
                }
                else {
-                       s.append(+1);
+                       s.push_back(1);
                }
        }
        if (all_zero) {
                return pow(log(y), x.nops()) / factorial(x.nops());
        }
-       return G_numeric(x, s, y);
+       std::vector<cln::cl_N> xv;
+       xv.reserve(x.nops());
+       for (lst::const_iterator it = x.begin(); it != x.end(); ++it)
+               xv.push_back(ex_to<numeric>(*it).to_cl_N());
+       cln::cl_N result = G_numeric(xv, s, ex_to<numeric>(y).to_cl_N());
+       return numeric(result);
 }
 
 
@@ -1247,7 +1269,8 @@ static ex G2_eval(const ex& x_, const ex& y)
        if (x.op(0) == y) {
                return G(x_, y).hold();
        }
-       lst s;
+       std::vector<int> s;
+       s.reserve(x.nops());
        bool all_zero = true;
        bool crational = true;
        for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
@@ -1261,10 +1284,10 @@ static ex G2_eval(const ex& x_, const ex& y)
                        all_zero = false;
                }
                if ( !ex_to<numeric>(*it).is_real() && ex_to<numeric>(*it).imag() < 0 ) {
-                       s.append(-1);
+                       s.push_back(-1);
                }
                else {
-                       s.append(+1);
+                       s.push_back(+1);
                }
        }
        if (all_zero) {
@@ -1276,7 +1299,12 @@ static ex G2_eval(const ex& x_, const ex& y)
        if (crational) {
                return G(x_, y).hold();
        }
-       return G_numeric(x, s, y);
+       std::vector<cln::cl_N> xv;
+       xv.reserve(x.nops());
+       for (lst::const_iterator it = x.begin(); it != x.end(); ++it)
+               xv.push_back(ex_to<numeric>(*it).to_cl_N());
+       cln::cl_N result = G_numeric(xv, s, ex_to<numeric>(y).to_cl_N());
+       return numeric(result);
 }
 
 
@@ -1306,7 +1334,8 @@ static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
        if (x.op(0) == y) {
                return G(x_, s_, y).hold();
        }
-       lst sn;
+       std::vector<int> sn;
+       sn.reserve(s.nops());
        bool all_zero = true;
        for (lst::const_iterator itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
                if (!(*itx).info(info_flags::numeric)) {
@@ -1320,25 +1349,30 @@ static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
                }
                if ( ex_to<numeric>(*itx).is_real() ) {
                        if ( *its >= 0 ) {
-                               sn.append(+1);
+                               sn.push_back(1);
                        }
                        else {
-                               sn.append(-1);
+                               sn.push_back(-1);
                        }
                }
                else {
                        if ( ex_to<numeric>(*itx).imag() > 0 ) {
-                               sn.append(+1);
+                               sn.push_back(1);
                        }
                        else {
-                               sn.append(-1);
+                               sn.push_back(-1);
                        }
                }
        }
        if (all_zero) {
                return pow(log(y), x.nops()) / factorial(x.nops());
        }
-       return G_numeric(x, sn, y);
+       std::vector<cln::cl_N> xn;
+       xn.reserve(x.nops());
+       for (lst::const_iterator it = x.begin(); it != x.end(); ++it)
+               xn.push_back(ex_to<numeric>(*it).to_cl_N());
+       cln::cl_N result = G_numeric(xn, sn, ex_to<numeric>(y).to_cl_N());
+       return numeric(result);
 }
 
 
@@ -1360,7 +1394,8 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
        if (x.op(0) == y) {
                return G(x_, s_, y).hold();
        }
-       lst sn;
+       std::vector<int> sn;
+       sn.reserve(s.nops());
        bool all_zero = true;
        bool crational = true;
        for (lst::const_iterator itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
@@ -1378,18 +1413,18 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
                }
                if ( ex_to<numeric>(*itx).is_real() ) {
                        if ( *its >= 0 ) {
-                               sn.append(+1);
+                               sn.push_back(1);
                        }
                        else {
-                               sn.append(-1);
+                               sn.push_back(-1);
                        }
                }
                else {
                        if ( ex_to<numeric>(*itx).imag() > 0 ) {
-                               sn.append(+1);
+                               sn.push_back(1);
                        }
                        else {
-                               sn.append(-1);
+                               sn.push_back(-1);
                        }
                }
        }
@@ -1402,7 +1437,12 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
        if (crational) {
                return G(x_, s_, y).hold();
        }
-       return G_numeric(x, sn, y);
+       std::vector<cln::cl_N> xn;
+       xn.reserve(x.nops());
+       for (lst::const_iterator it = x.begin(); it != x.end(); ++it)
+               xn.push_back(ex_to<numeric>(*it).to_cl_N());
+       cln::cl_N result = G_numeric(xn, sn, ex_to<numeric>(y).to_cl_N());
+       return numeric(result);
 }
 
 
@@ -2241,7 +2281,7 @@ bool convert_parameter_H_to_Li(const lst& l, lst& m, lst& s, ex& pf)
                }
        }
        if (has_negative_parameters) {
-               for (int i=0; i<m.nops(); i++) {
+               for (std::size_t i=0; i<m.nops(); i++) {
                        if (m.op(i) < 0) {
                                m.let_op(i) = -m.op(i);
                                s.append(-1);
@@ -2281,7 +2321,7 @@ struct map_trafo_H_convert_to_Li : public map_function
                                        s.let_op(0) = s.op(0) * arg;
                                        return pf * Li(m, s).hold();
                                } else {
-                                       for (int i=0; i<m.nops(); i++) {
+                                       for (std::size_t i=0; i<m.nops(); i++) {
                                                s.append(1);
                                        }
                                        s.let_op(0) = s.op(0) * arg;
@@ -2363,7 +2403,7 @@ struct map_trafo_H_reduce_trailing_zeros : public map_function
                                        
                                        //
                                        parameter.remove_last();
-                                       int lastentry = parameter.nops();
+                                       std::size_t lastentry = parameter.nops();
                                        while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
                                                lastentry--;
                                        }
@@ -2452,9 +2492,9 @@ ex trafo_H_mult(const ex& h1, const ex& h2)
                        hlong = h2.op(0).op(0);
                }
        }
-       for (int i=0; i<=hlong.nops(); i++) {
+       for (std::size_t i=0; i<=hlong.nops(); i++) {
                lst newparameter;
-               int j=0;
+               std::size_t j=0;
                for (; j<i; j++) {
                        newparameter.append(hlong[j]);
                }
@@ -2482,7 +2522,7 @@ struct map_trafo_H_mult : public map_function
                        ex result = 1;
                        ex firstH;
                        lst Hlst;
-                       for (int pos=0; pos<e.nops(); pos++) {
+                       for (std::size_t pos=0; pos<e.nops(); pos++) {
                                if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
                                        std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
                                        if (name == "H") {
@@ -2516,7 +2556,7 @@ struct map_trafo_H_mult : public map_function
                        if (Hlst.nops() > 0) {
                                ex buffer = trafo_H_mult(firstH, Hlst.op(0));
                                result *= buffer;
-                               for (int i=1; i<Hlst.nops(); i++) {
+                               for (std::size_t i=1; i<Hlst.nops(); i++) {
                                        result *= Hlst.op(i);
                                }
                                result = result.expand();
@@ -2544,7 +2584,7 @@ ex trafo_H_1tx_prepend_zero(const ex& e, const ex& arg)
        if (name == "H") {
                h = e;
        } else {
-               for (int i=0; i<e.nops(); i++) {
+               for (std::size_t i=0; i<e.nops(); i++) {
                        if (is_a<function>(e.op(i))) {
                                std::string name = ex_to<function>(e.op(i)).get_name();
                                if (name == "H") {
@@ -2576,7 +2616,7 @@ ex trafo_H_prepend_one(const ex& e, const ex& arg)
        if (name == "H") {
                h = e;
        } else {
-               for (int i=0; i<e.nops(); i++) {
+               for (std::size_t i=0; i<e.nops(); i++) {
                        if (is_a<function>(e.op(i))) {
                                std::string name = ex_to<function>(e.op(i)).get_name();
                                if (name == "H") {
@@ -2607,7 +2647,7 @@ ex trafo_H_1tx_prepend_minusone(const ex& e, const ex& arg)
        if (name == "H") {
                h = e;
        } else {
-               for (int i=0; i<e.nops(); i++) {
+               for (std::size_t i=0; i<e.nops(); i++) {
                        if (is_a<function>(e.op(i))) {
                                std::string name = ex_to<function>(e.op(i)).get_name();
                                if (name == "H") {
@@ -2640,7 +2680,7 @@ ex trafo_H_1mxt1px_prepend_minusone(const ex& e, const ex& arg)
        if (name == "H") {
                h = e;
        } else {
-               for (int i=0; i<e.nops(); i++) {
+               for (std::size_t i = 0; i < e.nops(); i++) {
                        if (is_a<function>(e.op(i))) {
                                std::string name = ex_to<function>(e.op(i)).get_name();
                                if (name == "H") {
@@ -2671,7 +2711,7 @@ ex trafo_H_1mxt1px_prepend_one(const ex& e, const ex& arg)
        if (name == "H") {
                h = e;
        } else {
-               for (int i=0; i<e.nops(); i++) {
+               for (std::size_t i = 0; i < e.nops(); i++) {
                        if (is_a<function>(e.op(i))) {
                                std::string name = ex_to<function>(e.op(i)).get_name();
                                if (name == "H") {
@@ -2709,7 +2749,7 @@ struct map_trafo_H_1mx : public map_function
                                // special cases if all parameters are either 0, 1 or -1
                                bool allthesame = true;
                                if (parameter.op(0) == 0) {
-                                       for (int i=1; i<parameter.nops(); i++) {
+                                       for (std::size_t i = 1; i < parameter.nops(); i++) {
                                                if (parameter.op(i) != 0) {
                                                        allthesame = false;
                                                        break;
@@ -2725,7 +2765,7 @@ struct map_trafo_H_1mx : public map_function
                                } else if (parameter.op(0) == -1) {
                                        throw std::runtime_error("map_trafo_H_1mx: cannot handle weights equal -1!");
                                } else {
-                                       for (int i=1; i<parameter.nops(); i++) {
+                                       for (std::size_t i = 1; i < parameter.nops(); i++) {
                                                if (parameter.op(i) != 1) {
                                                        allthesame = false;
                                                        break;
@@ -2751,7 +2791,7 @@ struct map_trafo_H_1mx : public map_function
                                        map_trafo_H_1mx recursion;
                                        ex buffer = recursion(H(newparameter, arg).hold());
                                        if (is_a<add>(buffer)) {
-                                               for (int i=0; i<buffer.nops(); i++) {
+                                               for (std::size_t i = 0; i < buffer.nops(); i++) {
                                                        res -= trafo_H_prepend_one(buffer.op(i), arg);
                                                }
                                        } else {
@@ -2765,13 +2805,13 @@ struct map_trafo_H_1mx : public map_function
                                        map_trafo_H_1mx recursion;
                                        map_trafo_H_mult unify;
                                        ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
-                                       int firstzero = 0;
+                                       std::size_t firstzero = 0;
                                        while (parameter.op(firstzero) == 1) {
                                                firstzero++;
                                        }
-                                       for (int i=firstzero-1; i<parameter.nops()-1; i++) {
+                                       for (std::size_t i = firstzero-1; i < parameter.nops()-1; i++) {
                                                lst newparameter;
-                                               int j=0;
+                                               std::size_t j=0;
                                                for (; j<=i; j++) {
                                                        newparameter.append(parameter[j+1]);
                                                }
@@ -2810,7 +2850,7 @@ struct map_trafo_H_1overx : public map_function
                                // special cases if all parameters are either 0, 1 or -1
                                bool allthesame = true;
                                if (parameter.op(0) == 0) {
-                                       for (int i=1; i<parameter.nops(); i++) {
+                                       for (std::size_t i = 1; i < parameter.nops(); i++) {
                                                if (parameter.op(i) != 0) {
                                                        allthesame = false;
                                                        break;
@@ -2820,7 +2860,7 @@ struct map_trafo_H_1overx : public map_function
                                                return pow(-1, parameter.nops()) * H(parameter, 1/arg).hold();
                                        }
                                } else if (parameter.op(0) == -1) {
-                                       for (int i=1; i<parameter.nops(); i++) {
+                                       for (std::size_t i = 1; i < parameter.nops(); i++) {
                                                if (parameter.op(i) != -1) {
                                                        allthesame = false;
                                                        break;
@@ -2832,7 +2872,7 @@ struct map_trafo_H_1overx : public map_function
                                                       / factorial(parameter.nops())).expand());
                                        }
                                } else {
-                                       for (int i=1; i<parameter.nops(); i++) {
+                                       for (std::size_t i = 1; i < parameter.nops(); i++) {
                                                if (parameter.op(i) != 1) {
                                                        allthesame = false;
                                                        break;
@@ -2855,7 +2895,7 @@ struct map_trafo_H_1overx : public map_function
                                        map_trafo_H_1overx recursion;
                                        ex buffer = recursion(H(newparameter, arg).hold());
                                        if (is_a<add>(buffer)) {
-                                               for (int i=0; i<buffer.nops(); i++) {
+                                               for (std::size_t i = 0; i < buffer.nops(); i++) {
                                                        res += trafo_H_1tx_prepend_zero(buffer.op(i), arg);
                                                }
                                        } else {
@@ -2870,7 +2910,7 @@ struct map_trafo_H_1overx : public map_function
                                        map_trafo_H_1overx recursion;
                                        ex buffer = recursion(H(newparameter, arg).hold());
                                        if (is_a<add>(buffer)) {
-                                               for (int i=0; i<buffer.nops(); i++) {
+                                               for (std::size_t i = 0; i < buffer.nops(); i++) {
                                                        res += trafo_H_1tx_prepend_zero(buffer.op(i), arg) - trafo_H_1tx_prepend_minusone(buffer.op(i), arg);
                                                }
                                        } else {
@@ -2884,13 +2924,13 @@ struct map_trafo_H_1overx : public map_function
                                        map_trafo_H_1overx recursion;
                                        map_trafo_H_mult unify;
                                        ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
-                                       int firstzero = 0;
+                                       std::size_t firstzero = 0;
                                        while (parameter.op(firstzero) == 1) {
                                                firstzero++;
                                        }
-                                       for (int i=firstzero-1; i<parameter.nops()-1; i++) {
+                                       for (std::size_t i = firstzero-1; i < parameter.nops() - 1; i++) {
                                                lst newparameter;
-                                               int j=0;
+                                               std::size_t j = 0;
                                                for (; j<=i; j++) {
                                                        newparameter.append(parameter[j+1]);
                                                }
@@ -2931,7 +2971,7 @@ struct map_trafo_H_1mxt1px : public map_function
                                // special cases if all parameters are either 0, 1 or -1
                                bool allthesame = true;
                                if (parameter.op(0) == 0) {
-                                       for (int i=1; i<parameter.nops(); i++) {
+                                       for (std::size_t i = 1; i < parameter.nops(); i++) {
                                                if (parameter.op(i) != 0) {
                                                        allthesame = false;
                                                        break;
@@ -2943,7 +2983,7 @@ struct map_trafo_H_1mxt1px : public map_function
                                                       / factorial(parameter.nops())).expand());
                                        }
                                } else if (parameter.op(0) == -1) {
-                                       for (int i=1; i<parameter.nops(); i++) {
+                                       for (std::size_t i = 1; i < parameter.nops(); i++) {
                                                if (parameter.op(i) != -1) {
                                                        allthesame = false;
                                                        break;
@@ -2955,7 +2995,7 @@ struct map_trafo_H_1mxt1px : public map_function
                                                       / factorial(parameter.nops())).expand());
                                        }
                                } else {
-                                       for (int i=1; i<parameter.nops(); i++) {
+                                       for (std::size_t i = 1; i < parameter.nops(); i++) {
                                                if (parameter.op(i) != 1) {
                                                        allthesame = false;
                                                        break;
@@ -2978,7 +3018,7 @@ struct map_trafo_H_1mxt1px : public map_function
                                        map_trafo_H_1mxt1px recursion;
                                        ex buffer = recursion(H(newparameter, arg).hold());
                                        if (is_a<add>(buffer)) {
-                                               for (int i=0; i<buffer.nops(); i++) {
+                                               for (std::size_t i = 0; i < buffer.nops(); i++) {
                                                        res -= trafo_H_1mxt1px_prepend_one(buffer.op(i), arg) + trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
                                                }
                                        } else {
@@ -2993,7 +3033,7 @@ struct map_trafo_H_1mxt1px : public map_function
                                        map_trafo_H_1mxt1px recursion;
                                        ex buffer = recursion(H(newparameter, arg).hold());
                                        if (is_a<add>(buffer)) {
-                                               for (int i=0; i<buffer.nops(); i++) {
+                                               for (std::size_t i = 0; i < buffer.nops(); i++) {
                                                        res -= trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
                                                }
                                        } else {
@@ -3007,13 +3047,13 @@ struct map_trafo_H_1mxt1px : public map_function
                                        map_trafo_H_1mxt1px recursion;
                                        map_trafo_H_mult unify;
                                        ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
-                                       int firstzero = 0;
+                                       std::size_t firstzero = 0;
                                        while (parameter.op(firstzero) == 1) {
                                                firstzero++;
                                        }
-                                       for (int i=firstzero-1; i<parameter.nops()-1; i++) {
+                                       for (std::size_t i = firstzero - 1; i < parameter.nops() - 1; i++) {
                                                lst newparameter;
-                                               int j=0;
+                                               std::size_t j=0;
                                                for (; j<=i; j++) {
                                                        newparameter.append(parameter[j+1]);
                                                }
@@ -3087,7 +3127,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
                        }
                }
 
-               for (int i=0; i<x1.nops(); i++) {
+               for (std::size_t i = 0; i < x1.nops(); i++) {
                        if (!x1.op(i).info(info_flags::integer)) {
                                return H(x1, x2).hold();
                        }
@@ -3159,7 +3199,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
                // ensure that the realpart of the argument is positive
                if (cln::realpart(x) < 0) {
                        x = -x;
-                       for (int i=0; i<m.nops(); i++) {
+                       for (std::size_t i = 0; i < m.nops(); i++) {
                                if (m.op(i) != 0) {
                                        m.let_op(i) = -m.op(i);
                                        res *= -1;
@@ -3668,7 +3708,7 @@ cln::cl_N zeta_do_Hoelder_convolution(const std::vector<int>& m_, const std::vec
        s_p[0] = s_p[0] * cln::cl_N("1/2");
        // convert notations
        int sig = 1;
-       for (int i=0; i<s_.size(); i++) {
+       for (std::size_t i = 0; i < s_.size(); i++) {
                if (s_[i] < 0) {
                        sig = -sig;
                        s_p[i] = -s_p[i];
diff --git a/ginac/input_lexer.h b/ginac/input_lexer.h
deleted file mode 100644 (file)
index 8cffe6b..0000000
+++ /dev/null
@@ -1,69 +0,0 @@
-/** @file input_lexer.h
- *
- *  Lexical analyzer definition for reading expressions.
- *  This file must be processed with flex. */
-
-/*
- *  GiNaC Copyright (C) 1999-2008 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
- */
-
-#ifndef __GINAC_INPUT_LEXER_H__
-#define __GINAC_INPUT_LEXER_H__
-
-extern "C" {
-#include <stdio.h>
-}
-       
-#include "config.h"
-
-// yacc stack type
-#define YYSTYPE ex
-
-// lex functions/variables
-extern int ginac_yyerror(char *s);
-extern int ginac_yylex(void);
-extern void ginac_yyrestart(FILE *f);
-extern char *ginac_yytext;
-
-namespace GiNaC {
-
-class ex;
-
-/** Set the input string to be parsed by ginac_yyparse() (used internally). */
-extern void set_lexer_string(const std::string &s);
-
-/** Get name of symbol/index (used internally). */
-extern std::string get_symbol_name(const ex & s);
-
-/** Set the list of predefined symbols for the lexer (used internally). */
-extern void set_lexer_symbols(ex l);
-
-/** Check whether lexer symbol was predefined (vs. created by the lexer, e.g. function names). */
-extern bool is_lexer_symbol_predefined(const ex &s);
-
-/** The expression parser function (used internally). */
-extern int ginac_yyparse(void);
-
-/** The expression returned by the parser (used internally). */
-extern ex parsed_ex;
-
-/** Get error message from the parser. */
-extern std::string get_parser_error(void);
-
-} // namespace GiNaC
-
-#endif // ndef __GINAC_INPUT_LEXER_H__
diff --git a/ginac/input_lexer.ll b/ginac/input_lexer.ll
deleted file mode 100644 (file)
index 0a7bbec..0000000
+++ /dev/null
@@ -1,211 +0,0 @@
-/** @file input_lexer.ll
- *
- *  Lexical analyzer definition for reading expressions.
- *  This file must be processed with flex. */
-
-/*
- *  GiNaC Copyright (C) 1999-2008 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
- */
-
-
-/*
- *  Definitions
- */
-
-%pointer
-
-%{
-#include <iostream>
-#include <string>
-#include <map>
-#include <stdexcept>
-
-#include "input_lexer.h"
-#include "ex.h"
-#include "constant.h"
-#include "fail.h"
-#include "numeric.h"
-#include "symbol.h"
-#include "lst.h"
-#include "idx.h"
-
-using namespace GiNaC;
-namespace GiNaC {
-
-#include "input_parser.h"
-
-} // namespace GiNaC
-
-// Table of all used symbols/indices
-struct sym_def {
-       sym_def() : predefined(false) {}
-       sym_def(const ex &s, bool predef) : sym(s), predefined(predef) {}
-       ~sym_def() {}
-
-       sym_def(const sym_def &other) {sym = other.sym; predefined = other.predefined;}
-       const sym_def &operator=(const sym_def &other)
-       {
-               if (this != &other) {
-                       sym = other.sym;
-                       predefined = other.predefined;
-               }
-               return *this;
-       }
-
-       ex sym;
-       bool predefined;        // true = user supplied symbol, false = lexer generated symbol
-};
-typedef std::map<std::string, sym_def> sym_tab;
-static sym_tab syms;
-
-// lex input function
-static int lexer_input(char *buf, int max_size);
-#define YY_INPUT(buf, result, max_size) (result = lexer_input(buf, max_size))
-%}
-
-       /* Abbreviations */
-D      [0-9]
-E      [elEL][-+]?{D}+
-A      [a-zA-Z_]
-AN     [0-9a-zA-Z_]
-
-
-/*
- *  Lexical rules
- */
-
-%%
-[ \t\n]+               /* skip whitespace */
-
-                       /* special values */
-Pi                     ginac_yylval = Pi; return T_LITERAL;
-Euler                  ginac_yylval = Euler; return T_LITERAL;
-Catalan                        ginac_yylval = Catalan; return T_LITERAL;
-FAIL                   ginac_yylval = *new fail(); return T_LITERAL;
-I                      ginac_yylval = I; return T_NUMBER;
-Digits                 ginac_yylval = (long)Digits; return T_DIGITS;
-
-                       /* comparison */
-"=="                   return T_EQUAL;
-"!="                   return T_NOTEQ;
-"<="                   return T_LESSEQ;
-">="                   return T_GREATEREQ;
-
-                       /* numbers */
-{D}+                   |
-"#"{D}+"R"{AN}+                |
-"#b"([01])+            |
-"#o"[0-7]+             |
-"#x"[0-9a-fA-F]+       |
-{D}+"."{D}*({E})?      |
-{D}*"."{D}+({E})?      |
-{D}+{E}                        ginac_yylval = numeric(yytext); return T_NUMBER;
-
-                       /* symbols */
-{A}{AN}*               {
-                               sym_tab::const_iterator i = syms.find(yytext);
-                               if (i == syms.end()) {
-                                       symbol tmp(yytext);
-                                       ginac_yylval = tmp;
-                                       syms[yytext] = sym_def(tmp, false);
-                               } else
-                                       ginac_yylval = (*i).second.sym;
-                               return T_SYMBOL;
-                       }
-
-                       /* end of input */
-<<EOF>>                        return T_EOF;
-
-                       /* everything else */
-.                      return *yytext;
-
-%%
-
-
-/*
- *  Routines
- */
-
-// The string from which we will read
-static std::string lexer_string;
-
-// The current position within the string
-static int curr_pos = 0;
-
-// Input function that reads from string
-static int lexer_input(char *buf, int max_size)
-{
-       int actual = lexer_string.length() - curr_pos;
-       if (actual > max_size)
-               actual = max_size;
-       if (actual <= 0)
-               return YY_NULL;
-       lexer_string.copy(buf, actual, curr_pos);
-       curr_pos += actual;
-       return actual;
-}
-
-// EOF encountered, terminate the scanner
-int ginac_yywrap()
-{
-       return 1;
-}
-
-namespace GiNaC {
-
-// Set the input string
-void set_lexer_string(const std::string &s)
-{
-       lexer_string = s;
-       curr_pos = 0;
-}
-
-// Get name of symbol/index
-std::string get_symbol_name(const ex & s)
-{
-       if (is_a<symbol>(s))
-               return ex_to<symbol>(s).get_name();
-       else if (is_a<idx>(s) && is_a<symbol>(s.op(0)))
-               return ex_to<symbol>(s.op(0)).get_name();
-       else
-               throw (std::runtime_error("get_symbol_name(): unexpected expression type"));
-}
-
-// Set the list of predefined symbols/indices
-void set_lexer_symbols(ex l)
-{
-       syms.clear();
-       if (!is_exactly_a<lst>(l))
-               return;
-       for (unsigned i=0; i<l.nops(); i++) {
-               const ex &o = l.op(i);
-               if (is_a<symbol>(o) || (is_a<idx>(o) && is_a<symbol>(o.op(0))))
-                       syms[get_symbol_name(o)] = sym_def(o, true);
-       }
-}
-
-// Check whether symbol/index was predefined
-bool is_lexer_symbol_predefined(const ex &s)
-{
-       sym_tab::const_iterator i = syms.find(get_symbol_name(s));
-       if (i == syms.end())
-               return false;
-       else
-               return (*i).second.predefined;
-}
-
-} // namespace GiNaC
diff --git a/ginac/input_parser.yy b/ginac/input_parser.yy
deleted file mode 100644 (file)
index b28e63a..0000000
+++ /dev/null
@@ -1,201 +0,0 @@
-/** @file input_parser.yy
- *
- *  Input grammar definition for reading expressions.
- *  This file must be processed with yacc/bison. */
-
-/*
- *  GiNaC Copyright (C) 1999-2008 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
- */
-
-
-/*
- *  Definitions
- */
-
-%{
-#include <stdexcept>
-
-#include "ex.h"
-#include "input_lexer.h"
-#include "relational.h"
-#include "operators.h"
-#include "symbol.h"
-#include "lst.h"
-#include "power.h"
-#include "exprseq.h"
-#include "idx.h"
-#include "indexed.h"
-#include "matrix.h"
-#include "inifcns.h"
-
-namespace GiNaC {
-
-#define YYERROR_VERBOSE 1
-
-// Parsed output expression
-ex parsed_ex;
-
-// Last error message returned by parser
-static std::string parser_error;
-
-// Prototypes
-ex attach_index(const ex & base, ex i, bool covariant);
-%}
-
-/* Tokens (T_LITERAL means a literal value returned by the parser, but not
-   of class numeric or symbol (e.g. a constant or the FAIL object)) */
-%token T_EOF T_NUMBER T_SYMBOL T_LITERAL T_DIGITS T_EQUAL T_NOTEQ T_LESSEQ T_GREATEREQ
-
-/* Operator precedence and associativity */
-%right '='
-%left T_EQUAL T_NOTEQ
-%left '<' '>' T_LESSEQ T_GREATEREQ
-%left '+' '-'
-%left '*' '/' '%'
-%nonassoc NEG
-%right '^'
-%left '.' '~'
-%nonassoc '!'
-
-%start input
-
-
-/*
- *  Grammar rules
- */
-
-%%
-input  : exp T_EOF {
-               try {
-                       parsed_ex = $1;
-                       YYACCEPT;
-               } catch (std::exception &err) {
-                       parser_error = err.what();
-                       YYERROR;
-               }
-       }
-       ;
-
-exp    : T_NUMBER              {$$ = $1;}
-       | T_SYMBOL {
-               if (is_lexer_symbol_predefined($1))
-                       $$ = $1.eval();
-               else
-                       throw (std::runtime_error("unknown symbol '" + get_symbol_name($1) + "'"));
-       }
-       | T_LITERAL             {$$ = $1;}
-       | T_DIGITS              {$$ = $1;}
-       | T_SYMBOL '(' exprseq ')' {
-               std::string n = get_symbol_name($1);
-               if (n == "sqrt") {
-                       if ($3.nops() != 1)
-                               throw (std::runtime_error("too many arguments to sqrt()"));
-                       $$ = sqrt($3.op(0));
-               } else if (n == "pow" || n == "power") {
-                 if ($3.nops() != 2) 
-                         throw std::invalid_argument("wrong number of arguments to pow()");
-                       $$ = power($3.op(0), $3.op(0));
-               } else {
-                       unsigned i = function::find_function(n, $3.nops());
-                       $$ = function(i, ex_to<exprseq>($3)).eval(1);
-               }
-       }
-       | exp T_EQUAL exp       {$$ = $1 == $3;}
-       | exp T_NOTEQ exp       {$$ = $1 != $3;}
-       | exp '<' exp           {$$ = $1 < $3;}
-       | exp T_LESSEQ exp      {$$ = $1 <= $3;}
-       | exp '>' exp           {$$ = $1 > $3;}
-       | exp T_GREATEREQ exp   {$$ = $1 >= $3;}
-       | exp '+' exp           {$$ = $1 + $3;}
-       | exp '-' exp           {$$ = $1 - $3;}
-       | exp '*' exp           {$$ = $1 * $3;}
-       | exp '/' exp           {$$ = $1 / $3;}
-       | '-' exp %prec NEG     {$$ = -$2;}
-       | '+' exp %prec NEG     {$$ = $2;}
-       | exp '^' exp           {$$ = pow($1, $3);}
-       | exp '.' exp           {$$ = attach_index($1, $3, true);}
-       | exp '~' exp           {$$ = attach_index($1, $3, false);}
-       | exp '!'               {$$ = factorial($1);}
-       | '(' exp ')'           {$$ = $2;}
-       | '{' list_or_empty '}' {$$ = $2;}
-       | '[' matrix ']'        {$$ = lst_to_matrix(ex_to<lst>($2));}
-       ;
-
-exprseq        : exp                   {$$ = exprseq($1);}
-       | exprseq ',' exp       {exprseq es(ex_to<exprseq>($1)); $$ = es.append($3);}
-       ;
-
-list_or_empty: /* empty */     {$$ = *new lst;}
-       | list                  {$$ = $1;}
-       ;
-
-list   : exp                   {$$ = lst($1);}
-       | list ',' exp          {lst l(ex_to<lst>($1)); $$ = l.append($3);}
-       ;
-
-matrix : '[' row ']'           {$$ = lst($2);}
-       | matrix ',' '[' row ']' {lst l(ex_to<lst>($1)); $$ = l.append($4);}
-       ;
-
-row    : exp                   {$$ = lst($1);}
-       | row ',' exp           {lst l(ex_to<lst>($1)); $$ = l.append($3);}
-       ;
-
-
-/*
- *  Routines
- */
-
-%%
-// Attach index to expression
-ex attach_index(const ex & base, ex i, bool covariant)
-{
-       // Toggle index variance if necessary
-       if (is_a<varidx>(i)) {
-               const varidx &vi = ex_to<varidx>(i);
-               if (vi.is_covariant() != covariant)
-                       i = vi.toggle_variance();
-       } else if (!covariant)
-               throw (std::runtime_error("index '" + get_symbol_name(i) + "' is not a varidx and cannot be contravariant"));
-
-       // Add index to an existing indexed object, or create a new indexed
-       // object if there are no indices yet
-       if (is_a<indexed>(base)) {
-               const ex &b = base.op(0);
-               exvector iv;
-               for (unsigned n=1; n<base.nops(); n++)
-                       iv.push_back(base.op(n));
-               iv.push_back(i);
-               return indexed(b, iv);
-       } else
-               return indexed(base, i);
-}
-
-// Get last error encountered by parser
-std::string get_parser_error(void)
-{
-       return parser_error;
-}
-
-} // namespace GiNaC
-
-// Error print routine (store error string in parser_error)
-int ginac_yyerror(char *s)
-{
-       GiNaC::parser_error = std::string(s) + " at " + std::string(ginac_yytext);
-       return 0;
-}
index 821e184d3acc0e8f014e67216f90aba91390778c..d108bd4c04be6bb315a3966c3f4a716d77a45504 100644 (file)
@@ -1418,7 +1418,7 @@ int matrix::fraction_free_elimination(const bool det)
                // Searching the first non-zero element in-place here instead of calling
                // pivot() allows us to do no more substitutions and back-substitutions
                // than are actually necessary.
-               int indx = r0;
+               unsigned indx = r0;
                while ((indx<m) &&
                       (tmp_n[indx*n+c0].subs(srl, subs_options::no_pattern).expand().is_zero()))
                        ++indx;
index 3d33c78a129eeaeaaf9074bc291e76624e403da8..63f92e1492d4fae3af2d8467f66793b0cc40bb68 100644 (file)
@@ -637,7 +637,7 @@ ex mul::eval_ncmul(const exvector & v) const
        return inherited::eval_ncmul(v);
 }
 
-bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatches, lst & repls)
+bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatches, exmap& repls)
 {      
        ex origbase;
        int origexponent;
@@ -669,7 +669,7 @@ bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatch
                patternexpsign = 1;
        }
 
-       lst saverepls = repls;
+       exmap saverepls = repls;
        if (origexponent < patternexponent || origexpsign != patternexpsign || !origbase.match(patternbase,saverepls))
                return false;
        repls = saverepls;
@@ -688,7 +688,7 @@ bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatch
   * that already have been replaced by previous substitutions and matched[i]
   * is true for factors that have been matched by the current match.
   */
-bool algebraic_match_mul_with_mul(const mul &e, const ex &pat, lst &repls,
+bool algebraic_match_mul_with_mul(const mul &e, const ex &pat, exmap& repls,
                int factor, int &nummatches, const std::vector<bool> &subsed,
                std::vector<bool> &matched)
 {
@@ -698,7 +698,7 @@ bool algebraic_match_mul_with_mul(const mul &e, const ex &pat, lst &repls,
        for (size_t i=0; i<e.nops(); ++i) {
                if(subsed[i] || matched[i])
                        continue;
-               lst newrepls = repls;
+               exmap newrepls = repls;
                int newnummatches = nummatches;
                if (tryfactsubs(e.op(i), pat.op(factor), newnummatches, newrepls)) {
                        matched[i] = true;
@@ -721,7 +721,7 @@ bool mul::has(const ex & pattern, unsigned options) const
        if(!(options&has_options::algebraic))
                return basic::has(pattern,options);
        if(is_a<mul>(pattern)) {
-               lst repls;
+               exmap repls;
                int nummatches = std::numeric_limits<int>::max();
                std::vector<bool> subsed(seq.size(), false);
                std::vector<bool> matched(seq.size(), false);
@@ -745,7 +745,7 @@ ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const
 retry1:
                        int nummatches = std::numeric_limits<int>::max();
                        std::vector<bool> currsubsed(seq.size(), false);
-                       lst repls;
+                       exmap repls;
                        
                        if(!algebraic_match_mul_with_mul(*this, it->first, repls, 0, nummatches, subsed, currsubsed))
                                continue;
@@ -754,10 +754,10 @@ retry1:
                                if (currsubsed[j])
                                        subsed[j] = true;
                        ex subsed_pattern
-                               = it->first.subs(ex(repls), subs_options::no_pattern);
+                               = it->first.subs(repls, subs_options::no_pattern);
                        divide_by *= power(subsed_pattern, nummatches);
                        ex subsed_result
-                               = it->second.subs(ex(repls), subs_options::no_pattern);
+                               = it->second.subs(repls, subs_options::no_pattern);
                        multiply_by *= power(subsed_result, nummatches);
                        goto retry1;
 
@@ -765,14 +765,14 @@ retry1:
 
                        for (size_t j=0; j<this->nops(); j++) {
                                int nummatches = std::numeric_limits<int>::max();
-                               lst repls;
+                               exmap repls;
                                if (!subsed[j] && tryfactsubs(op(j), it->first, nummatches, repls)){
                                        subsed[j] = true;
                                        ex subsed_pattern
-                                               = it->first.subs(ex(repls), subs_options::no_pattern);
+                                               = it->first.subs(repls, subs_options::no_pattern);
                                        divide_by *= power(subsed_pattern, nummatches);
                                        ex subsed_result
-                                               = it->second.subs(ex(repls), subs_options::no_pattern);
+                                               = it->second.subs(repls, subs_options::no_pattern);
                                        multiply_by *= power(subsed_result, nummatches);
                                }
                        }
index 5a7b94f8733adcebc7d1b83d7b5390ea23aef696..a5969f317d035bb48db1ce491c12491d491d3f60 100644 (file)
@@ -124,7 +124,7 @@ bool ncmul::info(unsigned inf) const
        return inherited::info(inf);
 }
 
-typedef std::vector<int> intvector;
+typedef std::vector<std::size_t> uintvector;
 
 ex ncmul::expand(unsigned options) const
 {
@@ -134,8 +134,8 @@ ex ncmul::expand(unsigned options) const
        
        // Now, look for all the factors that are sums and remember their
        // position and number of terms.
-       intvector positions_of_adds(expanded_seq.size());
-       intvector number_of_add_operands(expanded_seq.size());
+       uintvector positions_of_adds(expanded_seq.size());
+       uintvector number_of_add_operands(expanded_seq.size());
 
        size_t number_of_adds = 0;
        size_t number_of_expanded_terms = 1;
@@ -167,7 +167,7 @@ ex ncmul::expand(unsigned options) const
        exvector distrseq;
        distrseq.reserve(number_of_expanded_terms);
 
-       intvector k(number_of_adds);
+       uintvector k(number_of_adds);
 
        /* Rename indices in the static members of the product */
        exvector expanded_seq_mod;
index d8c49b8cf4b35da972b19074b506aa0380c3725e..7c3e6c84aecf5ac6d62db29452de77f981791a95 100644 (file)
@@ -35,13 +35,23 @@ int lexer::gettok()
                        return token_type::identifier;
        }
 
-       // Number: [0-9.]+
+       // Number: [0-9]+([.][0-9]*(eE[+-][0-9]+)*)*
        if (isdigit(c) || c == '.') {
                str = "";
                do {
                        str += c;
                        c = input->get();
                } while (isdigit(c) || c == '.');
+               if (c == 'E' || c == 'e') {
+                       str += 'E';
+                       c = input->get();
+                       if (isdigit(c))
+                               str += '+';
+                       do {
+                               str += c;
+                               c = input->get();
+                       } while (isdigit(c));
+               }
                return token_type::number;
        }
 
index e4956ba74d4510d234beeee023722b2c7240a2f9..1dd767778f9581cea6ad0b6e8790d400b38733c8 100644 (file)
@@ -4,23 +4,28 @@
 namespace GiNaC
 {
 
-const symbol&
+symbol
 find_or_insert_symbol(const std::string& name, symtab& syms, const bool strict)
 {
        symtab::const_iterator p = syms.find(name);
-       if (p != syms.end())
-               return p->second.first;
+       if (p != syms.end()) {
+               if (is_a<symbol>(p->second))
+                       return ex_to<symbol>(p->second);
+               else
+                       throw std::invalid_argument(
+                               std::string("find_or_insert_symbol: name \"")
+                               + name + "\" does not correspond to a symbol");
+       }
+
 
        if (strict)
                throw std::invalid_argument(
                                std::string("find_or_insert_symbol: symbol \"") 
                                + name + "\" not found");
 
-       // false means this symbol was created by parser 
-       const std::pair<symbol, bool> tmp = std::make_pair(symbol(name), false);
-
-       symtab::iterator i = syms.insert(symtab::value_type(name, tmp)).first;
-       return i->second.first;
+       const symbol sy(name);
+       syms[name] = sy;
+       return sy;
 }
 
 }
index d0a3d55870839b627a981ef113f97a51e82291a7..efcb6c620f309f0413f9f4e1279164eb4eda1c09 100644 (file)
@@ -11,14 +11,11 @@ namespace GiNaC
 {
 
 /**
- * Establishes correspondence between the strings and symbols.
+ * Establishes correspondence between the strings and expressions.
  * The parser will create missing symbols (if not instructed otherwise,
  * in which case it fails if the expression contains unknown symbols).
- * The .second element of pair helps to distinguish between the user
- * supplied symbols and parser generated ones. The .first is the symbol
- * itself
  */
-typedef std::map<std::string, std::pair<symbol, bool> > symtab;
+typedef std::map<std::string, ex> symtab;
 
 /**
  * Find the symbol with the @a name in the symbol table @a syms.
@@ -26,7 +23,7 @@ typedef std::map<std::string, std::pair<symbol, bool> > symtab;
  * If symbol is missing and @a strict = false, insert it, otherwise
  * throw an exception.
  */
-extern const symbol&
+extern symbol
 find_or_insert_symbol(const std::string& name, symtab& syms,
                      const bool strict);
 
index 1e905cf5ff1133e4d2074cf05dd5992a9d813467..6411ba1416e740a85705ceee3153815b33e9364f 100644 (file)
@@ -60,21 +60,25 @@ ex parser::parse_paren_expr()
 }
 
 extern numeric* _num_1_p;
+extern ex _ex0;
 
 /// unary_expr: [+-] expression
-ex parser::parse_unary_expr(const int s)
+ex parser::parse_unary_expr()
 {
-       // consume '-' (or '+')
-       get_next_tok();
-       ex v = parse_expression();
-       switch (s) {
-               case '-':
-                       return (new mul(v, *_num_1_p))->setflag(status_flags::dynallocated);
-               case '+':
-                       return v;
-               default:
-                       Parse_error_("invalid unary operator \"" << char(s) << "\"");
-       }
+       // Unlike most other parse_* method this one does NOT consume
+       // current token so parse_binop_rhs() knows what kind of operator
+       // is being parsed.
+       
+       // There are different kinds of expressions which need to be handled:
+       // -a+b 
+       // -(a) 
+       // +a
+       // +(a)
+       // Delegete the work to parse_binop_rhs(), otherwise we end up
+       // duplicating it here. 
+       ex lhs = _ex0; // silly trick
+       ex e = parse_binop_rhs(0, lhs);
+       return e;
 }
 
 /// primary: identifier_expr | number_expr | paren_expr | unary_expr 
@@ -88,9 +92,8 @@ ex parser::parse_primary()
                case '(': 
                         return parse_paren_expr();
                case '-':
-                        return parse_unary_expr('-');
                case '+':
-                        return parse_unary_expr('+');
+                        return parse_unary_expr();
                case lexer::token_type::literal:
                         return parse_literal_expr();
                case lexer::token_type::eof:
@@ -118,6 +121,7 @@ ex parser::parse_number_expr()
 /// literal_expr: 'I' | 'Pi' | 'Euler' | 'Catalan'
 ex parser::parse_literal_expr()
 {
+       get_next_tok(); // consume the literal
        if (scanner->str == "I")
                return I;
        else if (scanner->str == "Pi")
@@ -134,6 +138,13 @@ ex parser::operator()(std::istream& input)
        scanner->switch_input(&input);
        get_next_tok();
        ex ret = parse_expression();
+       // parse_expression() stops if it encounters an unknown token.
+       // This is not a bug: since the parser is recursive checking
+       // whether the next token is valid is responsibility of the caller.
+       // Hence make sure nothing is left in the stream:
+       if (token != lexer::token_type::eof)
+               Parse_error("expected EOF");
+
        return ret;
 }
 
@@ -150,9 +161,9 @@ int parser::get_next_tok()
        return token;
 }
 
-parser::parser(const symtab& syms_, const prototype_table& funcs_,
-              const bool strict_) : strict(strict_), funcs(funcs_),
-       syms(syms_)
+parser::parser(const symtab& syms_, const bool strict_,
+              const prototype_table& funcs_) : strict(strict_),
+       funcs(funcs_), syms(syms_)
 {
        scanner = new lexer();
 }
index 3e10b8fc016bacc9a685918c969a89d36c451841..b46ec6509ba2462509da836cd858a06681abc239 100644 (file)
@@ -1,4 +1,5 @@
 #ifndef GINAC_PARSER_HPP_
+#define GINAC_PARSER_HPP_
 
 #include "parse_context.hpp"
 #include <stdexcept>
@@ -48,7 +49,7 @@ class parser
        ex parse_number_expr();
 
        /// unary_expr: [+-] expression
-       ex parse_unary_expr(const int c);
+       ex parse_unary_expr();
 
        /// literal_expr: 'I' | 'Pi' | 'Euler' | 'Catalan'
        ex parse_literal_expr();
@@ -61,8 +62,8 @@ public:
         *        symbol is encountered.
         */
        parser(const symtab& syms_ = symtab(),
-               const prototype_table& funcs_ = get_default_reader(),
-               const bool strict_ = false);
+               const bool strict_ = false,
+               const prototype_table& funcs_ = get_default_reader());
        ~parser();
 
        /// parse the stream @a input
@@ -75,10 +76,15 @@ public:
        { 
                return syms; 
        }
+       /// read/write access to the symbol table
+       symtab& get_syms()
+       {
+               return syms;
+       }
 
-private:
        /// If true, throw an exception if an unknown symbol is encountered.
-       const bool strict;
+       bool strict;
+private:
        /**
         * Function/ctor table, maps a prototype (which is a name and number
         * arguments) to a C++ function. Used for parsing identifier_expr's
diff --git a/ginac/parser/parser_compat.cpp b/ginac/parser/parser_compat.cpp
new file mode 100644 (file)
index 0000000..d0723ee
--- /dev/null
@@ -0,0 +1,49 @@
+/// @file parser_compat.cpp Parser interface compatible with the old
+///       (bison/flex based) parser.
+#include "ex.h"
+#include "idx.h"
+#include "lst.h"
+#include "parser.hpp"
+#include <string>
+#include <iostream>
+
+namespace GiNaC
+{
+static symtab make_symtab(const ex& l);
+
+ptr<basic> ex::construct_from_string_and_lst(const std::string &s, const ex &l)
+{
+       static const bool strict = true;
+       symtab syms = make_symtab(l);
+       parser reader(syms, strict); 
+       ex parsed_ex = reader(s);
+       return parsed_ex.bp;
+}
+
+static std::string get_symbol_name(const ex & s);
+
+static symtab make_symtab(const ex& l)
+{
+       symtab syms;
+       if (is_exactly_a<lst>(l)) {
+               for (std::size_t i = 0; i < l.nops(); i++) {
+                       const ex &o = l.op(i);
+                       if (is_a<symbol>(o) || (is_a<idx>(o) && is_a<symbol>(o.op(0))))
+                               syms[get_symbol_name(o)] = o;
+               }
+       }
+       return syms;
+}
+
+static std::string get_symbol_name(const ex & s)
+{
+       if (is_a<symbol>(s))
+               return ex_to<symbol>(s).get_name();
+       else if (is_a<idx>(s) && is_a<symbol>(s.op(0)))
+               return ex_to<symbol>(s.op(0)).get_name();
+       else
+               throw (std::invalid_argument("get_symbol_name(): unexpected expression type"));
+}
+
+}
+
diff --git a/ginac/polynomial/cra_garner.cpp b/ginac/polynomial/cra_garner.cpp
new file mode 100644 (file)
index 0000000..b400adb
--- /dev/null
@@ -0,0 +1,88 @@
+#include <cln/integer.h>
+#include <cln/modinteger.h>
+#include <vector>
+#include <cstddef>
+#include "cra_garner.hpp"
+
+namespace cln
+{
+using std::vector;
+using std::size_t;
+
+static cl_I 
+retract_symm(const cl_MI& x, const cl_modint_ring& R,
+            const cl_I& modulus)
+{
+       cl_I result = R->retract(x);
+       if (result > (modulus >> 1))
+               result = result - modulus;
+       return result;
+}
+
+static void
+compute_recips(vector<cl_MI>& dst, 
+              const vector<cl_I>& moduli)
+{
+       for (size_t k = 1; k < moduli.size(); ++k) {
+               cl_modint_ring R = find_modint_ring(moduli[k]);
+               cl_MI product = R->canonhom(moduli[0]);
+               for (size_t i = 1; i < k; ++i)
+                       product = product*moduli[i];
+               dst[k-1] = recip(product);
+       }
+}
+
+static void
+compute_mix_radix_coeffs(vector<cl_I>& dst,
+                        const vector<cl_I>& residues,
+                        const vector<cl_I>& moduli,
+                        const vector<cl_MI>& recips)
+{
+       dst[0] = residues[0];
+
+       do {
+               cl_modint_ring R = find_modint_ring(moduli[1]);
+               cl_MI tmp = R->canonhom(residues[0]);
+               cl_MI next = (R->canonhom(residues[1]) - tmp)*recips[0];
+               dst[1] = retract_symm(next, R, moduli[1]);
+       } while (0);
+
+       for (size_t k = 2; k < residues.size(); ++k) {
+               cl_modint_ring R = find_modint_ring(moduli[k]);
+               cl_MI tmp = R->canonhom(dst[k-1]);
+
+               for (size_t j = k - 1 /* NOT k - 2 */; j-- != 0; )
+                       tmp = tmp*moduli[j] + R->canonhom(dst[j]);
+
+               cl_MI next = (R->canonhom(residues[k]) - tmp)*recips[k-1];
+               dst[k] = retract_symm(next, R, moduli[k]);
+       }
+}
+
+static cl_I
+mixed_radix_2_ordinary(const vector<cl_I>& mixed_radix_coeffs,
+                      const vector<cl_I>& moduli)
+{
+       size_t k = mixed_radix_coeffs.size() - 1;
+       cl_I u = mixed_radix_coeffs[k];
+       for (; k-- != 0; )
+               u = u*moduli[k] + mixed_radix_coeffs[k];
+       return u;
+}
+
+cl_I integer_cra(const vector<cl_I>& residues,
+                const vector<cl_I>& moduli)
+{
+
+       vector<cl_MI> recips(moduli.size() - 1);
+       compute_recips(recips, moduli);
+
+       vector<cl_I> coeffs(moduli.size());
+       compute_mix_radix_coeffs(coeffs, residues, moduli, recips);
+       cl_I result = mixed_radix_2_ordinary(coeffs, moduli);
+
+       return result;
+}
+
+} // namespace cln
+
diff --git a/ginac/polynomial/cra_garner.hpp b/ginac/polynomial/cra_garner.hpp
new file mode 100644 (file)
index 0000000..eff8e9b
--- /dev/null
@@ -0,0 +1,12 @@
+#ifndef CL_INTEGER_CRA
+#define CL_INTEGER_CRA
+#include <vector>
+#include <cln/integer.h>
+
+namespace cln
+{
+extern cl_I integer_cra(const std::vector<cl_I>& residues,
+                       const std::vector<cl_I>& moduli);
+}
+#endif // CL_INTEGER_CRA
+
diff --git a/ginac/polynomial/debug.hpp b/ginac/polynomial/debug.hpp
new file mode 100644 (file)
index 0000000..901ef87
--- /dev/null
@@ -0,0 +1,29 @@
+#ifndef GINAC_MOD_GCD_DEBUG_HPP
+#define GINAC_MOD_GCD_DEBUG_HPP
+#include <iostream>
+#include <string>
+#include <stdexcept>
+#include <sstream>
+
+#define DEBUG_PREFIX __func__ << ':' << __LINE__ << ": "
+#define EXCEPTION_PREFIX std::string(__func__) + std::string(": ") +
+
+#define Dout2(stream, msg) \
+do {                                                                   \
+       stream << DEBUG_PREFIX << msg << std::endl << std::flush;       \
+} while (0)
+#define Dout(msg) Dout2(std::cout, msg)
+
+#define bug3_on(condition, the_exception, msg)                         \
+do {                                                                   \
+       if (unlikely(condition)) {                                      \
+               std::ostringstream err_stream;                          \
+               Dout2(err_stream, "BUG: " << msg);                      \
+               throw the_exception(err_stream.str());                  \
+       }                                                               \
+} while (0)
+
+#define bug_on(condition, msg) bug3_on(condition, std::logic_error, msg)
+
+#endif // GINAC_MOD_GCD_DEBUG_HPP
+
diff --git a/ginac/polynomial/gcd_euclid.tcc b/ginac/polynomial/gcd_euclid.tcc
new file mode 100644 (file)
index 0000000..deae252
--- /dev/null
@@ -0,0 +1,45 @@
+#ifndef GINAC_GCD_EUCLID
+#define GINAC_GCD_EUCLID
+#include "upoly.hpp"
+#include "remainder.tcc"
+#include "normalize.tcc"
+#include "debug.hpp"
+#include "upoly_io.hpp"
+
+namespace GiNaC
+{
+
+static bool
+gcd_euclid(umodpoly& c, umodpoly /* passed by value */ a, umodpoly b)
+{
+       if (a.size() == 0) {
+               c.clear();
+               return true;
+       }
+       if (b.size() == 0) {
+               c.clear();
+               return true;
+       }
+       bug_on(a[0].ring()->modulus != b[0].ring()->modulus,
+               "different moduli");
+
+       normalize_in_field(a);
+       normalize_in_field(b);
+       if (degree(a) < degree(b))
+               std::swap(a, b);
+
+       umodpoly r;
+       while (b.size() != 0) {
+               remainder_in_field(r, a, b); 
+               a = b;
+               b = r;
+       }
+       normalize_in_field(a);
+       c = a;
+       return false;
+}
+
+} // namespace GiNaC
+
+#endif // GINAC_GCD_EUCLID
+
diff --git a/ginac/polynomial/mod_gcd.cpp b/ginac/polynomial/mod_gcd.cpp
new file mode 100644 (file)
index 0000000..8aa8480
--- /dev/null
@@ -0,0 +1,165 @@
+#include "upoly.hpp"
+#include "gcd_euclid.tcc"
+#include "cra_garner.hpp"
+#include <cln/random.h>
+#include <cln/numtheory.h>
+#include <stdexcept>
+#include "debug.hpp"
+
+namespace GiNaC
+{
+/**
+ * @brief Remove the integer content from univariate polynomials A and B.
+ *
+ * As a byproduct compute the GCD of contents.
+ */
+static void remove_content(upoly& A, upoly& B, upoly::value_type& c)
+{
+       // remove the integer
+       upoly::value_type acont, bcont;
+       normalize_in_ring(A, &acont);
+       normalize_in_ring(B, &bcont);
+       c = gcd(acont, bcont);
+};
+
+/// Check if @a candidate divides both @a A and @a B
+static bool
+do_division_check(const upoly& A, const upoly& B, const upoly& candidate)
+{
+       upoly r1;
+       remainder_in_ring(r1, A, candidate);
+       if (r1.size() != 0)
+               return false;
+
+       upoly r2;
+       remainder_in_ring(r2, B, candidate);
+       if (r2.size() != 0)
+               return false;
+
+       return true;
+}
+
+/**
+ * Given two GCD candidates H \in Z/q[x], and C \in Z/p[x] (where p is a prime)
+ * compute the candidate in Z/(q*p) with CRA (chinise remainder algorithm)
+ *
+ * @param H \in Z/q[x] GCD candidate, will be updated by this function
+ * @param q modulus of H, will NOT be updated by this function
+ * @param C \in Z/p[x] GCD candidate, will be updated by this function
+ * @param p modulus of C
+ */
+static void
+update_the_candidate(upoly& H, const upoly::value_type& q,
+                    const umodpoly& C,
+                    const upoly::value_type& p,
+                    const cln::cl_modint_ring& R)
+{
+       typedef upoly::value_type ring_t;
+       std::vector<ring_t> moduli(2);
+       moduli[0] = q;
+       moduli[1] = p;
+       if (H.size() < C.size())
+               H.resize(C.size());
+
+       for (std::size_t  i = C.size(); i-- != 0; ) {
+               std::vector<ring_t> coeffs(2);
+               coeffs[0] = H[i];
+               coeffs[1] = R->retract(C[i]); 
+               H[i] = integer_cra(coeffs, moduli);
+       }
+}
+
+/// Convert Z/p[x] -> Z[x]
+static void retract(upoly& p, const umodpoly& cp,
+                   const cln::cl_modint_ring& Rp)
+{
+       p.resize(cp.size());
+       for (std::size_t i = cp.size(); i-- != 0; )
+               p[i] = Rp->retract(cp[i]);
+}
+
+
+/// Find the prime which is > p, and does NOT divide g
+static void find_next_prime(cln::cl_I& p, const cln::cl_I& g)
+{
+       do {
+               ++p;
+               p = nextprobprime(p);
+       } while (zerop(mod(g, p)));
+}
+
+/// Compute the GCD of univariate polynomials A, B \in Z[x]
+void mod_gcd(upoly& result, upoly A, upoly B)
+{
+       typedef upoly::value_type ring_t;
+       ring_t content_gcd;
+       // remove the integer content
+       remove_content(A, B, content_gcd);
+
+       // compute the coefficient bound for GCD(a, b)
+       ring_t g = gcd(lcoeff(A), lcoeff(B));
+       std::size_t max_gcd_degree = std::min(degree(A), degree(B));
+       ring_t limit = (ring_t(1) << max_gcd_degree)*g*
+                      std::min(max_coeff(A), max_coeff(B));
+       ring_t q(0);
+       upoly H;
+
+       ring_t p(2);
+       while (true) {
+               find_next_prime(p, g);
+
+               // Map the polynomials onto Z/p[x]
+               cln::cl_modint_ring Rp = cln::find_modint_ring(p);
+               cln::cl_MI gp = Rp->canonhom(g);
+               umodpoly ap(A.size()), bp(B.size());
+               make_umodpoly(ap, A, Rp);
+               make_umodpoly(bp, B, Rp);
+
+               // Compute the GCD in Z/p[x]
+               umodpoly cp;
+               gcd_euclid(cp, ap, bp);
+               bug_on(cp.size() == 0, "gcd(ap, bp) = 0, with ap = " <<
+                                       ap << ", and bp = " << bp);
+
+
+               // Normalize the candidate so that its leading coefficient
+               // is g mod p
+               umodpoly::value_type norm_factor = gp*recip(lcoeff(cp));
+               bug_on(zerop(norm_factor), "division in a field give 0");
+
+               lcoeff(cp) = gp;
+               for (std::size_t k = cp.size() - 1; k-- != 0; )
+                       cp[k] = cp[k]*norm_factor;
+
+
+               // check for unlucky homomorphisms
+               if (degree(cp) < max_gcd_degree) {
+                       q = p;
+                       max_gcd_degree = degree(cp);
+                       retract(H, cp, Rp);
+               } else {
+                       update_the_candidate(H, q, cp, p, Rp);
+                       q = q*p;
+               }
+
+               if (q > limit) {
+                       result = H;
+                       normalize_in_ring(result);
+                       // if H divides both A and B it's a GCD
+                       if (do_division_check(A, B, result)) {
+                               result *= content_gcd;
+                               break; 
+                       }
+                       // H does not divide A and/or B, look for
+                       // another one
+               } else if (degree(cp) == 0) {
+                       // Polynomials are relatively prime
+                       result.resize(1);
+                       result[0] = content_gcd;
+                       break;
+               }
+       }
+}
+
+} // namespace GiNaC
+
diff --git a/ginac/polynomial/mod_gcd.hpp b/ginac/polynomial/mod_gcd.hpp
new file mode 100644 (file)
index 0000000..43e1e92
--- /dev/null
@@ -0,0 +1,11 @@
+#ifndef GINAC_MOD_GCD_HPP
+#define GINAC_MOD_GCD_HPP
+#include "upoly.hpp"
+
+namespace GiNaC
+{
+extern void mod_gcd(upoly& result, upoly A, upoly B);
+}
+
+#endif // GINAC_MOD_GCD_HPP
+
diff --git a/ginac/polynomial/normalize.tcc b/ginac/polynomial/normalize.tcc
new file mode 100644 (file)
index 0000000..e246307
--- /dev/null
@@ -0,0 +1,93 @@
+#ifndef GINAC_UPOLY_NORMALIZE_TCC
+#define GINAC_UPOLY_NORMALIZE_TCC
+#include "upoly.hpp"
+#include "ring_traits.hpp"
+#include "debug.hpp"
+
+namespace GiNaC
+{
+
+/// Make the univariate polynomial @a a \in F[x] unit normal.
+/// F should be a field.
+/// Returns true if the polynomial @x is already unit normal, and false
+/// otherwise.
+static bool normalize_in_field(umodpoly& a, cln::cl_MI* content_ = 0)
+{
+       if (a.size() == 0)
+               return true;
+       if (lcoeff(a) == the_one(a[0])) {
+               if (content_)
+                       *content_ = the_one(a[0]);
+               return true;
+       }
+
+       const cln::cl_MI lc_1 = recip(lcoeff(a));
+       for (std::size_t k = a.size(); k-- != 0; )
+               a[k] = a[k]*lc_1;
+       if (content_)
+               *content_ = lc_1;
+       return false;
+}
+
+/// Make the univariate polynomial @a x unit normal. This version is used
+/// for rings which are not fields. 
+/// Returns true if the polynomial @x is already unit normal, and false
+/// otherwise.
+template<typename T> bool
+normalize_in_ring(T& x, typename T::value_type* content_ = 0, int* unit_ = 0)
+{
+       typedef typename T::value_type ring_t;
+       static const ring_t one(1);
+       if (x.empty())
+               return true;
+
+       bool something_changed = false;
+       if (minusp(lcoeff(x))) {
+               something_changed = true;
+               if (unit_)
+                       *unit_ = -1;
+               for (std::size_t i = x.size(); i-- != 0; )
+                       x[i] = - x[i];
+       }
+
+       if (degree(x) == 0) {
+               if (content_)
+                       *content_ = x[0];
+               if (x[0] == one)
+                       return something_changed;
+               x[0] = one;
+               return false; // initial polynomial was unit normal
+       }
+               
+       // Compute the gcd of coefficients
+       ring_t content = lcoeff(x);
+       // We want this function to be fast when applied to unit normal
+       // polynomials. Hence we start from the leading coefficient.
+       for (std::size_t i = x.size() - 1; i-- != 0; ) {
+               if (content == one) {
+                       if (content_)
+                               *content_ = one;
+                       return something_changed;
+               }
+               content = gcd(x[i], content);
+       }
+
+       if (content == one) {
+               if (content_)
+                       *content_ = one;
+               return something_changed;
+       }
+       
+       lcoeff(x) = one;
+       for (std::size_t i = x.size() - 1; i-- != 0; )
+               x[i] = exquo(x[i], content);
+
+       if (content_)
+               *content_ = content;
+       return false; // initial polynomial was not unit normal
+}
+
+} // namespace GiNaC
+       
+#endif // GINAC_UPOLY_NORMALIZE_TCC
+
diff --git a/ginac/polynomial/remainder.tcc b/ginac/polynomial/remainder.tcc
new file mode 100644 (file)
index 0000000..4db2843
--- /dev/null
@@ -0,0 +1,116 @@
+#ifndef GINAC_UPOLY_REMAINDER_TCC
+#define GINAC_UPOLY_REMAINDER_TCC
+#include "upoly.hpp"
+#include "ring_traits.hpp"
+#include "upoly_io.hpp"
+#include "debug.hpp"
+
+namespace GiNaC
+{
+/**
+ * @brief Polynomial remainder for univariate polynomials over fields
+ *
+ * Given two univariate polynomials \f$a, b \in F[x]\f$, where F is
+ * a finite field (presumably Z/p) computes the remainder @a r, which is
+ * defined as \f$a = b q + r\f$. Returns true if the remainder is zero
+ * and false otherwise.
+ */
+static bool
+remainder_in_field(umodpoly& r, const umodpoly& a, const umodpoly& b)
+{
+       typedef cln::cl_MI field_t;
+
+       if (degree(a) < degree(b)) {
+               r = a;
+               return false;
+       }
+       // The coefficient ring is a field, so any 0 degree polynomial
+       // divides any other polynomial.
+       if (degree(b) == 0) {
+               r.clear();
+               return true;
+       }
+
+       r = a;
+       const field_t b_lcoeff = lcoeff(b);
+       for (std::size_t k = a.size(); k-- >= b.size(); ) {
+
+               // r -= r_k/b_n x^{k - n} b(x)
+               if (zerop(r[k]))
+                       continue;
+
+               field_t qk = div(r[k], b_lcoeff);
+               bug_on(zerop(qk), "division in a field yield zero: "
+                                  << r[k] << '/' << b_lcoeff);
+
+               // Why C++ is so off-by-one prone?
+               for (std::size_t j = k, i = b.size(); i-- != 0; --j) {
+                       if (zerop(b[i]))
+                               continue;
+                       r[j] = r[j] - qk*b[i];
+               }
+               bug_on(!zerop(r[k]), "polynomial division in field failed: " <<
+                                     "r[" << k << "] = " << r[k] << ", " <<
+                                     "r = " << r << ", b = " << b);
+
+       }
+
+       // Canonicalize the remainder: remove leading zeros. Give a hint
+       // to canonicalize(): we know degree(remainder) < degree(b) 
+       // (because the coefficient ring is a field), so 
+       // c_{degree(b)} \ldots c_{degree(a)} are definitely zero.
+       std::size_t from = degree(b) - 1;
+       canonicalize(r, from);
+       return r.empty();
+}
+
+/**
+ * @brief Polynomial remainder for univariate polynomials over a ring. 
+ *
+ * Given two univariate polynomials \f$a, b \in R[x]\f$, where R is
+ * a ring (presumably Z) computes the remainder @a r, which is
+ * defined as \f$a = b q + r\f$. Returns true if the remainder is zero
+ * and false otherwise.
+ */
+template<typename T>
+bool remainder_in_ring(T& r, const T& a, const T& b)
+{
+       typedef typename T::value_type ring_t;
+       if (degree(a) < degree(b)) {
+               r = a;
+               return false;
+       }
+       // N.B: don't bother to optimize division by constant
+
+       r = a;
+       const ring_t b_lcoeff = lcoeff(b);
+       for (std::size_t k = a.size(); k-- >= b.size(); ) {
+
+               // r -= r_k/b_n x^{k - n} b(x)
+               if (zerop(r[k]))
+                       continue;
+
+               const ring_t qk = div(r[k], b_lcoeff);
+
+               // Why C++ is so off-by-one prone?
+               for (std::size_t j = k, i = b.size(); i-- != 0; --j) {
+                       if (zerop(b[i]))
+                               continue;
+                       r[j] = r[j] - qk*b[i];
+               }
+
+               if (!zerop(r[k])) {
+                       // division failed, don't bother to continue
+                       break;
+               }
+       }
+
+       // Canonicalize the remainder: remove leading zeros. We can't say
+       // anything about the degree of the remainder here.
+       canonicalize(r);
+       return r.empty();
+}
+} // namespace GiNaC
+
+#endif // GINAC_UPOLY_REMAINDER_TCC
+
diff --git a/ginac/polynomial/ring_traits.hpp b/ginac/polynomial/ring_traits.hpp
new file mode 100644 (file)
index 0000000..10fb172
--- /dev/null
@@ -0,0 +1,32 @@
+#ifndef GINAC_RING_TRAITS_HPP
+#define GINAC_RING_TRAITS_HPP
+#include <cln/integer.h>
+#include <cln/modinteger.h>
+
+namespace cln
+{
+static inline cln::cl_I div(const cln::cl_I& x, const cln::cl_I& y)
+{
+       return cln::exquo(x, y);
+}
+
+static inline cln::cl_I get_ring_elt(const cln::cl_I& sample, const int val)
+{
+       return cln::cl_I(val);
+}
+
+static inline cln::cl_MI get_ring_elt(const cln::cl_MI& sample, const int val)
+{
+       return sample.ring()->canonhom(val);
+}
+
+template<typename T>
+static inline T the_one(const T& sample)
+{
+       return get_ring_elt(sample, 1);
+}
+
+} // namespace cln
+
+#endif // GINAC_RING_TRAITS_HPP
+
diff --git a/ginac/polynomial/upoly.hpp b/ginac/polynomial/upoly.hpp
new file mode 100644 (file)
index 0000000..3564ae0
--- /dev/null
@@ -0,0 +1,129 @@
+#ifndef GINAC_NEW_UPOLY_HPP
+#define GINAC_NEW_UPOLY_HPP
+#include <vector>
+#include <cstddef>
+#include <cln/integer.h>
+#include <cln/modinteger.h>
+#include <cln/integer_io.h>
+#include <stdexcept>
+#include <iterator>
+#include <limits>
+#include "debug.hpp"
+#include "compiler.h"
+
+namespace GiNaC
+{
+
+typedef std::vector<cln::cl_I> upoly;
+typedef std::vector<cln::cl_MI> umodpoly;
+
+template<typename T> static std::size_t degree(const T& p)
+{
+       return p.size() - 1;
+}
+
+template<typename T> static typename T::value_type lcoeff(const T& p)
+{
+       bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
+       return p[p.size() - 1];
+}
+
+template<typename T> static typename T::reference lcoeff(T& p)
+{
+       bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
+       return p[p.size() - 1];
+}
+
+template<typename T> static typename T::value_type max_coeff(const T& p)
+{
+       bug_on(p.empty(), "max_coeff of a zero polynomial is undefined");
+       typename T::value_type curr = p[p.size() - 1];
+       for (std::size_t i = p.size(); i-- != 0; ) {
+               if (p[i] > curr)
+                       curr = p[i];
+       }
+       return curr;
+}
+
+
+
+// Canonicalize the polynomial, i.e. remove leading zero coefficients
+template<typename T> static void
+canonicalize(T& p, const typename T::size_type hint =
+                  std::numeric_limits<typename T::size_type>::max())
+{
+       if (p.empty())
+               return;
+
+       std::size_t i = p.size() - 1;
+       // Be fast if the polynomial is already canonicalized
+       if (!zerop(p[i]))
+               return;
+
+       if (hint < p.size())
+               i = hint;
+
+       bool is_zero = false;
+       do {
+               if (!zerop(p[i])) {
+                       ++i;
+                       break;
+               }
+               if (i == 0) {
+                       is_zero = true;
+                       break;
+               }
+               --i;
+       } while (true);
+
+       if (is_zero) {
+               p.clear();
+               return;
+       }
+
+       bug_on(!zerop(p.at(i)), "p[" << i << "] = " << p[i] << " != 0 would be erased.");
+
+       typename T::const_iterator it = p.begin() + i;
+       for (std::size_t k = i; it != p.end(); ++it, ++k) {
+               bug_on(!zerop(*it), "p[" << k << "] = " << p[k] << " != 0 "
+                                  "would be erased.");
+       }
+
+       p.erase(p.begin() + i, p.end());
+
+       bug_on(!p.empty() && zerop(lcoeff(p)), "oops, lcoeff(p) = 0");
+}
+
+// Multiply a univariate polynomial p \in R[x] with a constant c \in R
+template<typename T> static T&
+operator*=(T& p, const typename T::value_type& c)
+{
+       if (p.empty())
+               return p;
+       if (zerop(c)) {
+               p.clear();
+               return p;
+       }
+       if (c == the_one(p[0]))
+               return p;
+
+       for (std::size_t i = p.size(); i-- != 0; )
+               p[i] = p[i]*c;
+       canonicalize(p);
+       return p;
+}
+
+// Convert Z[x] -> Z/p[x]
+
+static void
+make_umodpoly(umodpoly& up, const upoly& p, const cln::cl_modint_ring& R)
+{
+       for (std::size_t i = p.size(); i-- != 0; )
+               up[i] = R->canonhom(p[i]);
+       canonicalize(up);
+}
+
+} // namespace GiNaC
+
+#endif // GINAC_NEW_UPOLY_HPP
+
diff --git a/ginac/polynomial/upoly_io.cpp b/ginac/polynomial/upoly_io.cpp
new file mode 100644 (file)
index 0000000..35b623a
--- /dev/null
@@ -0,0 +1,53 @@
+#include <iostream>
+#include <string>
+#include "upoly.hpp"
+#include "upoly_io.hpp"
+#include <cln/integer_io.h>
+
+namespace GiNaC
+{
+using std::ostream;
+using std::string;
+
+template<typename T> static void 
+print(const T& p, ostream& os, const string& varname = string("x"))
+{
+       if (p.size() == 0)
+               os << '0';
+
+       bool seen_nonzero = false;
+
+       for (std::size_t i = p.size(); i-- != 0;  ) {
+               if (zerop(p[i])) {
+                       if (seen_nonzero)
+                               continue;
+                       os << "+ [WARNING: 0]*" << varname << "^" << i << "]";
+                       continue;
+               }
+               seen_nonzero = true;
+               os << "+ (" << p[i] << ")";
+               if (i != 0)
+                       os << "*" << varname;
+               if (i > 1)
+                       os << '^' << i;
+               os << " ";
+       }
+}
+
+#define DEFINE_OPERATOR_OUT(type)                                      \
+std::ostream& operator<<(std::ostream& os, const type& p)              \
+{                                                                      \
+       print(p, os);                                                   \
+       return os;                                                      \
+}                                                                      \
+void dbgprint(const type& p)                                           \
+{                                                                      \
+       print(p, std::cerr);                                            \
+}
+
+DEFINE_OPERATOR_OUT(upoly);
+DEFINE_OPERATOR_OUT(umodpoly);
+#undef DEFINE_OPERATOR_OUT
+
+} // namespace GiNaC
+
diff --git a/ginac/polynomial/upoly_io.hpp b/ginac/polynomial/upoly_io.hpp
new file mode 100644 (file)
index 0000000..a9fb259
--- /dev/null
@@ -0,0 +1,12 @@
+#ifndef GINAC_UPOLY_IO_HPP
+#define GINAC_UPOLY_IO_HPP
+#include <iostream>
+#include "upoly.hpp"
+
+namespace GiNaC
+{
+extern std::ostream& operator<<(std::ostream&, const upoly& );
+extern std::ostream& operator<<(std::ostream&, const umodpoly& );
+}
+
+#endif // GINAC_UPOLY_IO_HPP
index ec339b9a660afc6aac21fa9307b6bec5e26144ce..a4a33c576ca415e0a3baa7e6b089fd8dae22b81f 100644 (file)
@@ -635,7 +635,7 @@ bool power::has(const ex & other, unsigned options) const
 }
 
 // from mul.cpp
-extern bool tryfactsubs(const ex &, const ex &, int &, lst &);
+extern bool tryfactsubs(const ex &, const ex &, int &, exmap&);
 
 ex power::subs(const exmap & m, unsigned options) const
 {      
@@ -651,9 +651,13 @@ ex power::subs(const exmap & m, unsigned options) const
 
        for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
                int nummatches = std::numeric_limits<int>::max();
-               lst repls;
-               if (tryfactsubs(*this, it->first, nummatches, repls))
-                       return (ex_to<basic>((*this) * power(it->second.subs(ex(repls), subs_options::no_pattern) / it->first.subs(ex(repls), subs_options::no_pattern), nummatches))).subs_one_level(m, options);
+               exmap repls;
+               if (tryfactsubs(*this, it->first, nummatches, repls)) {
+                       ex anum = it->second.subs(repls, subs_options::no_pattern);
+                       ex aden = it->first.subs(repls, subs_options::no_pattern);
+                       ex result = (*this)*power(anum/aden, nummatches);
+                       return (ex_to<basic>(result)).subs_one_level(m, options);
+               }
        }
 
        return subs_one_level(m, options);
@@ -864,7 +868,6 @@ ex power::expand_add(const add & a, int n, unsigned options) const
        intvector k(m-1);
        intvector k_cum(m-1); // k_cum[l]:=sum(i=0,l,k[l]);
        intvector upper_limit(m-1);
-       int l;
 
        for (size_t l=0; l<m-1; ++l) {
                k[l] = 0;
@@ -875,7 +878,7 @@ ex power::expand_add(const add & a, int n, unsigned options) const
        while (true) {
                exvector term;
                term.reserve(m+1);
-               for (l=0; l<m-1; ++l) {
+               for (std::size_t l = 0; l < m - 1; ++l) {
                        const ex & b = a.op(l);
                        GINAC_ASSERT(!is_exactly_a<add>(b));
                        GINAC_ASSERT(!is_exactly_a<power>(b) ||
@@ -890,7 +893,7 @@ ex power::expand_add(const add & a, int n, unsigned options) const
                                term.push_back(power(b,k[l]));
                }
 
-               const ex & b = a.op(l);
+               const ex & b = a.op(m - 1);
                GINAC_ASSERT(!is_exactly_a<add>(b));
                GINAC_ASSERT(!is_exactly_a<power>(b) ||
                             !is_exactly_a<numeric>(ex_to<power>(b).exponent) ||
@@ -904,7 +907,7 @@ ex power::expand_add(const add & a, int n, unsigned options) const
                        term.push_back(power(b,n-k_cum[m-2]));
 
                numeric f = binomial(numeric(n),numeric(k[0]));
-               for (l=1; l<m-1; ++l)
+               for (std::size_t l = 1; l < m - 1; ++l)
                        f *= binomial(numeric(n-k_cum[l-1]),numeric(k[l]));
 
                term.push_back(f);
@@ -912,12 +915,19 @@ ex power::expand_add(const add & a, int n, unsigned options) const
                result.push_back(ex((new mul(term))->setflag(status_flags::dynallocated)).expand(options));
 
                // increment k[]
-               l = m-2;
-               while ((l>=0) && ((++k[l])>upper_limit[l])) {
+               bool done = false;
+               std::size_t l = m - 2;
+               while ((++k[l]) > upper_limit[l]) {
                        k[l] = 0;
-                       --l;
+                       if (l != 0)
+                               --l;
+                       else {
+                               done = true;
+                               break;
+                       }
                }
-               if (l<0) break;
+               if (done)
+                       break;
 
                // recalc k_cum[] and upper_limit[]
                k_cum[l] = (l==0 ? k[0] : k_cum[l-1]+k[l]);
index 26dff67beefd88ed12e38ba4a03b25f4e2b212e2..d00100b13ab60e156a10f1ca9c80b8c40fbcad10 100644 (file)
@@ -156,7 +156,7 @@ public:
 
        // pattern matching
        bool has(const ex & other, unsigned options = 0) const { return inherited::has(other, options); }
-       bool match(const ex & pattern, lst & repl_lst) const { return inherited::match(pattern, repl_lst); }
+       bool match(const ex & pattern, exmap& repl_lst) const { return inherited::match(pattern, repl_lst); }
 protected:
        bool match_same_type(const basic & other) const { return true; }
 public:
index ea2a7f7742948781585716b2484ba01f47599e90..b5abfa1463a18c8bef5bf4696ea1e8bc650fb3fe 100644 (file)
@@ -111,7 +111,7 @@ unsigned wildcard::calchash() const
        return hashvalue;
 }
 
-bool wildcard::match(const ex & pattern, lst & repl_lst) const
+bool wildcard::match(const ex & pattern, exmap& repl_lst) const
 {
        // Wildcards must match each other exactly (this is required for
        // subs() to work properly because in the final step it substitutes
index ba4141243638ba0c2b39d265d167192d996a6f22..9729237c97573a0aa038f3d3889b1be39c6d3852 100644 (file)
@@ -41,7 +41,7 @@ public:
 
        // functions overriding virtual functions from base classes
 public:
-       bool match(const ex & pattern, lst & repl_lst) const;
+       bool match(const ex & pattern, exmap& repl_lst) const;
 
 protected:
        unsigned calchash() const;
index b059172ef5166e3b94a6ad1b3690d3e9c9c1088d..6b2c93764a9c14286d2c24a6bea7259d7aa1f00a 100644 (file)
@@ -426,9 +426,12 @@ static ex f_evalf2(const exprseq &e)
 
 static ex f_find(const exprseq &e)
 {
-       lst found;
+       exset found;
        e[0].find(e[1], found);
-       return found;
+       lst l;
+       for (exset::const_iterator i = found.begin(); i != found.end(); ++i)
+               l.append(*i);
+       return l;
 }
 
 static ex f_fsolve(const exprseq &e)
@@ -478,11 +481,14 @@ static ex f_map(const exprseq &e)
 
 static ex f_match(const exprseq &e)
 {
-       lst repl_lst;
-       if (e[0].match(e[1], repl_lst))
+       exmap repls;
+       if (e[0].match(e[1], repls)) {
+               lst repl_lst;
+               for (exmap::const_iterator i = repls.begin(); i != repls.end(); ++i)
+                       repl_lst.append(relational(i->first, i->second, relational::equal));
                return repl_lst;
-       else
-               return fail();
+       }
+       throw std::runtime_error("FAIL");
 }
 
 static ex f_normal2(const exprseq &e)
index 52280a3f475ce916517a6ff655d42d60f2955102..824d1feb06bd7bea1bcea954b8cbedafc0588398 100644 (file)
@@ -6,7 +6,7 @@ viewgar_LDADD = ../ginac/libginac.la
 
 bin_SCRIPTS = ginac-excompiler
 
-INCLUDES = -I$(srcdir)/../ginac -I../ginac
+AM_CPPFLAGS = -I$(srcdir)/../ginac -I../ginac -DIN_GINAC
 
 man_MANS = viewgar.1