- GiNaC-list - ginac.de

Integration
by Jason Dusek 08 Feb '05

08 Feb '05

Hi, I read that someone has recently added symbolic integration to GiNaC. Can it do this integrand: (x - 1)/(x*sqrt(x^3 + 1) or q / ||r|| dV ? Is it good enough for use, or should I wait a bit? -- _jason

2 1

Fwd: Performance problem with series
by Frieder Kalisch 03 Feb '05

03 Feb '05

Hi, here is a very odd performance behaviour of ginac/series: The operation ex e = 1/(a+sqrt(b*b+x*x)); e = e.series(x==1, 21); is more than a factor of 10^4(!) faster then outputting the result using cout << e << endl; I have measured these times using clocks(). You find the source of the program at the end of this email. What does this mean? This means, that series is very efficient (both time- and space-wise), because it reuses intermediate results heavily. On the other hand the result is nearly useless, because it is actually _very_ complicated, which shows up in any subsequent usage. In this specific example, a hand-written taylor-series using derivatives was 2 orders of magnitude slower (for 60th order) to compute the series, but the series (20th order) still was printed 1000 times faster than the standard series. In the light of this finding, something needs to be done. I suggest to essentially revert the logic of series, namely to compute series by differentiation and to fall back to an explicit composition if it does not work. Mathematica does it that way according to its documentation. At least ginac should provide this mode as an option, either by an additional taylor() method for ex or by means of series_options. Do You agree with this analysis and my conclusions? Are there any other ways to solve this problem? Regards Frieder Kalisch #include <ginac/ginac.h> #include <iostream> #include <ctime> using namespace GiNaC; using namespace std; int main(int argc, char **argv) { time_t t0, t1, t2; t0 = clock(); symbol x("x"), a("a"), b("b"); ex e = 1/(a+sqrt(b*b+x*x)); e = e.series(x==1, 21); t1 = clock(); cout << e << endl; t2 = clock(); double res = CLOCKS_PER_SEC; cerr << (t1-t0)/res << endl << (t2-t1)/res << endl; return 0; } // int main(int argc, char **argv) -- Frieder Kalisch Institut für theoretische Physik kalisch(a)tphys.uni-heidelberg.de Philosophenweg 19 +49-6221-549-433 D-69120 Heidelberg

3 4

Include Guards
by Jason Dusek 22 Jan '05

22 Jan '05

In the FAQ, they discuss a problem with multiple inclusions of a header where global symbols are defined: ================================================ Q: I have various modules, some global variables and... A: ...and you are getting weird results with symbols not being identical? Imagine this situation: A file globals.h declares symbol x("x"). Two other modules #include "globals.h", do something with symbol x and return an expression ex containing the symbol to the caller (e.g. to main()). When the caller combines the two expressions there are really two different symbols both with the same print-name x. This may cause terms like x-x to stay around which will never be simplified to 0 because the system has no knowledge that the two symbols are supposed to be the same. How should it? After all the declarations really happend in different compilation units! ================================================ Can include guards solve this problem? Do I *have* to make that evil factory initializer thingy? -- _jason

4 8

pyginac: numeric vs. ex and a failing test case
by Matti Peltomäki 17 Jan '05

17 Jan '05

Hi, Although being mostly about PyGiNaC I send this message on this list since the problem itself has also to do with GiNaC and may be interesting. I've been playing around with PyGiNaC and I accidentally found out, that running bin/check_matricies.py sometimes shows a failing test case, which looks like below (the test case seems to be the same as the one in GiNaC's check/check_matrices.cpp with the same name) -- snip -- Running __main__.funny_matrix_determinants.__doc__ Trying: funny_matrix_determinants() Expecting: 'Funny' matrix determinant test passed. ***************************************************************** Failure in example: funny_matrix_determinants() from line #1 of __main__.funny_matrix_determinants Expected: 'Funny' matrix determinant test passed. Got: Determinant of 5 x 5 matrix [[0,25/9+100/81*a,0,-24/35-6/5*c,24/35+6/5*c],[-5/4*c**(-1)*a*(c+3*a)*(b+a),0,2,0,-5/4*c**(-1)*a*(c+3*a)*(b+a)],[8*b**3*(2.7182818284590452354+1/2*c)*a**(-1),8/49*c*a**(-1),0,0,8*b**3*(2.7182818284590452354+1/2*c)*a**(-1)],[0,0,0,c**(-1)*b**2,-c**(-1)*b**2],[0,0,c*b**(-2),0,0]] was not found to vanish! 'Funny' matrix determinant test failed. -- snip -- I traced this back to sparse_tree() generating sometimes (sub)expressions of the form cos(1), i.e. a function with an argument of type numeric. Trying this out gives >>> from ginac import *; >>> cos(numeric(1)) <cginac.numeric object at 0xb7dda4f4> >>> print cos(numeric(1)) 0.5403023058681397174 >>> so that numerically evaluating the cosine leads to truncation errors and the determinant is no more exactly zero. Now, this evaluation seems to be exactly what GiNaC does if a numeric is given as a parameter to a function. OTOH, if the parameter is an ex, no evaluation occurs, even if the exs wraps a numeric. The following code #include <iostream> #include <ginac/ginac.h> using namespace GiNaC; using namespace std; int main(void) { cout << cos(numeric(1)) << endl; cout << cos(ex(numeric(1))) << endl; return 0; } gives 0.5403023058681397174 cos(1) This does not cause the same failure in the GiNaC test case since sparse_tree returns an ex. However, in PyGiNaC the return value is unwrapped from ex to whatever type it is, and with a numeric, the problem above occurs. So, in the end, I'm left with the situation that I don't know how to symbolically present something like cos(1) in PyGiNaC. Is it possible in some way not obvious to me? Best regards, Matti -- Matti Peltomaki, Undergraduate Research Assistant Laboratory of Physics Helsinki University of Technology 040 7576977, room Y308B

2 1

Python bindings to GiNaC
by Ondrej Certik 13 Jan '05

13 Jan '05

Hello all, I wrote new Python bindings to ginac, because PyGiNaC (http://cens.ioc.ee/projects/pyginac/) seems to be dead and I wasn't able to even compile it. Currently I only support classes: basic, symbol, numeric, ex. But more can be done quite easily, unfortunately I probably won't have time to continue on it. More info at: http://ondrej.certik.cz/pyginac.php Does anyone of you would like to use ginac in python, or does everybody use c++/cint? Yours, Ondrej Certik

4 8

Re: [GiNaC-list] Python bindings to GiNaC
by CRAIGHS＠Nationwide.com 04 Jan '05

04 Jan '05

C++ Ondrej Certik <ondrej(a)certik.cz> T To: ginac-list(a)zino.physik.uni-mainz.de Sent by: cc: ginac-list-admin(a)zino.phys ik.uni-mainz.de bcc: Subject: [GiNaC-list] Python bindings to GiNaC 01/03/2005 08:53 AM Please respond to ginac-list Hello all, I wrote new Python bindings to ginac, because PyGiNaC (http://cens.ioc.ee/projects/pyginac/) seems to be dead and I wasn't able to even compile it. Currently I only support classes: basic, symbol, numeric, ex. But more can be done quite easily, unfortunately I probably won't have time to continue on it. More info at: http://ondrej.certik.cz/pyginac.php Does anyone of you would like to use ginac in python, or does everybody use c++/cint? Yours, Ondrej Certik _______________________________________________ GiNaC-list mailing list GiNaC-list(a)ginac.de http://thep.physik.uni-mainz.de/mailman/listinfo/ginac-list

1 0

indexed
by Andrea Bedini 09 Dec '04

09 Dec '04

Hi, I'm new to ginac, I'd like to know if is it possible to unroll indexed objects. As example, I would like to do something like this a_i a~i = (a_1)^2 + (a_2)^2 + (a_3)^2 + ... or just substituting a_i values in a_i * b~i a_i = =(1,0,0) a_i b~i -> b~1 Thanks a lot for your work.

2 1

Strange behaviour of archive
by Phil Stephens 02 Dec '04

02 Dec '04

I have found the archiving procedure to behave strangely and was wondering if anyone has also seen this or if they know the solution. I am using GiNaC-1.2.1 with cln-1.1.7. This is all compiled with cygwin. What happens is the following. I have a database that I wish to store both the index and the referenced expression in the archive. For small archives the numerical factors on my expressions remain rational and the sqrt's remain un-evalf()'d. Unfortunately, it seems that when my archive gets large the rational expressions and the sqrt()'s evaluate to real numbers. I have done the following to confirm all of this. My software is using the same algorithm with no code changed but using two different models. The models dictate the size of the database. When I save and load with the smaller model, all the rational and sqrts are returned properly. When I do so with the larger model, they all return as real numbers. I have used a viewgar-like program to confirm that they have indeed been saved as real values, rather than the rationals/sqrts for the larger model, but the expressions I give to archive_ex have not been evalf'd. Can anyone explain why this may happen or if this is intentional in the design (speed of archiving for example)? It seems to occur between creating the archive and saving it to a file. It isn't vital to my software, just an aesthetic thing that would be useful for users. Just as an idea (since I don't really know how to 'reproduce' this behaviour) the gar file for two models that work properly are 25KB and 205KB but for the one that doesn't work right it is 1.1 MB. Thanks for any help!

1 0

how to use ginac output in gsl at runtime
by Matei Stroila 21 Nov '04

21 Nov '04

Hello, Here is what I would like to do and I need a piece of advice on this. I want to use ginac to compute the gradient of a function so that I can use this gradient for the conjugate gradient method in gsl (the GNU scientific library). I succesfully computed the gradient and I can output its expression with "csrc_double" to a string. How can I use this expression in the definiton of a function at runtime? I guess other people had this problem too, to use a CAS library in conjuction with a numeric library, at runtime. Thanks a lot for any suggestion, Matei

1 1

code: extended gcd
by Ralf Stephan 20 Nov '04

20 Nov '04

Hello, please include the following snippet. Its application domain depends on that of quo() which is univariate, currently. It is tested but not heavily so. It is not optimized. Thanks, ralf Input: x,y rational polynomials Computes u,v,d such that ux+vy = gcd(x,y). /////////// Extended GCD ex xgcd (const ex& a, const ex& b, ex& u, ex& v, const symbol& s) { const ex& x = a.normal(); const ex& y = b.normal(); if (is_exactly_a<numeric>(x) && is_exactly_a<numeric>(y) && x.info(info_flags::integer) && y.info(info_flags::integer)) { cln::cl_I cu, cv, cr; cr = cln::xgcd (cln::the<cln::cl_I>(ex_to<numeric>(x).to_cl_N()), cln::the<cln::cl_I>(ex_to<numeric>(y).to_cl_N()), &cu, &cv); u = numeric(cu); v = numeric(cv); return numeric(cr); } if (is_exactly_a<numeric>(y) && y.info(info_flags::integer)) { u = 0; v = numeric(1)/y; return 1; } if (is_exactly_a<numeric>(x) && x.info(info_flags::integer)) { v = 0; u = numeric(1)/x; return 1; } if (x.info(info_flags::rational_polynomial) && y.info(info_flags::rational_polynomial)) { ex d, p, q; d = gcd (x, y); int dx = degree(x, s), dy = degree(y, s); if (dx < dy) { p = y; q = x; } else { p = x; q = y; } ex u1,u2,u3,v1,v2,v3,w1,w2,w3; u1 = p; u2 = 1; u3 = 0; v1 = q; v2 = 0; v3 = 1; while (!v1.is_zero()) { ex Q = quo (u1, v1, s); w1 = u1 - Q*v1; w2 = u2 - Q*v2; w3 = u3 - Q*v3; u1 = v1.expand(); u2 = v2.expand(); u3 = v3.expand(); v1 = w1.expand(); v2 = w2.expand(); v3 = w3.expand(); } u = (u2 / ((u1/d).normal())).expand(); v = (u3 / ((u1/d).normal())).expand(); //cerr << "CHK: " << (u*x+v*y).expand() << endl << flush; return d; } throw invalid_argument("xgcd(): not a rational polynomial"); }

4 11