February 2002 - PURRS-devel

Re: [Fwd: fi_lib++] (2nd attempt: first time sent on 06/02/2002)
by Roberto Bagnara 15 Feb '02

15 Feb '02

J.Wolff v. Gudenberg wrote: > Dear Roberto Bagnara, > > For what application do you need complex interval arithmetic ? > Do you need the extended mode ? Dear Wolff, our field of application is data-flow analysis and abstract interpretation. More specifically we now need complex interval arithmetic in a project that deals with the automatic approximate solution of difference equations. We need the extended mode; we would also need intervals with open boundaries, but we can live without that for the moment. > In the next future we are not planning to implement complex arithmetic > for filib++. > but we would appreciate, if your group can add this facility for the > library. OK, I have just sent a message to the authors of CoStLy to see if anyone is doing the same thing. We will start working on that immediately if the answer is negative. > We are investigating a way to use the same cvs archive behind our > firewall. We have already put filib++ on our CVS server so that we can easily track our and your changes. >>From the math viewpoint of complex arithmetic, we think that it is not > trivial to define > the extended mode. We currently see no ressources on our site for > research in this area. > But we are positively interested in the subject. OK. I have just fetched the new version of filib++ (the one of February 6) and I am having some problems with it. Firstly, the LIBDIR directory must be created by hand. In fact, if one tries to do `make libs' things end up with ar rv /home/roberto/filib++/interval/macro/lib/fi_lib.a Interval.o acos.o acosh.o acot.o acoth.o asin.o asinh.o atan.o atanh.o cos.o cosh.o cot.o coth.o exp.o exp10.o exp2.o expm1.o log.o log10.o log1p.o log2.o pow.o q_acos.o q_acot.o q_acsh.o q_acth.o q_asin.o q_asnh.o q_atan.o q_atn1.o q_atnh.o q_cos.o q_cos1.o q_cosh.o q_cot.o q_coth.o q_cth1.o q_ep1.o q_epm1.o q_errm.o q_ex10.o q_exp.o q_exp2.o q_expm.o q_glbl.o q_lg10.o q_log.o q_log1.o q_log2.o q_rtrg.o q_sin.o q_sin1.o q_sinh.o q_tan.o q_tanh.o sin.o sinh.o sqr.o sqrt.o tan.o tanh.o ar: /home/roberto/filib++/interval/macro/lib/fi_lib.a: No such file or directorymake[2]: *** [/home/roberto/filib++/interval/macro/lib/fi_lib.a] Error 1 make[2]: Leaving directory `/home/roberto/filib++/interval/macro/src' make[1]: *** [lib] Error 2 make[1]: Leaving directory `/home/roberto/filib++/interval/macro' Supply variable GMAKE please! Making primitive.so from primitive.cpp in ieee: done. Building shared prim support library libprim.so.0: done. Setting symbolic link from lib/libprim.so.0 to lib/libprim.so: done. So one has to create the directory .../interval/macro/lib by hand to complete the compilation of the library. May I also suggest to invoke `ar' with the options `cruv' instead of just `rv'? Note also the confusing error message "Supply variable GMAKE please!". Another problem is that the distribution seems somewhat mutilated: in fact, the `test' directory seems missing and this is responsible for `make check' to fail and for useful things like `make clean' to print not very reassuring error messages. That is all for the moment. Please let me know when you update the distribution (BTW: it would be very useful to identify archives for different versions like filib++-2.0.1.tar.gz, filib++-2.0.2.tar.gz and so forth.) All the best Roberto -- Prof. Roberto Bagnara Computer Science Group Department of Mathematics, University of Parma, Italy http://www.cs.unipr.it/~bagnara/ mailto:bagnara@cs.unipr.it

1 0

purrs ./acinclude.m4 ./config.h.in ./configure ...
by Roberto Bagnara 15 Feb '02

15 Feb '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: roberto(a)cs.unipr.it 2002-02-15 11:20:05 Modified files: . : acinclude.m4 config.h.in configure.ac tests : Makefile.am Added files: m4 : ac_check_use_readline.m4 Log message: Prepare for using GNU readline. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/acinclude.m4.diff?cvsroot=p… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/config.h.in.diff?cvsroot=pu… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/configure.ac.diff?cvsroot=p… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/m4/ac_check_use_readline.m4… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/tests/Makefile.am.diff?cvsr…

1 0

Re: [Fwd: fi_lib++]
by Roberto Bagnara 06 Feb '02

06 Feb '02

German Tischler wrote: > Dear Mr. Bagnara, > > thank you for your interest in filib++ and your patches, I have > merged them into our source repository. Please find our answers > in the quoted text below. > > >>Date: Fri, 25 Jan 2002 19:32:30 +0100 >>From: Roberto Bagnara <bagnara(a)cs.unipr.it> >>Organization: Department of Mathematics, University of Parma, Italy >>X-Mailer: Mozilla 4.78 [en] (X11; U; Linux 2.4.17 i686) >>X-Accept-Language: en >>To: wolff(a)informatik.uni-wuerzburg.de, >> werner.hofschuster(a)math.uni-karlsruhe.de, >> walter.kraemer(a)math.uni-karlsruhe.de >>CC: purrs-devel(a)cs.unipr.it >>Subject: fi_lib++ >> >>Dear Sirs, >> >>first of all I would like to thank you for fi_lib++. >>I have just finished reading the documentation and >>it turns out that fi_lib++ has almost all the features >>of the library I was writing myself. I have a few questions: >> >>1) What is the official web site for the library? >> > > If there is news concerning the library, we will put it up at > > http://www2.informatik.uni-wuerzburg.de/mitarbeiter/wvg/Public > > >>2) Is there a mailing list for fi_lib++ users and developers? >> > > So far the public interest on the library has not been high. > Therefore we didn't yet see a cause for having a mailing list. > > >>3) If not, to whom should I report bugs and improvements? >> > > You should be able to get the quickest response on bug reports > and improvements by sending an email to me (German Tischler > <tanis(a)informatik.uni-wuerzburg.de>), but please also send > a carbon copy to Prof. Wolff von Gudenberg > (<wolff(a)informatik.uni-wuerzburg.de>). > > If you are interested in a current development version, just tell me. Dear German, thanks for your message. I am certainly interested in the current development version. What I would like to avoid is my group and myself having to maintain another version of filib++. It is for this reason that it is important to have a web site where updated versions can be retrieved quickly and a mailing list so that the users can be kept up to date. Our first need is an extension of the library to deal with complex numbers. A good starting point would seem to be CoStLy (http://www.uni-karlsruhe.de/~Markus.Neher/CoStLy.html). Do you plan to develop such an extension in the immediate future? All the best Roberto -- Prof. Roberto Bagnara Computer Science Group Department of Mathematics, University of Parma, Italy http://www.cs.unipr.it/~bagnara/ mailto:bagnara@cs.unipr.it

1 0

Re: Roots of unity
by Roberto Bagnara 05 Feb '02

05 Feb '02

Bob McElrath wrote: > > ginsh: > > evalf((-1)^(1/3)); > 0.5+0.86602540378443864673*I > > While this answer is technically correct: > exp(I*pi/3) = -1 > > It's not what I was expecting, and the cause of some headache today. I was > trying to write a little routine like RootOf to solve a cubic polynomial, and > was getting three complex answers. (at least one of them must be real if you > start with a real polynomial) > > How should I go about enforcing that (-1)^(1/3) = -1? I think it would be unreasonable to require that (-1)^(1/3) = -1, since the other solutions are as good as -1. Thus I believe GiNaC does the right thing by leaving the symbolic expression (-1)^(1/3) untouched. On the other hand, a little experiment with ginsh shows that > 1^(1/2); 1 so that here GiNaC does choose one root. I believed it does so on the grounds that even roots of non-negative real numbers are usually defined that way. However, still with ginsh, we get also > (-1)^(1/2); I so it seems this tradition is pushed a little too far. I think these issues should be clarified with the GiNaC team. Perhaps you are criticizing the behavior of evalf() when applied to (-1)^(1/3) and I agree with you. Since odd roots of real numbers always have a real solution, why not defining evalf() so that it does "the Right Thing"? The point is perhaps: how much would it cost? GiNaC people? Please forgive me if I failed to understand your point. > I'm interested in > polynomials over the reals, so the complex solutions don't interest me. This > behaviour of evalf is preventing me from extracting the real solution. > Here's my little function if anyone is interested: > > // Find the roots of a polynomial. > // The polynomial can be at most 4th order. > lst RootOf(const ex& poly, const ex& x) { > ex p = poly; > ex D,Q,R,S,T; > p.expand(); > ex a0 = poly.coeff(x,0); > ex a1 = poly.coeff(x,1); > ex a2 = poly.coeff(x,2); > ex a3 = poly.coeff(x,3); > ex a4 = poly.coeff(x,4); > switch(p.degree(x)) { > case 0: > return lst(); > case 1: > return lst(-a0/a1); > case 2: > return lst((-a1 + sqrt(a1*a1-4*a2*a0))/(2*a2), (-a1 - sqrt(a1*a1-4*a2*a0))/(2*a2)); > case 3: > a0 /= a3; > a1 /= a3; > a2 /= a3; > Q = (3*a1-a2*a2)/9; > R = (9*a2*a1-27*a0-2*a2*a2*a2)/54; > D = pow(Q,3) + pow(R,2); > S = pow(R+sqrt(D),numeric(1,3)); > T = pow(R-sqrt(D),numeric(1,3)); > // These formulas come from: http://mathworld.wolfram.com/CubicEquation.html > return lst( > -numeric(1,3)*a2+S+T, > -numeric(1,3)*a2-numeric(1,2)*(S+T)+numeric(1,2)*I*sqrt(ex(3))*(S-T), > -numeric(1,3)*a2-numeric(1,2)*(S+T)-numeric(1,2)*I*sqrt(ex(3))*(S-T) > ); > //case 4: > // http://mathworld.wolfram.com/QuarticEquation.html > // It seems this is defined ambiguously. eq. 30 wants "a real root" but there > // may be as many as 3! Can there 12 roots!??! > default: > throw(std::domain_error("RootOf: Cannot take the root of a polynomial with degree > 4")); > } > return lst(); > } If that is your problem, I believe the piece of code below will solve it. Beware: it is extracted from a library we are writing so that it only received very little testing. I hope it helps. All the best Roberto -- Prof. Roberto Bagnara Computer Science Group Department of Mathematics, University of Parma, Italy http://www.cs.unipr.it/~bagnara/ mailto:bagnara@cs.unipr.it /* Definitions for the algebraic equation solver. Copyright (C) 2002 Roberto Bagnara <bagnara(a)cs.unipr.it> This file is part of the Parma University's Recurrence Relation Solver (PURRS). The PURRS is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. The PURRS is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. For the most up-to-date information see the PURRS site: http://www.cs.unipr.it/purrs/ . */ #include <ginac/ginac.h> #include <cln/complex.h> typedef GiNaC::ex GExpr; typedef GiNaC::symbol GSymbol; typedef GiNaC::lst GList; typedef GiNaC::numeric GNumber; using namespace GiNaC; static GExpr cubic_root(const GExpr& e) { static GExpr one_third = GExpr(1)/3; return pow(e, one_third); } /*! Solve the equation \f$x^3 + a_1 x^2 + a_2 x + a_3 = 0\f$ and return <CODE>true</CODE> if and only if all the solutions are real. \f$x_1\f$ is guaranteed to be a real solution. */ bool solve_equation_3(const GNumber& a1, const GNumber& a2, const GNumber& a3, GExpr& x1, GExpr& x2, GExpr& x3) { GExpr Q = (3*a2 - pow(a1, 2)) / 9; GExpr R = (9*a1*a2 - 27*a3 -2*pow(a1, 3)) / 54; GExpr d = pow(Q, 3) + pow(R, 2); GExpr a1_div_3 = a1/3; if (d < 0) { GExpr sqrt_minus_Q = sqrt(-Q); GExpr theta = acos(-R/Q*sqrt_minus_Q); x1 = -a1_div_3 + 2*sqrt_minus_Q*cos(theta/3); x2 = -a1_div_3 + 2*sqrt_minus_Q*cos((theta+2*Pi)/3); x3 = -a1_div_3 + 2*sqrt_minus_Q*cos((theta+4*Pi)/3); } else { GExpr sqrt_d = sqrt(d); GExpr S = cubic_root(R + sqrt_d); GExpr T = cubic_root(R - sqrt_d); GExpr S_plus_T = S + T; GExpr t1 = -S_plus_T/2 - a1_div_3; GExpr t2 = (S - T)*I*sqrt(GExpr(3))/2; x1 = S_plus_T - a1_div_3; x2 = t1 + t2; x3 = t1 - t2; } // The roots are all real if and only if d <= 0. return d <= 0; }

3 4

purrs/src alg_eq_solver.cc
by Roberto Bagnara 04 Feb '02

04 Feb '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: roberto(a)cs.unipr.it 2002-02-04 21:33:48 Modified files: src : alg_eq_solver.cc Log message: Support for nested polynomials added. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/src/alg_eq_solver.cc.diff?c…

1 0

purrs/demos purrs_eqs_cgi.cc
by Roberto Bagnara 04 Feb '02

04 Feb '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: roberto(a)cs.unipr.it 2002-02-04 07:06:05 Modified files: demos : purrs_eqs_cgi.cc Log message: Recall the user the polynomial whose roots were sought. Be more vague about the amount of CPU time consumed. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/demos/purrs_eqs_cgi.cc.diff…

1 0

purrs demos/Makefile.am src/Makefile.am src/al ...
by Roberto Bagnara 03 Feb '02

03 Feb '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: roberto(a)cs.unipr.it 2002-02-03 20:36:15 Modified files: demos : Makefile.am src : Makefile.am alg_eq_solver.cc tests : Makefile.am Added files: src : poly_factor.cc poly_factor.hh Log message: When all else fails, try to factorize the polynomial. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/demos/Makefile.am.diff?cvsr… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/src/poly_factor.cc.diff?cvs… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/src/poly_factor.hh.diff?cvs… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/src/Makefile.am.diff?cvsroo… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/src/alg_eq_solver.cc.diff?c… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/tests/Makefile.am.diff?cvsr…

1 0

purrs/src rr_solver.cc
by Roberto Bagnara 03 Feb '02

03 Feb '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: roberto(a)cs.unipr.it 2002-02-03 20:34:48 Modified files: src : rr_solver.cc Log message: Added an assertion and a FIXME. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/src/rr_solver.cc.diff?cvsro…

1 0

purrs ./config.h.in ./configure.ac demos/purrs ...
by Roberto Bagnara 03 Feb '02

03 Feb '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: roberto(a)cs.unipr.it 2002-02-03 07:55:21 Modified files: . : config.h.in configure.ac demos : purrs_eqs_cgi.cc Log message: Report on the amount of CPU time used, not on elapsed time. Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/config.h.in.diff?cvsroot=pu… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/configure.ac.diff?cvsroot=p… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/demos/purrs_eqs_cgi.cc.diff…

1 0

purrs ./config.h.in ./configure.ac demos/purrs ...
by Roberto Bagnara 02 Feb '02

02 Feb '02

CVSROOT: /cvs/purrs Module name: purrs Changes by: roberto(a)cs.unipr.it 2002-02-02 22:17:33 Modified files: . : config.h.in configure.ac demos : purrs_eqs_cgi.cc Log message: CGIs should be as nice as possible to the systems running them: renice the process so that its priority is minimized. Properly support systems without gettimeofday(). Patches: http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/config.h.in.diff?cvsroot=pu… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/configure.ac.diff?cvsroot=p… http://www.cs.unipr.it/cgi-bin/cvsweb.cgi/purrs/demos/purrs_eqs_cgi.cc.diff…

1 0