EXAMS = exam_paranoia \
exam_heur_gcd \
+ match_bug \
+ parser_bugs \
exam_numeric_archive \
exam_numeric \
exam_powerlaws \
exam_archive \
exam_structure \
exam_hashmap \
- exam_misc
+ exam_misc \
+ exam_mod_gcd \
+ exam_cra
TIMES = time_dennyfliegner \
time_gammaseries \
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
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
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
--- /dev/null
+#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
+
--- /dev/null
+#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 << "]";
+}
+
--- /dev/null
+#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;
+}
+
--- /dev/null
+/**
+ * @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;
+}
+
--- /dev/null
+/// @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;
+}
+
#include <vector>
#include <stdexcept>
#include "ginac.h"
-#include "parser/parser.hpp"
#include "timer.h"
extern void randomify_symbol_serials();
using namespace std;
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)
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;
}
## 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
--- /dev/null
+// 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;
+ }
+}
+
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)}
@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:
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
@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>
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
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
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
* 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++)
}
/** 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,
// 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 {
// 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))
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
}
}
// 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:
#include "power.h"
#include "lst.h"
#include "relational.h"
-#include "input_lexer.h"
#include "utils.h"
namespace GiNaC {
/** 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);
}
* 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;
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
//////////
// 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;
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)
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)
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
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
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 {
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;
#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;
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());
}
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;
}
-// 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);
+ 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)));
}
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)) {
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);
}
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) {
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) {
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);
}
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)) {
}
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);
}
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) {
}
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 (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);
}
}
}
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);
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;
//
parameter.remove_last();
- int lastentry = parameter.nops();
+ std::size_t lastentry = parameter.nops();
while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
lastentry--;
}
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]);
}
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") {
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();
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") {
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") {
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") {
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") {
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") {
// 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;
} 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;
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 {
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]);
}
// 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;
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;
/ 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;
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 {
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 {
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]);
}
// 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;
/ 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;
/ 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;
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 {
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 {
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]);
}
}
}
- 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();
}
// 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;
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];
+++ /dev/null
-/** @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__
+++ /dev/null
-/** @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
+++ /dev/null
-/** @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;
-}
// 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;
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;
patternexpsign = 1;
}
- lst saverepls = repls;
+ exmap saverepls = repls;
if (origexponent < patternexponent || origexpsign != patternexpsign || !origbase.match(patternbase,saverepls))
return false;
repls = saverepls;
* 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)
{
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;
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);
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;
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;
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);
}
}
return inherited::info(inf);
}
-typedef std::vector<int> intvector;
+typedef std::vector<std::size_t> uintvector;
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;
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;
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;
}
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;
}
}
{
/**
- * 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.
* 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);
}
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
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:
/// 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")
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;
}
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();
}
#ifndef GINAC_PARSER_HPP_
+#define GINAC_PARSER_HPP_
#include "parse_context.hpp"
#include <stdexcept>
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();
* 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
{
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
--- /dev/null
+/// @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"));
+}
+
+}
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
+
--- /dev/null
+#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
}
// 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
{
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);
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;
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) ||
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) ||
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);
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]);
// 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:
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
// 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;
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)
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)
bin_SCRIPTS = ginac-excompiler
-INCLUDES = -I$(srcdir)/../ginac -I../ginac
+AM_CPPFLAGS = -I$(srcdir)/../ginac -I../ginac -DIN_GINAC
man_MANS = viewgar.1
git://www.ginac.de / ginac.git / commitdiff