- GiNaC-list - ginac.de

small GiNaC query
by keith.briggs＠bt.com 27 Jan '04

27 Jan '04

How do I get the program below to compute 1/6*s.1^3+1/2*s.1*s.2+1/3*s.3 ? I don't care about ordering of summands, but I need the s.1*s.2 and s.2*s.1 terms to be merged and s.1*s.1*s.1 to be converted to s.1^3. Thanks, Keith // KMB 2004 Jan 26 // g++ -Wall tryginac1.cc -lginac && a.out #include <iostream> #include <stdexcept> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; ex Z(int n) { if (n==0) return 1; ex r=0; symbol s("s"),i_sym("i"); idx i(i_sym,n); for (int j=1; j<=n; j++) { ex si=indexed(s,i); r+=si.subs(i==j)*Z(n-j); } return r.expand().collect(s)/n; // -> 1/6*s.2*s.1+1/3*s.1*s.2+1/3*s.3+1/6*s.1*s.1*s.1 return simplify_indexed(r/n); // -> 1/6*s.2*s.1+1/3*s.1*s.2+1/3*s.3+1/6*s.1*s.1*s.1 return collect_common_factors(r/n); // -> 1/6*s.1*(s.2+s.1*s.1)+1/3*s.3+1/3*s.1*s.2 return r.expand()/n; // -> 1/6*s.1*s.1*s.1+1/6*s.1*s.2+1/3*s.3+1/3*s.1*s.2 } int main(void) { try { cout<<Z(3)<<endl; } catch (exception &p) { cerr << p.what() << endl; return 1; } return 0; } Dr. Keith M. Briggs Senior Mathematician, Complexity Research, BT Exact http://more.btexact.com/people/briggsk2/ (internet) http://research.btexact.com/teralab/keithbriggs.html (BT intranet) phone: +44(0)1473 work: 641 911 home: 610 517 fax: 642 161 mail: Keith Briggs, Polaris 134, Adastral Park, Martlesham, Suffolk IP5 3RE, UK

2 1

Physical Units
by Lee Harding 23 Jan '04

23 Jan '04

I'm an engineer of the physical ilk (MSME), but semi-literate in the art of C++ programming. After discovering GiNaC today, I thought "nice library ". I've created a few programs for optimizing things, doing simple FEA (finite element analysis) and other sorts of calculations. Invariably, I end up restricting the domain to a consistent set of physical units (e.g. metric meter, kg, sec) or writing a lot of code to compensate for varying units across the input. For general-purpose tools, it can be a real pain in the shorts to get right. I haven't thought a great deal about why math libraries so consistently ignore the problem of unit consistency and conversion. But, I have failed to do it successfully with other libraries. For example, the Gnu project citrus[1] is a unit conversion library that implements it's own expression concept. Units and expressions are inextricably tied, I can't do my job without the old reliable kg/s/m^2 (aka. Newtons, aka. energy) after all. And, what is sin(pi)? Is that degrees or radians? It seems worthwhile to explore integrating physical units to GiNaC, if only to make like easier on poor engineers like myself. Is that a reasonable possibility? [1] http://citrus.sourceforge.net/

5 4

Re: [Fwd: Re: Physical Units]
by Bernd Speiser 23 Jan '04

23 Jan '04

Kai Ludwig wrote: > > -------- Original Message -------- > Subject: Re: Physical Units > From: "Kai Ludwig" <kai.ludwig(a)uni-tuebingen.de> > Date: Fri, January 23, 2004 9:10 am > To: <ginac-list(a)zino.physik.uni-mainz.de> > > > I've created a few programs for optimizing things, doing simple FEA > > (finite element analysis) and other sorts of calculations. Invariably, > > I end up restricting the domain to a consistent set of physical units > > (e.g. metric meter, kg, sec) or writing a lot of code to compensate for > > varying units across the input. For general-purpose tools, it can be a > > real pain in the shorts to get right. Indeed, omitting units (and, for that matter, dimensions!) can lead to a lot of problems. By the way: in my opinion, sin(pi) is a number without dimension and unit (dimensionless number). However, it is interesting what unit the argument of the sin function, i.e. the pi, should have. Of course, this must be dimensionless. It could be a pure number, but it could equally well have a unit. It could result from some calculation with dimensioned quantities, which result in a dimensionless one. The Quantity library mentioned below does solve this problem by (a) checking the fact that the argument to a sin (and of course other such functions) to be dimensionless, (b) providing a type PureNumber, which can be used if a dimensionless quantity does not have any unit, and (c) if you want to work with an angle a type PlaneAngle (the name is according to IUPAC standards, the so-call Green Book), which can have units of Radian, Degree (not to be mixed up with temperature!!), Minute, Second, and Grade (all according to the Green Book). The Radian is the SI unit. Some general comment on Quantity follows below. > > Have a look at, > > http://sourceforge.net/projects/quantity > > this is a C++ template library designed for > unit conversions and unit type checking (already > at compile time). > As far as I know, there will be a new (redesigned) > release in due course. > The current version 0.6 of Quantity is still somehow preliminary. In particular, a lot of checking is still done at run time, resulting in considerable overhead. I am currently working on a thorough re-design where most of the functionality is shifted to compile time, and the type design is much more secure. I am not really finished, since I learned about static template metaprogramming (which is used for this re-design in many instances) only recently. This is not simple. I used some ideas of W. Brown (http://www.oonumerics.org/tmpw01/brown.pdf) and the possibilities of A. Alexandrescu's Typelist concept (A. Alexandrescu, Modern C++ Design. Generic Programming and Design Patterns Applied, Addison-Wesley, Boston, 2001). Unfortunately, the re-design also necessitates a slight re-design of the interface, so version 0.6 will not be fully compatible with the upcoming, more recent versions. If anyone is interested in the library, please let me know and/or check the announcements and news on the sourceforge page of the project. regards Bernd -- =================================================================== Bernd Speiser Institut f"ur Organische Chemie Auf der Morgenstelle 18 D-72076 T"ubingen Germany phone: +49-7071-2976205 (office) +49-7071-2976242 (laboratory) fax: +49-7071-295518 (new!) e-mail: bernd.speiser(a)uni-tuebingen.de Internet: http://www.uni-tuebingen.de/speiser ===================================================================

1 0

GiNaC 1.1.6 available!
by Jens Vollinga 22 Jan '04

22 Jan '04

Hello! GiNaC 1.1.6 "Sho Chiku Bai" is now available! It brings to you dummy functions and several bug fixes. In detail: * Added a function option "dummy()" which means "no options". This simplifies the implementation of symbolic functions which are not to be further evaluated. * Removed a bug in the numerical evaluation of Li() that caused the system to hang for certain parameter combinations. * Fixed a bug in the calculation of hash values for indices that could lead to wrong results or bogus error messages from simplify_indexed(). * Fixed a bug in the evaluation of harmonic polylogarithms for complex arguments with positive imaginary part. So we have GiNaC-Team 3 points, cimex lectularius 0 points. It's where you expect it to be: ftp://ftpthep.physik.uni-mainz.de/pub/GiNaC/ Bye, Jens

1 0

Symbols that are implicitly dependent on other symbols
by Johan Rutgeerts 16 Jan '04

16 Jan '04

Hi, Is there some concept in Ginac of "symbols which are function of some other symbol, but without making this explicit"? I'll try to explain what I mean: I wanted to: - define some function f(q_1, q_2, ... q_m) - calculate its derivatives df/d(q_j) symbolically and evaluate for specific values of q_j. But just calculating the derivatives and substituting values for q_j with subs(q_j == value) is far too slow for my application. However, I can rewrite f: f(q_1, q_2,... ) == g(r_1, r_2,.., r_n), in witch the r_i are functions of q_1, q_2,...,q_m. And I have numeric methods to calculate the value of r_i and its derivatives d(r_i)/d(q_j) for specific values of q_j. So I would like to calculate df/d(q_j) as: df/d(q_j) == sum( dg/d(r_i) * d(r_i)/d(q_j), for all i), using symbolic differentiation for the dg/d(r_i) and using the numerically calculated values of the r_i and d(r_i)/d(q_j) in the evaluation (which should go a lot faster since g(r_i) is a fairly simple function). So for this I'd need to define g as a function of r_1 to r_n, but then the r_i must somehow know that they are function of q_1 to q_m, otherwise the symbolic derivative dg/d(q_j) is zero. So my question is: is there an easy way to do this in Ginac? What I have an in mind as a possible solution is some kind of symbol (let's call it "MySymbol" to be original:) that has awareness of its numerical value and the value of its derivatives to certain variables. This would then allow me to take the derivative and would eliminate the use of subs(), making evaluation faster. The class MySymbol would have methods like: void setValue(double value); void setDerivative(symbol var, unsigned int order, double value); so that you can write: symbol q_j; MySymbol r_i; r_i.setValue(1.0); // the numeric value of r_i is 1.0 r_i.setDerivative(q_j, 1, 3.0); // the numeric value of d(r_i)/d(q_j) is 3.0 r_i.setDerivative(q_j, 2, 2.0); // the numeric value of d^2(r_i)/d(q_j)^2 is 2.0 // all other derivatives are standard zero then e.g.: r_i.diff(q_j, 1) returns a "MySymbol" with "setValue(3.0)" and "setDerivative(q_j, 1, 2.0)", r_i.diff(q_j, 2) returns a "numeric(2.0)", r_i.diff(q_j, i) with i > 2 returns a "numeric(0)", r_i.diff(l, 1) returns a "numeric(0)" cause we never did a "setDerivative()" with l as parameter. This should allow to write things like: MySymbol y = r_i + r_i.diff(q_j, 2) + r_i.diff(l, 1); double f = ex_to<double>( y ).to_double(); cout << f << endl; // should be "3.0" and: cout << sin(r_i).diff(q_j,1) << endl; // should be "cos(1.0)*3.0" but also: cout << ( sin(r_i).diff(q_j,1) ).diff(q_j,1) << endl; which should return the value of "cos(1.0) * 2.0 - sin(1.0) * (3.0)^2" Does this seem feasible? Anyone see any big flaws, has objections, knows a better name then "MySymbol"? ;) A remark: It might also be nice if one could do things like: symbol q_j; MySymbol r_i; r_i.setValue(1.0); r_i.setDerivative(q, 1, 3.0); MySymbol y = r_i + r_i.diff(q_j, 1); double f = ex_to<double>( y ).to_double(); r_i.setValue(5.0); double g = ex_to<double>( y ).to_double(); // without reassigning y cout << f << ", " << g << endl; // prints "4.0, 8.0" but this implies coupling between x and the object returned by e.g. x.diff(...). Don't know if this is good. Many greetings, Johan Rutgeerts

1 0

complex conjugation
by Andrius Kurtinaitis 12 Dec '03

12 Dec '03

Hello, it appears that the GiNaC expressions do not support the complex conjugation... :-( I need it, but I do not know yet how to implement it. Has anybody tried to do this? Clearly conjugation is a function ;-) Ok, i'll try to read the tutorial. Let's see if I get this working... Best regards Andrius Kurtinaitis

3 6

composite function
by Volker Vogel 12 Dec '03

12 Dec '03

Hello, I'am trying to define composite functions, for example this one: / 1 for all x > 1 f(x) =| x*x for all 1 >= x >= -1 \ 1 for all x <-1 Later I want to evaluate this function with different x. I tried to do this allready using symbol x("x"); string polystr="if (abs(x)>1) {1}else{x*x}"; ex poly; poly=ex(polystr,x); cout << poly; but I do not get it working. I assume that I defined my polystr incorrect, anybody knows how to do that ? Grettings, Volker Vogel ______________________________________________________________________________ Horoskop, Comics, VIPs, Wetter, Sport und Lotto im WEB.DE Screensaver1.2 Kostenlos downloaden: http://screensaver.web.de/?mc=021110

2 1

Error message during compilation of cln 1.1.5
by Chris Dams 05 Dec '03

05 Dec '03

Hi everybody, Today I installed gcc 3.3.2, so I had to recompile several things. One of the things that I compiled was cln 1.1.5. This, however, gave the error message below. It is about something in the tests directory, so this did not stop me from building GiNaC (the "consistency of numeric types" exam of GiNaC is passed, so I suppose everything is okay), but it still is a bit irritating. Perhaps the maintainer of CLN, who is believed to be reading this mailing list, knows what could be done about this. Have a nice day, Chris Dams make[1]: Entering directory `/home/chrisd/download/cln-1.1.5/tests' gcc -g -O2 -I../include -I../include -I. -I../src -I../src/base -I../src/float -I../src/float/lfloat -I../src/integer -I../src/base/random -I../src/base/digitseq -I../src/base/digit -I../src/base -I../src/base -c exam.s -o exam.o exam.s: Assembler messages: exam.s:246: Error: file number 70 already allocated make[1]: *** [exam.o] Error 1 make[1]: Leaving directory `/home/chrisd/download/cln-1.1.5/tests' make: *** [install] Error 2

1 0

Intervals in GiNaC
by Simon Langley 01 Dec '03

01 Dec '03

Carl, I have an interval class which I use in GiNaC. It isn't as tidy as I'd like but it works in the following sense: * I'm only interested in polynomials so the interval class doesn't have other functions like sin/exp/etc just field operations and powers. * It assumes that either rational or s/w floating point arithmetic is used. There is a problem with h/w floating point in GiNaC because underflow/overflow errors may not get noted. That said because of the way GiNaC works you can mix floating point and rational arithmetic. * When used in Ginac expressions, interval expressions aren't simplified. I.e. you can write ex z = interval(-1,1) + x + interval(1, 3)*4 + interval(1.3, 3.1) but z.simplify() won't add the intervals and z.subs(x == 2.3) won't reduce the expression to an interval. I've written a procedure to evaluate a Ginac expression in intervals to an interval (again assuming you only need polynomial expressions). If you want it let me know. By the way are you working with Dan Richardson or has interest in intervals at Bath spread? Simon ~~~~~ Simon Langley c/o CEMS, UWE, BS16 1QY

1 0

automatic differentiation and change expression ordering
by Andrius Kurtinaitis 20 Nov '03

20 Nov '03

Hello, I have two questions: 1. is it possible to prevent automatic differentiation of some expression? For example, when I do (with ginac 1.1.4): DECLARE_FUNCTION_1P( intensity ) static ex intensity_eval( const ex & r ) { return intensity(r).hold(); } // here i tried to hold it :-) REGISTER_FUNCTION( intensity, eval_func(intensity_eval) ) int main() { symbol r("r"); ex der1 = mul( diff( intensity(r), r ), r ).hold(); ex deriv = der1.diff(r); cout << deriv << endl; return 0; } I get the following output: D[0,0](intensity)(r)*r+D[0](intensity)(r) But I would like to get something like this: D[0]( D[0](intensity)(r)*r )(r) Is it possible with ginac? 2. What is the simplest way to change the order in multiplication so that the differential operators print after all other terms: r*D[0,0](intensity)(r) instead of D[0,0](intensity)(r)*r And the third one ;-) Is there maybe some work done to output fderivative in latex style? Thanks in advance Andrius

2 1