...which is a C++20 feature. This fixes MSVC build.
Reported by Igor Machado <igor.machado@gmail.com>.
-Contacing the developers
-------------------------
+Contacting the developers
+-------------------------
If you have found a bug, have a patch or a question, or would like to
make a suggestion please send email to one of our public mailing lists
Christian Bauer <Christian.Bauer@uni-mainz.de>
Chris Dams <Chris.Dams@mi.infn.it>
Matthias Dellweg <dellweg@tp1.uni-duesseldorf.de>
- Richard B. Kreckel <kreckel@in.terlu.de>
Do Hoang Son <dhson@phys.hcmuns.edu.vn>
+ Oleg Finkelshteyn <olegfink@gmail.com>
Alexander Frink <Alexander.Frink@uni-mainz.de>
Vladimir V. Kisil <kisilv@maths.leeds.ac.uk>
+ Richard B. Kreckel <kreckel@in.terlu.de>
Vitaly Magerya <vmagerya@gmail.com>
Markus Nullmeier <markus.nullmeier@urz.uni-heidelberg.de>
Bernard Parisse <Bernard.Parisse@ujf-grenoble.fr>
from <http://pkg-config.freedesktop.org/>. Also, Python 3 is required.
To build the GiNaC tutorial and reference manual the doxygen utility
-(it can be downloaded from http://www.stack.nl/~dimitri/doxygen) and
-TeX are necessary.
+(it can be downloaded from https://www.doxygen.nl/) and TeX are necessary.
Known to work with:
- Linux on x86 and x86_64 using
skipped if you don't intend to use ginsh (i.e. you need the GiNaC library
for compiling another piece of a software).
7. (optional) To build the GiNaC tutorial and reference manual the doxygen
- utility (it can be downloaded from http://www.stack.nl/~dimitri/doxygen)
- and TeX are necessary.
+ utility (it can be downloaded from https://www.doxygen.nl/) and TeX are
+ necessary.
INSTALLATION
============
cmake/modules/FindGiNaC.cmake \
cmake/modules/FindLibDL.cmake
-BUILD_HELPERS = scripts/yaptu.py scripts/fixupind.py
+BUILD_HELPERS = scripts/yaptu.py
# All the rest of the distributed files
EXTRA_DIST = ginac.pc GiNaC.spec $(BUILD_HELPERS) $(CMAKE_FILES)
This file records noteworthy changes.
+1.8.7 (12 August 2023)
+* Fix series expansion of polynomial(x)^n for small and large n.
+* Fix bugs in internal parser from strings.
+* Make ginsh evaluate line-by-line in non-interactive mode.
+* Several build fixes.
+
+1.8.6 (8 February 2023)
+* Fix wrong numeric info on transcendental functions.
+* Fix crash of evaluation of binomial(n, k) with negative integer n, k.
+
+1.8.5 (1 January 2023)
+* Speed up multivariate polynomial factorization; fix it in some rare
+ corner cases where it didn't previously terminate.
+
+1.8.4 (19 September 2022)
+* Add support for sqrfree_parfrac().
+* Add info methods for transcendental functions.
+
+1.8.3 (23 March 2022)
+* series_to_poly() can be used from ginsh.
+* Fix power::to_polynomial() for posint exponents.
+* Fix power::subs() in some special cases.
+
+1.8.2 (1 January 2022)
+* Fix elusive bug in comparing relational objects.
+* Ensure modular_form_kernel::series() includes an Order term.
+
+1.8.1 (9 August 2021)
+* Add method relational::canonical() and improve conversion of relational to
+ Boolean (it now works on many simple symbolic cases).
+* Improve normalization of negative exponents.
+* Fix indexing multiply referenced objects with ex::operator[].
+* Make functions evalf() their arguments before doing own evalf().
+* Fix bugs in H_evalf() and in evaluation of iterated integrals.
+* Several portability improvements and compiler warning fixes.
+
1.8.0 (14 October 2020)
* New routines for the numerical evaluation of iterated integrals like
elliptic multiple polylogarithms or iterated integrals of modular forms.
define(GINAC_GET_LTVERSION,
[GINAC_HEADER_GETVAL(GINAC_LT_$1,[ginac/version.h])])
-dnl Usage: GINAC_STD_CXX_HEADERS
-dnl Check for standard C++ headers, bail out if something is missing.
-AC_DEFUN([GINAC_STD_CXX_HEADERS], [
-AC_CACHE_CHECK([for standard C++ header files], [ginac_cv_std_cxx_headers], [
- ginac_cv_std_cxx_headers="no"
- AC_LANG_PUSH([C++])
- AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[
- #include <algorithm>
- #include <cstddef>
- #include <cstdio>
- #include <cstdlib>
- #include <cstring>
- #include <cstdint>
- #include <ctime>
- #include <fstream>
- #include <functional>
- #include <iomanip>
- #include <ios>
- #include <iosfwd>
- #include <iostream>
- #include <iterator>
- #include <limits>
- #include <list>
- #include <map>
- #include <memory>
- #include <numeric>
- #include <ostream>
- #include <set>
- #include <sstream>
- #include <stack>
- #include <stdexcept>
- #include <string>
- #include <typeinfo>
- #include <utility>
- #include <vector>
- ]])], [ginac_cv_std_cxx_headers="yes"])
- AC_LANG_POP([C++])])
-if test "${ginac_cv_std_cxx_headers}" != "yes"; then
- AC_MSG_ERROR([Standard ISO C++ headers are missing])
-fi
-])
-
dnl Usage: GINAC_LIBREADLINE
dnl
dnl Check if GNU readline library and headers are avialable.
AC_DEFUN([GINAC_HAVE_RUSAGE],
[AC_CACHE_CHECK([whether struct rusage is declared in <sys/resource.h>],
ac_cv_have_rusage,
- [AC_TRY_COMPILE([#include <sys/times.h>
- #include <sys/resource.h>],
- [struct rusage resUsage;
- getrusage(RUSAGE_SELF, &resUsage);
- return 0;],
- [ac_cv_have_rusage=yes],
- [ac_cv_have_rusage=no])
+ [AC_COMPILE_IFELSE([AC_LANG_PROGRAM([[#include <sys/times.h>
+ #include <sys/resource.h>]],
+ [[struct rusage resUsage;
+ getrusage(RUSAGE_SELF, &resUsage);
+ return 0;]])],
+ [ac_cv_have_rusage=yes],
+ [ac_cv_have_rusage=no])
])
CONFIG_RUSAGE="no"
if test "$ac_cv_have_rusage" = yes; then
exam.gar
exam_archive
exam_clifford
+exam_collect
exam_color
exam_differentiation
exam_factor
+exam_function_exvector
exam_heur_gcd
exam_indexed
exam_inifcns
exam_numeric
exam_paranoia
exam_polygcd
+exam_relational
exam_collect_common_factors
exam_powerlaws
exam_pseries
exam_match
exam_parser
exam_numeric
+ exam_relational
exam_powerlaws
+ exam_collect
exam_inifcns
exam_inifcns_nstdsums
exam_inifcns_elliptic
exam_pgcd
exam_mod_gcd
exam_real_imag
- exam_chinrem_gcd)
+ exam_chinrem_gcd
+ exam_function_exvector
+)
set(ginac_checks
check_numeric
exam_match \
exam_parser \
exam_numeric \
+ exam_relational \
exam_powerlaws \
+ exam_collect \
exam_inifcns \
exam_inifcns_nstdsums \
exam_inifcns_elliptic \
exam_pgcd \
exam_mod_gcd \
exam_chinrem_gcd \
+ exam_function_exvector \
exam_real_imag
CHECKS = check_numeric \
exam_numeric_SOURCES = exam_numeric.cpp
exam_numeric_LDADD = ../ginac/libginac.la
+exam_relational_SOURCES = exam_relational.cpp
+exam_relational_LDADD = ../ginac/libginac.la
+
exam_powerlaws_SOURCES = exam_powerlaws.cpp
exam_powerlaws_LDADD = ../ginac/libginac.la
+exam_collect_SOURCES = exam_collect.cpp
+exam_collect_LDADD = ../ginac/libginac.la
+
exam_inifcns_SOURCES = exam_inifcns.cpp
exam_inifcns_LDADD = ../ginac/libginac.la
exam_chinrem_gcd_SOURCES = exam_chinrem_gcd.cpp
exam_chinrem_gcd_LDADD = ../ginac/libginac.la
+exam_function_exvector_SOURCES = exam_function_exvector.cpp
+exam_function_exvector_LDADD = ../ginac/libginac.la
+
check_numeric_SOURCES = check_numeric.cpp
check_numeric_LDADD = ../ginac/libginac.la
* Test of Chinese remainder algorithm. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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 underlying symbolic manipulations. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* manipulations. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* tests on these numbers like is_integer() etc... */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Here we test GiNaC's archiving system. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* A small program exposing historical bug in poly_cra function. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Here we test GiNaC's Clifford algebra objects. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
--- /dev/null
+/** @file exam_collect.cpp
+ *
+ */
+
+/*
+ * GiNaC Copyright (C) 1999-2023 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 "ginac.h"
+using namespace GiNaC;
+
+#include <iostream>
+using namespace std;
+
+// Correctness of .collect(Z[x], x).
+static unsigned exam_collect_1()
+{
+ unsigned result = 0;
+ symbol x("x"), y("y");
+ ex a = (pow(x, 3) - 2*pow(x, 2) + 4*x) * pow(y, 2)
+ + (pow(x, 2) - x - 1) * y
+ + (x + 1);
+ ex b = pow(y, 2) * pow(x, 3)
+ + (y - 2*pow(y, 2)) * pow(x, 2)
+ + (4*pow(y, 2) - y + 1) * x
+ + (1 - y);
+
+ ex a_x = collect(a, x);
+ if (a_x != b) {
+ clog << "collect(" << a << ", " << x << ") erroneously returned "
+ << a_x << " instead of " << b << endl;
+ ++result;
+ }
+
+ ex b_y = collect(b, y);
+ if (b_y != a) {
+ clog << "collect(" << b << ", " << y << ") erroneously returned "
+ << b_y << " instead of " << a << endl;
+ ++result;
+ }
+
+ ex amb_x = collect(a - b, x);
+ if (amb_x != 0) {
+ clog << "collect(" << a - b << ", " << x << ") erroneously returned "
+ << amb_x << " instead of 0" << endl;
+ ++result;
+ }
+
+ ex amb_y = collect(a - b, y);
+ if (amb_y != 0) {
+ clog << "collect(" << a - b << ", " << y << ") erroneously returned "
+ << amb_y << " instead of 0" << endl;
+ ++result;
+ }
+
+ return result;
+}
+
+// Consistency of .collect(Z[x,y], {x,y}) with .coeff(x).
+static unsigned exam_collect_2()
+{
+ unsigned result = 0;
+ symbol x("x"), y("y"), p("p"), q("q");
+
+ ex e1 = x + y;
+ ex e2 = 1 + p + q;
+ ex a = expand(e1 * e2);
+
+ ex a_x = a.collect(x).coeff(x, 1);
+ if (a_x != e2) {
+ clog << "collect(" << a << ", " << x << ") erroneously returned "
+ << a_x << " as coefficient of " << x << endl;
+ ++result;
+ }
+
+ ex a_p = a.collect(p).coeff(p, 1);
+ if (a_p != e1) {
+ clog << "collect(" << a << ", " << p << ") erroneously returned "
+ << a_p << " as coefficient of " << p << endl;
+ ++result;
+ }
+
+ ex a_xy = a.collect(lst{x,y});
+ ex ref = e2*x + e2*y;
+ if (a_xy != ref) {
+ clog << "collect(" << a << ", {" << x << ", " << y << "}) erroneously returned "
+ << a_xy << " instead of " << ref << endl;
+ ++result;
+ }
+
+ return result;
+}
+
+// Consistency of .collect(Z[f(x)], f(x)) with .coeff(f(x)).
+static unsigned exam_collect_3()
+{
+ unsigned result = 0;
+ symbol x("x"), p("p"), q("q");
+
+ for (unsigned deg = 2; deg < 7; ++deg) {
+
+ ex a1 = expand(pow(p + q + x, deg));
+ a1 = a1.collect(x);
+
+ ex a2 = expand(pow(p + q + sin(x), deg));
+ a2 = a2.collect(sin(x));
+
+ for (unsigned i = 0; i < deg; ++i) {
+ ex a1_i = a1.coeff(x, i);
+ ex a2_i = a2.coeff(sin(x), i);
+ if (!expand(a1_i - a2_i).is_zero()) {
+ clog << "collect(" << a1 << ",sin(x)) inconsistent with "
+ "collect(" << a2 << ",x)" << endl;
+ ++result;
+ }
+ }
+ }
+
+ return result;
+}
+
+unsigned exam_collect()
+{
+ unsigned result = 0;
+
+ cout << "examining collect coefficients" << flush;
+
+ result += exam_collect_1(); cout << '.' << flush;
+ result += exam_collect_2(); cout << '.' << flush;
+ result += exam_collect_3(); cout << '.' << flush;
+
+ return result;
+}
+
+int main(int argc, char** argv)
+{
+ return exam_collect();
+}
* Here we test GiNaC's color objects (su(3) Lie algebra). */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Tests for symbolic differentiation, including various functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Factorization test suite. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
return result;
}
+static unsigned exam_factor_magerya()
+{
+ // In 2017, Vitaly Magerya reported a class of biviariate polynomials
+ // where Hensel lifting sometimes failed to terminate.
+ // https://www.ginac.de/pipermail/ginac-list/2017-December/002162.html
+ unsigned result = 0;
+ ex e;
+ symbol x("x"), y("y");
+
+ e = (1+2*x+y)*(1+2*x-y)*(2*x-y)*(2*x+y);
+ result += check_factor_expanded(e);
+
+ e = (7+4*x+y)*(-5+2*x-y)*(-6+6*x+y)*y*(10+2*x+y);
+ result += check_factor_expanded(e);
+
+ e = (8+6*x-y)*(-5+4*x-y)*(-5+6*x+y)*(-2+2*x-y)*(2+4*x+y);
+ result += check_factor_expanded(e);
+
+ e = -(-4+4*x+5*y)*(1+4*x+5*y)*(2+3*y)*(1+2*x-y)*(4+2*x+y);
+ result += check_factor_expanded(e);
+
+ e = (-3+y-2*x)*(4+3*y-4*x)*(3+3*y+2*x)*(-2+3*y+2*x)*(-1+4*y+3*x);
+ result += check_factor_expanded(e);
+
+ e = (-9+7*x+y)*(-5+6*x+y)*(4+2*x+y)*(5+2*x-y)*(7+9*x-y)*(8+6*x-y);
+ result += check_factor_expanded(e);
+
+ e = pow(2*x-y,2)*(-1+6*x-y)*(-1+3*x-y)*(-2+4*x-y)*(1+4*x-y)*(4*x-y)*(2+4*x-y);
+ result += check_factor_expanded(e);
+
+ e = (5+2*y-3*x)*(-4+4*y+3*x)*(-3+4*y-2*x)*(4+5*y-x)*(3*y+2*x)*(-1+3*y+5*x)*(5+3*y+4*x);
+ result += check_factor_expanded(e);
+
+ return result;
+}
+
unsigned exam_factor()
{
unsigned result = 0;
result += exam_factor3(); cout << '.' << flush;
result += exam_factor_content(); cout << '.' << flush;
result += exam_factor_wang(); cout << '.' << flush;
+ result += exam_factor_magerya(); cout << '.' << flush;
return result;
}
--- /dev/null
+#ifdef IN_GINAC
+#include "ginac.h"
+#else
+#include "ginac/ginac.h"
+#endif
+
+#include <vector>
+#include <iostream>
+
+using namespace GiNaC;
+using namespace std;
+
+DECLARE_FUNCTION_2P(foobar);
+
+static bool eval_called = false;
+static exvector eval_called_with = {};
+static bool evalf_called = false;
+static exvector evalf_called_with = {};
+
+static void reset() {
+ eval_called_with.clear();
+ evalf_called_with.clear();
+ evalf_called = false;
+ eval_called = false;
+}
+
+static ex foobar_eval(const exvector& args) {
+ eval_called = true;
+ for (auto const& v: args)
+ eval_called_with.push_back(v);
+
+ return foobar(args[0], args[1]).hold();
+}
+
+static ex foobar_evalf(const exvector& args) {
+ evalf_called = true;
+ for (auto const& v: args)
+ evalf_called_with.push_back(v);
+ return foobar(args[0], args[1]).hold();
+}
+
+
+REGISTER_FUNCTION(foobar, eval_func(foobar_eval).
+ evalf_func(foobar_evalf));
+
+static int check_exvector_eval() {
+ symbol x("x"), y("y");
+ int err = 1;
+
+ reset();
+ ex e = foobar(x, y);
+ if (!eval_called) {
+ clog << "*** Error: " << __func__ << ": foobar_eval hasn't been called" << endl;
+ err *= 2;
+ }
+ if (eval_called_with.size() != 2) {
+ clog << "*** Error: " << __func__ << ": fobar_eval: expected 2 arguments, got " <<
+ eval_called_with.size() << endl;
+ err *= 3;
+ }
+ if (eval_called_with[0] != x) {
+ clog << "*** Error: " << __func__ << ": fobar_eval: wrong 1st argument, "
+ "expected " << x << ", got " << eval_called_with[0] << endl;
+ err *= 5;
+ }
+ if (eval_called_with[1] != y) {
+ clog << "*** Error: " << __func__ << ": fobar_eval: wrong 1st argument, "
+ "expected " << y << ", got " << eval_called_with[1] << endl;
+ err *= 7;
+ }
+ return err - 1;
+}
+
+static int check_exvector_evalf() {
+ int err = 1;
+
+ reset();
+ ex e = foobar(Pi, Euler);
+ e = e.evalf();
+
+ if (!evalf_called) {
+ clog << "*** Error: " << __func__ << ": foobar_evalf hasn't been called" << endl;
+ err *= 2;
+ }
+ if (evalf_called_with.size() != 2) {
+ clog << __func__ << ": foobar_evalf: expected 2 arguments, got " <<
+ evalf_called_with.size() << endl;
+ err *= 3;
+ }
+ if (!is_a<numeric>(evalf_called_with[0])) {
+ clog << "*** Error: " << __func__ << ": wrong 1st argument of foobar_evalf: "
+ "expected a real number, got " << evalf_called_with[0] << endl;
+ err *= 5;
+ }
+ if (!is_a<numeric>(evalf_called_with[1])) {
+ clog << "*** Error: " << __func__ << ": wrong 1st argument of foobar_evalf: "
+ "expected a real number, got " << evalf_called_with[0] << endl;
+ err *= 7;
+ }
+ return err - 1;
+}
+
+int main(int argc, char** argv) {
+ int ret = 0;
+ auto err = check_exvector_evalf();
+ if (err) {
+ ret |= 1;
+ clog << "*** Error " << (err + 1) << " (check_exvector_evalf)" << endl;
+ }
+ err = check_exvector_eval();
+ if (err) {
+ ret |= 2;
+ clog << "*** Error " << (err + 1) << " (check_exvector_evalf)" << endl;
+ }
+ return ret;
+}
* endless loop or (even worse) wrong result. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Here we test manipulations on GiNaC's indexed objects. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
using GiNaC::log;
int digitsbuf = Digits;
Digits = 17;
- ex prec = 5 * pow(10, -(ex)Digits);
+ // 15.01.2022: prec set to 10*pow(10,-Digits) to avoid exam failure in sporadic cases
+ ex prec = 10 * pow(10, -(ex)Digits);
numeric almostone("0.999999999999999999");
unsigned result = 0;
* Comparison data for the test of polylogarithms. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* These exams test solving small linear systems of symbolic equations. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* https://www.ginac.de/pipermail/ginac-devel/2006-April/000942.html */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Here we examine manipulations on GiNaC's symbolic matrices. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
bool caught = false;
try {
m5 = inverse(m4);
- } catch (std::runtime_error err) {
+ } catch (std::runtime_error &err) {
caught = true;
}
if (!caught) {
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
++result;
}
+ // This used to fail in GiNaC 1.8.2 with subs_options::no_pattern
+ e1 = 1/x;
+ e2 = e1.subs(x == 1/x);
+ if (!e2.is_equal(x)) {
+ clog << "(1/x).subs(x==1/x) erroneously returned " << e2 << " instead of x" << endl;
+ ++result;
+ }
+ e2 = e1.subs(x == 1/x, subs_options::no_pattern);
+ if (!e2.is_equal(x)) {
+ clog << "(1/x).subs(x==1/x, subs_options::no_pattern) erroneously returned " << e2 << " instead of x" << endl;
+ ++result;
+ }
+ e2 = e1.subs(x == 1/x, subs_options::algebraic);
+ if (!e2.is_equal(x)) {
+ clog << "(1/x).subs(x==1/x, subs_options::algebraic) erroneously returned " << e2 << " instead of x" << endl;
+ ++result;
+ }
+
return result;
}
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Rational function normalization test suite. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
result += check_normal(e, den);
}
+ // Negative exponents
+ e = (exp(2*x)-exp(-2*x))/(exp(x)-exp(-x));
+ ex en = e.normal();
+ // Either exp(x) or exp(-x) can be viewed as a "symbol" during run-time
+ // thus two different forms of the result are possible
+ ex r1 = (exp(2*x)+1)/exp(x) ;
+ ex r2 = (exp(-2*x)+1)/exp(-x);
+
+ if (!en.is_equal(r1) && !en.is_equal(r2)) {
+ clog << "normal form of " << e << " erroneously returned "
+ << en << " (should be " << r1 << " or " << r2 << ")" << endl;
+ result += 1;
+ }
+
return result;
}
e /= 2*pow(b, y/2)-3*pow(b, z/2);
d = 2*pow(b, y/2)+3*pow(b, z/2);
result += check_normal(e, d);
+
+ // Negative powers
+ e = (b -pow(b,-1));
+ e /= (pow(b, numeric(1,2)) - pow(b, numeric(-1,2)));
+ d = (b+1)*pow(b, numeric(-1,2));
+ result += check_normal(e, d);
}
return result;
* tests on these numbers like is_integer() etc... */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* these oopses for good, so we run those stupid tests... */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
return result;
}
+// Bug in collect()
+// cf. https://www.ginac.de/pipermail/ginac-list/2021-March/002337.html
+static unsigned exam_collect_multiply_referenced_lst()
+{
+ unsigned result = 0;
+ symbol x("x"), y("y");
+ ex a = x + y;
+ ex l = lst{x, y};
+ ex l2 = l; // make l a multiply referenced object
+
+ try {
+ ex b = collect(a, l);
+ } catch (const std::runtime_error & e) {
+ clog << "collect(" << ", " << l << ") threw a runtime_error("
+ << e.what() << ")" << endl;
+ ++result;
+ }
+
+ return result;
+}
+
unsigned exam_paranoia()
{
unsigned result = 0;
result += exam_paranoia24(); cout << '.' << flush;
result += exam_paranoia25(); cout << '.' << flush;
result += exam_paranoia26(); cout << '.' << flush;
+ result += exam_collect_multiply_referenced_lst(); cout << '.' << flush;
return result;
}
* Check for some silly bugs in the parser. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
}
}
+// Check that two strings parse to equal expressions.
+static int check_eq(ostream &err_str, parser &reader, const char *expr1, const char *expr2)
+{
+ const string srep1(expr1);
+ const string srep2(expr2);
+ ex e1, e2;
+ try{ e1 = reader(srep1); } catch (const exception &e) {
+ err_str << "\"" << srep1 << "\" failed to parse: "
+ << e.what() << endl;
+ return 1;
+ }
+ try{ e2 = reader(srep2); } catch (const exception &e) {
+ err_str << "\"" << srep2 << "\" failed to parse: "
+ << e.what() << endl;
+ return 1;
+ }
+ if (!(e1-e2).expand().is_zero()) {
+ err_str << "\"" << srep1 << "\" was misparsed as \""
+ << e1 << "\"" << endl;
+ return 1;
+ }
+ return 0;
+}
+
+// Tests for the interaction of the '^' operator with
+// the unary '+' and the unary '-' operators
+static int check5(ostream& err_str)
+{
+ parser reader;
+ return
+ +check_eq(err_str, reader, "3^2+1", "10")
+ +check_eq(err_str, reader, "3^+2-1", "8")
+ +check_eq(err_str, reader, "3^-2+1", "10/9")
+ +check_eq(err_str, reader, "3^-2/5", "1/45")
+ +check_eq(err_str, reader, "3^-2*5", "5/9")
+ +check_eq(err_str, reader, "3^-2-5", "-44/9")
+ +check_eq(err_str, reader, "3^(-2)+1", "10/9")
+ +check_eq(err_str, reader, "(3)^(-2)+1", "10/9")
+ +check_eq(err_str, reader, "+3^2+1", "10")
+ +check_eq(err_str, reader, "+3^+2+1", "10")
+ +check_eq(err_str, reader, "+3^-2+1", "10/9")
+ +check_eq(err_str, reader, "+3^-2/5", "1/45")
+ +check_eq(err_str, reader, "+3^-2*5", "5/9")
+ +check_eq(err_str, reader, "+3^-2-5", "-44/9")
+ +check_eq(err_str, reader, "-3^2+1", "-8")
+ +check_eq(err_str, reader, "-3^+2+1", "-8")
+ +check_eq(err_str, reader, "-3^-2+1", "8/9")
+ +check_eq(err_str, reader, "-3^-2/3", "-1/27")
+ +check_eq(err_str, reader, "1+2^3^4", "1+(2^81)")
+ +check_eq(err_str, reader, "2^3^4+1", "1+(2^81)")
+ +check_eq(err_str, reader, "2^+3^4+1", "1+(2^81)")
+ +check_eq(err_str, reader, "2^3^+4+1", "1+(2^81)");
+}
+
+// Tests for the interaction of the '*' operator with
+// the unary '+' and the unary '-' operators
+static int check6(ostream& err_str)
+{
+ parser reader;
+ return
+ +check_eq(err_str, reader, "3*+2-1", "5")
+ +check_eq(err_str, reader, "3*2+1", "7")
+ +check_eq(err_str, reader, "3*+2+1", "7")
+ +check_eq(err_str, reader, "3*-2+1", "-5")
+ +check_eq(err_str, reader, "3*-2/5", "-6/5")
+ +check_eq(err_str, reader, "3*-2*5", "-30")
+ +check_eq(err_str, reader, "3*-2-5", "-11")
+ +check_eq(err_str, reader, "3*(-2)+1", "-5")
+ +check_eq(err_str, reader, "(3)*(-2)+1", "-5")
+ +check_eq(err_str, reader, "+3*2+1", "7")
+ +check_eq(err_str, reader, "+3*+2+1", "7")
+ +check_eq(err_str, reader, "+3*-2+1", "-5")
+ +check_eq(err_str, reader, "+3*-2/5", "-6/5")
+ +check_eq(err_str, reader, "+3*-2*5", "-30")
+ +check_eq(err_str, reader, "+3*-2-5", "-11")
+ +check_eq(err_str, reader, "-3*2+1", "-5")
+ +check_eq(err_str, reader, "-3*+2+1", "-5")
+ +check_eq(err_str, reader, "-3*-2+1", "7")
+ +check_eq(err_str, reader, "-3*-2/3", "2")
+ +check_eq(err_str, reader, "1+2*3*4", "25")
+ +check_eq(err_str, reader, "2*3*4+1", "25")
+ +check_eq(err_str, reader, "2*+3*4+1", "25")
+ +check_eq(err_str, reader, "2*3*+4+1", "25");
+}
+
+// Tests for nested unary + and unary -
+static int check7(ostream& err_str)
+{
+ parser reader;
+ return
+ +check_eq(err_str, reader, "+1", "1")
+ +check_eq(err_str, reader, "++1", "1")
+ +check_eq(err_str, reader, "+-+1", "-1")
+ +check_eq(err_str, reader, "+-+-1", "1")
+ +check_eq(err_str, reader, "+-+-+1", "1")
+ +check_eq(err_str, reader, "100++--+1+10", "111");
+}
+
int main(int argc, char** argv)
{
cout << "examining old parser bugs" << flush;
errors += check2(err_str); cout << '.' << flush;
errors += check3(err_str); cout << '.' << flush;
errors += check4(err_str); cout << '.' << flush;
+ errors += check5(err_str); cout << '.' << flush;
+ errors += check6(err_str); cout << '.' << flush;
+ errors += check7(err_str); cout << '.' << flush;
if (errors) {
cout << "Yes, unfortunately:" << endl;
cout << err_str.str();
* rational function normalization in normalization.cpp. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* this code, it is a sanity check rather deeply rooted in GiNaC's classes. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Series expansion test (Laurent and Taylor series). */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
return result;
}
+// Test expansion of powers of polynomials.
+static unsigned exam_series15()
+{
+ unsigned result = 0;
+
+ ex e = pow(x + pow(x,2), 2);
+
+ result += check_series(e, 0, Order(1), 0);
+ result += check_series(e, 0, Order(x), 1);
+ result += check_series(e, 0, Order(pow(x,2)), 2);
+ result += check_series(e, 0, pow(x,2) + Order(pow(x,3)), 3);
+ result += check_series(e, 0, pow(x,2) + 2*pow(x,3) + Order(pow(x,4)), 4);
+ result += check_series(e, 0, pow(x,2) + 2*pow(x,3) + pow(x,4), 5);
+ result += check_series(e, 0, pow(x,2) + 2*pow(x,3) + pow(x,4), 6);
+
+ return result;
+}
+
unsigned exam_pseries()
{
unsigned result = 0;
result += exam_series12(); cout << '.' << flush;
result += exam_series13(); cout << '.' << flush;
result += exam_series14(); cout << '.' << flush;
+ result += exam_series15(); cout << '.' << flush;
return result;
}
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
--- /dev/null
+/** @file exam_relational.cpp
+ *
+ * Small test for the relational objects and their conversion to bool. */
+
+/*
+ * GiNaC Copyright (C) 1999-2023 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 "ginac.h"
+using namespace GiNaC;
+
+#include <iostream>
+using namespace std;
+
+// Elementary relations should fall back to numeric comparisons.
+static unsigned exam_relational_elementary()
+{
+ unsigned result = 0;
+ ex one = 1, two = 2;
+
+ if (one == two) {
+ clog << "'" << one << " == " << two << "' was converted to true." << endl;
+ result += 1;
+ }
+ if (!(one != two)) {
+ clog << "'" << one << " != " << two << "' was not converted to true." << endl;
+ result += 1;
+ }
+ if (!(one < two)) {
+ clog << "'" << one << " < " << two << "' was not converted to true." << endl;
+ result += 1;
+ }
+ if (!(one <= two)) {
+ clog << "'" << one << " <= " << two << "' was not converted to true." << endl;
+ result += 1;
+ }
+ if (one > two) {
+ clog << "'" << one << " > " << two << "' was converted to true." << endl;
+ result += 1;
+ }
+ if (one >= two) {
+ clog << "'" << one << " >= " << two << "' was converted to true." << endl;
+ result += 1;
+ }
+
+ return result;
+}
+
+// These should fall back to looking up info flags.
+static unsigned exam_relational_possymbol()
+{
+ unsigned result = 0;
+ possymbol p("p");
+
+ if (p == 0) {
+ clog << "for positive p, 'p == 0' was converted to true." << endl;
+ result += 1;
+ }
+ if (!(p != 0)) {
+ clog << "for positive p, 'p != 0' was not converted to true." << endl;
+ result += 1;
+ }
+ if (p < 0) {
+ clog << "for positive p, 'p < 0' was converted to true." << endl;
+ result += 1;
+ }
+ if (p <= 0) {
+ clog << "for positive p, 'p <= 0' was converted to true." << endl;
+ result += 1;
+ }
+ if (!(p > 0)) {
+ clog << "for positive p, 'p > 0' was not converted to true." << endl;
+ result += 1;
+ }
+ if (!(p >= 0)) {
+ clog << "for positive p, 'p >= 0' was not converted to true." << endl;
+ result += 1;
+ }
+
+ return result;
+}
+
+// Very simple arithmetic should be supported, too.
+static unsigned exam_relational_arith()
+{
+ unsigned result = 0;
+ possymbol p("p");
+
+ if (!(p + 2 > p + 1)) {
+ clog << "for positive p, 'p + 2 > p + 1' was not converted to true." << endl;
+ result += 1;
+ }
+ if (!(p > -p)) {
+ clog << "for positive p, 'p > -p' was not converted to true." << endl;
+ result += 1;
+ }
+ if (!(2*p > p)) {
+ clog << "for positive p, '2*p > p' was not converted to true." << endl;
+ result += 1;
+ }
+
+ return result;
+}
+
+// Comparisons should maintain ordering invariants
+static unsigned exam_relational_order()
+{
+ unsigned result = 0;
+ numeric i = 1ll<<32, j = i+1;
+ symbol a;
+ relational x = i==a, y = j==a;
+ if (x.compare(y) != -y.compare(x)) {
+ clog << "comparison should be antisymmetric." << endl;
+ result += 1;
+ }
+
+ return result;
+}
+
+unsigned exam_relational()
+{
+ unsigned result = 0;
+
+ cout << "examining relational objects" << flush;
+
+ result += exam_relational_elementary(); cout << '.' << flush;
+ result += exam_relational_possymbol(); cout << '.' << flush;
+ result += exam_relational_arith(); cout << '.' << flush;
+ result += exam_relational_order(); cout << '.' << flush;
+
+ return result;
+}
+
+int main(int argc, char** argv)
+{
+ return exam_relational();
+}
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
return result;
}
+unsigned exam_sqrfree_parfrac()
+{
+ symbol x("x");
+ // (expression, number of terms after partial fraction decomposition)
+ vector<pair<ex, unsigned>> exams = {
+ {ex("(x - 1) / (x^2*(x^2 + 2))", lst{x}), 3},
+ {ex("(1 - x^10) / x", lst{x}), 2},
+ {ex("(2*x^3 + x + 3) / ((x^2 + 1)^2)", lst{x}), 2},
+ {ex("1 / (x * (x+1)^2 * (x+2)^3)", lst{x}), 6},
+ {ex("(x*x + 3*x - 1) / (x^2*(x^2 + 2)^3)", lst{x}), 5},
+ {ex("(1 - x^10) / (x + 2)", lst{x}), 11},
+ {ex("(1 - x + 3*x^2) / (x^3 * (2+x)^2)", lst{x}), 5},
+ {ex("(1 - x) / (x^4 * (x - 2)^3)", lst{x}), 6},
+ {ex("(1 - 2*x + x^9) / (x^5 * (1 - x + x^2)^6)", lst{x}), 11}
+ };
+ unsigned result = 0;
+
+ cout << "\n"
+ << "examining square-free partial fraction decomposition" << flush;
+ for (auto e: exams) {
+ ex e1 = e.first;
+ ex e2 = sqrfree_parfrac(e1, x);
+ if (e2.nops() != e.second ||
+ !is_a<add>(e2) ||
+ !normal(e1-e2).is_zero()) {
+ clog << "sqrfree_parfrac(" << e1 << ", " << x << ") erroneously returned "
+ << e2 << endl;
+ ++result;
+ }
+ cout << '.' << flush;
+ }
+
+ return result;
+}
+
int main(int argc, char** argv)
{
- return exam_sqrfree();
+ unsigned result = 0;
+
+ result += exam_sqrfree();
+ result += exam_sqrfree_parfrac();
+
+ return result;
}
* Small test for the structure<> template. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* input in the consistency checks. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Utility functions for benchmarking. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* This program is based on work by
* Isabella Bierenbaum <bierenbaum@thep.physik.uni-mainz.de> and
* Dirk Kreimer <dkreimer@bu.edu>.
- * For details, please see <http://www.arXiv.org/abs/hep-th/0111192>.
+ * For details, please see <https://www.arXiv.org/abs/hep-th/0111192>.
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
return typeid(*vert).before(typeid(*n.vert));
// Are the indices of the top-level nodes different?
if (!(*vert==*n.vert))
- return (vert<n.vert);
+ return (*vert<*n.vert);
// Are the sets of children different, one by one?
return (children<n.children);
}
* after which e should be just a1^2. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Some timings on series expansion of the Gamma function around a pole. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* by Robert H. Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Fermat-test). */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Lewis and Michael Wester. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Time the parser. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Time the different GCD algorithms. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* A simple stop watch class. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* A simple stop watch class. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
# Output a system dependent set of variables, describing how to set the
# run time search path of shared libraries in an executable.
#
-# Copyright 1996-2015 Free Software Foundation, Inc.
+# Copyright 1996-2023 Free Software Foundation, Inc.
# Taken from GNU libtool, 2001
# Originally by Gordon Matzigkeit <gord@gnu.ai.mit.edu>, 1996
#
hardcode_direct=yes
hardcode_minus_L=yes
;;
- freebsd* | dragonfly*)
+ freebsd* | dragonfly* | midnightbsd*)
hardcode_libdir_flag_spec='-R$libdir'
hardcode_direct=yes
;;
freebsd[23].*)
library_names_spec='$libname$shrext$versuffix'
;;
- freebsd* | dragonfly*)
+ freebsd* | dragonfly* | midnightbsd*)
library_names_spec='$libname$shrext'
;;
gnu*)
m4_define([ginac_lt_age], GINAC_GET_LTVERSION([AGE]))
m4_define([ginac_lt_revision], GINAC_GET_LTVERSION([REVISION]))
-AC_INIT([GiNaC], ginac_version, [ginac-list@ginac.de], [ginac], [https://www.ginac.de/])
-AC_PREREQ(2.59)
+AC_INIT([GiNaC], [ginac_version], [ginac-list@ginac.de], [ginac], [https://www.ginac.de/])
+AC_PREREQ([2.59])
AC_CONFIG_SRCDIR(ginac/basic.cpp)
AC_CONFIG_AUX_DIR([config])
AC_CONFIG_HEADERS([config/config.h])
AC_PROG_CXX
AC_PROG_CXXCPP
AC_PROG_INSTALL
-AM_PROG_LIBTOOL
-AC_PROG_LEX
+LT_INIT
+AC_PROG_LEX([yywrap])
AC_PROG_YACC
AC_PATH_PROG(YACCEXE, $YACC, "")
AS_IF([test "x$LEX" = "x:" -a ! -f $srcdir/ginsh/ginsh_lexer.cpp],
AC_LANG([C++])
AX_CXX_COMPILE_STDCXX([11])
-dnl Make sure all the necessary standard headers are installed on the system.
-GINAC_STD_CXX_HEADERS
-
dnl We need to have CLN installed.
PKG_CHECK_MODULES(CLN, cln >= 1.2.2)
AC_LIB_LINKFLAGS_FROM_LIBS([CLN_RPATH], [$CLN_LIBS])
AC_PATH_PROG(MAKEINDEX, makeindex, "")
AC_PATH_PROG(MAKEINFO, makeinfo, "")
AC_PATH_PROG(DVIPS, dvips, "")
-AM_CONDITIONAL(CONFIG_TEX, [test ! \( -z "$LATEX" -o -z $"PDFLATEX" -o -z "$MAKEINDEX" -o -z "$DVIPS" \)])
+AM_CONDITIONAL(CONFIG_TEX, [test ! \( -z "$LATEX" -o -z "$PDFLATEX" -o -z "$MAKEINDEX" -o -z "$DVIPS" \)])
AC_PATH_PROG(FIG2DEV, fig2dev, "")
AM_CONDITIONAL(CONFIG_FIG2DEV, [test ! -z "$FIG2DEV"])
AS_IF([test -z "$FIG2DEV" -o -z "$MAKEINFO"],
add_custom_target(${thename}_html DEPENDS ${${thename}_HTML})
add_dependencies(info ${thename}_info)
add_dependencies(html ${thename}_html)
- install(FILES ${${thename}_INFO} DESTINATION "${SHARE_INSTALL_PREFIX}/info")
+ install(FILES ${${thename}_INFO} DESTINATION "${CMAKE_INSTALL_PREFIX}/share/info")
endmacro()
macro(pdflatex_process texfile)
set(_idx ${_dirname}/${_basename}.idx)
set(_ind ${_dirname}/${_basename}.ind)
set(_pdf ${_dirname}/${_basename}.pdf)
- set(_fixupind ${CMAKE_SOURCE_DIR}/scripts/fixupind.py)
add_custom_command(
OUTPUT ${_idx}
COMMAND ${PDFLATEX_COMPILER} ${texfile}
add_custom_command(
OUTPUT ${_ind}
COMMAND ${MAKEINDEX_COMPILER} ${_idx}
- COMMAND ${PYTHON} ${_fixupind} ${_idx}
WORKING_DIRECTORY ${_dirname}
DEPENDS ${texfile} ${_idx}
COMMENT "MAKEINDEX ${_basename}.idx")
// First using compile_ex
{
- double result;
+ double result = 0.0;
double point = 0.2;
start = clock();
for (int i=0; i<100000; ++i) {
// Then without compile_ex
{
- ex result;
+ ex result = 0.0;
ex point = 0.2;
start = clock();
for (int i=0; i<100000; ++i) {
* The value of g is taken to be equal to the order N.
*
* More details can be found at Wikipedia:
- * http://en.wikipedia.org/wiki/Lanczos_approximation.
+ * https://en.wikipedia.org/wiki/Lanczos_approximation.
*
* (C) 2006 Chris Dams
*
<small><i>This page is part of the <b><a
href="https://www.ginac.de/">GiNaC</a></b>
developer's reference. It was generated automatically by <a
-href="http://www.stack.nl/~dimitri/doxygen/index.html">doxygen</a>. For
+href="https://www.doxygen.nl/">doxygen</a>. For
an introduction, see the <a href="../tutorial/">tutorial</a>.</i></small>
</body>
</html>
cd pdflatex; \
${PDFLATEX} reference.tex ;\
${MAKEINDEX} reference.idx ;\
- ${PYTHON} $(abs_top_srcdir)/scripts/fixupind.py reference.ind; \
${PDFLATEX} reference.tex
reference.dvi: latex latex/reference.dvi
This is a tutorial that documents GiNaC @value{VERSION}, an open
framework for symbolic computation within the C++ programming language.
-Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+Copyright (C) 1999-2023 Johannes Gutenberg University Mainz, Germany
Permission is granted to make and distribute verbatim copies of
this manual provided the copyright notice and this permission notice
@page
@vskip 0pt plus 1filll
-Copyright @copyright{} 1999-2020 Johannes Gutenberg University Mainz, Germany
+Copyright @copyright{} 1999-2023 Johannes Gutenberg University Mainz, Germany
@sp 2
Permission is granted to make and distribute verbatim copies of
this manual provided the copyright notice and this permission notice
@section License
The GiNaC framework for symbolic computation within the C++ programming
-language is Copyright @copyright{} 1999-2020 Johannes Gutenberg
+language is Copyright @copyright{} 1999-2023 Johannes Gutenberg
University Mainz, Germany.
This program is free software; you can redistribute it and/or
The classes therein serve as foundation classes for GiNaC. CLN stands
for Class Library for Numbers or alternatively for Common Lisp Numbers.
In order to find out more about CLN's internals, the reader is referred to
-the documentation of that library. @inforef{Introduction, , cln}, for
+the documentation of that library. @xref{Top,,, cln, The CLN Manual}, for
more information. Suffice to say that it is by itself build on top of
another library, the GNU Multiple Precision library GMP, which is an
extremely fast library for arbitrary long integers and rationals as well
method, where the left hand side of the relation specifies the variable
to expand in and the right hand side the expansion point. They can also
be used for creating systems of equations that are to be solved for
-unknown variables. But the most common usage of objects of this class
+unknown variables.
+
+But the most common usage of objects of this class
is rather inconspicuous in statements of the form @code{if
(expand(pow(a+b,2))==a*a+2*a*b+b*b) @{...@}}. Here, an implicit
conversion from @code{relational} to @code{bool} takes place. Note,
however, that @code{==} here does not perform any simplifications, hence
@code{expand()} must be called explicitly.
+Simplifications of
+relationals may be more efficient if preceded by a call to
+@example
+ex relational::canonical() const
+@end example
+which returns an equivalent relation with the zero
+right-hand side. For example:
+@example
+possymbol p("p");
+relational rel = (p >= (p*p-1)/p);
+if (ex_to<relational>(rel.canonical().normal()))
+ cout << "correct inequality" << endl;
+@end example
+However, a user shall not expect that any inequality can be fully
+resolved by GiNaC.
+
@node Integrals, Matrices, Relations, Basic concepts
@c node-name, next, previous, up
@section Integrals
One solution to this dilemma is the @dfn{Visitor} design pattern,
which is implemented in GiNaC (actually, Robert Martin's Acyclic Visitor
variation, described in detail in
-@uref{http://objectmentor.com/publications/acv.pdf}). Instead of adding
+@uref{https://condor.depaul.edu/dmumaugh/OOT/Design-Principles/acv.pdf}). Instead of adding
virtual functions to the class hierarchy to implement operations, GiNaC
provides a single "bouncing" method @code{accept()} that takes an instance
of a special @code{visitor} class and redirects execution to the one
Note also, how factors with the same exponents are not fully factorized
with this method.
+@subsection Square-free partial fraction decomposition
+@cindex square-free partial fraction decomposition
+@cindex partial fraction decomposition
+@cindex @code{sqrfree_parfrac()}
+
+GiNaC also supports square-free partial fraction decomposition of
+rational functions:
+@example
+ex sqrfree_parfrac(const ex & a, const symbol & x);
+@end example
+It is called square-free because it assumes a square-free
+factorization of the input's denominator:
+@example
+ ...
+ symbol x("x");
+
+ ex rat = (x-4)/(pow(x,2)*(x+2));
+ cout << sqrfree_parfrac(rat, x) << endl;
+ // -> -2*x^(-2)+3/2*x^(-1)-3/2*(2+x)^(-1)
+@end example
+
@subsection Polynomial factorization
@cindex factorization
@cindex polynomial factorization
ex machin_pi(int degr)
@{
symbol x;
- ex pi_expansion = series_to_poly(atan(x).series(x,degr));
+ ex pi_expansion = series_to_poly(atan(x).series(x==0,degr));
ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
-4*pi_expansion.subs(x==numeric(1,239));
return pi_approx;
happens. Knowing this will enable you to write much more efficient
code. If you still have an uncertain feeling with copy-on-write
semantics, we recommend you have a look at the
-@uref{http://www.parashift.com/c++-faq-lite/, C++-FAQ lite} by
-Marshall Cline. Chapter 16 covers this issue and presents an
-implementation which is pretty close to the one in GiNaC.
+@uref{https://isocpp.org/faq, C++-FAQ's} chapter on memory management.
+It covers this issue and presents an implementation which is pretty
+close to the one in GiNaC.
@node Internal representation of products and sums, Package tools, Expressions are reference counted, Internal structures
directory @var{prefix} one should set the @var{PKG_CONFIG_PATH}
environment variable to @var{prefix}/lib/pkgconfig for this to work.}
@example
-g++ -o simple `pkg-config --cflags --libs ginac` simple.cpp
+g++ -o simple simple.cpp `pkg-config --cflags --libs ginac`
@end example
This command line might expand to (for example):
@example
-g++ -o simple -lginac -lcln simple.cpp
+g++ -o simple simple.cpp -lginac -lcln
@end example
Not only is the form using @command{pkg-config} easier to type, it will
Name: GiNaC
Description: C++ library for symbolic calculations
Version: @VERSION@
-Requires: cln >= 1.1.6
+Requires: cln >= 1.2.2
Libs: -L${libdir} -lginac @GINACLIB_RPATH@
Cflags: -I${includedir}
* Implementation of GiNaC's sums of expressions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's sums of expressions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Archiving of GiNaC expressions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
/** Delete cached unarchived expressions in all archive_nodes (mainly for debugging). */
void archive::forget()
{
- for_each(nodes.begin(), nodes.end(), std::mem_fun_ref(&archive_node::forget));
+ for_each(nodes.begin(), nodes.end(), std::mem_fn(&archive_node::forget));
}
/** Delete cached unarchived expressions from node (for debugging). */
* Archiving of GiNaC expressions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Assertion macro definition. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's ABC. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's ABC. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Helper templates to provide per-class information for class hierarchies. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's clifford algebra (Dirac gamma) objects. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
ev.reserve(e1.nops());
cv.reserve(e1.nops());
// separate clifford and non-clifford entries
- for (int i= 0; i < e1.nops(); ++i) {
+ for (size_t i= 0; i < e1.nops(); ++i) {
if (is_a<clifford>(e1.op(i)) && is_a<cliffordunit>(e1.op(i).op(0)))
cv.push_back(e1.op(i));
else
bool found=false, same_value_found=false;
ex dummy_ind=0;
ev.reserve(e1.nops());
- for (int i=0; i < e1.nops();++i) {
+ for (size_t i=0; i < e1.nops(); ++i) {
// look for a Clifford unit with the same metric and representation label,
// if found remember its index
if (is_a<clifford>(e1.op(i)) && ex_to<clifford>(e1.op(i)).get_representation_label() == rl
* Interface to GiNaC's clifford algebra (Dirac gamma) objects. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* if e contains at least one, otherwise returns -1
*
* @param e Expression to be processed
- * @ignore_ONE defines if clifford_ONE should be ignored in the search*/
+ * @param ignore_ONE defines if clifford_ONE should be ignored in the search */
int clifford_max_label(const ex & e, bool ignore_ONE = false);
/** Calculation of the norm in the Clifford algebra. */
* @param mu Index (must be of class varidx or a derived class)
* @param metr Metric (should be indexed, tensmetric or a derived class, or a matrix)
* @param rl Representation label
- * @param e Clifford unit object
* @return Clifford vector with given components */
ex lst_to_clifford(const ex & v, const ex & mu, const ex & metr, unsigned char rl = 0);
+
+/** List or vector conversion into the Clifford vector.
+ *
+ * @param v List or vector of coordinates
+ * @param e Clifford unit object
+ * @return Clifford vector with given components */
ex lst_to_clifford(const ex & v, const ex & e);
/** An inverse function to lst_to_clifford(). For given Clifford vector extracts
* Implementation of GiNaC's color (SU(3) Lie algebra) objects. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's color (SU(3) Lie algebra) objects. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Definition of optimizing macros. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's constant types and some special constants. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's constant types and some special constants. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Wrapper template for making GiNaC classes out of STL containers. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's light-weight expression handles. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
bp->dbgprinttree();
}
+/** Expand an expression.
+ * @param options see GiNaC::expand_options */
ex ex::expand(unsigned options) const
{
if (options == 0 && (bp->flags & status_flags::expanded)) // The "expanded" flag only covers the standard options; someone might want to re-expand with different options
* Interface to GiNaC's light-weight expression handles. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* of this class in every object file that makes use of our flyweights in
* order to guarantee proper initialization. Hence we put it into this
* file which is included by every relevant file anyways. This is modeled
- * after section 27.4.2.1.6 of the C++ standard, where cout and friends are
+ * after section [ios::Init] of the C++ standard, where cout and friends are
* set up.
*
* @see utils.cpp */
// operand access
size_t nops() const { return bp->nops(); }
ex op(size_t i) const { return bp->op(i); }
- ex operator[](const ex & index) const { return (*bp)[index]; }
- ex operator[](size_t i) const { return (*bp)[i]; }
+ ex operator[](const ex & index) const { return (const_cast<const basic&>(*bp))[index]; }
+ ex operator[](size_t i) const { return (const_cast<const basic&>(*bp))[i]; }
ex & let_op(size_t i);
ex & operator[](const ex & index);
ex & operator[](size_t i);
// Iterators
-class const_iterator : public std::iterator<std::random_access_iterator_tag, ex, ptrdiff_t, const ex *, const ex &> {
+class const_iterator {
friend class ex;
friend class const_preorder_iterator;
friend class const_postorder_iterator;
-
public:
+ using iterator_category = std::random_access_iterator_tag;
+ using value_type = ex;
+ using difference_type = ptrdiff_t;
+ using pointer = const ex *;
+ using reference = const ex &;
+
const_iterator() noexcept {}
private:
} // namespace internal
-class const_preorder_iterator : public std::iterator<std::forward_iterator_tag, ex, ptrdiff_t, const ex *, const ex &> {
+class const_preorder_iterator {
public:
+ using iterator_category = std::forward_iterator_tag;
+ using value_type = ex;
+ using difference_type = ptrdiff_t;
+ using pointer = const ex *;
+ using reference = const ex &;
const_preorder_iterator() noexcept {}
const_preorder_iterator(const ex &e, size_t n)
}
};
-class const_postorder_iterator : public std::iterator<std::forward_iterator_tag, ex, ptrdiff_t, const ex *, const ex &> {
+class const_postorder_iterator {
public:
+ using iterator_category = std::forward_iterator_tag;
+ using value_type = ex;
+ using difference_type = ptrdiff_t;
+ using pointer = const ex *;
+ using reference = const ex &;
const_postorder_iterator() noexcept {}
const_postorder_iterator(const ex &e, size_t n)
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* fast numerical integration. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* C code equivalent in double precision. The function pointer has type FUNCP_2P.
*
* @param expr Expression to be compiled
- * @param sym Symbol from the expression to become the function parameter
+ * @param sym1 Symbol from the expression to become the first function parameter
+ * @param sym2 Symbol from the expression to become the second function parameter
* @param fp Returned function pointer
* @param filename Name of the intermediate source code and so-file. If
* supplied, these intermediate files will not be deleted
* Takes an expression and produces a function pointer to the compiled and linked
* C code equivalent in double precision. The function pointer has type FUNCP_CUBA.
*
- * @param expr Expression to be compiled
- * @param sym Symbol from the expression to become the function parameter
+ * @param exprs List of expression to be compiled
+ * @param syms Symbols from the expression to become the function parameters
* @param fp Returned function pointer
* @param filename Name of the intermediate source code and so-file. If
* supplied, these intermediate files will not be deleted
* Implementation of expression pairs (building blocks of expairseq). */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Definition of expression pairs (building blocks of expairseq). */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of sequences of expression pairs. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to sequences of expression pairs. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's exprseq. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Definition of GiNaC's exprseq. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
typedef container<std::vector> exprseq;
+/** Declaration of container::reg_info for exprseq. */
+#ifndef _MSC_VER // workaround error C2766: explicit specialization; 'reg_info' has already been defined
+template<> registered_class_info exprseq::reg_info;
+#endif
+
// defined in exprseq.cpp
template<> bool exprseq::info(unsigned inf) const;
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
#include "normal.h"
#include "add.h"
+#include <type_traits>
#include <algorithm>
#include <limits>
#include <list>
namespace GiNaC {
+// anonymous namespace to hide all utility functions
+namespace {
+
#ifdef DEBUGFACTOR
#define DCOUT(str) cout << #str << endl
#define DCOUTVAR(var) cout << #var << ": " << var << endl
#define DCOUT2(str,var)
#endif // def DEBUGFACTOR
-// anonymous namespace to hide all utility functions
-namespace {
-
////////////////////////////////////////////////////////////////////////////////
// modular univariate polynomial code
typedef std::vector<cln::cl_I> upoly;
typedef vector<umodpoly> upvec;
-// COPY FROM UPOLY.HPP
+
+// COPY FROM UPOLY.H
// CHANGED size_t -> int !!!
template<typename T> static int degree(const T& p)
return p[p.size() - 1];
}
+/** Make the polynomial unit normal (having unit normal leading coefficient).
+ *
+ * @param[in, out] a polynomial to make unit normal
+ * @return true if polynomial a was already unit normal, false otherwise
+ */
static bool normalize_in_field(umodpoly& a)
{
if (a.size() == 0)
return false;
}
+/** Remove leading zero coefficients from polynomial.
+ *
+ * @param[in, out] p polynomial from which the zero leading coefficients will be removed
+ * @param[in] hint assume all coefficients of order ≥ hint are zero
+ */
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;
+ std::size_t i = min(p.size(), hint);
- if (hint < p.size())
- i = hint;
+ while ( i-- && zerop(p[i]) ) { }
- 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;
- }
-
- p.erase(p.begin() + i, p.end());
-}
-
-// END COPY FROM UPOLY.HPP
-
-static void expt_pos(umodpoly& a, unsigned int q)
-{
- if ( a.empty() ) return;
- cl_MI zero = a[0].ring()->zero();
- int deg = degree(a);
- a.resize(degree(a)*q+1, zero);
- for ( int i=deg; i>0; --i ) {
- a[i*q] = a[i];
- a[i] = zero;
- }
+ p.erase(p.begin() + i + 1, p.end());
}
-template<bool COND, typename T = void> struct enable_if
-{
- typedef T type;
-};
-
-template<typename T> struct enable_if<false, T> { /* empty */ };
+// END COPY FROM UPOLY.H
template<typename T> struct uvar_poly_p
{
return e;
}
-/** Divides all coefficients of the polynomial a by the integer x.
+/** Divides all coefficients of the polynomial a by the positive integer x.
* All coefficients are supposed to be divisible by x. If they are not, the
- * the<cl_I> cast will raise an exception.
+ * division will raise an exception.
*
* @param[in,out] a polynomial of which the coefficients will be reduced by x
- * @param[in] x integer that divides the coefficients
+ * @param[in] x positive integer that divides the coefficients
*/
static void reduce_coeff(umodpoly& a, const cl_I& x)
{
for (auto & i : a) {
// cln cannot perform this division in the modular field
cl_I c = R->retract(i);
- i = cl_MI(R, the<cl_I>(c / x));
+ i = cl_MI(R, exquopos(c, x));
}
}
} while ( k-- );
fill(r.begin()+n, r.end(), a[0].ring()->zero());
- canonicalize(r);
+ canonicalize(r, n);
}
/** Calculates quotient of a/b.
static bool unequal_one(const umodpoly& a)
{
- if ( a.empty() ) return true;
return ( a.size() != 1 || a[0] != a[0].ring()->one() );
}
return equal_one(c);
}
+/** Computes w^q mod a.
+ * Uses theorem 2.1 from A.K.Lenstra's PhD thesis; see exercise 8.13 in [GCL].
+ *
+ * @param[in] w polynomial
+ * @param[in] a modulus polynomial
+ * @param[in] q common modulus of w and a
+ * @param[out] r result
+ */
+static void expt_pos_Q(const umodpoly& w, const umodpoly& a, unsigned int q, umodpoly& r)
+{
+ if ( w.empty() ) return;
+ cl_MI zero = w[0].ring()->zero();
+ int deg = degree(w);
+ umodpoly buf(deg*q+1, zero);
+ for ( size_t i=0; i<=deg; ++i ) {
+ buf[i*q] = w[i];
+ }
+ rem(buf, a, r);
+}
+
// END modular univariate polynomial code
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
/** Calculates the Q matrix for a polynomial. Used by Berlekamp's algorithm.
+ *
+ * The implementation follows algorithm 8.5 of [GCL].
*
* @param[in] a_ modular polynomial
* @param[out] Q Q matrix
/** Berlekamp's modular factorization.
*
- * The implementation follows the algorithm in chapter 8 of [GCL].
+ * The implementation follows algorithm 8.4 of [GCL].
*
* @param[in] a modular polynomial
* @param[out] upv vector containing modular factors. if upv was not empty the
/** Calculates a^(1/prime).
*
- * @param[in] a polynomial
- * @param[in] prime prime number -> exponent 1/prime
- * @param[in] ap resulting polynomial
+ * @param[in] a polynomial
+ * @param[in] prime prime number -> exponent 1/prime
+ * @param[out] ap resulting polynomial
*/
static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap)
{
/** Distinct degree factorization (DDF).
*
- * The implementation follows the algorithm in chapter 8 of [GCL].
+ * The implementation follows algorithm 8.8 of [GCL].
*
* @param[in] a_ modular polynomial
* @param[out] degrees vector containing the degrees of the factors of the
int nhalf = degree(a)/2;
int i = 1;
- umodpoly w(2);
- w[0] = R->zero();
- w[1] = R->one();
+ umodpoly w = {R->zero(), R->one()};
umodpoly x = w;
while ( i <= nhalf ) {
- expt_pos(w, q);
umodpoly buf;
- rem(w, a, buf);
+ expt_pos_Q(w, a, q, buf);
w = buf;
- umodpoly wx = w - x;
- gcd(a, wx, buf);
+ gcd(a, w - x, buf);
if ( unequal_one(buf) ) {
degrees.push_back(i);
ddfactors.push_back(buf);
- }
- if ( unequal_one(buf) ) {
umodpoly buf2;
div(a, buf, buf2);
a = buf2;
return r;
}
-/** Calculates the bound for the modulus.
- * See [Mig].
+/** Calculates bound for the product of absolute values (modulus) of the roots.
+ * Uses Landau's inequality, see [Mig].
*/
-static inline cl_I calc_bound(const ex& a, const ex& x, int maxdeg)
+static inline cl_I calc_bound(const ex& a, const ex& x)
{
- cl_I maxcoeff = 0;
- cl_R coeff = 0;
+ cl_R radicand = 0;
for ( int i=a.degree(x); i>=a.ldegree(x); --i ) {
cl_I aa = abs(the<cl_I>(ex_to<numeric>(a.coeff(x, i)).to_cl_N()));
- if ( aa > maxcoeff ) maxcoeff = aa;
- coeff = coeff + square(aa);
+ radicand = radicand + square(aa);
}
- cl_I coeffnorm = ceiling1(the<cl_R>(cln::sqrt(coeff)));
- cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg));
- return ( B > maxcoeff ) ? B : maxcoeff;
+ return ceiling1(the<cl_R>(cln::sqrt(radicand)));
}
-/** Calculates the bound for the modulus.
- * See [Mig].
+/** Calculates bound for the product of absolute values (modulus) of the roots.
+ * Uses Landau's inequality, see [Mig].
*/
-static inline cl_I calc_bound(const upoly& a, int maxdeg)
+static inline cl_I calc_bound(const upoly& a)
{
- cl_I maxcoeff = 0;
- cl_R coeff = 0;
+ cl_R radicand = 0;
for ( int i=degree(a); i>=0; --i ) {
cl_I aa = abs(a[i]);
- if ( aa > maxcoeff ) maxcoeff = aa;
- coeff = coeff + square(aa);
+ radicand = radicand + square(aa);
}
- cl_I coeffnorm = ceiling1(the<cl_R>(cln::sqrt(coeff)));
- cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg));
- return ( B > maxcoeff ) ? B : maxcoeff;
+ return ceiling1(the<cl_R>(cln::sqrt(radicand)));
}
/** Hensel lifting as used by factor_univariate().
*
- * The implementation follows the algorithm in chapter 6 of [GCL].
+ * The implementation follows algorithm 6.1 of [GCL].
*
* @param[in] a_ primitive univariate polynomials
* @param[in] p prime number that does not divide lcoeff(a)
// calc bound B
int maxdeg = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_);
- cl_I maxmodulus = 2*calc_bound(a, maxdeg);
+ cl_I maxmodulus = ash(calc_bound(a), maxdeg+1); // = 2 * calc_bound(a) * 2^maxdeg
// step 1
cl_I alpha = lcoeff(a);
}
}
-/** Returns a new prime number.
+/** Returns a new small prime number.
*
- * @param[in] p prime number
- * @return next prime number after p
+ * @param[in] n an integer
+ * @return smallest prime greater than n
*/
-static unsigned int next_prime(unsigned int p)
-{
- static vector<unsigned int> primes;
- if (primes.empty()) {
- primes = {3, 5, 7};
- }
- if ( p >= primes.back() ) {
- unsigned int candidate = primes.back() + 2;
- while ( true ) {
- size_t n = primes.size()/2;
- for ( size_t i=0; i<n; ++i ) {
- if (candidate % primes[i])
- continue;
- candidate += 2;
- i=-1;
- }
- primes.push_back(candidate);
- if (candidate > p)
+static unsigned int next_prime(unsigned int n)
+{
+ static vector<unsigned int> primes = {2, 3, 5, 7};
+ unsigned int candidate = primes.back();
+ while (primes.back() <= n) {
+ candidate += 2;
+ bool is_prime = true;
+ for (size_t i=1; primes[i]*primes[i]<=candidate; ++i) {
+ if (candidate % primes[i] == 0) {
+ is_prime = false;
break;
+ }
}
- return candidate;
+ if (is_prime)
+ primes.push_back(candidate);
}
for (auto & it : primes) {
- if ( it > p ) {
+ if ( it > n ) {
return it;
}
}
throw logic_error("next_prime: should not reach this point!");
}
-/** Manages the splitting a vector of of modular factors into two partitions.
+/** Manages the splitting of a vector of modular factors into two partitions.
*/
class factor_partition
{
* with deg(s_i) < deg(a_i)
* and with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r
*
- * The implementation follows the algorithm in chapter 6 of [GCL].
+ * The implementation follows algorithm 6.3 of [GCL].
*
* @param[in] a vector of modular univariate polynomials
* @param[in] x symbol
*
* Solves s*a + t*b == 1 mod p^k given a,b.
*
- * The implementation follows the algorithm in chapter 6 of [GCL].
+ * The implementation follows algorithm 6.3 of [GCL].
*
* @param[in] a polynomial
* @param[in] b polynomial
* s_1*b_1 + ... + s_r*b_r == x^m mod p^k
* with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r
*
- * The implementation follows the algorithm in chapter 6 of [GCL].
+ * The implementation follows algorithm 6.3 of [GCL].
*
* @param a vector with univariate polynomials mod p^k
* @param x symbol
* s_1*b_1 + ... + s_r*b_r == c mod <I^(d+1),p^k>
* with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r
*
- * The implementation follows the algorithm in chapter 6 of [GCL].
+ * The implementation follows algorithm 6.2 of [GCL].
*
* @param a_ vector of multivariate factors mod p^k
* @param x symbol (equiv. x_1 in [GCL])
}
/** Multivariate Hensel lifting.
- * The implementation follows the algorithm in chapter 6 of [GCL].
+ * The implementation follows algorithm 6.4 of [GCL].
* Since we don't have a data type for modular multivariate polynomials, the
* respective operations are done in a GiNaC::ex and the function
* make_modular() is then called to make the coefficient modular p^l.
}
}
-/** Takes a factorized expression and puts the factors in a lst. The exponents
+/** Takes a factorized expression and puts the factors in a vector. The exponents
* of the factors are discarded, e.g. 7*x^2*(y+1)^4 --> {7,x,y+1}. The first
- * element of the list is always the numeric coefficient.
+ * element of the result is always the numeric coefficient.
*/
-static ex put_factors_into_lst(const ex& e)
+static exvector put_factors_into_vec(const ex& e)
{
- lst result;
+ exvector result;
if ( is_a<numeric>(e) ) {
- result.append(e);
+ result.push_back(e);
return result;
}
if ( is_a<power>(e) ) {
- result.append(1);
- result.append(e.op(0));
+ result.push_back(1);
+ result.push_back(e.op(0));
return result;
}
if ( is_a<symbol>(e) || is_a<add>(e) ) {
ex icont(e.integer_content());
- result.append(icont);
- result.append(e/icont);
+ result.push_back(icont);
+ result.push_back(e/icont);
return result;
}
if ( is_a<mul>(e) ) {
ex nfac = 1;
+ result.push_back(nfac);
for ( size_t i=0; i<e.nops(); ++i ) {
ex op = e.op(i);
if ( is_a<numeric>(op) ) {
nfac = op;
}
if ( is_a<power>(op) ) {
- result.append(op.op(0));
+ result.push_back(op.op(0));
}
if ( is_a<symbol>(op) || is_a<add>(op) ) {
- result.append(op);
+ result.push_back(op);
}
}
- result.prepend(nfac);
+ result[0] = nfac;
return result;
}
- throw runtime_error("put_factors_into_lst: bad term.");
+ throw runtime_error("put_factors_into_vec: bad term.");
}
/** Checks a set of numbers for whether each number has a unique prime factor.
*
- * @param[in] f list of numbers to check
+ * @param[in] f numbers to check
* @return true: if number set is bad, false: if set is okay (has unique
* prime factors)
*/
-static bool checkdivisors(const lst& f)
+static bool checkdivisors(const exvector& f)
{
- const int k = f.nops();
+ const int k = f.size();
numeric q, r;
vector<numeric> d(k);
- d[0] = ex_to<numeric>(abs(f.op(0)));
+ d[0] = ex_to<numeric>(abs(f[0]));
for ( int i=1; i<k; ++i ) {
- q = ex_to<numeric>(abs(f.op(i)));
+ q = ex_to<numeric>(abs(f[i]));
for ( int j=i-1; j>=0; --j ) {
r = d[j];
do {
*
* @param[in] u multivariate polynomial to be factored
* @param[in] vn leading coefficient of u in x (x==first symbol in syms)
- * @param[in] syms set of symbols that appear in u
- * @param[in] f lst containing the factors of the leading coefficient vn
+ * @param[in] x first symbol that appears in u
+ * @param[in] syms_wox remaining symbols that appear in u
+ * @param[in] f vector containing the factors of the leading coefficient vn
* @param[in,out] modulus integer modulus for random number generation (i.e. |a_i| < modulus)
* @param[out] u0 returns the evaluated (univariate) polynomial
* @param[out] a returns the valid evaluation points. must have initial size equal
* number of symbols-1 before calling generate_set
*/
-static void generate_set(const ex& u, const ex& vn, const exset& syms, const lst& f,
+static void generate_set(const ex& u, const ex& vn, const ex& x, const exset& syms_wox, const exvector& f,
numeric& modulus, ex& u0, vector<numeric>& a)
{
- const ex& x = *syms.begin();
while ( true ) {
++modulus;
// generate a set of integers ...
u0 = u;
ex vna = vn;
ex vnatry;
- exset::const_iterator s = syms.begin();
- ++s;
+ auto s = syms_wox.begin();
for ( size_t i=0; i<a.size(); ++i ) {
do {
a[i] = mod(numeric(rand()), 2*modulus) - modulus;
}
if ( !is_a<numeric>(vn) ) {
// ... and for which the evaluated factors have each an unique prime factor
- lst fnum = f;
- fnum.let_op(0) = fnum.op(0) * u0.content(x);
- for ( size_t i=1; i<fnum.nops(); ++i ) {
- if ( !is_a<numeric>(fnum.op(i)) ) {
- s = syms.begin();
- ++s;
+ exvector fnum = f;
+ fnum[0] = fnum[0] * u0.content(x);
+ for ( size_t i=1; i<fnum.size(); ++i ) {
+ if ( !is_a<numeric>(fnum[i]) ) {
+ s = syms_wox.begin();
for ( size_t j=0; j<a.size(); ++j, ++s ) {
- fnum.let_op(i) = fnum.op(i).subs(*s == a[j]);
+ fnum[i] = fnum[i].subs(*s == a[j]);
}
}
}
// forward declaration
static ex factor_sqrfree(const ex& poly);
-/** Multivariate factorization.
- *
- * The implementation is based on the algorithm described in [Wan].
- * An evaluation homomorphism (a set of integers) is determined that fulfills
- * certain criteria. The evaluated polynomial is univariate and is factorized
- * by factor_univariate(). The main work then is to find the correct leading
- * coefficients of the univariate factors. They have to correspond to the
- * factors of the (multivariate) leading coefficient of the input polynomial
- * (as defined for a specific variable x). After that the Hensel lifting can be
- * performed.
- *
- * @param[in] poly expanded, square free polynomial
- * @param[in] syms contains the symbols in the polynomial
- * @return factorized polynomial
+/** Used by factor_multivariate().
*/
-static ex factor_multivariate(const ex& poly, const exset& syms)
-{
- const ex& x = *syms.begin();
-
- // make polynomial primitive
- ex unit, cont, pp;
- poly.unitcontprim(x, unit, cont, pp);
- if ( !is_a<numeric>(cont) ) {
- return unit * factor_sqrfree(cont) * factor_sqrfree(pp);
- }
-
- // factor leading coefficient
- ex vn = pp.collect(x).lcoeff(x);
- ex vnlst;
- if ( is_a<numeric>(vn) ) {
- vnlst = lst{vn};
- }
- else {
- ex vnfactors = factor(vn);
- vnlst = put_factors_into_lst(vnfactors);
- }
-
- const unsigned int maxtrials = 3;
- numeric modulus = (vnlst.nops() > 3) ? vnlst.nops() : 3;
- vector<numeric> a(syms.size()-1, 0);
-
- // try now to factorize until we are successful
- while ( true ) {
+struct factorization_ctx {
+ const ex poly, x; // polynomial, first symbol x...
+ const exset syms_wox; // ...remaining symbols w/o x
+ ex unit, cont, pp; // unit * cont * pp == poly
+ ex vn; exvector vnlst; // leading coeff, factors of leading coeff
+ numeric modulus; // incremented each time we try
+ /** returns factors or empty if it did not succeed */
+ ex try_next_evaluation_homomorphism()
+ {
+ constexpr unsigned maxtrials = 3;
+ vector<numeric> a(syms_wox.size(), 0);
unsigned int trialcount = 0;
unsigned int prime;
int factor_count = 0;
int min_factor_count = -1;
ex u, delta;
- ex ufac, ufaclst;
+ ex ufac;
+ exvector ufaclst;
// try several evaluation points to reduce the number of factors
while ( trialcount < maxtrials ) {
// generate a set of valid evaluation points
- generate_set(pp, vn, syms, ex_to<lst>(vnlst), modulus, u, a);
+ generate_set(pp, vn, x, syms_wox, vnlst, modulus, u, a);
ufac = factor_univariate(u, x, prime);
- ufaclst = put_factors_into_lst(ufac);
- factor_count = ufaclst.nops()-1;
- delta = ufaclst.op(0);
+ ufaclst = put_factors_into_vec(ufac);
+ factor_count = ufaclst.size()-1;
+ delta = ufaclst[0];
if ( factor_count <= 1 ) {
// irreducible
- return poly;
+ return lst{pp};
}
if ( min_factor_count < 0 ) {
// first time here
vector<ex> C(factor_count);
if ( is_a<numeric>(vn) ) {
// easy case
- for ( size_t i=1; i<ufaclst.nops(); ++i ) {
- C[i-1] = ufaclst.op(i).lcoeff(x);
+ for ( size_t i=1; i<ufaclst.size(); ++i ) {
+ C[i-1] = ufaclst[i].lcoeff(x);
}
} else {
// difficult case.
// we use the property of the ftilde having a unique prime factor.
// details can be found in [Wan].
// calculate ftilde
- vector<numeric> ftilde(vnlst.nops()-1);
+ vector<numeric> ftilde(vnlst.size()-1);
for ( size_t i=0; i<ftilde.size(); ++i ) {
- ex ft = vnlst.op(i+1);
- auto s = syms.begin();
- ++s;
+ ex ft = vnlst[i+1];
+ auto s = syms_wox.begin();
for ( size_t j=0; j<a.size(); ++j ) {
ft = ft.subs(*s == a[j]);
++s;
vector<ex> D(factor_count, 1);
if ( delta == 1 ) {
for ( int i=0; i<factor_count; ++i ) {
- numeric prefac = ex_to<numeric>(ufaclst.op(i+1).lcoeff(x));
+ numeric prefac = ex_to<numeric>(ufaclst[i+1].lcoeff(x));
for ( int j=ftilde.size()-1; j>=0; --j ) {
int count = 0;
- while ( irem(prefac, ftilde[j]) == 0 ) {
- prefac = iquo(prefac, ftilde[j]);
+ numeric q;
+ while ( irem(prefac, ftilde[j], q) == 0 ) {
+ prefac = q;
++count;
}
if ( count ) {
used_flag[j] = true;
- D[i] = D[i] * pow(vnlst.op(j+1), count);
+ D[i] = D[i] * pow(vnlst[j+1], count);
}
}
C[i] = D[i] * prefac;
}
} else {
for ( int i=0; i<factor_count; ++i ) {
- numeric prefac = ex_to<numeric>(ufaclst.op(i+1).lcoeff(x));
+ numeric prefac = ex_to<numeric>(ufaclst[i+1].lcoeff(x));
for ( int j=ftilde.size()-1; j>=0; --j ) {
int count = 0;
- while ( irem(prefac, ftilde[j]) == 0 ) {
- prefac = iquo(prefac, ftilde[j]);
+ numeric q;
+ while ( irem(prefac, ftilde[j], q) == 0 ) {
+ prefac = q;
++count;
}
while ( irem(ex_to<numeric>(delta)*prefac, ftilde[j]) == 0 ) {
numeric g = gcd(prefac, ex_to<numeric>(ftilde[j]));
prefac = iquo(prefac, g);
delta = delta / (ftilde[j]/g);
- ufaclst.let_op(i+1) = ufaclst.op(i+1) * (ftilde[j]/g);
+ ufaclst[i+1] = ufaclst[i+1] * (ftilde[j]/g);
++count;
}
if ( count ) {
used_flag[j] = true;
- D[i] = D[i] * pow(vnlst.op(j+1), count);
+ D[i] = D[i] * pow(vnlst[j+1], count);
}
}
C[i] = D[i] * prefac;
}
}
if ( some_factor_unused ) {
- continue;
+ return lst{}; // next try
}
}
-
+
// multiply the remaining content of the univariate polynomial into the
// first factor
if ( delta != 1 ) {
C[0] = C[0] * delta;
- ufaclst.let_op(1) = ufaclst.op(1) * delta;
+ ufaclst[1] = ufaclst[1] * delta;
}
// set up evaluation points
- EvalPoint ep;
vector<EvalPoint> epv;
- auto s = syms.begin();
- ++s;
+ auto s = syms_wox.begin();
for ( size_t i=0; i<a.size(); ++i ) {
- ep.x = *s++;
- ep.evalpoint = a[i].to_int();
- epv.push_back(ep);
+ epv.emplace_back(EvalPoint{*s++, a[i].to_int()});
}
// calc bound p^l
int maxdeg = 0;
for ( int i=1; i<=factor_count; ++i ) {
- if ( ufaclst.op(i).degree(x) > maxdeg ) {
+ if ( ufaclst[i].degree(x) > maxdeg ) {
maxdeg = ufaclst[i].degree(x);
}
}
- cl_I B = 2*calc_bound(u, x, maxdeg);
+ cl_I B = ash(calc_bound(u, x), maxdeg+1); // = 2 * calc_bound(u,x) * 2^maxdeg
cl_I l = 1;
cl_I pl = prime;
while ( pl < B ) {
l = l + 1;
pl = pl * prime;
}
-
+
// set up modular factors (mod p^l)
- cl_modint_ring R = find_modint_ring(expt_pos(cl_I(prime),l));
- upvec modfactors(ufaclst.nops()-1);
- for ( size_t i=1; i<ufaclst.nops(); ++i ) {
- umodpoly_from_ex(modfactors[i-1], ufaclst.op(i), x, R);
+ cl_modint_ring R = find_modint_ring(pl);
+ upvec modfactors(ufaclst.size()-1);
+ for ( size_t i=1; i<ufaclst.size(); ++i ) {
+ umodpoly_from_ex(modfactors[i-1], ufaclst[i], x, R);
}
// try Hensel lifting
- ex res = hensel_multivar(pp, x, epv, prime, l, modfactors, C);
+ return hensel_multivar(pp, x, epv, prime, l, modfactors, C);
+ }
+};
+
+/** Multivariate factorization.
+ *
+ * The implementation is based on the algorithm described in [Wan].
+ * An evaluation homomorphism (a set of integers) is determined that fulfills
+ * certain criteria. The evaluated polynomial is univariate and is factorized
+ * by factor_univariate(). The main work then is to find the correct leading
+ * coefficients of the univariate factors. They have to correspond to the
+ * factors of the (multivariate) leading coefficient of the input polynomial
+ * (as defined for a specific variable x). After that the Hensel lifting can be
+ * performed. This is done in round-robin for each x in syms until success.
+ *
+ * @param[in] poly expanded, square free polynomial
+ * @param[in] syms contains the symbols in the polynomial
+ * @return factorized polynomial
+ */
+static ex factor_multivariate(const ex& poly, const exset& syms)
+{
+ // set up one factorization context for each symbol
+ vector<factorization_ctx> ctx_in_x;
+ for (auto x : syms) {
+ exset syms_wox; // remaining syms w/o x
+ copy_if(syms.begin(), syms.end(),
+ inserter(syms_wox, syms_wox.end()), [x](const ex& y){ return y != x; });
+
+ factorization_ctx ctx{poly, x, syms_wox};
+
+ // make polynomial primitive
+ poly.unitcontprim(x, ctx.unit, ctx.cont, ctx.pp);
+ if ( !is_a<numeric>(ctx.cont) ) {
+ // content is a polynomial in one or more of remaining syms, let's start over
+ return ctx.unit * factor_sqrfree(ctx.cont) * factor_sqrfree(ctx.pp);
+ }
+
+ // find factors of leading coefficient
+ ctx.vn = ctx.pp.collect(x).lcoeff(x);
+ ctx.vnlst = put_factors_into_vec(factor(ctx.vn));
+
+ ctx.modulus = (ctx.vnlst.size() > 3) ? ctx.vnlst.size() : numeric(3);
+
+ ctx_in_x.push_back(ctx);
+ }
+
+ // try an evaluation homomorphism for each context in round-robin
+ auto ctx = ctx_in_x.begin();
+ while ( true ) {
+
+ ex res = ctx->try_next_evaluation_homomorphism();
+
if ( res != lst{} ) {
- ex result = cont * unit;
+ // found the factors
+ ex result = ctx->cont * ctx->unit;
for ( size_t i=0; i<res.nops(); ++i ) {
- result *= res.op(i).content(x) * res.op(i).unit(x);
- result *= res.op(i).primpart(x);
+ ex unit, cont, pp;
+ res.op(i).unitcontprim(ctx->x, unit, cont, pp);
+ result *= unit * cont * pp;
}
return result;
}
+
+ // switch context for next symbol
+ if (++ctx == ctx_in_x.end()) {
+ ctx = ctx_in_x.begin();
+ }
}
}
* Polynomial factorization. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* functions or polynomials inside function arguments.
*
* @param[in] poly expression to factorize
- * @param[in] option options to influence the factorization
+ * @param[in] options see GiNaC::factor_options
* @return factorized expression
*/
extern ex factor(const ex& poly, unsigned options = 0);
* somewhat obsolete (most of this can be replaced by exceptions). */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* somewhat obsolete (most of this can be replaced by exceptions). */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of abstract derivatives of functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to abstract derivatives of functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Collection of all flags used through the GiNaC framework. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Please do not modify it directly, edit function.cppy instead!
* function.py options: maxargs=@maxargs@
*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
}
current_serial = serial;
if (opt.evalf_use_exvector_args)
- return ((evalf_funcp_exvector)(opt.evalf_f))(seq);
+ return ((evalf_funcp_exvector)(opt.evalf_f))(eseq);
switch (opt.nparams) {
// the following lines have been generated for max. @maxargs@ parameters
+++ for N in range(1, maxargs + 1):
* This file was generated automatically from function.hppy.
* Please do not modify it directly, edit function.hppy instead!
*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* This include file includes all other public GiNaC headers. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Replacement for map<> using hash tables. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's indices. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's indices. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
/** Construct index with given value and dimension.
*
* @param v Value of index (numeric or symbolic)
- * @param dim Dimension of index space (numeric or symbolic)
- * @return newly constructed index */
+ * @param dim Dimension of index space (numeric or symbolic) */
explicit idx(const ex & v, const ex & dim);
// functions overriding virtual functions from base classes
*
* @param v Value of index (numeric or symbolic)
* @param dim Dimension of index space (numeric or symbolic)
- * @param covariant Make covariant index (default is contravariant)
- * @return newly constructed index */
+ * @param covariant Make covariant index (default is contravariant) */
varidx(const ex & v, const ex & dim, bool covariant = false);
// functions overriding virtual functions from base classes
* @param v Value of index (numeric or symbolic)
* @param dim Dimension of index space (numeric or symbolic)
* @param covariant Make covariant index (default is contravariant)
- * @param dotted Make covariant dotted (default is undotted)
- * @return newly constructed index */
+ * @param dotted Make covariant dotted (default is undotted) */
spinidx(const ex & v, const ex & dim = 2, bool covariant = false, bool dotted = false);
// functions overriding virtual functions from base classes
* Implementation of GiNaC's indexed expressions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's indexed expressions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
public:
/** Construct indexed object with no index.
*
- * @param b Base expression
- * @return newly constructed indexed object */
+ * @param b Base expression */
indexed(const ex & b);
/** Construct indexed object with one index. The index must be of class idx.
*
* @param b Base expression
- * @param i1 The index
- * @return newly constructed indexed object */
+ * @param i1 The index */
indexed(const ex & b, const ex & i1);
/** Construct indexed object with two indices. The indices must be of class idx.
*
* @param b Base expression
* @param i1 First index
- * @param i2 Second index
- * @return newly constructed indexed object */
+ * @param i2 Second index */
indexed(const ex & b, const ex & i1, const ex & i2);
/** Construct indexed object with three indices. The indices must be of class idx.
* @param b Base expression
* @param i1 First index
* @param i2 Second index
- * @param i3 Third index
- * @return newly constructed indexed object */
+ * @param i3 Third index */
indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3);
/** Construct indexed object with four indices. The indices must be of class idx.
* @param i1 First index
* @param i2 Second index
* @param i3 Third index
- * @param i4 Fourth index
- * @return newly constructed indexed object */
+ * @param i4 Fourth index */
indexed(const ex & b, const ex & i1, const ex & i2, const ex & i3, const ex & i4);
/** Construct indexed object with two indices and a specified symmetry. The
* @param b Base expression
* @param symm Symmetry of indices
* @param i1 First index
- * @param i2 Second index
- * @return newly constructed indexed object */
+ * @param i2 Second index */
indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2);
/** Construct indexed object with three indices and a specified symmetry.
* @param symm Symmetry of indices
* @param i1 First index
* @param i2 Second index
- * @param i3 Third index
- * @return newly constructed indexed object */
+ * @param i3 Third index */
indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3);
/** Construct indexed object with four indices and a specified symmetry. The
* @param i1 First index
* @param i2 Second index
* @param i3 Third index
- * @param i4 Fourth index
- * @return newly constructed indexed object */
+ * @param i4 Fourth index */
indexed(const ex & b, const symmetry & symm, const ex & i1, const ex & i2, const ex & i3, const ex & i4);
/** Construct indexed object with a specified vector of indices. The indices
* must be of class idx.
*
* @param b Base expression
- * @param iv Vector of indices
- * @return newly constructed indexed object */
+ * @param iv Vector of indices */
indexed(const ex & b, const exvector & iv);
/** Construct indexed object with a specified vector of indices and
*
* @param b Base expression
* @param symm Symmetry of indices
- * @param iv Vector of indices
- * @return newly constructed indexed object */
+ * @param iv Vector of indices */
indexed(const ex & b, const symmetry & symm, const exvector & iv);
// internal constructors
* Implementation of GiNaC's initially known functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
return Order(x).hold();
}
+static ex Order_power(const ex & x, const ex & e)
+{
+ // Order(x)^e -> Order(x^e) for positive integer e
+ if (is_exactly_a<numeric>(e) && e.info(info_flags::posint))
+ return Order(pow(x, e));
+ // NB: For negative exponents, the above could be wrong.
+ // This is because series() produces Order(x^n) to denote the order where
+ // it gave up. So, Order(x^n) can also be an x^(n+1) term if the x^n term
+ // vanishes. In this situation, 1/Order(x^n) can also be a x^(-n-1) term.
+ // Transforming it to Order(x^-n) would miss that.
+
+ return power(Order(x), e).hold();
+}
+
static ex Order_expl_derivative(const ex & arg, const symbol & s)
{
return Order(arg.diff(s));
expl_derivative_func(Order_expl_derivative).
conjugate_func(Order_conjugate).
real_part_func(Order_real_part).
- imag_part_func(Order_imag_part));
+ imag_part_func(Order_imag_part).
+ power_func(Order_power));
//////////
// Solve linear system
* Interface to GiNaC's initially known functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* some related stuff. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
*/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
for (std::size_t i = 0; i < size; ++i)
x[i] = x[i]/y;
+ // 24.03.2021: this block can be outside the loop over r
+ cln::cl_RA p(2);
+ bool adjustp;
+ do {
+ adjustp = false;
+ for (std::size_t i = 0; i < size; ++i) {
+ // 24.03.2021: replaced (x[i] == cln::cl_RA(1)/p) by (cln::zerop(x[i] - cln::cl_RA(1)/p)
+ // in the case where we compare a float with a rational, CLN behaves differently in the two versions
+ if (cln::zerop(x[i] - cln::cl_RA(1)/p) ) {
+ p = p/2 + cln::cl_RA(3)/2;
+ adjustp = true;
+ continue;
+ }
+ }
+ } while (adjustp);
+ cln::cl_RA q = p/(p-1);
+
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;
- }
- }
- } 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) {
// x -> 1-x
if (has_minus_one) {
map_trafo_H_convert_to_Li filter;
- return filter(H(m, numeric(x)).hold()).evalf();
+ // 09.06.2021: bug fix: don't forget a possible minus sign from the case realpart(x) < 0
+ res *= filter(H(m, numeric(x)).hold()).evalf();
+ return res;
}
map_trafo_H_1mx trafo;
res *= trafo(H(m, xtemp).hold());
* functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
return power(exp(x), a).hold();
}
+static bool exp_info(const ex & x, unsigned inf)
+{
+ switch (inf) {
+ case info_flags::expanded:
+ case info_flags::real:
+ return x.info(inf);
+ case info_flags::positive:
+ case info_flags::nonnegative:
+ return x.info(info_flags::real);
+ default:
+ return false;
+ }
+}
+
REGISTER_FUNCTION(exp, eval_func(exp_eval).
evalf_func(exp_evalf).
+ info_func(exp_info).
expand_func(exp_expand).
derivative_func(exp_deriv).
real_part_func(exp_real_part).
return conjugate_function(log(x)).hold();
}
+static bool log_info(const ex & x, unsigned inf)
+{
+ switch (inf) {
+ case info_flags::expanded:
+ return x.info(inf);
+ case info_flags::real:
+ return x.info(info_flags::positive);
+ default:
+ return false;
+ }
+}
+
REGISTER_FUNCTION(log, eval_func(log_eval).
evalf_func(log_evalf).
+ info_func(log_info).
expand_func(log_expand).
derivative_func(log_deriv).
series_func(log_series).
return sin(x.conjugate());
}
+static bool trig_info(const ex & x, unsigned inf)
+{
+ switch (inf) {
+ case info_flags::expanded:
+ case info_flags::real:
+ return x.info(inf);
+ default:
+ return false;
+ }
+}
+
REGISTER_FUNCTION(sin, eval_func(sin_eval).
evalf_func(sin_evalf).
+ info_func(trig_info).
derivative_func(sin_deriv).
real_part_func(sin_real_part).
imag_part_func(sin_imag_part).
}
REGISTER_FUNCTION(cos, eval_func(cos_eval).
+ info_func(trig_info).
evalf_func(cos_evalf).
derivative_func(cos_deriv).
real_part_func(cos_real_part).
REGISTER_FUNCTION(tan, eval_func(tan_eval).
evalf_func(tan_evalf).
+ info_func(trig_info).
derivative_func(tan_deriv).
series_func(tan_series).
real_part_func(tan_real_part).
return conjugate_function(asin(x)).hold();
}
+static bool asin_info(const ex & x, unsigned inf)
+{
+ switch (inf) {
+ case info_flags::expanded:
+ return x.info(inf);
+ default:
+ return false;
+ }
+}
+
REGISTER_FUNCTION(asin, eval_func(asin_eval).
evalf_func(asin_evalf).
+ info_func(asin_info).
derivative_func(asin_deriv).
conjugate_func(asin_conjugate).
latex_name("\\arcsin"));
REGISTER_FUNCTION(acos, eval_func(acos_eval).
evalf_func(acos_evalf).
+ info_func(asin_info). // Flags of acos are shared with asin functions
derivative_func(acos_deriv).
conjugate_func(acos_conjugate).
latex_name("\\arccos"));
return conjugate_function(atan(x)).hold();
}
+static bool atan_info(const ex & x, unsigned inf)
+{
+ switch (inf) {
+ case info_flags::expanded:
+ case info_flags::real:
+ return x.info(inf);
+ case info_flags::positive:
+ case info_flags::negative:
+ case info_flags::nonnegative:
+ return x.info(info_flags::real) && x.info(inf);
+ default:
+ return false;
+ }
+}
+
REGISTER_FUNCTION(atan, eval_func(atan_eval).
evalf_func(atan_evalf).
+ info_func(atan_info).
derivative_func(atan_deriv).
series_func(atan_series).
conjugate_func(atan_conjugate).
return -y*power(power(x,_ex2)+power(y,_ex2),_ex_1);
}
+static bool atan2_info(const ex & y, const ex & x, unsigned inf)
+{
+ switch (inf) {
+ case info_flags::expanded:
+ case info_flags::real:
+ return y.info(inf) && x.info(inf);
+ case info_flags::positive:
+ case info_flags::negative:
+ case info_flags::nonnegative:
+ return y.info(info_flags::real) && x.info(info_flags::real)
+ && y.info(inf);
+ default:
+ return false;
+ }
+}
+
REGISTER_FUNCTION(atan2, eval_func(atan2_eval).
+ evalf_func(atan2_evalf).
+ info_func(atan2_info).
evalf_func(atan2_evalf).
derivative_func(atan2_deriv));
REGISTER_FUNCTION(sinh, eval_func(sinh_eval).
evalf_func(sinh_evalf).
+ info_func(atan_info). // Flags of sinh are shared with atan functions
derivative_func(sinh_deriv).
real_part_func(sinh_real_part).
imag_part_func(sinh_imag_part).
REGISTER_FUNCTION(cosh, eval_func(cosh_eval).
evalf_func(cosh_evalf).
+ info_func(exp_info). // Flags of cosh are shared with exp functions
derivative_func(cosh_deriv).
real_part_func(cosh_real_part).
imag_part_func(cosh_imag_part).
REGISTER_FUNCTION(tanh, eval_func(tanh_eval).
evalf_func(tanh_evalf).
+ info_func(atan_info). // Flags of tanh are shared with atan functions
derivative_func(tanh_deriv).
series_func(tanh_series).
real_part_func(tanh_real_part).
REGISTER_FUNCTION(asinh, eval_func(asinh_eval).
evalf_func(asinh_evalf).
+ info_func(atan_info). // Flags of asinh are shared with atan functions
derivative_func(asinh_deriv).
conjugate_func(asinh_conjugate));
REGISTER_FUNCTION(acosh, eval_func(acosh_eval).
evalf_func(acosh_evalf).
+ info_func(asin_info). // Flags of acosh are shared with asin functions
derivative_func(acosh_deriv).
conjugate_func(acosh_conjugate));
REGISTER_FUNCTION(atanh, eval_func(atanh_eval).
evalf_func(atanh_evalf).
+ info_func(asin_info). // Flags of atanh are shared with asin functions
derivative_func(atanh_deriv).
series_func(atanh_series).
conjugate_func(atanh_conjugate));
* Implementation of GiNaC's symbolic integral. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's symbolic integral. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's integration kernels for iterated integrals. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
subs_q_expansion do_subs_q_expansion(qbar, order);
ex res = do_subs_q_expansion(P).series(qbar,order);
+ res += Order(pow(qbar,order));
+ res = res.series(qbar,order);
return res;
}
symbol qbar("qbar");
// test with a random number and random expansion
- return series_to_poly(P.series(qbar,18)).subs(qbar==numeric(1,937)).evalf().info(info_flags::numeric);
+ return series_to_poly(q_expansion_modular_form(qbar,18)).subs(qbar==numeric(1,937)).evalf().info(info_flags::numeric);
}
ex modular_form_kernel::Laurent_series(const ex & qbar, int order) const
{
- ex res = series_to_poly(P.series(qbar,order+1));
+ ex res = series_to_poly(q_expansion_modular_form(qbar,order+1));
res = C_norm * res/qbar;
res = res.series(qbar,order);
return res;
* Interface to GiNaC's integration kernels for iterated integrals. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
#include "basic.h"
#include "archive.h"
#include "numeric.h"
+#include "lst.h"
#include <cln/complex.h>
#include <vector>
* Implementation of GiNaC's lst. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
return inherited::info(inf);
}
+template bool container<std::list>::info(unsigned) const;
GINAC_BIND_UNARCHIVER(lst);
} // namespace GiNaC
* Definition of GiNaC's lst. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
typedef container<std::list> lst;
+/** Declaration of container::reg_info for lst. */
+#ifndef _MSC_VER // workaround error C2766: explicit specialization; 'reg_info' has already been defined
+template<> registered_class_info lst::reg_info;
+#endif
+
/** Specialization of container::get_default_flags() for lst. */
template<> inline unsigned lst::get_default_flags() { return status_flags::not_shareable; }
// defined in lst.cpp
template<> bool lst::info(unsigned inf) const;
+extern template bool container<std::list>::info(unsigned) const;
GINAC_DECLARE_UNARCHIVER(lst);
} // namespace GiNaC
* Implementation of symbolic matrices */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to symbolic matrices */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's products of expressions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's products of expressions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's non-commutative products of expressions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's non-commutative products of expressions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* computation, square-free factorization and rational function normalization. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* @param cb pointer to expression that will receive the cofactor of b, or nullptr
* @param check_args check whether a and b are polynomials with rational
* coefficients (defaults to "true")
+ * @param options see GiNaC::gcd_options
* @return the GCD as a new expression */
ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned options)
{
// Factorize denominator and compute cofactors
epvector yun = sqrfree_yun(denom, x);
- size_t yun_max_exponent = yun.empty() ? 0 : ex_to<numeric>(yun.back().coeff).to_int();
exvector factor, cofac;
+ size_t dim = 0;
for (size_t i=0; i<yun.size(); i++) {
numeric i_exponent = ex_to<numeric>(yun[i].coeff);
for (size_t j=0; j<i_exponent; j++) {
factor.push_back(pow(yun[i].rest, j+1));
+ dim += degree(yun[i].rest, x);
ex prod = _ex1;
for (size_t k=0; k<yun.size(); k++) {
if (yun[k].coeff == i_exponent)
cofac.push_back(prod.expand());
}
}
- size_t num_factors = factor.size();
//std::clog << "factors : " << exprseq(factor) << std::endl;
//std::clog << "cofactors: " << exprseq(cofac) << std::endl;
- // Construct coefficient matrix for decomposition
- int max_denom_deg = denom.degree(x);
- matrix sys(max_denom_deg + 1, num_factors);
- matrix rhs(max_denom_deg + 1, 1);
- for (int i=0; i<=max_denom_deg; i++) {
- for (size_t j=0; j<num_factors; j++)
- sys(i, j) = cofac[j].coeff(x, i);
- rhs(i, 0) = red_numer.coeff(x, i);
+ // Construct linear system for decomposition
+ matrix sys(dim, dim);
+ matrix rhs(dim, 1);
+ matrix vars(dim, 1);
+ for (size_t i=0, n=0, f=0; i<yun.size(); i++) {
+ size_t i_expo = to_int(ex_to<numeric>(yun[i].coeff));
+ for (size_t j=0; j<i_expo; j++) {
+ for (size_t k=0; k<size_t(degree(yun[i].rest, x)); k++) {
+ GINAC_ASSERT(n < dim && f < factor.size());
+
+ // column n of coefficient matrix
+ for (size_t r=0; r+k<dim; r++) {
+ sys(r+k, n) = cofac[f].coeff(x, r);
+ }
+
+ // element n of right hand side vector
+ rhs(n, 0) = red_numer.coeff(x, n);
+
+ // element n of free variables vector
+ vars(n, 0) = symbol();
+
+ n++;
+ }
+ f++;
+ }
}
//std::clog << "coeffs: " << sys << std::endl;
//std::clog << "rhs : " << rhs << std::endl;
- // Solve resulting linear system
- matrix vars(num_factors, 1);
- for (size_t i=0; i<num_factors; i++)
- vars(i, 0) = symbol();
+ // Solve resulting linear system and sum up decomposed fractions
matrix sol = sys.solve(vars, rhs);
+//std::clog << "sol : " << sol << std::endl;
+ ex sum = red_poly;
+ for (size_t i=0, n=0, f=0; i<yun.size(); i++) {
+ size_t i_expo = to_int(ex_to<numeric>(yun[i].coeff));
+ for (size_t j=0; j<i_expo; j++) {
+ ex frac_numer = 0;
+ for (size_t k=0; k<size_t(degree(yun[i].rest, x)); k++) {
+ GINAC_ASSERT(n < dim && f < factor.size());
+ frac_numer += sol(n, 0) * pow(x, k);
+ n++;
+ }
+ sum += frac_numer / factor[f];
- // Sum up decomposed fractions
- ex sum = 0;
- for (size_t i=0; i<num_factors; i++)
- sum += sol(i, 0) / factor[i];
+ f++;
+ }
+ }
- return red_poly + sum;
+ return sum;
}
// they can be rationalised more efficiently
if (is_a<function>(e_replaced) && is_ex_the_function(e_replaced, exp)) {
for (auto & it : repl) {
- if (is_a<function>(it.second) && is_ex_the_function(e_replaced, exp)) {
+ if (is_a<function>(it.second) && is_ex_the_function(it.second, exp)) {
ex ratio = normal(e_replaced.op(0) / it.second.op(0));
if (is_a<numeric>(ratio) && ex_to<numeric>(ratio).is_rational()) {
// Different exponents can be treated as powers of the same basic equation
return dynallocate<lst>({replace_with_symbol(*this, repl, rev_lookup, modifier), _ex1});
normal_map_function map_normal;
- int nmod = modifier.nops(); // To watch new modifiers to the replacement list
- lst result = dynallocate<lst>({replace_with_symbol(map(map_normal), repl, rev_lookup, modifier), _ex1});
- for (int imod = nmod; imod < modifier.nops(); ++imod) {
+ size_t nmod = modifier.nops(); // To watch new modifiers to the replacement list
+ ex result = replace_with_symbol(map(map_normal), repl, rev_lookup, modifier);
+ for (size_t imod = nmod; imod < modifier.nops(); ++imod) {
exmap this_repl;
this_repl.insert(std::make_pair(modifier.op(imod).op(0), modifier.op(imod).op(1)));
- result = ex_to<lst>(result.subs(this_repl, subs_options::no_pattern));
+ result = result.subs(this_repl, subs_options::no_pattern);
}
- return result;
+ // Sometimes we may obtain negative powers, they need to be placed to denominator
+ if (is_a<power>(result) && result.op(1).info(info_flags::negative))
+ return dynallocate<lst>({_ex1, power(result.op(0), -result.op(1))});
+ else
+ return dynallocate<lst>({result, _ex1});
}
exvector nums, dens;
nums.reserve(seq.size()+1);
dens.reserve(seq.size()+1);
- int nmod = modifier.nops(); // To watch new modifiers to the replacement list
+ size_t nmod = modifier.nops(); // To watch new modifiers to the replacement list
for (auto & it : seq) {
ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, modifier);
nums.push_back(n.op(0));
auto num_it = nums.begin(), num_itend = nums.end();
auto den_it = dens.begin(), den_itend = dens.end();
//std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl;
- for (int imod = nmod; imod < modifier.nops(); ++imod) {
+ for (size_t imod = nmod; imod < modifier.nops(); ++imod) {
while (num_it != num_itend) {
*num_it = num_it->subs(modifier.op(imod), subs_options::no_pattern);
++num_it;
exvector num; num.reserve(seq.size());
exvector den; den.reserve(seq.size());
ex n;
- int nmod = modifier.nops(); // To watch new modifiers to the replacement list
+ size_t nmod = modifier.nops(); // To watch new modifiers to the replacement list
for (auto & it : seq) {
n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, modifier);
num.push_back(n.op(0));
num.push_back(n.op(0));
den.push_back(n.op(1));
auto num_it = num.begin(), num_itend = num.end();
- auto den_it = den.begin(), den_itend = den.end();
- for (int imod = nmod; imod < modifier.nops(); ++imod) {
+ auto den_it = den.begin();
+ for (size_t imod = nmod; imod < modifier.nops(); ++imod) {
while (num_it != num_itend) {
*num_it = num_it->subs(modifier.op(imod), subs_options::no_pattern);
++num_it;
ex power::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
{
// Normalize basis and exponent (exponent gets reassembled)
- int nmod = modifier.nops(); // To watch new modifiers to the replacement list
+ size_t nmod = modifier.nops(); // To watch new modifiers to the replacement list
ex n_basis = ex_to<basic>(basis).normal(repl, rev_lookup, modifier);
- for (int imod = nmod; imod < modifier.nops(); ++imod)
+ for (size_t imod = nmod; imod < modifier.nops(); ++imod)
n_basis = n_basis.subs(modifier.op(imod), subs_options::no_pattern);
nmod = modifier.nops();
ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup, modifier);
- for (int imod = nmod; imod < modifier.nops(); ++imod)
+ for (size_t imod = nmod; imod < modifier.nops(); ++imod)
n_exponent = n_exponent.subs(modifier.op(imod), subs_options::no_pattern);
n_exponent = n_exponent.op(0) / n_exponent.op(1);
// Re-insert replaced symbols
if (!repl.empty()) {
- for(int i=0; i < modifier.nops(); ++i)
+ for(size_t i=0; i < modifier.nops(); ++i)
e = e.subs(modifier.op(i), subs_options::no_pattern);
e = e.subs(repl, subs_options::no_pattern);
}
if (repl.empty())
return e.op(0);
else {
- for(int i=0; i < modifier.nops(); ++i)
+ for(size_t i=0; i < modifier.nops(); ++i)
e = e.subs(modifier.op(i), subs_options::no_pattern);
return e.op(0).subs(repl, subs_options::no_pattern);
if (repl.empty())
return e.op(1);
else {
- for(int i=0; i < modifier.nops(); ++i)
+ for(size_t i=0; i < modifier.nops(); ++i)
e = e.subs(modifier.op(i), subs_options::no_pattern);
return e.op(1).subs(repl, subs_options::no_pattern);
if (repl.empty())
return e;
else {
- for(int i=0; i < modifier.nops(); ++i)
+ for(size_t i=0; i < modifier.nops(); ++i)
e = e.subs(modifier.op(i), subs_options::no_pattern);
return e.subs(repl, subs_options::no_pattern);
ex power::to_polynomial(exmap & repl) const
{
if (exponent.info(info_flags::posint))
- return pow(basis.to_rational(repl), exponent);
+ return pow(basis.to_polynomial(repl), exponent);
else if (exponent.info(info_flags::negint))
{
ex basis_pref = collect_common_factors(basis);
* computation, square-free factorization and rational function normalization. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* of special functions or implement the interface to the bignum package. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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 Binomial coefficients. It computes the binomial coefficients. For
* integer n and k and positive n this is the number of ways of choosing k
- * objects from n distinct objects. If n is negative, the formula
- * binomial(n,k) == (-1)^k*binomial(k-n-1,k) is used to compute the result. */
+ * objects from n distinct objects. If n is a negative integer, the formula
+ * binomial(n,k) == (-1)^k*binomial(k-n-1,k) (if k>=0)
+ * binomial(n,k) == (-1)^(n-k)*binomial(-k-1,n-k) (otherwise)
+ * is used to compute the result. */
const numeric binomial(const numeric &n, const numeric &k)
{
if (n.is_integer() && k.is_integer()) {
else
return *_num0_p;
} else {
- return _num_1_p->power(k)*binomial(k-n-(*_num1_p),k);
+ if (k.is_nonneg_integer())
+ return _num_1_p->power(k)*binomial(k-n-(*_num1_p), k);
+ else
+ return _num_1_p->power(n-k)*binomial(-k-(*_num1_p), n-k);
}
}
* Makes the interface to the underlying bignum package available. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's overloaded operators. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's overloaded operators. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Debugging facilities for parser. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
**/
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's lexer. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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 lexer. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Code to deal with binary operators. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
}
}
+/// unary_expr: [+-] expression
+ex parser::parse_unary_expr()
+{
+ // Parse a binary expression with the priority of exponentiation
+ // or higher. Ignore the overall sign, because parse_primary()
+ // handles it for us.
+ get_next_tok(); // Skip [+-]
+ ex lhs = parse_primary();
+ ex e = parse_binop_rhs(get_tok_prec('^'), lhs);
+ return e;
+}
+
extern const numeric* _num_1_p;
static ex make_minus_expr(const exvector& args)
return dynallocate<mul>(args[0], rest);
}
+static ex make_power_expr(const exvector& args)
+{
+ size_t n = args.size();
+ ex p = pow(args[n - 2], args[n - 1]);
+ for (size_t i = n - 2; i > 0; i--) {
+ p = pow(args[i - 1], p);
+ }
+ return p;
+}
+
static ex make_binop_expr(const int binop, const exvector& args)
{
switch (binop) {
case '/':
return make_divide_expr(args);
case '^':
- if (args.size() != 2)
- throw std::invalid_argument(
- std::string(__func__)
- + ": power should have exactly 2 operands");
- return pow(args[0], args[1]);
+ return make_power_expr(args);
default:
throw std::invalid_argument(
std::string(__func__)
* Implementation of the parser context. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to parser context. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's parser. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
#include "mul.h"
#include "constant.h"
#include "function.h"
+#include "operators.h"
#include <cstdint> // for uintptr_t
#include <sstream>
return list;
}
-extern const ex _ex0;
-
-/// unary_expr: [+-] expression
-ex parser::parse_unary_expr()
-{
- // 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)
- // Delegate 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
ex parser::parse_primary()
{
case '{':
return parse_lst_expr();
case '-':
+ return -parse_unary_expr();
case '+':
return parse_unary_expr();
case lexer::token_type::literal:
* Interface to the parser. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Parser interface compatible with the old (bison/flex based) parser. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Chinese remainder algorithm. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GCD functions using Chinese remainder algorithm. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Utility functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to utility functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Garner's algorithm. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to Garner's algorithm. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Utility macros and functions for debugging. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of polynomial division in Z/Zp. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to polynomial division in Z/Zp. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Euclidean GCD and supporting functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Functions for finding an evaluation point. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Numerical evaluation of univariate polynomials. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* GCD using Euclidean algorithm. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Several GCD algorithms. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Heuristic GCD code. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Utility function for interpolation. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Chinese remainder algorithm. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of modular GCD. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to modular GCD. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Functions for Newton interpolation. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Functions to normalize polynomials in a field. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Functions to normalize polynomials in a field. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Functions to optimize the choice of variable for multivariate GCD. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* computations. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* GCD for polynomials in prime field. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GCD functions for polynomials over prime fields. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Chinese remainder algorithm. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Function to calculate the pseudo-remainder. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Factory for prime numbers. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Function to find primitive part of a multivariate polynomial. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Functions calculating remainders. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Functions calculating remainders. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Functions for polynomial ring arithmetic. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Utility functions for modular arithmetic. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* GCD function for univariate polynomials using PRS method. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to polynomials with integer and modular coefficients. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Input/Output function for univariate polynomials. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Output operators for polynomials. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's symbolic exponentiation (basis^exponent). */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
static void print_sym_pow(const print_context & c, const symbol &x, int exp)
{
// Optimal output of integer powers of symbols to aid compiler CSE.
- // C.f. ISO/IEC 14882:2011, section 1.9 [intro execution], paragraph 15
+ // C.f. ISO/IEC 14882:2011, section 1.9 [intro.execution], paragraph 15
// to learn why such a parenthesation is really necessary.
if (exp == 1) {
x.print(c);
if (!are_ex_trivially_equal(basis, subsed_basis)
|| !are_ex_trivially_equal(exponent, subsed_exponent))
- return power(subsed_basis, subsed_exponent).subs_one_level(m, options);
+ return dynallocate<power>(subsed_basis, subsed_exponent);
if (!(options & subs_options::algebraic))
return subs_one_level(m, options);
* Interface to GiNaC's symbolic exponentiation (basis^exponent). */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of helper classes for expression output. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Definition of helper classes for expression output. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* methods for series expansion. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* non-terminating series.
*
* @param rel_ expansion variable and point (must hold a relational)
- * @param ops_ vector of {coefficient, power} pairs (coefficient must not be zero)
- * @return newly constructed pseries */
+ * @param ops_ vector of {coefficient, power} pairs (coefficient must not be zero) */
pseries::pseries(const ex &rel_, const epvector &ops_)
: seq(ops_)
{
// which can easily be solved given the starting value c_0 = (a_0)^p.
// For the more general case where the leading coefficient of A(x) is not
// a constant, just consider A2(x) = A(x)*x^m, with some integer m and
- // repeat the above derivation. The leading power of C2(x) = A2(x)^2 is
- // then of course x^(p*m) but the recurrence formula still holds.
+ // repeat the above derivation. The leading power of C2(x) = A2(x)^p is
+ // then of course a_0^p*x^(p*m) but the recurrence formula still holds.
if (seq.empty()) {
// as a special case, handle the empty (zero) series honoring the
else
return *this;
}
-
- const int ldeg = ldegree(var);
- if (!(p*ldeg).is_integer())
+
+ const int base_ldeg = ldegree(var);
+ if (!(p*base_ldeg).is_integer())
throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
+ int new_ldeg = (p*base_ldeg).to_int();
+
+ const int base_deg = degree(var);
+ int new_deg = deg;
+ if (p.is_pos_integer()) {
+ // No need to compute beyond p*base_deg.
+ new_deg = std::min((p*base_deg).to_int(), deg);
+ }
// adjust number of coefficients
- int numcoeff = deg - (p*ldeg).to_int();
+ int numcoeff = new_deg - new_ldeg;
+ if (new_deg < deg)
+ numcoeff += 1;
+
if (numcoeff <= 0) {
- epvector epv { expair(Order(_ex1), deg) };
- return dynallocate<pseries>(relational(var,point), std::move(epv));
+ return dynallocate<pseries>(relational(var, point),
+ epvector{{Order(_ex1), deg}});
}
// O(x^n)^(-m) is undefined
// Compute coefficients of the powered series
exvector co;
co.reserve(numcoeff);
- co.push_back(pow(coeff(var, ldeg), p));
+ co.push_back(pow(coeff(var, base_ldeg), p));
for (int i=1; i<numcoeff; ++i) {
ex sum = _ex0;
for (int j=1; j<=i; ++j) {
- ex c = coeff(var, j + ldeg);
+ ex c = coeff(var, j + base_ldeg);
if (is_order_function(c)) {
co.push_back(Order(_ex1));
break;
} else
sum += (p * j - (i - j)) * co[i - j] * c;
}
- co.push_back(sum / coeff(var, ldeg) / i);
+ co.push_back(sum / coeff(var, base_ldeg) / i);
}
// Construct new series (of non-zero coefficients)
epvector new_seq;
bool higher_order = false;
for (int i=0; i<numcoeff; ++i) {
- if (!co[i].is_zero())
- new_seq.emplace_back(expair(co[i], p * ldeg + i));
+ if (!co[i].is_zero()) {
+ new_seq.emplace_back(expair(co[i], new_ldeg + i));
+ }
if (is_order_function(co[i])) {
higher_order = true;
break;
}
}
- if (!higher_order)
- new_seq.emplace_back(expair(Order(_ex1), p * ldeg + numcoeff));
+ if (!higher_order && new_deg == deg) {
+ new_seq.emplace_back(expair{Order(_ex1), new_deg});
+ }
return pseries(relational(var,point), std::move(new_seq));
}
if (!(real_ldegree*numexp).is_integer())
throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
- ex e = basis.series(r, (order + real_ldegree*(1-numexp)).to_int(), options);
+ int extra_terms = (real_ldegree*(1-numexp)).to_int();
+ ex e = basis.series(r, order + std::max(0, extra_terms), options);
ex result;
try {
* Interface to class for extended truncated power series. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Reference-counted pointer template. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* GiNaC's class registrar (for class basic and all classes derived from it). */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* GiNaC's class registrar (for class basic and all classes derived from it). */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of relations between expressions */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
return (o < oth.o) ? -1 : 1;
break;
}
- const int lcmpval = lh.compare(oth.rh);
- return (lcmpval!=0) ? lcmpval : rh.compare(oth.lh);
+ const int lcmpval = lh.compare(oth.lh);
+ return (lcmpval!=0) ? lcmpval : rh.compare(oth.rh);
}
bool relational::match_same_type(const basic & other) const
* unequal or undecidable). */
relational::operator relational::safe_bool() const
{
- const ex df = lh-rh;
- if (!is_exactly_a<numeric>(df))
- // cannot decide on non-numerical results
- return o==not_equal ? make_safe_bool(true) : make_safe_bool(false);
-
- switch (o) {
- case equal:
- return make_safe_bool(ex_to<numeric>(df).is_zero());
- case not_equal:
- return make_safe_bool(!ex_to<numeric>(df).is_zero());
- case less:
- return make_safe_bool(ex_to<numeric>(df)<(*_num0_p));
- case less_or_equal:
- return make_safe_bool(ex_to<numeric>(df)<=(*_num0_p));
- case greater:
- return make_safe_bool(ex_to<numeric>(df)>(*_num0_p));
- case greater_or_equal:
- return make_safe_bool(ex_to<numeric>(df)>=(*_num0_p));
- default:
- throw(std::logic_error("invalid relational operator"));
+ const ex df = lh-rh; // like ::canonical() method
+ // We treat numeric and symbolic expression differently
+ if (is_exactly_a<numeric>(df)) {
+ switch (o) {
+ case equal:
+ return make_safe_bool(ex_to<numeric>(df).is_zero());
+ case not_equal:
+ return make_safe_bool(!ex_to<numeric>(df).is_zero());
+ case less:
+ return make_safe_bool(ex_to<numeric>(df)<(*_num0_p));
+ case less_or_equal:
+ return make_safe_bool(ex_to<numeric>(df)<=(*_num0_p));
+ case greater:
+ return make_safe_bool(ex_to<numeric>(df)>(*_num0_p));
+ case greater_or_equal:
+ return make_safe_bool(ex_to<numeric>(df)>=(*_num0_p));
+ default:
+ throw(std::logic_error("invalid relational operator"));
+ }
+ } else {
+ // The conversion for symbolic expressions is based on the info flags
+ switch (o) {
+ case equal:
+ return make_safe_bool(df.is_zero());
+ case not_equal:
+ return make_safe_bool(! df.is_zero());
+ case less:
+ return make_safe_bool(df.info(info_flags::negative));
+ case less_or_equal:
+ return make_safe_bool((-df).info(info_flags::nonnegative));
+ case greater:
+ return make_safe_bool(df.info(info_flags::positive));
+ case greater_or_equal:
+ return make_safe_bool(df.info(info_flags::nonnegative));
+ default:
+ throw(std::logic_error("invalid relational operator"));
+ }
}
}
+/** Returns an equivalent relational with zero right-hand side.
+ */
+ex relational::canonical() const
+{
+ return relational(lh-rh, _ex0, o);
+}
+
} // namespace GiNaC
* Interface to relations between expressions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
void archive(archive_node& n) const override;
/** Read (a.k.a. deserialize) object from archive. */
void read_archive(const archive_node& n, lst& syms) override;
+ ex canonical() const;
+
protected:
ex eval_ncmul(const exvector & v) const override;
bool match_same_type(const basic & other) const override;
* in GiNaC functions */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* in GiNaC functions */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Wrapper template for making GiNaC classes out of C++ structures. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's symbolic objects. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's symbolic objects. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's symmetry definitions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's symmetry definitions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Implementation of GiNaC's special tensors. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's special tensors. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* but not of any interest to the user of the library. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
return factorial(n).div(d);
}
-//////////
-// flyweight chest of numbers is initialized here:
-//////////
/** How many static objects were created? Only the first one must create
* the static flyweights on the heap. */
int library_init::count = 0;
-// static numeric -120
-const numeric *_num_120_p;
-const ex _ex_120 = _ex_120;
-
-// static numeric -60
-const numeric *_num_60_p;
-const ex _ex_60 = _ex_60;
-
-// static numeric -48
-const numeric *_num_48_p;
-const ex _ex_48 = _ex_48;
-
-// static numeric -30
-const numeric *_num_30_p;
-const ex _ex_30 = _ex_30;
-
-// static numeric -25
-const numeric *_num_25_p;
-const ex _ex_25 = _ex_25;
-
-// static numeric -24
-const numeric *_num_24_p;
-const ex _ex_24 = _ex_24;
-
-// static numeric -20
-const numeric *_num_20_p;
-const ex _ex_20 = _ex_20;
-
-// static numeric -18
-const numeric *_num_18_p;
-const ex _ex_18 = _ex_18;
-
-// static numeric -15
-const numeric *_num_15_p;
-const ex _ex_15 = _ex_15;
-
-// static numeric -12
-const numeric *_num_12_p;
-const ex _ex_12 = _ex_12;
-
-// static numeric -11
-const numeric *_num_11_p;
-const ex _ex_11 = _ex_11;
-
-// static numeric -10
-const numeric *_num_10_p;
-const ex _ex_10 = _ex_10;
-
-// static numeric -9
-const numeric *_num_9_p;
-const ex _ex_9 = _ex_9;
-
-// static numeric -8
-const numeric *_num_8_p;
-const ex _ex_8 = _ex_8;
-
-// static numeric -7
-const numeric *_num_7_p;
-const ex _ex_7 = _ex_7;
-
-// static numeric -6
-const numeric *_num_6_p;
-const ex _ex_6 = _ex_6;
-
-// static numeric -5
-const numeric *_num_5_p;
-const ex _ex_5 = _ex_5;
-
-// static numeric -4
-const numeric *_num_4_p;
-const ex _ex_4 = _ex_4;
-
-// static numeric -3
-const numeric *_num_3_p;
-const ex _ex_3 = _ex_3;
-
-// static numeric -2
-const numeric *_num_2_p;
-const ex _ex_2 = _ex_2;
-
-// static numeric -1
-const numeric *_num_1_p;
-const ex _ex_1 = _ex_1;
-
-// static numeric -1/2
-const numeric *_num_1_2_p;
-const ex _ex_1_2= _ex_1_2;
-
-// static numeric -1/3
-const numeric *_num_1_3_p;
-const ex _ex_1_3= _ex_1_3;
-
-// static numeric -1/4
-const numeric *_num_1_4_p;
-const ex _ex_1_4= _ex_1_4;
-
-// static numeric 0
-const numeric *_num0_p;
-const basic *_num0_bp;
-const ex _ex0 = _ex0;
-
-// static numeric 1/4
-const numeric *_num1_4_p;
-const ex _ex1_4 = _ex1_4;
-
-// static numeric 1/3
-const numeric *_num1_3_p;
-const ex _ex1_3 = _ex1_3;
-
-// static numeric 1/2
-const numeric *_num1_2_p;
-const ex _ex1_2 = _ex1_2;
-
-// static numeric 1
-const numeric *_num1_p;
-const ex _ex1 = _ex1;
-
-// static numeric 2
-const numeric *_num2_p;
-const ex _ex2 = _ex2;
-
-// static numeric 3
-const numeric *_num3_p;
-const ex _ex3 = _ex3;
-
-// static numeric 4
-const numeric *_num4_p;
-const ex _ex4 = _ex4;
-
-// static numeric 5
-const numeric *_num5_p;
-const ex _ex5 = _ex5;
-
-// static numeric 6
-const numeric *_num6_p;
-const ex _ex6 = _ex6;
-
-// static numeric 7
-const numeric *_num7_p;
-const ex _ex7 = _ex7;
-
-// static numeric 8
-const numeric *_num8_p;
-const ex _ex8 = _ex8;
-
-// static numeric 9
-const numeric *_num9_p;
-const ex _ex9 = _ex9;
-
-// static numeric 10
-const numeric *_num10_p;
-const ex _ex10 = _ex10;
-
-// static numeric 11
-const numeric *_num11_p;
-const ex _ex11 = _ex11;
-
-// static numeric 12
-const numeric *_num12_p;
-const ex _ex12 = _ex12;
-
-// static numeric 15
-const numeric *_num15_p;
-const ex _ex15 = _ex15;
-
-// static numeric 18
-const numeric *_num18_p;
-const ex _ex18 = _ex18;
-
-// static numeric 20
-const numeric *_num20_p;
-const ex _ex20 = _ex20;
-
-// static numeric 24
-const numeric *_num24_p;
-const ex _ex24 = _ex24;
-
-// static numeric 25
-const numeric *_num25_p;
-const ex _ex25 = _ex25;
-
-// static numeric 30
-const numeric *_num30_p;
-const ex _ex30 = _ex30;
-
-// static numeric 48
-const numeric *_num48_p;
-const ex _ex48 = _ex48;
-
-// static numeric 60
-const numeric *_num60_p;
-const ex _ex60 = _ex60;
-
-// static numeric 120
-const numeric *_num120_p;
-const ex _ex120 = _ex120;
-
/** Ctor of static initialization helpers. The fist call to this is going
* to initialize the library, the others do nothing. */
library_init::library_init()
void library_init::init_unarchivers() { }
+
+//////////
+// Flyweight chest of numbers is re-initialized here. Note that this works
+// because the numeric* have been dynallocated by the library_init ctor before
+// (with the first module that has a static library_init object), so the
+// assignments here only increment their refcounts.
+//////////
+
+// static numeric -120
+const numeric *_num_120_p;
+const ex _ex_120 = ex(*_num_120_p);
+
+// static numeric -60
+const numeric *_num_60_p;
+const ex _ex_60 = ex(*_num_60_p);
+
+// static numeric -48
+const numeric *_num_48_p;
+const ex _ex_48 = ex(*_num_48_p);
+
+// static numeric -30
+const numeric *_num_30_p;
+const ex _ex_30 = ex(*_num_30_p);
+
+// static numeric -25
+const numeric *_num_25_p;
+const ex _ex_25 = ex(*_num_25_p);
+
+// static numeric -24
+const numeric *_num_24_p;
+const ex _ex_24 = ex(*_num_24_p);
+
+// static numeric -20
+const numeric *_num_20_p;
+const ex _ex_20 = ex(*_num_20_p);
+
+// static numeric -18
+const numeric *_num_18_p;
+const ex _ex_18 = ex(*_num_18_p);
+
+// static numeric -15
+const numeric *_num_15_p;
+const ex _ex_15 = ex(*_num_15_p);
+
+// static numeric -12
+const numeric *_num_12_p;
+const ex _ex_12 = ex(*_num_12_p);
+
+// static numeric -11
+const numeric *_num_11_p;
+const ex _ex_11 = ex(*_num_11_p);
+
+// static numeric -10
+const numeric *_num_10_p;
+const ex _ex_10 = ex(*_num_10_p);
+
+// static numeric -9
+const numeric *_num_9_p;
+const ex _ex_9 = ex(*_num_9_p);
+
+// static numeric -8
+const numeric *_num_8_p;
+const ex _ex_8 = ex(*_num_8_p);
+
+// static numeric -7
+const numeric *_num_7_p;
+const ex _ex_7 = ex(*_num_7_p);
+
+// static numeric -6
+const numeric *_num_6_p;
+const ex _ex_6 = ex(*_num_6_p);
+
+// static numeric -5
+const numeric *_num_5_p;
+const ex _ex_5 = ex(*_num_5_p);
+
+// static numeric -4
+const numeric *_num_4_p;
+const ex _ex_4 = ex(*_num_4_p);
+
+// static numeric -3
+const numeric *_num_3_p;
+const ex _ex_3 = ex(*_num_3_p);
+
+// static numeric -2
+const numeric *_num_2_p;
+const ex _ex_2 = ex(*_num_2_p);
+
+// static numeric -1
+const numeric *_num_1_p;
+const ex _ex_1 = ex(*_num_1_p);
+
+// static numeric -1/2
+const numeric *_num_1_2_p;
+const ex _ex_1_2= ex(*_num_1_2_p);
+
+// static numeric -1/3
+const numeric *_num_1_3_p;
+const ex _ex_1_3= ex(*_num_1_3_p);
+
+// static numeric -1/4
+const numeric *_num_1_4_p;
+const ex _ex_1_4= ex(*_num_1_4_p);
+
+// static numeric 0
+const numeric *_num0_p;
+const basic *_num0_bp;
+const ex _ex0 = ex(*_num0_p);
+
+// static numeric 1/4
+const numeric *_num1_4_p;
+const ex _ex1_4 = ex(*_num1_4_p);
+
+// static numeric 1/3
+const numeric *_num1_3_p;
+const ex _ex1_3 = ex(*_num1_3_p);
+
+// static numeric 1/2
+const numeric *_num1_2_p;
+const ex _ex1_2 = ex(*_num1_2_p);
+
+// static numeric 1
+const numeric *_num1_p;
+const ex _ex1 = ex(*_num1_p);
+
+// static numeric 2
+const numeric *_num2_p;
+const ex _ex2 = ex(*_num2_p);
+
+// static numeric 3
+const numeric *_num3_p;
+const ex _ex3 = ex(*_num3_p);
+
+// static numeric 4
+const numeric *_num4_p;
+const ex _ex4 = ex(*_num4_p);
+
+// static numeric 5
+const numeric *_num5_p;
+const ex _ex5 = ex(*_num5_p);
+
+// static numeric 6
+const numeric *_num6_p;
+const ex _ex6 = ex(*_num6_p);
+
+// static numeric 7
+const numeric *_num7_p;
+const ex _ex7 = ex(*_num7_p);
+
+// static numeric 8
+const numeric *_num8_p;
+const ex _ex8 = ex(*_num8_p);
+
+// static numeric 9
+const numeric *_num9_p;
+const ex _ex9 = ex(*_num9_p);
+
+// static numeric 10
+const numeric *_num10_p;
+const ex _ex10 = ex(*_num10_p);
+
+// static numeric 11
+const numeric *_num11_p;
+const ex _ex11 = ex(*_num11_p);
+
+// static numeric 12
+const numeric *_num12_p;
+const ex _ex12 = ex(*_num12_p);
+
+// static numeric 15
+const numeric *_num15_p;
+const ex _ex15 = ex(*_num15_p);
+
+// static numeric 18
+const numeric *_num18_p;
+const ex _ex18 = ex(*_num18_p);
+
+// static numeric 20
+const numeric *_num20_p;
+const ex _ex20 = ex(*_num20_p);
+
+// static numeric 24
+const numeric *_num24_p;
+const ex _ex24 = ex(*_num24_p);
+
+// static numeric 25
+const numeric *_num25_p;
+const ex _ex25 = ex(*_num25_p);
+
+// static numeric 30
+const numeric *_num30_p;
+const ex _ex30 = ex(*_num30_p);
+
+// static numeric 48
+const numeric *_num48_p;
+const ex _ex48 = ex(*_num48_p);
+
+// static numeric 60
+const numeric *_num60_p;
+const ex _ex60 = ex(*_num60_p);
+
+// static numeric 120
+const numeric *_num120_p;
+const ex _ex120 = ex(*_num120_p);
+
// comment skeleton for header files
* of any interest to the user of the library. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
// Generates all distinct permutations of a multiset.
// (Based on Aaron Williams' algorithm 1 from "Loopless Generation of
// Multiset Permutations using a Constant Number of Variables by Prefix
- // Shifts." <http://webhome.csc.uvic.ca/~haron/CoolMulti.pdf>)
+ // Shifts." <https://dl.acm.org/doi/pdf/10.5555/1496770.1496877>)
struct coolmulti {
// element of singly linked list
struct element {
* Utilities for summing over multiple indices */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* GiNaC library version information. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
#define GINACLIB_MINOR_VERSION 8
/* Micro version of GiNaC */
-#define GINACLIB_MICRO_VERSION 0
+#define GINACLIB_MICRO_VERSION 7
// GiNaC library version information. It has very little to do with GiNaC
// version number. In particular, library version is OS dependent.
// increasing. This doesn't matter, though: there is not incurred cost
// for numbers that are omitted, except for shrinking the available space
// of leftover numbers. Not something we need to worry about yet. ;-)
-// TODO, when setting GINAC_LT_REVISION to 0:
+//
+// On Linux, the SONAME is libginac.so.$(GINAC_LT_CURRENT)-$(GINAC_LT_AGE).
+//
+// TODO, when breaking the SONAME:
// * change matrix inverse to use default argument (twice)
-// * remove interfaces marked as deprecated
-#define GINAC_LT_CURRENT 11
-#define GINAC_LT_REVISION 0
-#define GINAC_LT_AGE 0
+// * check for interfaces marked as deprecated
+#define GINAC_LT_CURRENT 12
+#define GINAC_LT_REVISION 6
+#define GINAC_LT_AGE 1
/*
* GiNaC archive file version information.
* Implementation of GiNaC's wildcard objects. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* Interface to GiNaC's wildcard objects. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
.BI series( expression ", " relation-or-symbol ", " order )
\- series expansion
.br
+.BI series_to_poly( series )
+\- convert a series into a polynomial by dropping the Order() term
+.br
.BI sprem( expression ", " expression ", " symbol )
\- sparse pseudo-remainder of polynomials
.br
.BI sqrfree( "expression [" ", " symbol-list] )
\- square-free factorization of a polynomial
.br
+.BI sqrfree_parfrac( expression ", " symbol )
+\- square-free partial fraction decomposition of rational function
+.br
.BI sqrt( expression )
\- square root
.br
.PP
CLN \- A Class Library for Numbers, Bruno Haible
.SH COPYRIGHT
-Copyright \(co 1999-2020 Johannes Gutenberg Universit\(:at Mainz, Germany
+Copyright \(co 1999-2023 Johannes Gutenberg Universit\(:at 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
*
* Global definitions for ginsh.
*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* functions. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
* This file must be processed with flex. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
YY_FATAL_ERROR("input in flex scanner failed");
result = n;
#endif
- } else if (((result = fread(buf, 1, max_size, yyin)) == 0) && ferror(yyin))
- YY_FATAL_ERROR("input in flex scanner failed");
+ } else {
+ int c = '*', n;
+ for (n = 0; n < max_size && (c = getc(yyin)) != EOF && c != '\n'; ++n)
+ buf[n] = (char)c;
+ if (c == '\n')
+ buf[n++] = (char)c;
+ if (c == EOF && ferror(yyin))
+ YY_FATAL_ERROR("input in flex scanner failed");
+ result = n;
+ }
return result;
}
* This file must be processed with yacc/bison. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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
return e[0].series(e[1], ex_to<numeric>(e[2]).to_int());
}
+static ex f_series_to_poly(const exprseq &e)
+{
+ CHECK_ARG(0, pseries, series_to_poly);
+ return series_to_poly(ex_to<pseries>(e[0]));
+}
+
static ex f_sprem(const exprseq &e)
{
return sprem(e[0], e[1], e[2]);
return sqrfree(e[0], ex_to<lst>(e[1]));
}
+static ex f_sqrfree_parfrac(const exprseq &e)
+{
+ return sqrfree_parfrac(e[0], ex_to<symbol>(e[1]));
+}
+
static ex f_subs3(const exprseq &e)
{
CHECK_ARG(1, lst, subs);
{"rem", f_rem, 3},
{"resultant", f_resultant, 3},
{"series", f_series, 3},
+ {"series_to_poly", f_series_to_poly, 1},
{"sprem", f_sprem, 3},
{"sqrfree", f_sqrfree1, 1},
{"sqrfree", f_sqrfree2, 2},
+ {"sqrfree_parfrac", f_sqrfree_parfrac, 2},
{"sqrt", f_sqrt, 1},
{"subs", f_subs2, 2},
{"subs", f_subs3, 3},
void greeting(void)
{
cout << "ginsh - GiNaC Interactive Shell (GiNaC V" << GINACLIB_VERSION << ")" << endl;
- cout << " __, _______ Copyright (C) 1999-2020 Johannes Gutenberg University Mainz,\n"
+ cout << " __, _______ Copyright (C) 1999-2023 Johannes Gutenberg University Mainz,\n"
<< " (__) * | Germany. This is free software with ABSOLUTELY NO WARRANTY.\n"
<< " ._) i N a C | You are welcome to redistribute it under certain conditions.\n"
<< "<-------------' For details type `warranty;'.\n" << endl;
#
# Check for baseline language coverage in the compiler for the specified
# version of the C++ standard. If necessary, add switches to CXX and
-# CXXCPP to enable support. VERSION may be '11' (for the C++11 standard)
-# or '14' (for the C++14 standard).
+# CXXCPP to enable support. VERSION may be '11', '14', '17', or '20' for
+# the respective C++ standard version.
#
# The second argument, if specified, indicates whether you insist on an
# extended mode (e.g. -std=gnu++11) or a strict conformance mode (e.g.
# -std=c++11). If neither is specified, you get whatever works, with
-# preference for an extended mode.
+# preference for no added switch, and then for an extended mode.
#
# The third argument, if specified 'mandatory' or if left unspecified,
# indicates that baseline support for the specified C++ standard is
# Copyright (c) 2015 Moritz Klammler <moritz@klammler.eu>
# Copyright (c) 2016, 2018 Krzesimir Nowak <qdlacz@gmail.com>
# Copyright (c) 2019 Enji Cooper <yaneurabeya@gmail.com>
+# Copyright (c) 2020 Jason Merrill <jason@redhat.com>
+# Copyright (c) 2021 Jörn Heusipp <osmanx@problemloesungsmaschine.de>
#
# Copying and distribution of this file, with or without modification, are
# permitted in any medium without royalty provided the copyright notice
# and this notice are preserved. This file is offered as-is, without any
# warranty.
-#serial 11
+#serial 18
dnl This macro is based on the code from the AX_CXX_COMPILE_STDCXX_11 macro
dnl (serial version number 13).
m4_if([$1], [11], [ax_cxx_compile_alternatives="11 0x"],
[$1], [14], [ax_cxx_compile_alternatives="14 1y"],
[$1], [17], [ax_cxx_compile_alternatives="17 1z"],
+ [$1], [20], [ax_cxx_compile_alternatives="20"],
[m4_fatal([invalid first argument `$1' to AX_CXX_COMPILE_STDCXX])])dnl
m4_if([$2], [], [],
[$2], [ext], [],
AC_LANG_PUSH([C++])dnl
ac_success=no
+ m4_if([$2], [], [dnl
+ AC_CACHE_CHECK(whether $CXX supports C++$1 features by default,
+ ax_cv_cxx_compile_cxx$1,
+ [AC_COMPILE_IFELSE([AC_LANG_SOURCE([_AX_CXX_COMPILE_STDCXX_testbody_$1])],
+ [ax_cv_cxx_compile_cxx$1=yes],
+ [ax_cv_cxx_compile_cxx$1=no])])
+ if test x$ax_cv_cxx_compile_cxx$1 = xyes; then
+ ac_success=yes
+ fi])
+
m4_if([$2], [noext], [], [dnl
if test x$ac_success = xno; then
for alternative in ${ax_cxx_compile_alternatives}; do
dnl HP's aCC needs +std=c++11 according to:
dnl http://h21007.www2.hp.com/portal/download/files/unprot/aCxx/PDF_Release_Notes/769149-001.pdf
dnl Cray's crayCC needs "-h std=c++11"
+ dnl MSVC needs -std:c++NN for C++17 and later (default is C++14)
for alternative in ${ax_cxx_compile_alternatives}; do
- for switch in -std=c++${alternative} +std=c++${alternative} "-h std=c++${alternative}"; do
- cachevar=AS_TR_SH([ax_cv_cxx_compile_cxx$1_$switch])
+ for switch in -std=c++${alternative} +std=c++${alternative} "-h std=c++${alternative}" MSVC; do
+ if test x"$switch" = xMSVC; then
+ dnl AS_TR_SH maps both `:` and `=` to `_` so -std:c++17 would collide
+ dnl with -std=c++17. We suffix the cache variable name with _MSVC to
+ dnl avoid this.
+ switch=-std:c++${alternative}
+ cachevar=AS_TR_SH([ax_cv_cxx_compile_cxx$1_${switch}_MSVC])
+ else
+ cachevar=AS_TR_SH([ax_cv_cxx_compile_cxx$1_$switch])
+ fi
AC_CACHE_CHECK(whether $CXX supports C++$1 features with $switch,
$cachevar,
[ac_save_CXX="$CXX"
_AX_CXX_COMPILE_STDCXX_testbody_new_in_11
)
-
dnl Test body for checking C++14 support
m4_define([_AX_CXX_COMPILE_STDCXX_testbody_14],
_AX_CXX_COMPILE_STDCXX_testbody_new_in_14
)
+dnl Test body for checking C++17 support
+
m4_define([_AX_CXX_COMPILE_STDCXX_testbody_17],
_AX_CXX_COMPILE_STDCXX_testbody_new_in_11
_AX_CXX_COMPILE_STDCXX_testbody_new_in_14
_AX_CXX_COMPILE_STDCXX_testbody_new_in_17
)
+dnl Test body for checking C++20 support
+
+m4_define([_AX_CXX_COMPILE_STDCXX_testbody_20],
+ _AX_CXX_COMPILE_STDCXX_testbody_new_in_11
+ _AX_CXX_COMPILE_STDCXX_testbody_new_in_14
+ _AX_CXX_COMPILE_STDCXX_testbody_new_in_17
+ _AX_CXX_COMPILE_STDCXX_testbody_new_in_20
+)
+
+
dnl Tests for new features in C++11
m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_11], [[
#error "This is not a C++ compiler"
-#elif __cplusplus < 201103L
+// MSVC always sets __cplusplus to 199711L in older versions; newer versions
+// only set it correctly if /Zc:__cplusplus is specified as well as a
+// /std:c++NN switch:
+// https://devblogs.microsoft.com/cppblog/msvc-now-correctly-reports-__cplusplus/
+#elif __cplusplus < 201103L && !defined _MSC_VER
#error "This is not a C++11 compiler"
#error "This is not a C++ compiler"
-#elif __cplusplus < 201402L
+#elif __cplusplus < 201402L && !defined _MSC_VER
#error "This is not a C++14 compiler"
#error "This is not a C++ compiler"
-#elif __cplusplus < 201703L
+#elif __cplusplus < 201703L && !defined _MSC_VER
#error "This is not a C++17 compiler"
} // namespace cxx17
-#endif // __cplusplus < 201703L
+#endif // __cplusplus < 201703L && !defined _MSC_VER
+
+]])
+
+
+dnl Tests for new features in C++20
+
+m4_define([_AX_CXX_COMPILE_STDCXX_testbody_new_in_20], [[
+
+#ifndef __cplusplus
+
+#error "This is not a C++ compiler"
+
+#elif __cplusplus < 202002L && !defined _MSC_VER
+
+#error "This is not a C++20 compiler"
+
+#else
+
+#include <version>
+
+namespace cxx20
+{
+
+// As C++20 supports feature test macros in the standard, there is no
+// immediate need to actually test for feature availability on the
+// Autoconf side.
+
+} // namespace cxx20
+
+#endif // __cplusplus < 202002L && !defined _MSC_VER
]])
-# host-cpu-c-abi.m4 serial 12
-dnl Copyright (C) 2002-2019 Free Software Foundation, Inc.
+# host-cpu-c-abi.m4 serial 15
+dnl Copyright (C) 2002-2023 Free Software Foundation, Inc.
dnl This file is free software; the Free Software Foundation
dnl gives unlimited permission to copy and/or distribute it,
dnl with or without modifications, as long as this notice is preserved.
# be generating 64-bit code.
AC_COMPILE_IFELSE(
[AC_LANG_SOURCE(
- [[#if defined __powerpc64__ || defined _ARCH_PPC64
+ [[#if defined __powerpc64__ || defined __LP64__
int ok;
#else
error fail
#ifndef __ia64__
#undef __ia64__
#endif
+#ifndef __loongarch64__
+#undef __loongarch64__
+#endif
#ifndef __m68k__
#undef __m68k__
#endif
dnl Sets the HOST_CPU_C_ABI_32BIT variable to 'yes' if the C language ABI
-dnl (application binary interface) is a 32-bit one, or to 'no' otherwise.
+dnl (application binary interface) is a 32-bit one, to 'no' if it is a 64-bit
+dnl one, or to 'unknown' if unknown.
dnl This is a simplified variant of gl_HOST_CPU_C_ABI.
AC_DEFUN([gl_HOST_CPU_C_ABI_32BIT],
[
case "$gl_cv_host_cpu_c_abi" in
i386 | x86_64-x32 | arm | armhf | arm64-ilp32 | hppa | ia64-ilp32 | mips | mipsn32 | powerpc | riscv*-ilp32* | s390 | sparc)
gl_cv_host_cpu_c_abi_32bit=yes ;;
- *)
+ x86_64 | alpha | arm64 | hppa64 | ia64 | mips64 | powerpc64 | powerpc64-elfv2 | riscv*-lp64* | s390x | sparc64 )
gl_cv_host_cpu_c_abi_32bit=no ;;
+ *)
+ gl_cv_host_cpu_c_abi_32bit=unknown ;;
esac
else
case "$host_cpu" in
gl_cv_host_cpu_c_abi_32bit=yes
;;
+ # CPUs that only support a 64-bit ABI.
+changequote(,)dnl
+ alpha | alphaev[4-8] | alphaev56 | alphapca5[67] | alphaev6[78] \
+ | mmix )
+changequote([,])dnl
+ gl_cv_host_cpu_c_abi_32bit=no
+ ;;
+
changequote(,)dnl
i[34567]86 )
changequote([,])dnl
# be generating 64-bit code.
AC_COMPILE_IFELSE(
[AC_LANG_SOURCE(
- [[#if defined __powerpc64__ || defined _ARCH_PPC64
+ [[#if defined __powerpc64__ || defined __LP64__
int ok;
#else
error fail
;;
*)
- gl_cv_host_cpu_c_abi_32bit=no
+ gl_cv_host_cpu_c_abi_32bit=unknown
;;
esac
fi
-# lib-ld.m4 serial 9
-dnl Copyright (C) 1996-2003, 2009-2019 Free Software Foundation, Inc.
+# lib-ld.m4 serial 10
+dnl Copyright (C) 1996-2003, 2009-2023 Free Software Foundation, Inc.
dnl This file is free software; the Free Software Foundation
dnl gives unlimited permission to copy and/or distribute it,
dnl with or without modifications, as long as this notice is preserved.
*-*-aix*)
AC_COMPILE_IFELSE(
[AC_LANG_SOURCE(
- [[#if defined __powerpc64__ || defined _ARCH_PPC64
+ [[#if defined __powerpc64__ || defined __LP64__
int ok;
#else
error fail
-# lib-link.m4 serial 28
-dnl Copyright (C) 2001-2019 Free Software Foundation, Inc.
+# lib-link.m4 serial 33
+dnl Copyright (C) 2001-2023 Free Software Foundation, Inc.
dnl This file is free software; the Free Software Foundation
dnl gives unlimited permission to copy and/or distribute it,
dnl with or without modifications, as long as this notice is preserved.
AC_LIB_WITH_FINAL_PREFIX([
eval additional_includedir=\"$includedir\"
eval additional_libdir=\"$libdir\"
+ eval additional_libdir2=\"$exec_prefix/$acl_libdirstem2\"
+ eval additional_libdir3=\"$exec_prefix/$acl_libdirstem3\"
])
AC_ARG_WITH(PACK[-prefix],
-[[ --with-]]PACK[[-prefix[=DIR] search for ]PACKLIBS[ in DIR/include and DIR/lib
- --without-]]PACK[[-prefix don't search for ]PACKLIBS[ in includedir and libdir]],
+[[ --with-]]PACK[[-prefix[=DIR] search for ]]PACKLIBS[[ in DIR/include and DIR/lib
+ --without-]]PACK[[-prefix don't search for ]]PACKLIBS[[ in includedir and libdir]],
[
if test "X$withval" = "Xno"; then
use_additional=no
AC_LIB_WITH_FINAL_PREFIX([
eval additional_includedir=\"$includedir\"
eval additional_libdir=\"$libdir\"
+ eval additional_libdir2=\"$exec_prefix/$acl_libdirstem2\"
+ eval additional_libdir3=\"$exec_prefix/$acl_libdirstem3\"
])
else
additional_includedir="$withval/include"
additional_libdir="$withval/$acl_libdirstem"
- if test "$acl_libdirstem2" != "$acl_libdirstem" \
- && test ! -d "$withval/$acl_libdirstem"; then
- additional_libdir="$withval/$acl_libdirstem2"
- fi
+ additional_libdir2="$withval/$acl_libdirstem2"
+ additional_libdir3="$withval/$acl_libdirstem3"
fi
fi
])
+ if test "X$additional_libdir2" = "X$additional_libdir"; then
+ additional_libdir2=
+ fi
+ if test "X$additional_libdir3" = "X$additional_libdir"; then
+ additional_libdir3=
+ fi
dnl Search the library and its dependencies in $additional_libdir and
- dnl $LDFLAGS. Using breadth-first-seach.
+ dnl $LDFLAGS. Use breadth-first search.
LIB[]NAME=
LTLIB[]NAME=
INC[]NAME=
shrext=
fi
if test $use_additional = yes; then
- dir="$additional_libdir"
- dnl The same code as in the loop below:
- dnl First look for a shared library.
- if test -n "$acl_shlibext"; then
- if test -f "$dir/$libname$shrext"; then
- found_dir="$dir"
- found_so="$dir/$libname$shrext"
- else
- if test "$acl_library_names_spec" = '$libname$shrext$versuffix'; then
- ver=`(cd "$dir" && \
- for f in "$libname$shrext".*; do echo "$f"; done \
- | sed -e "s,^$libname$shrext\\\\.,," \
- | sort -t '.' -n -r -k1,1 -k2,2 -k3,3 -k4,4 -k5,5 \
- | sed 1q ) 2>/dev/null`
- if test -n "$ver" && test -f "$dir/$libname$shrext.$ver"; then
- found_dir="$dir"
- found_so="$dir/$libname$shrext.$ver"
+ for additional_libdir_variable in additional_libdir additional_libdir2 additional_libdir3; do
+ if test "X$found_dir" = "X"; then
+ eval dir=\$$additional_libdir_variable
+ if test -n "$dir"; then
+ dnl The same code as in the loop below:
+ dnl First look for a shared library.
+ if test -n "$acl_shlibext"; then
+ if test -f "$dir/$libname$shrext" && acl_is_expected_elfclass < "$dir/$libname$shrext"; then
+ found_dir="$dir"
+ found_so="$dir/$libname$shrext"
+ else
+ if test "$acl_library_names_spec" = '$libname$shrext$versuffix'; then
+ ver=`(cd "$dir" && \
+ for f in "$libname$shrext".*; do echo "$f"; done \
+ | sed -e "s,^$libname$shrext\\\\.,," \
+ | sort -t '.' -n -r -k1,1 -k2,2 -k3,3 -k4,4 -k5,5 \
+ | sed 1q ) 2>/dev/null`
+ if test -n "$ver" && test -f "$dir/$libname$shrext.$ver" && acl_is_expected_elfclass < "$dir/$libname$shrext.$ver"; then
+ found_dir="$dir"
+ found_so="$dir/$libname$shrext.$ver"
+ fi
+ else
+ eval library_names=\"$acl_library_names_spec\"
+ for f in $library_names; do
+ if test -f "$dir/$f" && acl_is_expected_elfclass < "$dir/$f"; then
+ found_dir="$dir"
+ found_so="$dir/$f"
+ break
+ fi
+ done
+ fi
+ fi
fi
- else
- eval library_names=\"$acl_library_names_spec\"
- for f in $library_names; do
- if test -f "$dir/$f"; then
+ dnl Then look for a static library.
+ if test "X$found_dir" = "X"; then
+ if test -f "$dir/$libname.$acl_libext" && ${AR-ar} -p "$dir/$libname.$acl_libext" | acl_is_expected_elfclass; then
found_dir="$dir"
- found_so="$dir/$f"
- break
+ found_a="$dir/$libname.$acl_libext"
fi
- done
+ fi
+ if test "X$found_dir" != "X"; then
+ if test -f "$dir/$libname.la"; then
+ found_la="$dir/$libname.la"
+ fi
+ fi
fi
fi
- fi
- dnl Then look for a static library.
- if test "X$found_dir" = "X"; then
- if test -f "$dir/$libname.$acl_libext"; then
- found_dir="$dir"
- found_a="$dir/$libname.$acl_libext"
- fi
- fi
- if test "X$found_dir" != "X"; then
- if test -f "$dir/$libname.la"; then
- found_la="$dir/$libname.la"
- fi
- fi
+ done
fi
if test "X$found_dir" = "X"; then
for x in $LDFLAGS $LTLIB[]NAME; do
dir=`echo "X$x" | sed -e 's/^X-L//'`
dnl First look for a shared library.
if test -n "$acl_shlibext"; then
- if test -f "$dir/$libname$shrext"; then
+ if test -f "$dir/$libname$shrext" && acl_is_expected_elfclass < "$dir/$libname$shrext"; then
found_dir="$dir"
found_so="$dir/$libname$shrext"
else
| sed -e "s,^$libname$shrext\\\\.,," \
| sort -t '.' -n -r -k1,1 -k2,2 -k3,3 -k4,4 -k5,5 \
| sed 1q ) 2>/dev/null`
- if test -n "$ver" && test -f "$dir/$libname$shrext.$ver"; then
+ if test -n "$ver" && test -f "$dir/$libname$shrext.$ver" && acl_is_expected_elfclass < "$dir/$libname$shrext.$ver"; then
found_dir="$dir"
found_so="$dir/$libname$shrext.$ver"
fi
else
eval library_names=\"$acl_library_names_spec\"
for f in $library_names; do
- if test -f "$dir/$f"; then
+ if test -f "$dir/$f" && acl_is_expected_elfclass < "$dir/$f"; then
found_dir="$dir"
found_so="$dir/$f"
break
fi
dnl Then look for a static library.
if test "X$found_dir" = "X"; then
- if test -f "$dir/$libname.$acl_libext"; then
+ if test -f "$dir/$libname.$acl_libext" && ${AR-ar} -p "$dir/$libname.$acl_libext" | acl_is_expected_elfclass; then
found_dir="$dir"
found_a="$dir/$libname.$acl_libext"
fi
dnl standard /usr/lib.
if test "$enable_rpath" = no \
|| test "X$found_dir" = "X/usr/$acl_libdirstem" \
- || test "X$found_dir" = "X/usr/$acl_libdirstem2"; then
+ || test "X$found_dir" = "X/usr/$acl_libdirstem2" \
+ || test "X$found_dir" = "X/usr/$acl_libdirstem3"; then
dnl No hardcoding is needed.
LIB[]NAME="${LIB[]NAME}${LIB[]NAME:+ }$found_so"
else
fi
additional_includedir="$basedir/include"
;;
+ */$acl_libdirstem3 | */$acl_libdirstem3/)
+ basedir=`echo "X$found_dir" | sed -e 's,^X,,' -e "s,/$acl_libdirstem3/"'*$,,'`
+ if test "$name" = '$1'; then
+ LIB[]NAME[]_PREFIX="$basedir"
+ fi
+ additional_includedir="$basedir/include"
+ ;;
esac
if test "X$additional_includedir" != "X"; then
dnl Potentially add $additional_includedir to $INCNAME.
for dep in $dependency_libs; do
case "$dep" in
-L*)
- additional_libdir=`echo "X$dep" | sed -e 's/^X-L//'`
- dnl Potentially add $additional_libdir to $LIBNAME and $LTLIBNAME.
+ dependency_libdir=`echo "X$dep" | sed -e 's/^X-L//'`
+ dnl Potentially add $dependency_libdir to $LIBNAME and $LTLIBNAME.
dnl But don't add it
dnl 1. if it's the standard /usr/lib,
dnl 2. if it's /usr/local/lib and we are using GCC on Linux,
dnl 3. if it's already present in $LDFLAGS or the already
dnl constructed $LIBNAME,
dnl 4. if it doesn't exist as a directory.
- if test "X$additional_libdir" != "X/usr/$acl_libdirstem" \
- && test "X$additional_libdir" != "X/usr/$acl_libdirstem2"; then
+ if test "X$dependency_libdir" != "X/usr/$acl_libdirstem" \
+ && test "X$dependency_libdir" != "X/usr/$acl_libdirstem2" \
+ && test "X$dependency_libdir" != "X/usr/$acl_libdirstem3"; then
haveit=
- if test "X$additional_libdir" = "X/usr/local/$acl_libdirstem" \
- || test "X$additional_libdir" = "X/usr/local/$acl_libdirstem2"; then
+ if test "X$dependency_libdir" = "X/usr/local/$acl_libdirstem" \
+ || test "X$dependency_libdir" = "X/usr/local/$acl_libdirstem2" \
+ || test "X$dependency_libdir" = "X/usr/local/$acl_libdirstem3"; then
if test -n "$GCC"; then
case $host_os in
linux* | gnu* | k*bsd*-gnu) haveit=yes;;
haveit=
for x in $LDFLAGS $LIB[]NAME; do
AC_LIB_WITH_FINAL_PREFIX([eval x=\"$x\"])
- if test "X$x" = "X-L$additional_libdir"; then
+ if test "X$x" = "X-L$dependency_libdir"; then
haveit=yes
break
fi
done
if test -z "$haveit"; then
- if test -d "$additional_libdir"; then
- dnl Really add $additional_libdir to $LIBNAME.
- LIB[]NAME="${LIB[]NAME}${LIB[]NAME:+ }-L$additional_libdir"
+ if test -d "$dependency_libdir"; then
+ dnl Really add $dependency_libdir to $LIBNAME.
+ LIB[]NAME="${LIB[]NAME}${LIB[]NAME:+ }-L$dependency_libdir"
fi
fi
haveit=
for x in $LDFLAGS $LTLIB[]NAME; do
AC_LIB_WITH_FINAL_PREFIX([eval x=\"$x\"])
- if test "X$x" = "X-L$additional_libdir"; then
+ if test "X$x" = "X-L$dependency_libdir"; then
haveit=yes
break
fi
done
if test -z "$haveit"; then
- if test -d "$additional_libdir"; then
- dnl Really add $additional_libdir to $LTLIBNAME.
- LTLIB[]NAME="${LTLIB[]NAME}${LTLIB[]NAME:+ }-L$additional_libdir"
+ if test -d "$dependency_libdir"; then
+ dnl Really add $dependency_libdir to $LTLIBNAME.
+ LTLIB[]NAME="${LTLIB[]NAME}${LTLIB[]NAME:+ }-L$dependency_libdir"
fi
fi
fi
;;
-l*)
dnl Handle this in the next round.
- names_next_round="$names_next_round "`echo "X$dep" | sed -e 's/^X-l//'`
+ dnl But on GNU systems, ignore -lc options, because
+ dnl - linking with libc is the default anyway,
+ dnl - linking with libc.a may produce an error
+ dnl "/usr/bin/ld: dynamic STT_GNU_IFUNC symbol `strcmp' with pointer equality in `/usr/lib/libc.a(strcmp.o)' can not be used when making an executable; recompile with -fPIE and relink with -pie"
+ dnl or may produce an executable that always crashes, see
+ dnl <https://lists.gnu.org/archive/html/grep-devel/2020-09/msg00052.html>.
+ dep=`echo "X$dep" | sed -e 's/^X-l//'`
+ if test "X$dep" != Xc \
+ || case $host_os in
+ linux* | gnu* | k*bsd*-gnu) false ;;
+ *) true ;;
+ esac; then
+ names_next_round="$names_next_round $dep"
+ fi
;;
*.la)
dnl Handle this in the next round. Throw away the .la's
dir="$next"
dnl No need to hardcode the standard /usr/lib.
if test "X$dir" != "X/usr/$acl_libdirstem" \
- && test "X$dir" != "X/usr/$acl_libdirstem2"; then
+ && test "X$dir" != "X/usr/$acl_libdirstem2" \
+ && test "X$dir" != "X/usr/$acl_libdirstem3"; then
rpathdirs="$rpathdirs $dir"
fi
next=
-L*) dir=`echo "X$opt" | sed -e 's,^X-L,,'`
dnl No need to hardcode the standard /usr/lib.
if test "X$dir" != "X/usr/$acl_libdirstem" \
- && test "X$dir" != "X/usr/$acl_libdirstem2"; then
+ && test "X$dir" != "X/usr/$acl_libdirstem2" \
+ && test "X$dir" != "X/usr/$acl_libdirstem3"; then
rpathdirs="$rpathdirs $dir"
fi
next= ;;
-# lib-prefix.m4 serial 14
-dnl Copyright (C) 2001-2005, 2008-2019 Free Software Foundation, Inc.
+# lib-prefix.m4 serial 20
+dnl Copyright (C) 2001-2005, 2008-2023 Free Software Foundation, Inc.
dnl This file is free software; the Free Software Foundation
dnl gives unlimited permission to copy and/or distribute it,
dnl with or without modifications, as long as this notice is preserved.
])
dnl AC_LIB_PREPARE_MULTILIB creates
-dnl - a variable acl_libdirstem, containing the basename of the libdir, either
-dnl "lib" or "lib64" or "lib/64",
-dnl - a variable acl_libdirstem2, as a secondary possible value for
-dnl acl_libdirstem, either the same as acl_libdirstem or "lib/sparcv9" or
-dnl "lib/amd64".
+dnl - a function acl_is_expected_elfclass, that tests whether standard input
+dn; has a 32-bit or 64-bit ELF header, depending on the host CPU ABI,
+dnl - 3 variables acl_libdirstem, acl_libdirstem2, acl_libdirstem3, containing
+dnl the basename of the libdir to try in turn, either "lib" or "lib64" or
+dnl "lib/64" or "lib32" or "lib/sparcv9" or "lib/amd64" or similar.
AC_DEFUN([AC_LIB_PREPARE_MULTILIB],
[
- dnl There is no formal standard regarding lib and lib64.
- dnl On glibc systems, the current practice is that on a system supporting
+ dnl There is no formal standard regarding lib, lib32, and lib64.
+ dnl On most glibc systems, the current practice is that on a system supporting
dnl 32-bit and 64-bit instruction sets or ABIs, 64-bit libraries go under
- dnl $prefix/lib64 and 32-bit libraries go under $prefix/lib. We determine
- dnl the compiler's default mode by looking at the compiler's library search
- dnl path. If at least one of its elements ends in /lib64 or points to a
- dnl directory whose absolute pathname ends in /lib64, we assume a 64-bit ABI.
- dnl Otherwise we use the default, namely "lib".
+ dnl $prefix/lib64 and 32-bit libraries go under $prefix/lib. However, on
+ dnl Arch Linux based distributions, it's the opposite: 32-bit libraries go
+ dnl under $prefix/lib32 and 64-bit libraries go under $prefix/lib.
+ dnl We determine the compiler's default mode by looking at the compiler's
+ dnl library search path. If at least one of its elements ends in /lib64 or
+ dnl points to a directory whose absolute pathname ends in /lib64, we use that
+ dnl for 64-bit ABIs. Similarly for 32-bit ABIs. Otherwise we use the default,
+ dnl namely "lib".
dnl On Solaris systems, the current practice is that on a system supporting
dnl 32-bit and 64-bit instruction sets or ABIs, 64-bit libraries go under
dnl $prefix/lib/64 (which is a symlink to either $prefix/lib/sparcv9 or
AC_REQUIRE([AC_CANONICAL_HOST])
AC_REQUIRE([gl_HOST_CPU_C_ABI_32BIT])
- case "$host_os" in
- solaris*)
- AC_CACHE_CHECK([for 64-bit host], [gl_cv_solaris_64bit],
- [AC_COMPILE_IFELSE(
- [AC_LANG_SOURCE(
- [[#ifdef _LP64
- int ok;
- #else
- error fail
- #endif
- ]])],
- [gl_cv_solaris_64bit=yes],
- [gl_cv_solaris_64bit=no])
- ]);;
- esac
+ AC_CACHE_CHECK([for ELF binary format], [gl_cv_elf],
+ [AC_EGREP_CPP([Extensible Linking Format],
+ [#if defined __ELF__ || (defined __linux__ && defined __EDG__)
+ Extensible Linking Format
+ #endif
+ ],
+ [gl_cv_elf=yes],
+ [gl_cv_elf=no])
+ ])
+ if test $gl_cv_elf = yes; then
+ # Extract the ELF class of a file (5th byte) in decimal.
+ # Cf. https://en.wikipedia.org/wiki/Executable_and_Linkable_Format#File_header
+ if od -A x < /dev/null >/dev/null 2>/dev/null; then
+ # Use POSIX od.
+ func_elfclass ()
+ {
+ od -A n -t d1 -j 4 -N 1
+ }
+ else
+ # Use BSD hexdump.
+ func_elfclass ()
+ {
+ dd bs=1 count=1 skip=4 2>/dev/null | hexdump -e '1/1 "%3d "'
+ echo
+ }
+ fi
+ # Use 'expr', not 'test', to compare the values of func_elfclass, because on
+ # Solaris 11 OpenIndiana and Solaris 11 OmniOS, the result is 001 or 002,
+ # not 1 or 2.
+changequote(,)dnl
+ case $HOST_CPU_C_ABI_32BIT in
+ yes)
+ # 32-bit ABI.
+ acl_is_expected_elfclass ()
+ {
+ expr "`func_elfclass | sed -e 's/[ ]//g'`" = 1 > /dev/null
+ }
+ ;;
+ no)
+ # 64-bit ABI.
+ acl_is_expected_elfclass ()
+ {
+ expr "`func_elfclass | sed -e 's/[ ]//g'`" = 2 > /dev/null
+ }
+ ;;
+ *)
+ # Unknown.
+ acl_is_expected_elfclass ()
+ {
+ :
+ }
+ ;;
+ esac
+changequote([,])dnl
+ else
+ acl_is_expected_elfclass ()
+ {
+ :
+ }
+ fi
dnl Allow the user to override the result by setting acl_cv_libdirstems.
AC_CACHE_CHECK([for the common suffixes of directories in the library search path],
[acl_cv_libdirstems],
- [acl_libdirstem=lib
+ [dnl Try 'lib' first, because that's the default for libdir in GNU, see
+ dnl <https://www.gnu.org/prep/standards/html_node/Directory-Variables.html>.
+ acl_libdirstem=lib
acl_libdirstem2=
+ acl_libdirstem3=
case "$host_os" in
solaris*)
dnl See Solaris 10 Software Developer Collection > Solaris 64-bit Developer's Guide > The Development Environment
dnl "Portable Makefiles should refer to any library directories using the 64 symbolic link."
dnl But we want to recognize the sparcv9 or amd64 subdirectory also if the
dnl symlink is missing, so we set acl_libdirstem2 too.
- if test $gl_cv_solaris_64bit = yes; then
- acl_libdirstem=lib/64
+ if test $HOST_CPU_C_ABI_32BIT = no; then
+ acl_libdirstem2=lib/64
case "$host_cpu" in
- sparc*) acl_libdirstem2=lib/sparcv9 ;;
- i*86 | x86_64) acl_libdirstem2=lib/amd64 ;;
+ sparc*) acl_libdirstem3=lib/sparcv9 ;;
+ i*86 | x86_64) acl_libdirstem3=lib/amd64 ;;
esac
fi
;;
*)
dnl If $CC generates code for a 32-bit ABI, the libraries are
- dnl surely under $prefix/lib, not $prefix/lib64.
- if test "$HOST_CPU_C_ABI_32BIT" != yes; then
- dnl The result is a property of the system. However, non-system
- dnl compilers sometimes have odd library search paths. Therefore
- dnl prefer asking /usr/bin/gcc, if available, rather than $CC.
- searchpath=`(if test -f /usr/bin/gcc \
- && LC_ALL=C /usr/bin/gcc -print-search-dirs >/dev/null 2>/dev/null; then \
- LC_ALL=C /usr/bin/gcc -print-search-dirs; \
- else \
- LC_ALL=C $CC -print-search-dirs; \
- fi) 2>/dev/null \
- | sed -n -e 's,^libraries: ,,p' | sed -e 's,^=,,'`
- if test -n "$searchpath"; then
- acl_save_IFS="${IFS= }"; IFS=":"
- for searchdir in $searchpath; do
- if test -d "$searchdir"; then
- case "$searchdir" in
- */lib64/ | */lib64 ) acl_libdirstem=lib64 ;;
- */../ | */.. )
- # Better ignore directories of this form. They are misleading.
- ;;
- *) searchdir=`cd "$searchdir" && pwd`
- case "$searchdir" in
- */lib64 ) acl_libdirstem=lib64 ;;
- esac ;;
- esac
- fi
- done
- IFS="$acl_save_IFS"
+ dnl surely under $prefix/lib or $prefix/lib32, not $prefix/lib64.
+ dnl Similarly, if $CC generates code for a 64-bit ABI, the libraries
+ dnl are surely under $prefix/lib or $prefix/lib64, not $prefix/lib32.
+ dnl Find the compiler's search path. However, non-system compilers
+ dnl sometimes have odd library search paths. But we can't simply invoke
+ dnl '/usr/bin/gcc -print-search-dirs' because that would not take into
+ dnl account the -m32/-m31 or -m64 options from the $CC or $CFLAGS.
+ searchpath=`(LC_ALL=C $CC $CPPFLAGS $CFLAGS -print-search-dirs) 2>/dev/null \
+ | sed -n -e 's,^libraries: ,,p' | sed -e 's,^=,,'`
+ if test $HOST_CPU_C_ABI_32BIT != no; then
+ # 32-bit or unknown ABI.
+ if test -d /usr/lib32; then
+ acl_libdirstem2=lib32
+ fi
+ fi
+ if test $HOST_CPU_C_ABI_32BIT != yes; then
+ # 64-bit or unknown ABI.
+ if test -d /usr/lib64; then
+ acl_libdirstem3=lib64
+ fi
+ fi
+ if test -n "$searchpath"; then
+ acl_save_IFS="${IFS= }"; IFS=":"
+ for searchdir in $searchpath; do
+ if test -d "$searchdir"; then
+ case "$searchdir" in
+ */lib32/ | */lib32 ) acl_libdirstem2=lib32 ;;
+ */lib64/ | */lib64 ) acl_libdirstem3=lib64 ;;
+ */../ | */.. )
+ # Better ignore directories of this form. They are misleading.
+ ;;
+ *) searchdir=`cd "$searchdir" && pwd`
+ case "$searchdir" in
+ */lib32 ) acl_libdirstem2=lib32 ;;
+ */lib64 ) acl_libdirstem3=lib64 ;;
+ esac ;;
+ esac
+ fi
+ done
+ IFS="$acl_save_IFS"
+ if test $HOST_CPU_C_ABI_32BIT = yes; then
+ # 32-bit ABI.
+ acl_libdirstem3=
+ fi
+ if test $HOST_CPU_C_ABI_32BIT = no; then
+ # 64-bit ABI.
+ acl_libdirstem2=
fi
fi
;;
esac
test -n "$acl_libdirstem2" || acl_libdirstem2="$acl_libdirstem"
- acl_cv_libdirstems="$acl_libdirstem,$acl_libdirstem2"
+ test -n "$acl_libdirstem3" || acl_libdirstem3="$acl_libdirstem"
+ acl_cv_libdirstems="$acl_libdirstem,$acl_libdirstem2,$acl_libdirstem3"
])
- # Decompose acl_cv_libdirstems into acl_libdirstem and acl_libdirstem2.
+ dnl Decompose acl_cv_libdirstems into acl_libdirstem, acl_libdirstem2, and
+ dnl acl_libdirstem3.
+changequote(,)dnl
acl_libdirstem=`echo "$acl_cv_libdirstems" | sed -e 's/,.*//'`
- acl_libdirstem2=`echo "$acl_cv_libdirstems" | sed -e '/,/s/.*,//'`
+ acl_libdirstem2=`echo "$acl_cv_libdirstems" | sed -e 's/^[^,]*,//' -e 's/,.*//'`
+ acl_libdirstem3=`echo "$acl_cv_libdirstems" | sed -e 's/^[^,]*,[^,]*,//' -e 's/,.*//'`
+changequote([,])dnl
])
+++ /dev/null
-
-# encoding: utf-8
-# TeX Live 2012 seems to dislike files produces by doxygen (1.8.x.y)
-# In particular, makeindex(1) program creates invalid index entries like
-# \hyperpage{NNN_}
-# (note the trailing underscore in the page number). This breaks automatic
-# builds and is very annoying. Hence this script. It replaces (broken)
-# \hyperpage{NNN_} with \hyperpage{NNN}.
-# Note: this is an ugly work around, a proper fix is welcome.
-import sys, os, re
-
-def fixupind(fname):
- """ Fix \\hyperpage{NNN_} entries in the ind file @var{fname} """
- tmpout = fname + '.tmp'
- inp = open(fname)
- out = open(tmpout, 'wt')
- rx = re.compile('(hyperpage)[{]([0-9]+)[_][}]')
- for line in inp:
- out.write(re.sub(rx, '\\1{\\2}', line))
- out.flush()
- out.close()
- inp.close()
- os.rename(tmpout, fname)
-
-if __name__ == '__main__':
- if len(sys.argv) <= 1:
- sys.exit(1)
- fixupind(sys.argv[1])
- sys.exit(0)
-
GiNaC Tutorial \- An open framework for symbolic computation within the
C++ programming language
.SH COPYRIGHT
-Copyright \(co 1999-2020 Johannes Gutenberg Universit\(:at Mainz, Germany
+Copyright \(co 1999-2023 Johannes Gutenberg Universit\(:at 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
* GiNaC archive file viewer. */
/*
- * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2023 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