cout << "." << flush;
}
+ Digits = digitsbuf;
+
return result;
}
static unsigned inifcns_test_legacy()
{
+ int digitsbuf = Digits;
Digits = 17;
ex prec = 5 * pow(10, -(ex)Digits);
result++;
}
+ Digits = digitsbuf;
+
return result;
}
polynomial/mod_gcd.cpp \
polynomial/cra_garner.cpp \
polynomial/gcd_euclid.h \
+polynomial/remainder.cpp \
polynomial/remainder.h \
+polynomial/normalize.cpp \
polynomial/normalize.h \
polynomial/upoly.h \
polynomial/ring_traits.h \
GINAC_ASSERT(serial<registered_functions().size());
const function_options &opt = registered_functions()[serial];
- // No derivative defined? Then return abstract derivative object
- if (opt.derivative_f == nullptr)
- return fderivative(serial, diff_param, seq);
-
- current_serial = serial;
- if (opt.derivative_use_exvector_args)
- return ((derivative_funcp_exvector)(opt.derivative_f))(seq, diff_param);
- switch (opt.nparams) {
- // the following lines have been generated for max. @maxargs@ parameters
+ if (opt.derivative_f) {
+ // Invoke the defined derivative function.
+ current_serial = serial;
+ if (opt.derivative_use_exvector_args)
+ return ((derivative_funcp_exvector)(opt.derivative_f))(seq, diff_param);
+ switch (opt.nparams) {
+ // the following lines have been generated for max. @maxargs@ parameters
+++ for N in range(1, maxargs + 1):
- case @N@:
- return ((derivative_funcp_@N@)(opt.derivative_f))(@seq('seq[%(n)d]', N, 0)@, diff_param);
+ case @N@:
+ return ((derivative_funcp_@N@)(opt.derivative_f))(@seq('seq[%(n)d]', N, 0)@, diff_param);
---
- // end of generated lines
+ // end of generated lines
+ }
}
- throw(std::logic_error("function::pderivative(): no diff function defined"));
+ // No derivative defined? Fall back to abstract derivative object.
+ return fderivative(serial, diff_param, seq);
}
ex function::expl_derivative(const symbol & s) const // explicit differentiation
GINAC_ASSERT(serial<registered_functions().size());
const function_options &opt = registered_functions()[serial];
- // No explicit derivative defined? Then this function shall not be called!
- if (opt.expl_derivative_f == nullptr)
- throw(std::logic_error("function::expl_derivative(): explicit derivation is called, but no such function defined"));
-
- current_serial = serial;
- if (opt.expl_derivative_use_exvector_args)
- return ((expl_derivative_funcp_exvector)(opt.expl_derivative_f))(seq, s);
- switch (opt.nparams) {
- // the following lines have been generated for max. @maxargs@ parameters
+ if (opt.expl_derivative_f) {
+ // Invoke the defined explicit derivative function.
+ current_serial = serial;
+ if (opt.expl_derivative_use_exvector_args)
+ return ((expl_derivative_funcp_exvector)(opt.expl_derivative_f))(seq, s);
+ switch (opt.nparams) {
+ // the following lines have been generated for max. @maxargs@ parameters
+++ for N in range(1, maxargs + 1):
- case @N@:
- return ((expl_derivative_funcp_@N@)(opt.expl_derivative_f))(@seq('seq[%(n)d]', N, 0)@, s);
+ case @N@:
+ return ((expl_derivative_funcp_@N@)(opt.expl_derivative_f))(@seq('seq[%(n)d]', N, 0)@, s);
---
- // end of generated lines
+ // end of generated lines
+ }
}
+ // There is no fallback for explicit deriviative.
+ throw(std::logic_error("function::expl_derivative(): explicit derivation is called, but no such function defined"));
}
ex function::power(const ex & power_param) const // power of function
GINAC_ASSERT(serial<registered_functions().size());
const function_options &opt = registered_functions()[serial];
- if (opt.power_f == nullptr)
- return (new GiNaC::power(*this, power_param))->setflag(status_flags::dynallocated |
- status_flags::evaluated);
-
- current_serial = serial;
- if (opt.power_use_exvector_args)
- return ((power_funcp_exvector)(opt.power_f))(seq, power_param);
- switch (opt.nparams) {
- // the following lines have been generated for max. @maxargs@ parameters
+ if (opt.power_f) {
+ // Invoke the defined power function.
+ current_serial = serial;
+ if (opt.power_use_exvector_args)
+ return ((power_funcp_exvector)(opt.power_f))(seq, power_param);
+ switch (opt.nparams) {
+ // the following lines have been generated for max. @maxargs@ parameters
+++ for N in range(1, maxargs + 1):
- case @N@:
- return ((power_funcp_@N@)(opt.power_f))(@seq('seq[%(n)d]', N, 0)@, power_param);
+ case @N@:
+ return ((power_funcp_@N@)(opt.power_f))(@seq('seq[%(n)d]', N, 0)@, power_param);
---
- // end of generated lines
+ // end of generated lines
+ }
}
- throw(std::logic_error("function::power(): no power function defined"));
+ // No power function defined? Fall back to returning a power object.
+ return (new GiNaC::power(*this, power_param))->setflag(status_flags::dynallocated |
+ status_flags::evaluated);
}
ex function::expand(unsigned options) const
GINAC_ASSERT(serial<registered_functions().size());
const function_options &opt = registered_functions()[serial];
- // No expand defined? Then return the same function with expanded arguments (if required)
- if (opt.expand_f == nullptr) {
- // Only expand arguments when asked to do so
- if (options & expand_options::expand_function_args)
- return inherited::expand(options);
- else
- return (options == 0) ? setflag(status_flags::expanded) : *this;
- }
-
- current_serial = serial;
- if (opt.expand_use_exvector_args)
- return ((expand_funcp_exvector)(opt.expand_f))(seq, options);
- switch (opt.nparams) {
- // the following lines have been generated for max. @maxargs@ parameters
+ if (opt.expand_f) {
+ // Invoke the defined expand function.
+ current_serial = serial;
+ if (opt.expand_use_exvector_args)
+ return ((expand_funcp_exvector)(opt.expand_f))(seq, options);
+ switch (opt.nparams) {
+ // the following lines have been generated for max. @maxargs@ parameters
+++ for N in range(1, maxargs + 1):
- case @N@:
- return ((expand_funcp_@N@)(opt.expand_f))(@seq('seq[%(n)d]', N, 0)@, options);
+ case @N@:
+ return ((expand_funcp_@N@)(opt.expand_f))(@seq('seq[%(n)d]', N, 0)@, options);
---
- // end of generated lines
+ // end of generated lines
+ }
}
- throw(std::logic_error("function::expand(): no expand of function defined"));
+ // No expand function defined? Return the same function with expanded arguments (if required)
+ if (options & expand_options::expand_function_args)
+ return inherited::expand(options);
+ else
+ return (options == 0) ? setflag(status_flags::expanded) : *this;
}
std::vector<function_options> & function::registered_functions()
*
* @see get_symbol_stats */
struct sym_desc {
+ /** Initialize symbol, leave other variables uninitialized */
+ sym_desc(const ex& s)
+ : sym(s), deg_a(0), deg_b(0), ldeg_a(0), ldeg_b(0), max_deg(0), max_lcnops(0)
+ { }
+
/** Reference to symbol */
ex sym;
if (it.sym.is_equal(s)) // If it's already in there, don't add it a second time
return;
- sym_desc d;
- d.sym = s;
- v.push_back(d);
+ v.push_back(sym_desc(s));
}
// Collect all symbols of an expression (used internally by get_symbol_stats())
--- /dev/null
+/** @file normalize.h
+ *
+ * Functions to normalize polynomials in a field. */
+
+/*
+ * GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include "normalize.h"
+
+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.
+bool normalize_in_field(umodpoly& a, cln::cl_MI* content_)
+{
+ 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;
+}
+
+} // namespace GiNaC
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;
-}
+bool normalize_in_field(umodpoly& a, cln::cl_MI* content_ = 0);
/// Make the univariate polynomial @a x unit normal. This version is used
/// for rings which are not fields.
* @see get_symbol_stats */
struct sym_desc
{
+ /** Initialize symbol, leave other variables uninitialized */
+ sym_desc(const ex& s)
+ : sym(s), deg_a(0), deg_b(0), ldeg_a(0), ldeg_b(0), max_deg(0), max_lcnops(0)
+ { }
+
/** Reference to symbol */
ex sym;
if (it.sym.is_equal(s)) // If it's already in there, don't add it a second time
return;
}
- sym_desc d;
- d.sym = s;
- v.push_back(d);
+ v.push_back(sym_desc(s));
}
// Collect all symbols of an expression (used internally by get_symbol_stats())
--- /dev/null
+/** @file remainder.h
+ *
+ * Functions calculating remainders. */
+
+/*
+ * GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include "remainder.h"
+#include "ring_traits.h"
+#include "upoly_io.h"
+#include "debug.h"
+
+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.
+ */
+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();
+}
+
+} // namespace GiNaC
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();
-}
+bool remainder_in_field(umodpoly& r, const umodpoly& a, const umodpoly& b);
/**
* @brief Polynomial remainder for univariate polynomials over a ring.
// Convert Z[x] -> Z/p[x]
-static void
+static inline void
make_umodpoly(umodpoly& up, const upoly& p, const cln::cl_modint_ring& R)
{
for (std::size_t i = p.size(); i-- != 0; )
* Definitions
*/
+%option nounput
+
%pointer
%{
git://www.ginac.de / ginac.git / commitdiff