- GiNaC-list - ginac.de

[patch] Allow exprseq in the arguments to lsolve()
by Vitaly Magerya 31 May '18

31 May '18

Hi, folks. Currently lsolve(eqns, vars) takes two lists as arguments, but then only accesses them via op(), so essentially any other container should work equally well. Here's a small patch to make it support exprseq in addition to lst. The motivation is that std::vector<ex>, upon which exprseq is based, can be more convenient to work with than std::list<ex>. Note that there's one additional places where lst is required, but exprseq would work too: ex::subs(). It would make sense to support exprseq there too. Also note that testing for lst is done via two different patterns in the source code: sometimes it is 'info(info_flags::list)', other times it is 'is_a<lst>'. It's not clear to me why this difference exists, or what should be the preferred way.

2 1

[patch] Add solve_algo::markowitz for large sparse systems
by Vitaly Magerya 29 May '18

29 May '18

Hi, folks. I've been trying to apply GiNaC to solving systems of sparse linear equations: it is currently too slow, and I have an improvement. The setup of the issue is this: I have a system of 2279 equations in 749 variables, which you can find in 'sys749.txt' at [1]. When I give this system directly to lsolve(), like this: // compile with: // c++ -std=c++11 -o sys749 sys749.cpp -lginac -lcln -ldl #include <ginac/ginac.h> #include <ginac/parser.h> #include <fstream> using namespace GiNaC; using namespace std; int main(int argc, char *argv[]) { symtab table; parser reader(table); ifstream i("sys749.txt"); ex data = reader(i); lst eqns; for (auto &&e : data.op(0)) { if (!e.is_zero()) eqns.append(expand(e) == 0); } cout << lsolve(eqns, data.op(1), solve_algo::gauss); } ... after ~5 hours of time and 1.5GB of RAM, I conclude that this will never terminate. In contrast, it takes ~2.5 minutes for Maple to solve this system, so I know we can do better. The main slowdown for GiNaC comes from the classical 'fill-in' problem during Gaussian elimination. Consider a sparse matrix like this: [* * * * *] [* * 0 0 0] [* 0 * 0 0] [* 0 0 * 0] [* 0 0 0 *] When you apply naive Gaussian elimination to this matrix, on the first step you will subtract the first row from the rest of the matrix -- this will make all the zero entries non-zero (this is the "fill-in" of zeros), thus making the matrix dense, and killing all the performance because each subsequent row elimination step will require N normal() calls. On the other hand, if you where to permute the matrix like so: [* 0 0 0 *] [0 * 0 0 *] [0 0 * 0 *] [0 0 0 * *] [* * * * *] ... then naive Gaussian elimination would fully preserve the sparsity structure, and at each step most of the operations you'll need to perform would amount to subtracting zeros, and could be optimized away (so, nothing to normal()ize). This sort of permutation to reduce fill-in is an essential part of solving sparse linear systems, and is the main reason for Maple's performance, I think. In fact, if you'll just sort the list of equations by the number of terms before passing that list into lsolve(), this already allows GiNaC to solve sys749.txt, in under an hour (don't remember the exact timing). We can do better than this, though. There's an overview of the current state of the art in this sort of techniques at [2]. Out of all the tricks described there, the most straightforward to me is the Markowitz ordering, which is a heuristic that says: at each step of Gaussian elimination, choose the pivot (i,j) that minimizes the amount of fill-in you'll need to perform in this step, approximated as (r[i]-1)*(c[j]-1), where r[i] is the number of non-zero entries in row i, and c[j] is the number of non-zero entries in column j. This heuristic is expensive to calculate, but it is well worth it for sys749. The patch I'm attaching implements Markowitz-ordered Gaussian elimination as a new solve_algo::markowitz. An additional part of the slowness is the part of lsolve() that constructs the coefficient matrix: it does so by calling coeff() repeatedly for every variable -- this alone takes about 10 minutes for sys749. To optimize this part I first obtain the full set of all symbols in an equation, and then call coeff() only if the symbol in question belongs to that set. This changes 10 minutes into 10 seconds. Still, I think there's more room for improvement here: it should be possible to collect the whole coefficient row in one pass over the equation, without calling coeff() multiple times. Since it takes about 1.3 seconds to read the whole expression from disk, it should take no more than that to determine the full coefficient matrix, I think. If you know how to achieve that though, I'd be interested in learning. In summary, with the proposed optimizations GiNaC can solve sys749 in 60-90 seconds (differs every time; at one point I was getting a consistent 45 seconds, but I've probably changed something since then), which is actually faster than Maple, at least for this particular system. I would appreciate some correctness testing on this new algorithm, because Markowitz ordering requires full pivoting, and that is trickier to get right than the partial pivoting which GiNaC has had used until now, especially when 'vars' in matrix::solve() is wider than n*1. I've modified the check_lsolve test to also test this algorithm, but it doesn't look comprehensive. Style changes or algorithmic suggestions are also welcome. [1] http://tx97.net/~magv/sys749.txt [2] http://faculty.cse.tamu.edu/davis/publications_files/survey_tech_report.pdf

1 0

psi (digamma) solving
by Matej Usaj 06 Mar '18

06 Mar '18

Hi All, I'm trying to solve function that contains a call to psi (digamma) but my call to fsolve fails every time, ie: > f=psi(x/2)==0; psi(1/2*x)==0 > fsolve(f,x,2.7,3.1); fsolve(): function does not evaluate numerically but I'd expect the result to be x~=2.923... Am I missing something here? Thanks for your help! Cheers, Matej

2 1

previous release archive
by Eric Deveaud 22 Feb '18

22 Feb '18

Hello, we've setup few days ago some tools using ginac-1.73 as I have to reproduce this using the same version, I would like to know if you provide release archive for old versions or should I go grabbing the git snapshot tag for the release 1.7.3 regards Eric

2 1

Re: [GiNaC-list] Problem of installing ginac on mac os
by Goutam Das 02 Feb '18

02 Feb '18

@Vladimir, thanks but I do not prefer to install Homebrew or MacPorts as most of the time they affects due to broken packages and also much bulky. So I install whatever needed for my work from the source itself. In this case I have installed cln successfully but failed in installing Ginac. @Isuru, after configure, in the Makefile I can see the following: CPPFLAGS = /Users/gkdgoutam/Softwares/HEP_Softwares/Install/Cln134/include: CXXFLAGS = -g -O2

1 3

Release 1.7.3
by Richard B. Kreckel 01 Feb '18

01 Feb '18

Hi, A new GiNaC version has just been released. These are the changes: * Speed up some corner cases of polynomial factorization. * Speed up some corner cases of gcd computation. * Speed up factorization of pre-factored polynomials. * Make matrix::solve() work with hidden (non-normal) zeros. * Fix bug in evaluation of real and imaginary parts of powers. * Fix some LaTeX output, add new function symbol::get_TeX_name(). The following people have contributed patches to this release: - Vladimir V. Kisil - Richard B. Kreckel - Vitaly Magerya Happy hacking! -richy.

1 0

matrix::solve()-related problems
by Vitaly Magerya 01 Feb '18

01 Feb '18

Hi, folks. I've run into a number of problems with matrix solution routines in GiNaC, and I'd like to sort them out one by one. The basic overview of the problems is this: 1) matrix::gauss_elimination checks for zero entries via a simple .is_zero() call without a .normal() call before, so it misses zeros in non-normal form. For example it will not detect that this expression: 1/(x - 1) + 1/(x + 1) - 2*x/(x*x - 1) ... is zero. This leads to "division by zero" errors. Test case: #include <ginac/ginac.h> using namespace GiNaC; int main() { symbol x("x"); matrix mx { {1,0,0}, {0,1/(x+1)-1/(x-1)+2/(x*x-1),1}, {0,1,2}, }; matrix zero(mx.rows(), 1); matrix tmp(mx.cols(), 1); for (unsigned i = 0; i < tmp.nops(); i++) { tmp.let_op(i) = symbol(); } mx.solve(tmp, zero, solve_algo::gauss); return 0; } 2) matrix matrix::solve does the same erorr, so if gauss_elimination returned a matrix with un-normal zeros, it too may error out with "division by zero". Test case: #include <ginac/ginac.h> using namespace GiNaC; int main() { symbol x("x"); matrix mx { {1,0,0}, {0,1/(x+1)-1/(x-1)+2/(x*x-1),1}, {0,0,0}, }; matrix zero(mx.rows(), 1); matrix tmp(mx.cols(), 1); for (unsigned i = 0; i < tmp.nops(); i++) { tmp.let_op(i) = symbol(); } mx.solve(tmp, zero, solve_algo::gauss); return 0; } 3) The default echelon form algorithm, Bareiss elimination, (matrix::fraction_free_elimination) is unbearably slow for sparse matrices, and should never be used with them. This is because of explosive common factor growth during its evaluation. For an illustration, using matrix::fraction_free_elimination on a matrix like this: {{a11,a12,a13,a14}, { 0,a22,a23,a24}, { 0, 0,a33,a34}, { 0, 0, 0,a44}} ... will give you a result like this: {{a11, a12, a13, a14}, { 0,a11*a22, a11*a23, a11*a24}, { 0, 0,a11*a22*a33, a11*a22*a34}, { 0, 0, 0,a11*a22*a33*a44}} Now imagine a similar 100x100 matrix with slightly more complex diagonal elements. Now, there are multiple ways of dealing with these problems. The practical workaround for the moment is to: a) use solve_algo::gauss for all sparse-ish matrices b) normal()ize the whole matrix before passing it into solve() There are places where this workaround is not available. One of such places is matrix::inverse(): it doesn't provide an 'algo' argument, and always defaults to solve_algo::automatic (which means fraction_free_elimination). So for the start I propose to add that argument (see the attached patch). The reservation I have with this patch is other functions (solve() and determinant()) have their own enumerations for the 'algo' argument, while inverse() will be reusing the one from solve(). Not the cleanest thing in the world from API usage point of view. What do you think of it?

4 17

Problem of installing ginac on mac os
by Goutam Das 29 Jan '18

29 Jan '18

I am trying to install ginac on mac os 10.13.2. The Cln1.3.4 is installed without error. I have tried to installed both way: 1. ./configure works with showing === The following minor problems have been detected by configure. === Please check the messages below before running "make". === (see the section 'Common Problems' in the INSTALL file) == ginsh will not compile, because readline headers could not be found. == Documentation may not build, because fig2dev/texinfo is missing. == Reference manual will not build, because doxygen is missing. Configuration of GiNaC 1.7.2 done. Now type "make". These above warning I am ignoring. Next, while doing make I am getting the following error: Making all in ginac /Applications/Xcode.app/Contents/Developer/usr/bin/make all-am /bin/sh ../libtool --tag=CXX --mode=compile g++ -std=gnu++11 -DHAVE_CONFIG_H -I. -I../config -DLIBEXECDIR='"/usr/local/libexec/"' /Users/gkdgoutam/Softwares/HEP_Softwares/Install/Cln134/include/cln: -g -O2 -MT libginac_la-add.lo -MD -MP -MF .deps/libginac_la-add.Tpo -c -o libginac_la-add.lo `test -f 'add.cpp' || echo './'`add.cpp libtool: compile: g++ -std=gnu++11 -DHAVE_CONFIG_H -I. -I../config -DLIBEXECDIR=\"/usr/local/libexec/\" /Users/gkdgoutam/Softwares/HEP_Softwares/Install/Cln134/include/cln: -g -O2 -MT libginac_la-add.lo -MD -MP -MF .deps/libginac_la-add.Tpo -c add.cpp -fno-common -DPIC -o .libs/libginac_la-add.o clang: *error: *no such file or directory: '/Users/gkdgoutam/Softwares/HEP_Softwares/Install/Cln134/include/cln:' make[2]: *** [libginac_la-add.lo] Error 1 make[1]: *** [all] Error 2 make: *** [all-recursive] Error 1 2. I tried to compile with cmake, however there I am getting the following problem: cmake ../ginac-1.7.2 -- Checking for module 'cln' -- No package 'cln' found -- found CLN [/Users/gkdgoutam/Softwares/HEP_Softwares/Install/Cln134/include/cln/version.h], version 1.3.4 CMake Error at /Applications/CMake.app/Contents/share/cmake-3.8/Modules/FindPackageHandleStandardArgs.cmake:137 (message): Could NOT find CLN (missing: CLN_LIBRARIES) (found suitable version "1.3.4", minimum required is "1.2.2") Call Stack (most recent call first): /Applications/CMake.app/Contents/share/cmake-3.8/Modules/FindPackageHandleStandardArgs.cmake:377 (_FPHSA_FAILURE_MESSAGE) cmake/modules/FindCLN.cmake:116 (FIND_PACKAGE_HANDLE_STANDARD_ARGS) CMakeLists.txt:30 (find_package) -- Configuring incomplete, errors occurred! See also "/Users/gkdgoutam/Softwares/HEP_Softwares/Source/ginac_build/CMakeFiles/CMakeOutput.log". In both cases, the problem seems to come from cln, but could not able to fix. Is it possible to get some help? I am using : gcc --version Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.13.sdk/usr/include/c++/4.2.1 Apple LLVM version 9.0.0 (clang-900.0.39.2) Target: x86_64-apple-darwin17.3.0 Thread model: posix InstalledDir: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin

3 3

Bug? Segmentation while inverting a matrix
by Patrick Schulz 09 Jan '18

09 Jan '18

Hi, i'm getting a segmentation fault while inverting a matrix (the matrix itself is invertable, i checked with the symbolic package for octave). Can anyone give me some pointers why this is happening? Is the matrix ill-formed in some way or is this a bug? Many thanks in advance! I'm running archlinux (Kernel 4.14.12), my compiler is g++ 7.2.1, the GiNaC version i'm using is 1.7.2. Attached is a sample listing which produces the error (at least on my system). Greetings, Patrick

4 5

factor() hangs on a simple expression, randomly
by Vitaly Magerya 28 Dec '17

28 Dec '17

Hi, folks. I have this problem: $ ginsh ginsh - GiNaC Interactive Shell (GiNaC V1.7.2) __, _______ Copyright (C) 1999-2017 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. > factor((2*x-l)^2*(-1+6*x-l)*(-1+3*x-l)*(-2+4*x-l)*(1+4*x-l)*(4*x-l)*(2+4*x-l)); At this point ginsh either hangs (80% of cases), or gives the correct result (20% of cases). As you can see the expression is already factored, so a hang is pretty unexpected here. Note that if you'll drop two of the factors from that expression, GiNaC won't hang anymore. Here are the cumulative timings you'll typically see via 'perf top' during this hang: 91.54% libginac.so.6.0.2 [.] GiNaC::(anonymous namespace)::factor_multivariate 81.32% libginac.so.6.0.2 [.] GiNaC::(anonymous namespace)::hensel_multivar 79.67% libginac.so.6.0.2 [.] GiNaC::(anonymous namespace)::multivar_diophant 79.67% libginac.so.6.0.2 [.] GiNaC::(anonymous namespace)::univar_diophant 79.67% libginac.so.6.0.2 [.] GiNaC::(anonymous namespace)::multiterm_eea_lift 72.77% libginac.so.6.0.2 [.] GiNaC::(anonymous namespace)::factor_sqrfree 66.18% libginac.so.6.0.2 [.] GiNaC::(anonymous namespace)::eea_lift 64.28% libginac.so.6.0.2 [.] GiNaC::factor 26.99% libginac.so.6.0.2 [.] GiNaC::(anonymous namespace)::operator* 23.67% ginsh [.] yyparse 23.67% ginsh [.] f_factor 16.91% libginac.so.6.0.2 [.] cln::cl_MI::~cl_MI 15.26% libginac.so.6.0.2 [.] cln::cl_I::~cl_I 13.66% libginac.so.6.0.2 [.] cln::cl_MI::cl_MI 13.55% libginac.so.6.0.2 [.] cln::cl_RA::~cl_RA So, any ideas on the nature of this problem, and how to fix it?

2 9