- GiNaC-list - ginac.de

FORTRAN output
by Warren Weckesser 28 Mar '07

28 Mar '07

I am working on a program that uses ginac to generate source code. One part of the program generates fortran-90 code. Currently I am experimenting with using the csrc output format, but with the print function of the power class redefined as in section 6.3.1 of the online tutorial. Is there a better way? Has anyone implemented a full-blown fortran output format? Best regards, Warren

2 1

Bug with evalm()?
by Martin Sandve Alnæs 26 Mar '07

26 Mar '07

I use ginac 1.3.4. I believe there is a bug in the implementation of evalm(). Its behaviour should be to do nothing if it is not applicable to the ex it is called on, right? The minimal example below throws an exception: #include <ginac/ginac.h> using namespace std; using namespace GiNaC; int main(int argc, char **argv) { ex x = symbol("x"); ex e = sin(x); ex s = e.series(x == 0, 3); cout << "This will crash:" << endl; ex m = s.evalm(); return 0; } The error is: terminate called after throwing an instance of 'std::range_error' what(): basic::let_op(): pseries has no operands Aborted (core dumped) The call stack from gdb is: (gdb) where #0 0xffffe410 in __kernel_vsyscall () #1 0xb7b12770 in raise () from /lib/tls/i686/cmov/libc.so.6 #2 0xb7b13ef3 in abort () from /lib/tls/i686/cmov/libc.so.6 #3 0xb7d02520 in __gnu_cxx::__verbose_terminate_handler () from /usr/lib/libstdc++.so.6 #4 0xb7cfff55 in std::set_unexpected () from /usr/lib/libstdc++.so.6 #5 0xb7cfff92 in std::terminate () from /usr/lib/libstdc++.so.6 #6 0xb7d000ca in __cxa_throw () from /usr/lib/libstdc++.so.6 #7 0xb7d9a4a6 in GiNaC::basic::let_op (this=0x80597d8, i=1) at basic.cpp:256 #8 0xb7d96dfd in GiNaC::basic::map (this=0x80595e0, f=@0xb7f8e824) at basic.cpp:316 #9 0xb7d962a9 in GiNaC::basic::evalm (this=0x80595e0) at basic.cpp:482 #10 0x0804b47a in GiNaC::ex::evalm (this=0xbfd820d4) at /usr/include/ginac/ex.h:118 #11 0x0804b214 in main () at test.cpp:12 Of course, it doesn't make much sense to call evalm() on this expression, but it is a problem in an application where we want to call evalm() "just in case" on all ex objects that pass a certain point in the application. Thus a "do nothing" behaviour for evalm() is of great importance. martin

6 20

Bug with matrix subtraction
by Martin Sandve Alnæs 26 Mar '07

26 Mar '07

If I subtract two equal matrices, I get the scalar 0. This messes up later calculations, since a scalar doesn't have op(i), transpose(), etc. If I do a+(-b) instead of a-b, I get a zero matrix as wanted. Note that the code is in Python using the swiginac bindings, which now calls evalm() automatically between operations. In [32]: m = matrix(2,2,[x,0,0,0]) In [33]: m-m Out[33]: 0 In [34]: m+(-m) Out[34]: [[0,0],[0,0]] So we have the silly relation that "a-b != a+(-b)", which can only be called a bug. -- --- Martin Sandve Alnæs PhD student, Cardiac Computations Simula Research Laboratory, Norway

3 4

multi index storage of expressions
by Charlls Quarra 21 Feb '07

21 Feb '07

Hi, I wanted to define a simple matrix of derivatives (i.e: a Jacobian) first as an indexed object: ex f = a*x + b*y + c*x*y + ... ex g = u*pow(x,2) + v*pow(y,2) + w*x,y + .... idx i(symbol("is") , 2) , j(symbol("js") , 2); ex Jacobian = diff( f*delta_tensor(i , 0) + g*delta_tensor(i,1) , x )*delta_tensor( 0 , j ) + diff( f*delta_tensor(i , 0) + g*delta_tensor(i,1) , y )*delta_tensor( 1 , j ); but then this failed on me (because 0 and 1 are not acceptable arguments where an idx should be expected) Why am i trying to do this as an indexed object and not with a matrix? well, because later, i would also want to make computations with an object like D( f_{i} ) / dx_{j} dx_{k} but afaik there isnt a "three indexed" matrix object in ginac, or how they are usually called, cubic arrays. correct? __________________________________________________ Preguntá. Respondé. Descubrí. Todo lo que querías saber, y lo que ni imaginabas, está en Yahoo! Respuestas (Beta). ¡Probalo ya! http://www.yahoo.com.ar/respuestas

4 14

Re: [GiNaC-list] table and exhashmap
by Jerome BENOIT 19 Feb '07

19 Feb '07

Hi, Sheplyakov Alexei wrote: > On Mon, Feb 19, 2007 at 03:45:03PM +0800, Jerome BENOIT wrote: > >> Let be more specific. >> >> During some kind of computation I have to compute (multi-indexed) >> coefficients: >> 1] these coefficients are long to compute; >> 2] they are computed by sub set through (consuming) recurrence formulae; >> 3] these coefficients are meant to be used randomly and very often. > >> Clearly the remember machinery can be use to overcome point [1] and [3]: >> for GiNaC function this is transparent for the user. >> >> Point [2] causes me some difficulty: >> I want to compute each subset only once. > >> Hence the idea to compute each subset through a GiNaC function >> (with the remember option) which return tables, > > Why do you need _GiNaC_ function to compute coefficients? What's wrong > with C++ function, e.g. > > ex compute(const ex& data) > { > static exmap remember; > exmap::iterator i = remember.find(data); > if (i != remember.end()) > return i->second; > > ex result; > // actually compute it, making use of compute() itself > > remember.insert(exmap::value_type(data, result)); > return result; > } > >> Right now, my guess is that I have to introduce a super class >> with a exhashmap<SOMETHING> as (principal) member: I thought >> that there was a better GiNaC way to do so, but apparently not. > > Unless I'm missing something fundamental, this approach seems to be > unnecessary complicated. I am more than agree: hence my question. I guess that I still too Maple oriented. Thanks you very much, Jerome

1 0

table and exhashmap
by Jerome BENOIT 19 Feb '07

19 Feb '07

Hello List, I am new in GiNaC and in C++ (but very familiar with another CAS and C). I would like to implement a function which returns table. Unless I missed something, there is no ready to use table object in GiNaC. My first thought was to use exhashmap<T> , but I am no more sure about it as I do not see how to proceed. Any hint is welcomes. Thanks in advance, Jerome -- Jerome BENOIT jgmbenoit_at_mailsnare_dot_net

2 5

list of list
by Jerome BENOIT 17 Feb '07

17 Feb '07

Hello List, I am rather a newbie in GiNaC and in C++ (but I am very familiar with another CAS and with C): I have just noticed that: lst list; list = lst(lst(1,1),lst(2,4),lst(3,9)); gives a list of list as expected, whereas lst list; list = lst(1,1),lst(2,4),lst(3,9); gives the first list (namely lst(1,1)) is it a bug ? or is it a wrong thing to do ? Thanks in advance, Jerome

3 6

Release of GiNaC 1.3.7
by Jens Vollinga 05 Feb '07

05 Feb '07

Hello everybody, GiNaC 1.3.7 is out and available. The changes are: - Fixed bug in expansion of power. - Fixed bugs in functions S() (infinite loop), H() and zeta() (wrong results). - Rebuilt bison related files with bison version 2.3 (fixes parse error bugs). - Improved power::eval(). - Improved libreadline detection. As always, this release can be downloaded from http://www.ginac.de/Download.html Best wishes, Jens

1 0

Teaching material for GiNaC
by Chris Dams 02 Feb '07

02 Feb '07

Dear all, For the C++ course that has been held at my university, I have done two of the lectures and talked about GiNaC. For use by others in similar circumstances, you find my teaching materials attached to this course. Labeling the sections as "First hour" and "Second hour" turned out to be a bit optimistic. I used two hours for both the first and second section. Mainly because the students needed quite a bit of time for the final exercise in each section. Best wishes, Chris

1 0

collecting on similar series expansion
by Charlls Quarra 01 Feb '07

01 Feb '07

Hi , I'm just trying out ginac for the first time. I have an expression with sumands of this type (in latex format): 1-\frac{1}{3} \frac{ p D[0](F)(q) D[0,0](F)(q) \nu \epsilon^{2} dt^{5} \delta}{m^{3}}-\frac{1}{9} \frac{ p D[0](F)(q) D[0,0](F)(q)^{2} \alpha^{2} \nu^{2} \epsilon^{2} dt^{9} \F(q) \delta \gamma}{m^{5}}-\frac{1}{18} \frac{ p^{2} D[0,0](F)(q) D[0,0,0](F)(q) \beta \alpha^{2} \epsilon^{2} dt^{8} \F(q) \delta \gamma}{m^{5}}+\frac{1}{2} \frac{ D[0](F)(q) D[0,0](F)(q) \alpha \nu \epsilon dt^{6} \F(q)}{m^{3}} .... so i want to collect terms that have the same coefficients of the form D[0..i times] (F)(q) * D[ 0.. j times] (F)(q) .... but the polynomial collect won't work. Initially i thought about substituting each D[ ] (F)(q) with a different symbol, then collect on the multivariate polynomials of these symbols and then finally substitute back the D[] expressions, however, I don't know how the collect call is supposed to receive the symbols that it must collect, the code looks like this: ex jacobian = diff( q + dq_new*dt , q )*diff( p + dp_new*dt , p ) - diff( q + dq_new*dt , p )*diff( p + dp_new*dt , q ); cout << " the jacobian of the integrator is " << expand(jacobian) << endl; exmap diff_map , inv_diff_map; symbol foo[10]; diff_map[ F(q) ] = foo[0]; inv_diff_map[ foo[0] ] = F(q); diff_map[ F(q).diff(q) ] = foo[1]; inv_diff_map[ foo[1] ] = F(q).diff(q); diff_map[ F(q).diff(q,2) ] = foo[2]; inv_diff_map[ foo[2] ] = F(q).diff(q,2); diff_map[ F(q).diff(q,3) ] = foo[3]; inv_diff_map[ foo[3] ] = F(q).diff(q,3); diff_map[ F(q).diff(q,4) ] = foo[4]; inv_diff_map[ foo[4] ] = F(q).diff(q,4); ex rord = collect( jacobian.subs( diff_map ) , { foo[0] , foo[1] , foo[2] , foo[3] , foo[4] } ); // i know this is an error in C++, but this is the idea rord = rord.subs( inv_diff_map ); cout << "reordered we obtain " << rord << endl; any thoughts on this? __________________________________________________ Preguntá. Respondé. Descubrí. Todo lo que querías saber, y lo que ni imaginabas, está en Yahoo! Respuestas (Beta). ¡Probalo ya! http://www.yahoo.com.ar/respuestas

2 4