January 2002 - GiNaC-devel

Re: Alternative bernoulli(const numeric &)
by Markus Nullmeier 01 Feb '02

01 Feb '02

From: Richard B. Kreckel <kreckel(a)thep.physik.uni-mainz.de> > Hmm, I think a simple cast to long would get the job done, in the same > manner as the numeric ctors from builtin types work. Yes, a change to `long' would push the limit for even arguments to 65534. However I think the real problem lies in the "philosophical" nature of my question. I guess nobody will want to calculate Bernoulli numbers this big (the limit of the 1.0.3 code seems to be 8190). Thus I think things could be left as they are, since the CLN manual hints that conversions from `long' are less efficient. By altogether abolishing the theoretical limit and letting CLN calculate (p3-i)*(p2-i), we would slow the procedure down by some per cent without any real gain. But I suppose this road should be taken if 23168 did become an issue, like if (p < 23168) { normal_inner_loop; } else { slow_inner_loop_with_CLN; } If you like this better I can make a patch. Markus

3 11

GiNaC 1.0.5 released
by Christian Bauer 30 Jan '02

30 Jan '02

Hi! You probably won't see it on national TV but it's new, improved, and finger-licking good. :-) GiNaC 1.0.5 is binary compatible to 1.0.x and has - a slightly more versatile degree()/coeff()/collect() facility - a bugfix for the output of imaginary numbers As usual, you can get it from ftp://ftpthep.physik.uni-mainz.de/pub/GiNaC Bye, Christian -- / Coding on PowerPC and proud of it \/ http://www.uni-mainz.de/~bauec002/

3 3

Parenthesization bug
by Roberto Bagnara 29 Jan '02

29 Jan '02

Dear all, there is a bug somewhere in GiNaC whereby some non-redundant parentheses are not output. The problem is shown below by means of simple `ginsh' session, but it can be reproduced equally well with C++ code. $ ginsh ginsh - GiNaC Interactive Shell (GiNaC V1.0.4) __, _______ Copyright (C) 1999-2002 Johannes Gutenberg University Mainz, (__) * | Germany. This is free software with ABSOLUTELY NO WARRANTY. ._) i N a C | You are welcome to redistribute it under certain conditions. <-------------' For details type `warranty;'. Type ?? for a list of help topics. > 2*I^(1/3); 2*I^(1/3) > (2*I)^(1/3); 2*I^(1/3) > 2*(I^(1/3)); 2*I^(1/3) > If this is confirmed to be a bug, we would need to develop a fix or a workaround quite urgently. Thanks a lot Roberto -- Prof. Roberto Bagnara Computer Science Group Department of Mathematics, University of Parma, Italy http://www.cs.unipr.it/~bagnara/ mailto:bagnara@cs.unipr.it

3 5

1.0.4 unleashed
by Richard B. Kreckel 24 Jan '02

24 Jan '02

Like all other 1.0.n releases, this one claims binary compatibility to 1.0.0. Cebix has opened a branch in CVS now where ABI-breaking patches may go in (but is too shy to make the announcement). Anyway, the news for 1.0.4 are: * Speedup in expand(). * Faster Bernoulli numbers (thanks to Markus Nullmeier). * Some minor bugfixes and documentation updates. .rpm's already on FTP, .deb's to follow soon thru all woody-mirrors. Enjoy -richy.

1 0

PATCH: remove_all() method added to the containers
by Roberto Bagnara 20 Jan '02

20 Jan '02

Dear all, of course there are several ways in the current version of GiNaC to erase all the elements from a GiNaC container. However, as far as I can tell, they are either obscure or inefficient or both (removing one element at a time is inefficient, assigning from an empty container is obscure and still measurably inefficient and so forth). Removing all the elements of a container is a very common operation and forcing users to resort to kludges is, I believe, not a very good idea. Moreover, not providing the right method for a common operation reduces the implementor's latitude for future choices (what if you would like to use a different representation for containers? What if this new representation is such that the user's kludges for erasing all the elements are seriously inefficient?) You don't want your users to ask you whether _today_ it is best to clear all the elements using the technique A or by the tecnique B, right? It is not by chance that all the STL containers have a clear() method. Enough for advocacy, please find below a patch adding a remove_all() method to the containers. I feel that the name "remove_all" fits well with the others. The patch allows `make check' to complete successfully and the documentation looks good. I have been careful to follow your coding style. Please let me know if I haven't been successful and I will submit a revised patch. All the best Roberto -- Prof. Roberto Bagnara Computer Science Group Department of Mathematics, University of Parma, Italy http://www.cs.unipr.it/~bagnara/ mailto:bagnara@cs.unipr.it diff -rcp GiNaC-1.0.3.orig/doc/tutorial/ginac.texi GiNaC-1.0.3/doc/tutorial/ginac.texi *** GiNaC-1.0.3.orig/doc/tutorial/ginac.texi Fri Dec 21 12:38:03 2001 --- GiNaC-1.0.3/doc/tutorial/ginac.texi Sat Jan 19 18:22:28 2002 *************** canonical form. *** 1193,1198 **** --- 1193,1199 ---- @cindex @code{prepend()} @cindex @code{remove_first()} @cindex @code{remove_last()} + @cindex @code{remove_all()} The GiNaC class @code{lst} serves for holding a @dfn{list} of arbitrary expressions. These are sometimes used to supply a variable number of *************** and @code{prepend()} methods: *** 1230,1242 **** // ... @end example ! Finally you can remove the first or last element of a list with @code{remove_first()} and @code{remove_last()}: @example // ... l.remove_first(); // l is now @{x, 2, y, x+y, 4*x@} l.remove_last(); // l is now @{x, 2, y, x+y@} @} @end example --- 1231,1251 ---- // ... @end example ! You can remove the first or last element of a list with @code{remove_first()} and @code{remove_last()}: @example // ... l.remove_first(); // l is now @{x, 2, y, x+y, 4*x@} l.remove_last(); // l is now @{x, 2, y, x+y@} + @end example + + Finally, you can remove all the elements of a list with + @code{remove_all()}: + + @example + // ... + l.remove_all(); // l is now empty @} @end example diff -rcp GiNaC-1.0.3.orig/ginac/container.pl GiNaC-1.0.3/ginac/container.pl *** GiNaC-1.0.3.orig/ginac/container.pl Wed Dec 19 11:54:04 2001 --- GiNaC-1.0.3/ginac/container.pl Sat Jan 19 18:19:14 2002 *************** protected: *** 250,255 **** --- 250,256 ---- public: virtual ${CONTAINER} & append(const ex & b); virtual ${CONTAINER} & remove_last(void); + virtual ${CONTAINER} & remove_all(void); ${PREPEND_INTERFACE} ${SORT_INTERFACE} protected: *************** ${CONTAINER} & ${CONTAINER}::remove_last *** 547,552 **** --- 548,560 ---- return *this; } + ${CONTAINER} & ${CONTAINER}::remove_all(void) + { + ensure_if_modifiable(); + seq.clear(); + return *this; + } + ${PREPEND_IMPLEMENTATION} ${SORT_IMPLEMENTATION}

2 1

Alternative bernoulli(const numeric &)
by Markus Nullmeier 18 Jan '02

18 Jan '02

Hello, being unware of GiNaC's implemetation, I wrote my own Bernoulli number function after reading Peter Stone's straightforward description in his Maple Worksheet http://minyos.its.rmit.edu.au/~e05254/MAPLE/ENGMATH/TAYLOR/BERNUM.MWS Later I adapted my function to bernoulli() in numeric.cpp, as it turned out to be about twice as fast as GiNaC 1.0.3 (gcc 2.95, -O2). Thus maybe you'll find the enclosed patch useful in some way or the other. Hopefully no hidden bugs remain. Markus --- GiNaC-1.0.3/ginac/numeric.cpp Wed Dec 19 12:11:46 2001 +++ GiNaC-1.0.3_new_bernoulli/ginac/numeric.cpp Fri Jan 18 14:44:16 2002 @@ -1512,7 +1512,7 @@ // But if somebody works with the n'th Bernoulli number she is likely to // also need all previous Bernoulli numbers. So we need a complete remember // table and above divide and conquer algorithm is not suited to build one - // up. The code below is adapted from Pari's function bernvec(). + // up. // // (There is an interesting relation with the tangent polynomials described // in `Concrete Mathematics', which leads to a program twice as fast as our @@ -1520,38 +1520,38 @@ // addition to the remember table. This doubles the memory footprint so // we don't use it.) + unsigned n = nn.to_int(); + // the special cases not covered by the algorithm below - if (nn.is_equal(_num1)) - return _num_1_2; - if (nn.is_odd()) - return _num0; - + if (n & 1) + return (n == 1) ? _num_1_2 : _num0; + if (!n) + return _num1; + // store nonvanishing Bernoulli numbers here static std::vector< cln::cl_RA > results; - static int highest_result = 0; - // algorithm not applicable to B(0), so just store it - if (results.empty()) - results.push_back(cln::cl_RA(1)); - - int n = nn.to_long(); - for (int i=highest_result; i<n/2; ++i) { - cln::cl_RA B = 0; - long n = 8; - long m = 5; - long d1 = i; - long d2 = 2*i-1; - for (int j=i; j>0; --j) { - B = cln::cl_I(n*m) * (B+results[j]) / (d1*d2); - n += 4; - m += 2; - d1 -= 1; - d2 -= 2; - } - B = (1 - ((B+1)/(2*i+3))) / (cln::cl_I(1)<<(2*i+2)); - results.push_back(B); - ++highest_result; - } - return results[n/2]; + static unsigned next_r = 0; + + // algorithm not applicable to B(2), so just store it + if (!next_r) { + results.push_back(cln::recip(cln::cl_RA(6))); + next_r = 4; + } + for (unsigned p = next_r; p <= n; p += 2) { + cln::cl_I c = 1; + cln::cl_RA b = cln::cl_RA(1-p)/2; + unsigned p3 = p+3; + unsigned p2 = p+2; + unsigned pm = p-2; + unsigned i, k; + for (i=2, k=0; i <= pm; i += 2, k++) { + c = cln::exquo(c * ((p3 - i)*(p2 - i)), (i - 1)*i); + b = b + c * results[k]; + } + results.push_back(- b / (p + 1)); + next_r += 2; + } + return results[n/2 - 1]; }

2 3

End of a nightmare: patch for GiNaC's acinclude.m4
by Roberto Bagnara 15 Jan '02

15 Jan '02

A student of mine started having troubles with her Linux machine at home soon after she started playing with GiNaC. The machine seemed to become instable after a few days of work and she reinstalled GNU/Linux several times because of that. A few days ago, we discovered that /dev/null was no longer a character device: it was an ordinary file and that caused the system not to boot properly. We recreated the /dev/null device and everything seemed to go well until this morning, when the device disappeared again. BUT this time Tatiana provided the relevant bit of information; here it is: ... config.status: creating doc/reference/Makefile config.status: creating config.h **** The following problems have been detected by configure. **** Please check the messages below before running "make". **** (see the section 'Common Problems' in the INSTALL file) ** No suitable installed version of CLN could be found. deleting cache /dev/null [root@crystal GiNaC-1.0.3]# Aaargh!!! It was GiNaC configuration erasing /dev/null! Yeah, right, there is no need to run `configure' while being root. However, I feel in this case the price that had to be payed was a bit too high. That is why I propose the following patch is applied to `acinclude.m4': if `cache_file' is `/dev/null' do not delete it. All the best, Roberto -- Prof. Roberto Bagnara Computer Science Group Department of Mathematics, University of Parma, Italy http://www.cs.unipr.it/~bagnara/ mailto:bagnara@cs.unipr.it diff -rcp GiNaC-1.0.3.orig/acinclude.m4 GiNaC-1.0.3/acinclude.m4 *** GiNaC-1.0.3.orig/acinclude.m4 Tue Nov 20 18:42:05 2001 --- GiNaC-1.0.3/acinclude.m4 Tue Jan 15 09:56:33 2002 *************** if test "x${ginac_error}" = "xyes"; then *** 86,93 **** if test "x${ginac_warning_txt}" != "x"; then echo "${ginac_warning_txt}" fi ! echo "deleting cache ${cache_file}" ! rm -f $cache_file else if test x$ginac_warning = xyes; then echo "=== The following minor problems have been detected by configure." --- 86,95 ---- if test "x${ginac_warning_txt}" != "x"; then echo "${ginac_warning_txt}" fi ! if test "x$cache_file" != "x/dev/null" ! echo "deleting cache ${cache_file}" ! rm -f $cache_file ! fi else if test x$ginac_warning = xyes; then echo "=== The following minor problems have been detected by configure."

2 1

News from the gTybalt corner
by Stefan Weinzierl 15 Jan '02

15 Jan '02

Dear all, I wrote a library based on GiNaC which can handle the expansion of transcendental functions in a small parameter, like for example 3F2(-2eps,-2eps;1-eps,1-2eps,1-2eps;x) = 1 + 4 Li2(x) eps^2 + O(eps3) This is a domain where commercial CAS usually have no clue at all. The library is available under GPL from http://www.fis.unipr.it/~stefanw/nestedsums I should also say that at the moment you need GiNaC version 0.8.3, a migration to the actual GiNaC version is planned for the future. The rest of the mail is more for the developpers of GiNaC: First of all my thanks to the developpers of GiNaC, GiNaC proved to be quite a solid piece of software. For a part of my program I could do a benchmark test with a corresponding program written in FORM, and it turned out that the GiNaC/C++ code was roughly a factor two faster (as long as all expressions fitted into the available RAM). Since FORM is mainly known for its speed, this is quite an achievement. To the more technical points: The algorithms for the expansion of transcendental functions relies heavily on various algebras, e.g. you have two elements a1 and a2 in an algebra A and you can multiply them to get a third element: a1 * a2 = a3 Although the algebras are all commutative, I implemented them as non-commutative objects, then GiNaC automatically groups elements of the same algebra together and I only have to supply the actual multiplication routine in the method simplify_ncmul I think, this is handled elegantly and effciently in GiNaC. However, the expand function does not handle cases if not all types are the same. I can even blame myself for making the suggestion that the class add should throw an exception if it encounters an expression like 1 + A, where A is non-commutative. Right now one is supposed to write ONE_in_A + A, where ONE_in_A is the unit in the algebra of A. With more than one algebra this becomes rather ugly and inefficient. My suggestion would be therefore to think about a class mixed_type_add in GiNaC, which should have a similar relation with class add as class ncmul has with class mul: If all terms in a sum are of the same type, they will end up in class add, otherwise in this new class. I think, that can be implemented efficiently, such that users who do not care about non-commutative objects will not suffer any severe penalty. In addition pedantics can check at run-time their expressions and start a panic attack, if they encounter a class mixed_type_add. In short, one would have a container, where you can put an "apple" and a "potato" in. I had a look at the GiNaC source code, how one would do it, but since it involves quite a bit of cross links, I would not directly volunteer for it. What do you think ? Best wishes, Stefan

2 2

The specification of sqrfree()
by Roberto Bagnara 09 Jan '02

09 Jan '02

Hi there, we are developing a GiNaC-based recurrence relations solver. During this work we have found that GiNaC's documentation is not very precise about what a "square-free factorization" is. Below you find what we believe is a sensible definition (which also seems to be compatible with the current implementation). Please, check if that is also consistent with the specification of GiNaC (we would like to avoid relying on non-features that may disappear on a subsequent release). A polynomial p(X) in Q[X] is said <EM>square-free</EM> if, whenever two polynomials q(X) and r(X) in Q[X] are such that p(X) = q(X)^2*r(X), q(X) is constant. The sqrfree function computes polynomials p1(X), ..., pk(X) and positive integers e1, ..., ek such that p(X) = p1(X)^e1 * ... * pk(X)^ek and p1(X) * ... * pk(X) is square-free. If you agree with the above definition we are willing to provide Doxygen and TeXinfo patches against the current sources. All the best, the PURRS team -- Prof. Roberto Bagnara Computer Science Group Department of Mathematics, University of Parma, Italy http://www.cs.unipr.it/~bagnara/ mailto:bagnara@cs.unipr.it

2 1