Dear All, In my GiNaC code I oftenly use some objects which have special properties with respect to power. For example, if values of a symbol s are assumed to be only -1, 0, or 1, then s^3=s. Another example is provided by "imaginary units" for double and dual numbers. It will be good for such objects to use such simplifying identities on the fly in power_eval() method rather then by user calls to subs() afterwards. I propose a patch to GiNaC which allows to supply such simplification rules from a user-defined classes. It requires only two new virtual methods to basic class and should not involve much overheads to the general GiNaC performance. I give an example below (which should work with the patch) for a symbol which is assumed to take only values -1 or 1, thus satisfying to s^2=1 rule. If general feeling toward this proposition will be positive I would send a more complete patch including documentation. Best wishes, Vladimir PS Somehow I was not able to use macros DEFAULT_UNARCHIVE() and DEFAULT_COMPARE() for my class, since compiler produces the error: error: expected constructor, destructor, or type conversion before 'int' What may be the reason for this? -- Vladimir V. Kisil email: kisilv@maths.leeds.ac.uk -- www: http://maths.leeds.ac.uk/~kisilv/ #include <iostream> #include <stdexcept> #include <ginac/ginac.h> using namespace std; using namespace GiNaC; class myrealsymb : public symbol { GINAC_DECLARE_REGISTERED_CLASS(myrealsymb, symbol) public: myrealsymb(const std::string & initname); virtual ex eval_power_basis(const ex & exp) const; virtual bool info(unsigned inf) const; }; GINAC_IMPLEMENT_REGISTERED_CLASS(myrealsymb,symbol) myrealsymb::myrealsymb(): symbol(){}; myrealsymb::myrealsymb(const std::string & initname) : symbol(initname) {;} ex myrealsymb::eval_power_basis(const ex & exp) const { if (is_a<numeric>(exp)) { if (ex_to<numeric>(exp).is_even()) return numeric(1); if (ex_to<numeric>(exp).is_odd()) return *this; } return power(*this,exp).hold(); } bool myrealsymb::info(unsigned inf) const { return (inf == info_flags::is_power_basis) || symbol::info(inf); } myrealsymb::myrealsymb(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) { } void myrealsymb::archive(archive_node &n) const { symbol::archive(n); } ex myrealsymb::unarchive(const archive_node &n, lst &sym_lst) { return (new myrealsymb(n, sym_lst))->symbol::setflag(status_flags::dynallocated); } int myrealsymb::compare_same_type(const basic & other) const { return symbol::compare_same_type(other); } int main(){ try { ex e=myrealsymb("s"); cout << pow(e,2) << endl; // -> 1 cout << pow(1+e,2).expand() << endl; // -> 2+2*s } catch (exception &p) { cerr << "Got problem: " << p.what() << endl; } return 0; }
Hi Vladimir, You wrote:
In my GiNaC code I oftenly use some objects which have special properties with respect to power. For example, if values of a symbol s are assumed to be only -1, 0, or 1, then s^3=s. Another example is provided by "imaginary units" for double and dual numbers.
Thanks a lot for your input, but this time I am not convinced it would be good to include your contribution. About five years ago, I was doing something very similar. I had polynomials of theta functions of different (symbolic) arguments and it was convenient enough to write tree traversers and convert the theta functions to individual symbols and back again to theta functions. Just replacing s^2 by s was not hard. But it also was hardly enough! This was, because theta functions offer a number of other, very useful identities: 1) If the argument of theta factorizes, one can transform theta(P*Q) -> theta(P)*theta(Q)+theta(-P)*theta(-Q). 2) One can find arguments that differ in sign only and transform theta(-P) -> 1-theta(P). 3) One can easily proof that square free factorization of polynomials of idempotent things is equal to the full factorization and do fine things with GiNaC's sqrfree function. n) ... My point is that with objects as specialized as theta functions or, equivalently, step functions, many specialized rules suddenly hold and it does not appear to make much sense adding virtual functions to the basic class in order to support just a few of them. I suggest doing these things in specialized code. Best wishes -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Dear Richard and Alexei, You are very preoccupied with slowdown of GiNac with new patch. I am including (under the signature) two outputs from "make check" with standard GiNaC and a patched one: *no obvious difference* in timing.
"RK" == Richard B Kreckel <kreckel@ginac.de> writes:
RK> My point is that with objects as specialized as theta functions RK> or, equivalently, step functions, many specialized rules RK> suddenly hold and it does not appear to make much sense adding RK> virtual functions to the basic class in order to support just a RK> few of them. It looks like you missed the point: no new rules are hard-coded into core GiNaC at all. The two virtual functions allow a user to include seamless automatic evaluation for *his own specialised* object in user-defined classes. It is up to the end user to decide which rule if any he wish to be applied to which object. The price for this freedom in term of speed is virtually nil. Best wishes, Vladimir -- Vladimir V. Kisil email: kisilv@maths.leeds.ac.uk -- www: http://maths.leeds.ac.uk/~kisilv/ GiNaC will now run through some basic timings: timing commutative expansion and substitution.... passed size: 100 200 400 800 time/s: 0.18 0.72 3.3 15. timing Laurent series expansion of Gamma function.... passed order: 20 25 30 35 time/s: 1.6 1.9 6.7 21. timing determinant of univariate symbolic Vandermonde matrices.... passed dim: 8x8 10x10 12x12 14x14 time/s: 0.16 0.40 3.5 30. timing determinant of polyvariate symbolic Toeplitz matrices.... passed dim: 7x7 8x8 9x9 10x10 time/s: 0.17 0.69 3.0 12. timing hash map operations.... passed size: 10000 50000 100000 500000 insert/s: 0.016 0.10 0.23 1.1 find/s: 0.012 0.052 0.11 0.62 erase/s: 0.0080 0.044 0.10 0.52 timing Lewis-Wester test A (divide factorials). passed 0.58s timing Lewis-Wester test B (sum of rational numbers). passed 0.0034s timing Lewis-Wester test C (gcd of big integers). passed 0.044s timing Lewis-Wester test D (normalized sum of rational fcns). passed 0.13s timing Lewis-Wester test E (normalized sum of rational fcns). passed 0.068s timing Lewis-Wester test F (gcd of 2-var polys). passed 0.017s timing Lewis-Wester test G (gcd of 3-var polys). passed 0.51s timing Lewis-Wester test H (det of 80x80 Hilbert). passed 1.5s timing Lewis-Wester test I (invert rank 40 Hilbert). passed 0.43s timing Lewis-Wester test J (check rank 40 Hilbert). passed 0.26s timing Lewis-Wester test K (invert rank 70 Hilbert). passed 2.9s timing Lewis-Wester test L (check rank 70 Hilbert). passed 1.5s timing Lewis-Wester test M1 (26x26 sparse, det). passed 0.058s timing Lewis-Wester test M2 (101x101 sparse, det) disabled timing Lewis-Wester test N (poly at rational fcns) disabled timing Lewis-Wester test O1 (three 15x15 dets)... passed 11.s (average) timing Lewis-Wester test O2 (Resultant) disabled timing Lewis-Wester test P (det of sparse rank 101). passed 0.38s timing Lewis-Wester test P' (det of less sparse rank 101). passed 0.98s timing Lewis-Wester test Q (charpoly(P)). passed 24.s timing Lewis-Wester test Q' (charpoly(P')). passed 45.s timing computation of antipodes in Yukawa theory...... passed 15.s (total) timing Fateman's polynomial expand benchmark. passed 83.s PASS: run_times ============================================================================ GiNaC will now run through some basic timings: timing commutative expansion and substitution.... passed size: 100 200 400 800 time/s: 0.18 0.73 3.3 15. timing Laurent series expansion of Gamma function.... passed order: 20 25 30 35 time/s: 1.6 1.9 6.8 22. timing determinant of univariate symbolic Vandermonde matrices.... passed dim: 8x8 10x10 12x12 14x14 time/s: 0.16 0.39 3.5 29. timing determinant of polyvariate symbolic Toeplitz matrices.... passed dim: 7x7 8x8 9x9 10x10 time/s: 0.16 0.71 3.0 12. timing hash map operations.... passed size: 10000 50000 100000 500000 insert/s: 0.020 0.10 0.22 1.1 find/s: 0.012 0.052 0.11 0.60 erase/s: 0.0080 0.040 0.10 0.50 timing Lewis-Wester test A (divide factorials). passed 0.58s timing Lewis-Wester test B (sum of rational numbers). passed 0.0036s timing Lewis-Wester test C (gcd of big integers). passed 0.041s timing Lewis-Wester test D (normalized sum of rational fcns). passed 0.14s timing Lewis-Wester test E (normalized sum of rational fcns). passed 0.076s timing Lewis-Wester test F (gcd of 2-var polys). passed 0.016s timing Lewis-Wester test G (gcd of 3-var polys). passed 0.54s timing Lewis-Wester test H (det of 80x80 Hilbert). passed 1.5s timing Lewis-Wester test I (invert rank 40 Hilbert). passed 0.44s timing Lewis-Wester test J (check rank 40 Hilbert). passed 0.25s timing Lewis-Wester test K (invert rank 70 Hilbert). passed 2.9s timing Lewis-Wester test L (check rank 70 Hilbert). passed 1.5s timing Lewis-Wester test M1 (26x26 sparse, det). passed 0.064s timing Lewis-Wester test M2 (101x101 sparse, det) disabled timing Lewis-Wester test N (poly at rational fcns) disabled timing Lewis-Wester test O1 (three 15x15 dets)... passed 11.s (average) timing Lewis-Wester test O2 (Resultant) disabled timing Lewis-Wester test P (det of sparse rank 101). passed 0.38s timing Lewis-Wester test P' (det of less sparse rank 101). passed 1.00s timing Lewis-Wester test Q (charpoly(P)). passed 23.s timing Lewis-Wester test Q' (charpoly(P')). passed 44.s timing computation of antipodes in Yukawa theory...... passed 14.s (total) timing Fateman's polynomial expand benchmark. passed 86.s PASS: run_times
Dear Vladimir, Vladimir Kisil wrote:
Dear Richard and Alexei,
You are very preoccupied with slowdown of GiNac with new patch. I am including (under the signature) two outputs from "make check" with standard GiNaC and a patched one: *no obvious difference* in timing.
No, I would have been surprised if your patch made a noticable difference. My concerns are rather that rooting something as specialized inside the class hierarchy will not be useful enough but open the floodgates for much more specialized ideas and make the code even harder to understand for novices. This is comparable to the saying "if your function has 13 parameters, it probably has one too... few." Have you tried using subs or custom tree traversers for your purpose? (Along the ideas sketched in the FAQ, I mean.) -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>
Dear Richy,
My concerns are rather that rooting something as specialized inside the class hierarchy will not be useful enough but open the floodgates for much more specialized ideas and make the code even harder to understand for novices.
Once more: nothing specialised is intended for the GiNaC core classes. This is only a pair of wires leading outside of the kernel where an end user can connect his own devices to.
Have you tried using subs or custom tree traversers for your purpose? (Along the ideas sketched in the FAQ, I mean.)
Imagine that any calculation with complex numbers will require manual substitution of I^2=-1, will it be convenient? The purpose of the extension is just to let a user a freedom to make some such simple rules work automatically. It is the responsibility of the end user to decide himself what kind of rules are suitable for his task. I think that this is a nice way to keep GiNaC kernel small and tidy and yet suitable for very wide set of questions. Best wishes, Vladimir -- Vladimir V. Kisil email: kisilv@maths.leeds.ac.uk http://maths.leeds.ac.uk/~kisilv
Hello, On Sat, Oct 28, 2006 at 08:57:27PM +0100, kisilv@maths.leeds.ac.uk wrote:
My concerns are rather that rooting something as specialized inside the class hierarchy will not be useful enough but open the floodgates for much more specialized ideas and make the code even harder to understand for novices.
Once more: nothing specialised is intended for the GiNaC core classes. This is only a pair of wires leading outside of the kernel where an end user can connect his own devices to.
What if someone else wants to insert another "pair of wires" to implement another weird feature? And than one more... making the whole thing more cumbersome and slow (yes, two extra useless method calls do not really improve the performance). Please stop adding more methods to core classes (unless there is a _really_ good reason for this), there are already much of them (evalm, evalf, eval_integ, eval_indexed, eval_ncmul). Pretty please.
Have you tried using subs or custom tree traversers for your purpose? (Along the ideas sketched in the FAQ, I mean.)
Imagine that any calculation with complex numbers will require manual substitution of I^2=-1, will it be convenient?
It depends. E.g., for someone who hardly ever uses complex numbers this is not a big deal. On the other hand,
+It is not unusual to use some objects which have special properties with +respect to exponentiation. For example, if a symbol @code{s} may take +only values @math{-1} or @math{1} then @math{s^2=1}.
This sounds like modular arithmetics (in Z_3, to be more specific). It would be nice if GiNaC had support of such a thing. CLN has modular integers, so _probably_ it is possible to do such calculations in more natural way (so that simplification of this kind are handled by CLN, just as I^2 -> -1).
The purpose of the extension is just to let a user a freedom to make some such simple rules work automatically.
IMHO, GiNaC's purpose is somewhat inverse: not to give as much freedom as possible, but to restrict the freedom to do some specific things effectively. BTW, other tools used by [particle] physicists (e.g. FORM) are even more restrictive, and there are very few complaints (if any) about that.
It is the responsibility of the end user to decide himself what kind of rules are suitable for his task.
The main point is not validity of those rules. Your implementation complicates things even for those who are not going to use any of such rules.
I think that this is a nice way to keep GiNaC kernel small and tidy and yet suitable for very wide set of questions.
IMHO, this is exaclty the way to make it more cumbersome without any positive side effects. Best regards, Alexei -- All science is either physics or stamp collecting.
Dear Alexei, It is a pity that this discussion is going round by round over the same path and we all reiterate the previous arguments.
"ASh" == Sheplyakov Alexei <varg@theor.jinr.ru> writes:
ASh> What if someone else wants to insert another "pair of wires" to ASh> implement another weird feature? Once more: the purpose of this wires *to avoid* implementation of weird features in the GiNaC kernel. For any such a request to add anything to the core GiNaC with special property with respect to power::eval() will follow a polite answer: "You can do it yourself in your own derived class." ASh> It depends. E.g., for someone who hardly ever uses complex ASh> numbers this is not a big deal. On the other hand, Anyone who hardly use any differential geometry will vote for exclusion of all "indexed" part from GiNaC, so.... ASh> This sounds like modular arithmetics (in Z_3, to be more ASh> specific). It would be nice if GiNaC had support of such a ASh> thing. CLN has modular integers, so _probably_ it is possible ASh> to do such calculations in more natural way (so that ASh> simplification of this kind are handled by CLN, just as I^2 -> ASh> -1). There are also symbols taking only values -1, 0, 1, with the property s^3=s, there are imaginary units for double and dual number (e^2=1, \epsilon^2=0), etc. Adding all of them separately "in a more natural way" either to GiNaC or CLN will really do it very cumbersome. ASh> IMHO, GiNaC's purpose is somewhat inverse: not to give as much ASh> freedom as possible, but to restrict the freedom to do some ASh> specific things effectively. BTW, other tools used by ASh> [particle] physicists (e.g. FORM) are even more restrictive, ASh> and there are very few complaints (if any) about that. Probably my real problem is that I am not a "[particle] physicist" ;-) I was also mislead by the phrase from the GiNaC tutorial: "we...tried to fill the gap by writing GiNaC. But of course its applications are in no way restricted to theoretical physics." ASh> The main point is not validity of those rules. Your ASh> implementation complicates things even for those who are not ASh> going to use any of such rules. Untrue. Anyone who do not need them may safely assume that they do not exist at all. Virtual methods basic::eval_power_base() and basic::eval_power_exponent() are not ever called from the core GiNaC. And they will not be ever called by any user code unless the user will take specific steps to call them. ASh> Please stop adding more methods to core classes Alright, assuming no second opinion exists. Best wishes, Vladimir -- Vladimir V. Kisil email: kisilv@maths.leeds.ac.uk -- www: http://maths.leeds.ac.uk/~kisilv/
Hello, On Mon, Oct 30, 2006 at 06:05:52PM +0000, Vladimir Kisil wrote:
It is a pity that this discussion is going round by round over the same path and we all reiterate the previous arguments.
So, let me first give a short summary: 1) I think adding more virtual methods to the `basic' class should be considered evil and used only as a last resort (if ever). 2) Calling newly added methods from the eval() methods of core classes is even more evil. 3) I think the proper way to implement the feature you need is to provide support of polynomials (rational functions, algebraic expressions, etc.) over modular integers and and [finite] algebraic extensions. 4) Until this is done, it is better to stick with subs(), map_function and visitors. As a "lower acceptable limit" it is possible to add new virtual method (sigh!) to the `basic' class (something like eval_pow) which *never ever* gets called by power::eval (or any other core classes) and does whatever you want. Note that such a method must be called explicitly (like evalm, evalf, eval_integ, simplify_indexed and so on). However, this is still dirty hack. It is just not that harmful. Now I'll try to elaborate, so (hopefully) my point will be more clear.
ASh> What if someone else wants to insert another "pair of wires" to ASh> implement another weird feature?
Once more: the purpose of this wires *to avoid* implementation of weird features in the GiNaC kernel.
It does not really matter if the actual code implementing $FEATURE resides in GiNaC or there exist some hook in the core code. In general, explicitly calling some method of the *derived class* from the parent one sounds very anti-OO. Unfortunately, there are already a number of such hooks in GiNaC (evalm, eval_integ, simplify_indexed, to name few). I have no idea how to fix this, but at least we should stop adding new ones. To re-iterate: clever overloading of existing virtual method[s] in a derived class is a good way to implement new features. Adding more eval_foo virtual methods to the `basic' class is not.
For any such a request to add anything to the core GiNaC with special property with respect to power::eval() will follow a polite answer: "You can do it yourself in your own derived class."
[Just kidding] So my polite answer is: "You can do it yourself with subs() without modifying power::eval at all". Seriously, at some moment you will _certainly_ want to hack mul::eval and/or add:eval (for more complete emulation of modular arithmetics) and ask for more virtual methods in basic.
ASh> This sounds like modular arithmetics (in Z_3, to be more ASh> specific). It would be nice if GiNaC had support of such a ASh> thing. CLN has modular integers, so _probably_ it is possible ASh> to do such calculations in more natural way (so that ASh> simplification of this kind are handled by CLN, just as I^2 -> ASh> -1).
There are also symbols taking only values -1, 0, 1,
Yes, polynomials and rational functions over such rings (or fields, as in your example) need to be handled too.
with the property s^3=s,
Modular arithmetics has in fact more generic property: for any x \in Z_p (with p being a prime number) x^p = x. It would be nice if power::eval handled such things (*without* evil hacks). For proper support of modular polynomials/rational functions it is also necessary to modify (at least) mul::eval and add::eval. Once again, it would be really nice to have such a feature[s] (e.g., to have reasonably fast _ordinary_ polynomial operations).
there are imaginary units for double and dual number (e^2=1, \epsilon^2=0), etc. Adding all of them separately "in a more natural way" either to GiNaC or CLN will really do it very cumbersome.
CLN _already_ supports modular integers. However, [finite] algebraic extensions (except ordinary complex numbers) are not supported. Polynomials over arbitrary rings are supported by CLN too, but AFAIK it implements primitive operations only (no GCD, no factorization, etc.). Thus, there are indeed a lot of non-trivial work to be done. Let me repeat once again: it is complicated, but _not_ cumbersome. BTW, I think at some moment you will implement most (if not all) of these in order to work with your "symbols which take only 3 values". So it is better to do this properly instead of adding random hacks.
ASh> It depends. E.g., for someone who hardly ever uses complex ASh> numbers this is not a big deal. On the other hand,
Anyone who hardly use any differential geometry will vote for exclusion of all "indexed" part from GiNaC, so....
They are free to patch away all indexed classes and related methods. However, doing so does not buy much: AFAIK, if the expression does not contain indexed objects, neither indexed::eval, nor eval_indexed & friends get called implicitly, so the only price is extra memory for vtable entries and disk space for library code (some operating systems won't even load to the RAM "unnecessary" portion of code). On the other hand, your patch adds two unnecessary calls for _every_ power object (see more detailed explanation below). I think this is absolutely unacceptable.
Probably my real problem is that I am not a "[particle] physicist" ;-)
I don't think this is problem :)
I was also mislead by the phrase from the GiNaC tutorial:
"we...tried to fill the gap by writing GiNaC. But of course its applications are in no way restricted to theoretical physics."
Tutorial also warns: "First of all, GiNaC's name must be read literally." And introductory paper (http://www.ginac.de/csSC-0004015.pdf) clearly states that minimalistic design of GiNaC (as opposed to "can do everything equally bad") _is a feature_. However, this does not prevent people from using GiNaC for other tasks. For example, Ecco (http://sourceforge.net/projects/echempp) has very little to do with high-energy physics.
ASh> The main point is not validity of those rules. Your ASh> implementation complicates things even for those who are not ASh> going to use any of such rules.
Untrue. Anyone who do not need them may safely assume that they do not exist at all.
Only modulo two unnecessary virtual method calls in power::eval. To be more precise, I mean this hunk:
+ // user defined classes may provide their own virtual methods for power evaluation + if (ebasis.info(info_flags::is_power_basis)) + return ebasis.eval_power_basis(eexponent); + + if (eexponent.info(info_flags::is_power_exponent)) + return eexponent.eval_power_exponent(ebasis);
E.g. evaluation of a simple x^2 now takes several cycles longer (to call symbol::info and numeric::info methods and check that no custom evaluation rules exist). Evaluation of monomial like x^2*y^3*z^4 now takes a little bit longer. Evaluation of polynomial like x^99 + y x^98 + y^2 x^97 + ... + y^99 now takes... oops, somewhat longer. Given the fact that large and ugly intermediate expressions are quite common in [multivariate] GCD calculation, this slowdown is not appreciated at all. "Things I do not use should not harm me".
Virtual methods basic::eval_power_base() and basic::eval_power_exponent() are not ever called from the core GiNaC.
First of all, it does not really matter if you call info() or directly eval_power_{base,exponent} -- in any case there are two extra calls [in the fast code path]. Secondly, checking if there is non-trivial implementation of virtual method is a bad thing on its own. It is best to provide sensible (e.g., no-op) default implementation and overload the method in the derived class. This way there is no need in bogus extra entries in info_flags (IMHO these should provide some mathematical details, not implementation ones), and users need not to implement custom info method.
And they will not be ever called by any user code unless the user will take specific steps to call them.
Unfortunately, this is not completely correct (see my explanation above). Best regards, Alexei -- All science is either physics or stamp collecting.
participants (4)
-
kisilv＠maths.leeds.ac.uk
-
Richard B. Kreckel
-
varg＠theor.jinr.ru
-
Vladimir Kisil