- GiNaC-devel - ginac.de

Upcoming release
by Jens Vollinga 20 Nov '07

20 Nov '07

Hi, I think I have done all patches now. For tomorrow evening I plan to do the releases 1.3.9 and 1.4.1. Regards, Jens

1 0

bugfix-only release?
by Alexei Sheplyakov 20 Nov '07

20 Nov '07

Hi! Since two serious bugs have been fixed, it would be nice to make a bugfix-only release (1.3.9 and 1.4.1). If anyone has objections, or some patches, or problems/issues to be resolved first, please let me know. Best regards, Alexei -- All science is either physics or stamp collecting.

2 3

[PATCH] Constructing expression from a string does not leak memory any more.
by Alexei Sheplyakov 18 Nov '07

18 Nov '07

Thanks to valgrind for spotting this. --- check/parser_memleak.cpp | 36 ++++++++++++++++++++++++++++++++++++ ginac/input_lexer.ll | 4 +++- 2 files changed, 39 insertions(+), 1 deletions(-) create mode 100644 check/parser_memleak.cpp diff --git a/check/parser_memleak.cpp b/check/parser_memleak.cpp new file mode 100644 index 0000000..f347ce9 --- /dev/null +++ b/check/parser_memleak.cpp @@ -0,0 +1,36 @@ +/** + * This small program exhibits the memory leak in the ginac_yylex(). + * Run it as + * + * valgrind --leak-check=yes ./parser_memleak + * + * or simply + * + * ulimit -v `expr 64 \* 1024` ./parser_memleak + */ +#include <iostream> +#include <stdexcept> +#include <ginac/ginac.h> +using namespace std; +using namespace GiNaC; + +int main(int argc, char** argv) { + const symbol x("x"), y("y"); + const lst syms(x, y); + // parser-generated symbol => memory leak. + static const char* str[] = { "x^2+2*x*y + cos(x)", "Li2(x/y) + log(y/x)" }; + + // depends on the amount of the available VM, compiler options, etc. + const unsigned N_max = 500000; + unsigned N=0; + ex e; + try { + for (; N < N_max; N++) { + e = ex(str[N & 1], syms); + } + } catch (std::bad_alloc) { + cerr << "N = " << N << endl; + return 1; + } + return 0; +} diff --git a/ginac/input_lexer.ll b/ginac/input_lexer.ll index f931394..bc568ce 100644 --- a/ginac/input_lexer.ll +++ b/ginac/input_lexer.ll @@ -119,7 +119,9 @@ Digits ginac_yylval = (long)Digits; return T_DIGITS; {A}{AN}* { sym_tab::const_iterator i = syms.find(yytext); if (i == syms.end()) { - syms[yytext] = sym_def(ginac_yylval = *(new symbol(yytext)), false); + symbol tmp(yytext); + ginac_yylval = tmp; + syms[yytext] = sym_def(tmp, false); } else ginac_yylval = (*i).second.sym; return T_SYMBOL; -- 1.5.3.2 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH] Fix memory leaks in mul::eval() and power::eval() (GiNaC 1.4, HEAD)
by Alexei Sheplyakov 17 Nov '07

17 Nov '07

The bug was introduced by the patch commit f418c6ee4d558c852e1fb95533af07a3ae43f409 Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Wed Jul 11 14:34:42 2007 +0400 (it was commited into the official branch as commit a602d34c225dceb3e53742a7b3e19a4b5e280485 Author: Jens Vollinga <vollinga(a)thep.physik.uni-mainz.de> Date: Wed Jul 11 21:07:40 2007 +0000) While at it, delete bogus try -- catch block in the ex::series(), since it papers over some serious issues (like this one) and does not seem to do anything useful. Thanks to valgrind (http://valgrind.kde.org) developers for a wonderful piece of software. Debugging this problem would take way too long without it. --- check/mul_eval_memleak.cpp | 119 ++++++++++++++++++++++++++++++++++++++++++++ ginac/mul.cpp | 4 +- ginac/power.cpp | 4 +- ginac/pseries.cpp | 6 +-- 4 files changed, 124 insertions(+), 9 deletions(-) create mode 100644 check/mul_eval_memleak.cpp diff --git a/check/mul_eval_memleak.cpp b/check/mul_eval_memleak.cpp new file mode 100644 index 0000000..78e9312 --- /dev/null +++ b/check/mul_eval_memleak.cpp @@ -0,0 +1,119 @@ +/** + * @file mul_eval_memleak.cpp Test for memory leak in {mul,power}::eval + * + * The bug was introduced by + * + * commit f418c6ee4d558c852e1fb95533af07a3ae43f409 + * Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> + * Date: Wed Jul 11 14:34:42 2007 +0400 + * (it was commited into the official branch as + * commit a602d34c225dceb3e53742a7b3e19a4b5e280485 + * Author: Jens Vollinga <vollinga(a)thep.physik.uni-mainz.de> + * Date: Wed Jul 11 21:07:40 2007 +0000) + */ + +/* + * GiNaC Copyright (C) 1999-2007 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 <iostream> +#include <ginac/ginac.h> +#include <string> +#include <stdexcept> +using namespace std; +using namespace GiNaC; + +const unsigned check_mul_eval_memleak(const unsigned N) +{ + // Expanding this expression into a Laurent series triggers the bug. + static const string e_str("\ +1/605927415293858601*tgamma(3-eps)^(-1)*tgamma(2-eps)*(5013234896802\ +*(-1+2*eps)*((-2539001/2)*eps^(-2)+(-7617003/2+1692800*log(920*scale\ +^(-1))+3385202*log(1301*scale^(-1)))*eps^(-1)+F)+(846400+2539001*eps\ +)*(2539001*(-1+2*eps)*((-2539001/2)*eps^(-2)+(-7617003/2+1692800*log\ +(920*scale^(-1))+3385202*log(1301*scale^(-1)))*eps^(-1)+F)+286523497\ +2800*(-1+eps)^(-1)*eps^(-2)*(920*scale^(-1))^(-2*eps)*(1301*scale^(-\ +1))^(-2*eps)*tgamma(1+eps)^2)+6061411748045832000*(-1+eps)^(-1)*eps^\ +(-2)*(920*scale^(-1))^(-2*eps)*(1301*scale^(-1))^(-2*eps)*tgamma(1+e\ +ps)^2+716056132401*((-2539001/2)*eps^(-2)+(-7617003/2+1692800*log(92\ +0*scale^(-1))+3385202*log(1301*scale^(-1)))*eps^(-1)+F))+71656139360\ +0/716056132401*tgamma(1-eps)*tgamma(2*eps)^(-1)*tgamma(eps)^2*tgamma\ +(3-eps)^(-1)*(920*scale^(-1))^(-4*eps)*tgamma(-1+2*eps)-2/6059274152\ +93858601*tgamma(3-eps)^(-1)*(2149010446400*(-1+2*eps)*((-2539001/2)*\ +eps^(-2)+(-7617003/2+1692800*log(920*scale^(-1))+3385202*log(1301*sc\ +ale^(-1)))*eps^(-1)+F)+2425134880977920000*(-1+eps)^(-1)*eps^(-2)*(9\ +20*scale^(-1))^(-2*eps)*(1301*scale^(-1))^(-2*eps)*tgamma(1+eps)^2-7\ +16056132401*((-2539001/2)*eps^(-2)+(-7617003/2+1692800*log(920*scale\ +^(-1))+3385202*log(1301*scale^(-1)))*eps^(-1)+F)+1692601*(-1+2*eps)*\ +(2539001*(-1+2*eps)*((-2539001/2)*eps^(-2)+(-7617003/2+1692800*log(9\ +20*scale^(-1))+3385202*log(1301*scale^(-1)))*eps^(-1)+F)+28652349728\ +00*(-1+eps)^(-1)*eps^(-2)*(920*scale^(-1))^(-2*eps)*(1301*scale^(-1)\ +)^(-2*eps)*tgamma(1+eps)^2))*tgamma(2-eps)+(1/716056132401*I)*tgamma\ +(-1+eps)*mb^2*(mb*scale^(-1))^(-2*eps)*((716392960000*I)*mb^(-2)*(92\ +0*scale^(-1))^(-2*eps)*tgamma(-2+eps)-(2864898145201*I)*mb^(-2)*(130\ +1*scale^(-1))^(-2*eps)*tgamma(-2+eps)-(716224526400*I)*tgamma(-1+eps\ +)*mb^(-2)*(920*scale^(-1))^(-2*eps))-3385202/605927415293858601*tgam\ +ma(3-eps)^(-1)*tgamma(2-eps)*(2539001*(-1+2*eps)*((-2539001/2)*eps^(\ +-2)+(-7617003/2+1692800*log(920*scale^(-1))+3385202*log(1301*scale^(\ +-1)))*eps^(-1)+F)+2865234972800*(-1+eps)^(-1)*eps^(-2)*(920*scale^(-\ +1))^(-2*eps)*(1301*scale^(-1))^(-2*eps)*tgamma(1+eps)^2+846201*((-25\ +39001/2)*eps^(-2)+(-7617003/2+1692800*log(920*scale^(-1))+3385202*lo\ +g(1301*scale^(-1)))*eps^(-1)+F))\ +"); + const symbol eps("eps"), scale("scale"), mb("mb"), F("F"); + const lst syms(eps, scale, mb, F); + const ex e0(e_str, syms); + + unsigned i = 0; + unsigned n_failures = 0; + + ex e; + try { + for (; i < N; i++) + e = e0.series(eps, 1).subs(Euler==0).expand(); + } catch (std::bad_alloc) { + return i; + } + return 0; +} + +#include <sys/types.h> +#include <sys/time.h> +#include <sys/resource.h> +#include <unistd.h> + +static void set_VM_limit(const unsigned long MB) { + const unsigned mem_lim_kb = MB*1024*1024; + struct rlimit lim; + lim.rlim_cur = mem_lim_kb; + lim.rlim_max = mem_lim_kb; + setrlimit(RLIMIT_AS, &lim); +} + +int main(int argc, char** argv) { + static const unsigned max_mem = 32; // megabytes + // otherwise one need wait for a long[er] time. + set_VM_limit(max_mem); + static const unsigned n_of_tests = 10000; + const unsigned n_loops = check_mul_eval_memleak(n_of_tests); + if (n_loops) { + cerr << "memory exhausted after " << n_loops << " loops" << endl; + return 1; + } + return 0; +} + diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 4e95a63..9e1ba48 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -469,8 +469,8 @@ ex mul::eval(int level) const // XXX: What is the best way to check if the polynomial is a primitive? numeric c = i->rest.integer_content(); - const numeric& lead_coeff = - ex_to<numeric>(ex_to<add>(i->rest).seq.begin()->coeff).div_dyn(c); + const numeric lead_coeff = + ex_to<numeric>(ex_to<add>(i->rest).seq.begin()->coeff).div(c); const bool canonicalizable = lead_coeff.is_integer(); // XXX: The main variable is chosen in a random way, so this code diff --git a/ginac/power.cpp b/ginac/power.cpp index 7d909e8..cbf196f 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -505,8 +505,8 @@ ex power::eval(int level) const // (2*x + 6*y)^(-4) -> 1/16*(x + 3*y)^(-4) if (num_exponent->is_integer() && is_exactly_a<add>(ebasis)) { numeric icont = ebasis.integer_content(); - const numeric& lead_coeff = - ex_to<numeric>(ex_to<add>(ebasis).seq.begin()->coeff).div_dyn(icont); + const numeric lead_coeff = + ex_to<numeric>(ex_to<add>(ebasis).seq.begin()->coeff).div(icont); const bool canonicalizable = lead_coeff.is_integer(); const bool unit_normal = lead_coeff.is_pos_integer(); diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index 0da6ae8..f4571ed 100644 --- a/ginac/pseries.cpp +++ b/ginac/pseries.cpp @@ -1280,11 +1280,7 @@ ex ex::series(const ex & r, int order, unsigned options) const else throw (std::logic_error("ex::series(): expansion point has unknown type")); - try { - e = bp->series(rel_, order, options); - } catch (std::exception &x) { - throw (std::logic_error(std::string("unable to compute series (") + x.what() + ")")); - } + e = bp->series(rel_, order, options); return e; } -- 1.5.3.2 -- All science is either physics or stamp collecting.

1 0

[PATCH] Fix memory leaks in mul::eval() and power::eval() (GiNaC 1.3)
by Alexei Sheplyakov 17 Nov '07

17 Nov '07

The bug was introduced by the patch commit f418c6ee4d558c852e1fb95533af07a3ae43f409 Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Wed Jul 11 14:34:42 2007 +0400 (it was commited into the official branch as commit a602d34c225dceb3e53742a7b3e19a4b5e280485 Author: Jens Vollinga <vollinga(a)thep.physik.uni-mainz.de> Date: Wed Jul 11 21:07:40 2007 +0000) While at it, delete bogus try -- catch block in the ex::series(), since it papers over some serious issues (like this one) and does not seem to do anything useful. Thanks to valgrind (http://valgrind.kde.org) developers for a wonderful piece of software. Debugging this problem would take way too long without it. --- check/mul_eval_memleak.cpp | 119 ++++++++++++++++++++++++++++++++++++++++++++ ginac/mul.cpp | 4 +- ginac/power.cpp | 4 +- ginac/pseries.cpp | 6 +-- 4 files changed, 124 insertions(+), 9 deletions(-) create mode 100644 check/mul_eval_memleak.cpp diff --git a/check/mul_eval_memleak.cpp b/check/mul_eval_memleak.cpp new file mode 100644 index 0000000..78e9312 --- /dev/null +++ b/check/mul_eval_memleak.cpp @@ -0,0 +1,119 @@ +/** + * @file mul_eval_memleak.cpp Test for memory leak in {mul,power}::eval + * + * The bug was introduced by + * + * commit f418c6ee4d558c852e1fb95533af07a3ae43f409 + * Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> + * Date: Wed Jul 11 14:34:42 2007 +0400 + * (it was commited into the official branch as + * commit a602d34c225dceb3e53742a7b3e19a4b5e280485 + * Author: Jens Vollinga <vollinga(a)thep.physik.uni-mainz.de> + * Date: Wed Jul 11 21:07:40 2007 +0000) + */ + +/* + * GiNaC Copyright (C) 1999-2007 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 <iostream> +#include <ginac/ginac.h> +#include <string> +#include <stdexcept> +using namespace std; +using namespace GiNaC; + +const unsigned check_mul_eval_memleak(const unsigned N) +{ + // Expanding this expression into a Laurent series triggers the bug. + static const string e_str("\ +1/605927415293858601*tgamma(3-eps)^(-1)*tgamma(2-eps)*(5013234896802\ +*(-1+2*eps)*((-2539001/2)*eps^(-2)+(-7617003/2+1692800*log(920*scale\ +^(-1))+3385202*log(1301*scale^(-1)))*eps^(-1)+F)+(846400+2539001*eps\ +)*(2539001*(-1+2*eps)*((-2539001/2)*eps^(-2)+(-7617003/2+1692800*log\ +(920*scale^(-1))+3385202*log(1301*scale^(-1)))*eps^(-1)+F)+286523497\ +2800*(-1+eps)^(-1)*eps^(-2)*(920*scale^(-1))^(-2*eps)*(1301*scale^(-\ +1))^(-2*eps)*tgamma(1+eps)^2)+6061411748045832000*(-1+eps)^(-1)*eps^\ +(-2)*(920*scale^(-1))^(-2*eps)*(1301*scale^(-1))^(-2*eps)*tgamma(1+e\ +ps)^2+716056132401*((-2539001/2)*eps^(-2)+(-7617003/2+1692800*log(92\ +0*scale^(-1))+3385202*log(1301*scale^(-1)))*eps^(-1)+F))+71656139360\ +0/716056132401*tgamma(1-eps)*tgamma(2*eps)^(-1)*tgamma(eps)^2*tgamma\ +(3-eps)^(-1)*(920*scale^(-1))^(-4*eps)*tgamma(-1+2*eps)-2/6059274152\ +93858601*tgamma(3-eps)^(-1)*(2149010446400*(-1+2*eps)*((-2539001/2)*\ +eps^(-2)+(-7617003/2+1692800*log(920*scale^(-1))+3385202*log(1301*sc\ +ale^(-1)))*eps^(-1)+F)+2425134880977920000*(-1+eps)^(-1)*eps^(-2)*(9\ +20*scale^(-1))^(-2*eps)*(1301*scale^(-1))^(-2*eps)*tgamma(1+eps)^2-7\ +16056132401*((-2539001/2)*eps^(-2)+(-7617003/2+1692800*log(920*scale\ +^(-1))+3385202*log(1301*scale^(-1)))*eps^(-1)+F)+1692601*(-1+2*eps)*\ +(2539001*(-1+2*eps)*((-2539001/2)*eps^(-2)+(-7617003/2+1692800*log(9\ +20*scale^(-1))+3385202*log(1301*scale^(-1)))*eps^(-1)+F)+28652349728\ +00*(-1+eps)^(-1)*eps^(-2)*(920*scale^(-1))^(-2*eps)*(1301*scale^(-1)\ +)^(-2*eps)*tgamma(1+eps)^2))*tgamma(2-eps)+(1/716056132401*I)*tgamma\ +(-1+eps)*mb^2*(mb*scale^(-1))^(-2*eps)*((716392960000*I)*mb^(-2)*(92\ +0*scale^(-1))^(-2*eps)*tgamma(-2+eps)-(2864898145201*I)*mb^(-2)*(130\ +1*scale^(-1))^(-2*eps)*tgamma(-2+eps)-(716224526400*I)*tgamma(-1+eps\ +)*mb^(-2)*(920*scale^(-1))^(-2*eps))-3385202/605927415293858601*tgam\ +ma(3-eps)^(-1)*tgamma(2-eps)*(2539001*(-1+2*eps)*((-2539001/2)*eps^(\ +-2)+(-7617003/2+1692800*log(920*scale^(-1))+3385202*log(1301*scale^(\ +-1)))*eps^(-1)+F)+2865234972800*(-1+eps)^(-1)*eps^(-2)*(920*scale^(-\ +1))^(-2*eps)*(1301*scale^(-1))^(-2*eps)*tgamma(1+eps)^2+846201*((-25\ +39001/2)*eps^(-2)+(-7617003/2+1692800*log(920*scale^(-1))+3385202*lo\ +g(1301*scale^(-1)))*eps^(-1)+F))\ +"); + const symbol eps("eps"), scale("scale"), mb("mb"), F("F"); + const lst syms(eps, scale, mb, F); + const ex e0(e_str, syms); + + unsigned i = 0; + unsigned n_failures = 0; + + ex e; + try { + for (; i < N; i++) + e = e0.series(eps, 1).subs(Euler==0).expand(); + } catch (std::bad_alloc) { + return i; + } + return 0; +} + +#include <sys/types.h> +#include <sys/time.h> +#include <sys/resource.h> +#include <unistd.h> + +static void set_VM_limit(const unsigned long MB) { + const unsigned mem_lim_kb = MB*1024*1024; + struct rlimit lim; + lim.rlim_cur = mem_lim_kb; + lim.rlim_max = mem_lim_kb; + setrlimit(RLIMIT_AS, &lim); +} + +int main(int argc, char** argv) { + static const unsigned max_mem = 32; // megabytes + // otherwise one need wait for a long[er] time. + set_VM_limit(max_mem); + static const unsigned n_of_tests = 10000; + const unsigned n_loops = check_mul_eval_memleak(n_of_tests); + if (n_loops) { + cerr << "memory exhausted after " << n_loops << " loops" << endl; + return 1; + } + return 0; +} + diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 99d811f..ab1a042 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -461,8 +461,8 @@ ex mul::eval(int level) const // XXX: What is the best way to check if the polynomial is a primitive? numeric c = i->rest.integer_content(); - const numeric& lead_coeff = - ex_to<numeric>(ex_to<add>(i->rest).seq.begin()->coeff).div_dyn(c); + const numeric lead_coeff = + ex_to<numeric>(ex_to<add>(i->rest).seq.begin()->coeff).div(c); const bool canonicalizable = lead_coeff.is_integer(); // XXX: The main variable is chosen in a random way, so this code diff --git a/ginac/power.cpp b/ginac/power.cpp index e839220..099f935 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -491,8 +491,8 @@ ex power::eval(int level) const // (2*x + 6*y)^(-4) -> 1/16*(x + 3*y)^(-4) if (num_exponent->is_integer() && is_exactly_a<add>(ebasis)) { numeric icont = ebasis.integer_content(); - const numeric& lead_coeff = - ex_to<numeric>(ex_to<add>(ebasis).seq.begin()->coeff).div_dyn(icont); + const numeric lead_coeff = + ex_to<numeric>(ex_to<add>(ebasis).seq.begin()->coeff).div(icont); const bool canonicalizable = lead_coeff.is_integer(); const bool unit_normal = lead_coeff.is_pos_integer(); diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index 2b948a7..6ff1d96 100644 --- a/ginac/pseries.cpp +++ b/ginac/pseries.cpp @@ -1215,11 +1215,7 @@ ex ex::series(const ex & r, int order, unsigned options) const else throw (std::logic_error("ex::series(): expansion point has unknown type")); - try { - e = bp->series(rel_, order, options); - } catch (std::exception &x) { - throw (std::logic_error(std::string("unable to compute series (") + x.what() + ")")); - } + e = bp->series(rel_, order, options); return e; } -- 1.5.3.2 -- All science is either physics or stamp collecting.

1 0

[PATCH] fix broken perprocessor instruction, so GiNaC compiles when ! HAVE_LIBDL.
by Alexei Sheplyakov 06 Nov '07

06 Nov '07

Hi! I've tried to build GiNaC 1.4 on my ReactOS box, it failed due to broken preprocessor instruction in excomplier.cpp. The patch below fixes it. --- ginac/excompiler.cpp | 4 +++- 1 files changed, 3 insertions(+), 1 deletions(-) diff --git a/ginac/excompiler.cpp b/ginac/excompiler.cpp index fe96123..4cf2841 100644 --- a/ginac/excompiler.cpp +++ b/ginac/excompiler.cpp @@ -32,7 +32,9 @@ #include <string> #include <vector> +#ifdef HAVE_CONFIG_H #include "config.h" +#endif #ifdef HAVE_LIBDL #include <dlfcn.h> @@ -302,7 +304,7 @@ void unlink_ex(const std::string filename) global_excompiler.unlink(filename); } -#elif // def HAVE_LIBDL +#else // def HAVE_LIBDL /* * In case no working libdl has been found by configure, the following function -- 1.5.3.2 -- All science is either physics or stamp collecting.

1 0

Re: [GiNaC-cvs] ginac flags.h mul.cpp power.cpp
by Richard B. Kreckel 05 Nov '07

05 Nov '07

CVS Service wrote: > Update of /home/cvs/GiNaC/ginac by vollinga > > Modified Files: > Tag: ginac_1-4 > flags.h mul.cpp power.cpp > Log Message: > Synced to HEAD: > - This helps mul::expand() and friends to recognize objects which have > no indices at all [Sheplyakov]. Chris, We are curious what can be done about the performance regression one of your patches for dummy index renaming introduced. The seemingly innocuous patch has been identified by Jens and Alexei: <http://www.ginac.de/pipermail/ginac-list/2007-September/001215.html>. At first, the above log message suggested to me that the problem has been dealt with. But that is not so: The program is still extremely slow. In fact, Alexei's patches are entirely unrelated. Since neither an explanation of what's happening nor a fix has been suggested, we would really appreciate some expert insight by you. Thanks in advance -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>

4 5

Restore the speed of polynomial operations.
by Alexei Sheplyakov 18 Sep '07

18 Sep '07

Hello! These patches commit 0633ed8082961673eedc092689e06fa39d6bc322 Author: Jens Vollinga <vollinga(a)thep.physik.uni-mainz.de> Date: Thu May 19 18:10:40 2005 +0000 Fixed bug in expanding expressions containing dummy indices. [V.Kisil] (see <http://www.ginac.de/pipermail/ginac-cvs/2005-May/001586.html>) commit 6b1b11264fe4872132fef8cb8b22cb983ba01736 Author: Chris Dams <Chris.Dams(a)mi.infn.it> Date: Thu Jul 13 20:02:38 2006 +0000 Needed to call get_all_dummy_indices_safely more often. (see http://www.ginac.de/pipermail/ginac-cvs/2006-July/001780.html) commit 4cf53676d7823ec3c55a90df7723f7f1492a95b2 Author: Chris Dams <Chris.Dams(a)mi.infn.it> Date: Mon Jul 24 16:47:43 2006 +0000 Improved dummy index renaming. (see http://www.ginac.de/pipermail/ginac-cvs/2006-July/001782.html) made polynomial operations (in particular expand()) substantially slower. For example, Fateman's polynomial expanding benchark degraded from 32sec. to 80sec. My own code which uses multivariate GCD (which in turn uses .expand() extensively) become 4x slower. I'll post 4 patches which address the issue in reply to this message. Comments, questions, suggestions, etc. are wellcome. Best regards, Alexei -- All science is either physics or stamp collecting.

3 12

Volunteer to continue the work on NTL factorisation bindings
by Remco Bloemen 06 Sep '07

06 Sep '07

I have created a simple factorisation for polynomials in Q[x] with real roots (source attached, based on [2]) but it's really slow. Creating a reasonable implementation would require much work. An NTL based implementation would be very powerfull and feasible, as said in [1]. I would like to volunteer to continue this work. I have the following questions: 1) Where can I obtain the source of this binding? (or should I start from the code in gTybalt, or from scratch?) 2) Is there any chance such a patch will be included? (this would imply dependecy on NTL) 2a) And if it is compile-time optional? 2b) And if it also included integration of rational functions? 2c) And if it also included integration of trancedental functions? 3) Do you have any related advice, hints or pointers? [1] http://www.ginac.de/pipermail/ginac-devel/2004-August/000695.html [2] http://yacas.sourceforge.net/Algochapter3.html

3 3

[PATCH] improved CLN C++ code output
by Alexei Sheplyakov 05 Sep '07

05 Sep '07

* small integers are printed as int (not as cl_I) * no more ugly x+-y output * powers always printed using cln::expt function. Please note: this patch is for development branch only, since it breaks binary compatibility (due to power::do_print_csrc_cl_N). I'll post patches for GiNaC 1.[34] later. --- ginac/add.cpp | 34 +++++++++++++++------------------- ginac/numeric.cpp | 51 ++++++++++++++++++++++++++++++++++++++++++++++----- ginac/power.cpp | 35 +++++++++++++++++++++-------------- ginac/power.h | 1 + 4 files changed, 83 insertions(+), 38 deletions(-) diff --git a/ginac/add.cpp b/ginac/add.cpp index b37ab79..98e8c79 100644 --- a/ginac/add.cpp +++ b/ginac/add.cpp @@ -23,6 +23,7 @@ #include <iostream> #include <stdexcept> #include <limits> +#include <string> #include "add.h" #include "mul.h" @@ -172,37 +173,32 @@ void add::do_print_csrc(const print_csrc & c, unsigned level) const if (precedence() <= level) c.s << "("; - // Print arguments, separated by "+" + // Print arguments, separated by "+" or "-" epvector::const_iterator it = seq.begin(), itend = seq.end(); + char separator = ' '; while (it != itend) { - // If the coefficient is -1, it is replaced by a single minus sign - if (it->coeff.is_equal(_ex1)) { + // If the coefficient is negative, separator is "-" + if (it->coeff.is_equal(_ex_1) || + ex_to<numeric>(it->coeff).numer().is_equal(*_num_1_p)) + separator = '-'; + c.s << separator; + if (it->coeff.is_equal(_ex1) || it->coeff.is_equal(_ex_1)) { it->rest.print(c, precedence()); - } else if (it->coeff.is_equal(_ex_1)) { - c.s << "-"; + } else if (ex_to<numeric>(it->coeff).numer().is_equal(*_num1_p) || + ex_to<numeric>(it->coeff).numer().is_equal(*_num_1_p)) + { it->rest.print(c, precedence()); - } else if (ex_to<numeric>(it->coeff).numer().is_equal(*_num1_p)) { - it->rest.print(c, precedence()); - c.s << "/"; - ex_to<numeric>(it->coeff).denom().print(c, precedence()); - } else if (ex_to<numeric>(it->coeff).numer().is_equal(*_num_1_p)) { - c.s << "-"; - it->rest.print(c, precedence()); - c.s << "/"; + c.s << '/'; ex_to<numeric>(it->coeff).denom().print(c, precedence()); } else { it->coeff.print(c, precedence()); - c.s << "*"; + c.s << '*'; it->rest.print(c, precedence()); } - // Separator is "+", except if the following expression would have a leading minus sign or the sign is sitting in parenthesis (as in a ctor) ++it; - if (it != itend - && (is_a<print_csrc_cl_N>(c) || !it->coeff.info(info_flags::real) // sign inside ctor arguments - || !(it->coeff.info(info_flags::negative) || (it->coeff.is_equal(*_num1_p) && is_exactly_a<numeric>(it->rest) && it->rest.info(info_flags::negative))))) - c.s << "+"; + separator = '+'; } if (!overall_coeff.is_zero()) { diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp index ef61ce7..74f0e50 100644 --- a/ginac/numeric.cpp +++ b/ginac/numeric.cpp @@ -399,6 +399,38 @@ static void print_real_csrc(const print_context & c, const cln::cl_R & x) } } +template<typename T1, typename T2> +static inline bool coerce(T1& dst, const T2& arg); + +/** + * @brief Check if CLN integer can be converted into int + * + * @sa http://www.ginac.de/pipermail/cln-list/2006-October/000248.html + */ +template<> static inline bool +coerce<int, cln::cl_I>(int& dst, const cln::cl_I& arg) { + static const cln::cl_I cl_max_int = + (cln::cl_I)(long)(std::numeric_limits<int>::max()); + static const cln::cl_I cl_min_int = + (cln::cl_I)(long)(std::numeric_limits<int>::min()); + if ((arg >= cl_min_int) && (arg <= cl_max_int)) { + dst = cl_I_to_int(arg); + return true; + } + return false; +} + +template<> static inline bool +coerce<unsigned int, cln::cl_I>(unsigned int& dst, const cln::cl_I& arg) { + static const cln::cl_I cl_max_uint = + (cln::cl_I)(unsigned long)(std::numeric_limits<unsigned int>::max()); + if ((! minusp(arg)) && (arg <= cl_max_uint)) { + dst = cl_I_to_uint(arg); + return true; + } + return false; +} + /** Helper function to print real number in C++ source format using cl_N types. * * @see numeric::print() */ @@ -406,11 +438,20 @@ static void print_real_cl_N(const print_context & c, const cln::cl_R & x) { if (cln::instanceof(x, cln::cl_I_ring)) { - // Integer number - c.s << "cln::cl_I(\""; - print_real_number(c, x); - c.s << "\")"; - + int dst; + // fixnum + if (coerce(dst, cln::the<cln::cl_I>(x))) { + // can be converted to native int + if (dst < 0) + c.s << "(-" << dst << ")"; + else + c.s << dst; + } else { + // bignum + c.s << "cln::cl_I(\""; + print_real_number(c, x); + c.s << "\")"; + } } else if (cln::instanceof(x, cln::cl_RA_ring)) { // Rational number diff --git a/ginac/power.cpp b/ginac/power.cpp index faa0d05..68bec71 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -50,7 +50,8 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(power, basic, print_func<print_latex>(&power::do_print_latex). print_func<print_csrc>(&power::do_print_csrc). print_func<print_python>(&power::do_print_python). - print_func<print_python_repr>(&power::do_print_python_repr)) + print_func<print_python_repr>(&power::do_print_python_repr). + print_func<print_csrc_cl_N>(&power::do_print_csrc_cl_N)) typedef std::vector<int> intvector; @@ -162,6 +163,21 @@ static void print_sym_pow(const print_context & c, const symbol &x, int exp) } } +void power::do_print_csrc_cl_N(const print_csrc_cl_N& c, unsigned level) const +{ + if (exponent.is_equal(_ex_1)) { + c.s << "recip("; + basis.print(c); + c.s << ')'; + return; + } + c.s << "expt("; + basis.print(c); + c.s << ", "; + exponent.print(c); + c.s << ')'; +} + void power::do_print_csrc(const print_csrc & c, unsigned level) const { // Integer powers of symbols are printed in a special, optimized way @@ -172,29 +188,20 @@ void power::do_print_csrc(const print_csrc & c, unsigned level) const c.s << '('; else { exp = -exp; - if (is_a<print_csrc_cl_N>(c)) - c.s << "recip("; - else - c.s << "1.0/("; + c.s << "1.0/("; } print_sym_pow(c, ex_to<symbol>(basis), exp); c.s << ')'; // <expr>^-1 is printed as "1.0/<expr>" or with the recip() function of CLN } else if (exponent.is_equal(_ex_1)) { - if (is_a<print_csrc_cl_N>(c)) - c.s << "recip("; - else - c.s << "1.0/("; + c.s << "1.0/("; basis.print(c); c.s << ')'; - // Otherwise, use the pow() or expt() (CLN) functions + // Otherwise, use the pow() function } else { - if (is_a<print_csrc_cl_N>(c)) - c.s << "expt("; - else - c.s << "pow("; + c.s << "pow("; basis.print(c); c.s << ','; exponent.print(c); diff --git a/ginac/power.h b/ginac/power.h index f5e65d6..3375040 100644 --- a/ginac/power.h +++ b/ginac/power.h @@ -88,6 +88,7 @@ protected: void do_print_csrc(const print_csrc & c, unsigned level) const; void do_print_python(const print_python & c, unsigned level) const; void do_print_python_repr(const print_python_repr & c, unsigned level) const; + void do_print_csrc_cl_N(const print_csrc_cl_N & c, unsigned level) const; ex expand_add(const add & a, int n, unsigned options) const; ex expand_add_2(const add & a, unsigned options) const; -- 1.4.4.4 Best regards, Alexei -- All science is either physics or stamp collecting.

1 2