Markus Neher wrote:
Dear colleagues,
we have finished a beta version of CoStLy 0.3 that should take care of the problems you experienced with CoStLy 0.2.
Concerning the power function z^p, the library has been redesigned and enhanced. The details of the new implementation can be found in the preliminary draft of a documentation which is also attached.
We have extensively tested the library and eliminated all errors that we detected. Please report any bugs that occur with your computations.
Best regards,
Ingo and Markus
Dear Ingo and Markus, thanks for the beta. In order to be able to use it we had to make the corrections you will find in the attached patch. They fall into 5 categories: 1) namespace errors whereby some std:: qualifications were missing; 2) namespace errors when compiling (as we have to) with FILIB_NAMESPACES defined; 3) incorrect uses of abs() (an int -> int function) instead of fabs() (double -> double); 4) undefined type `Complex'; 5) unused variables (causing annoying warnings). We are now using the revised beta version. However, we could not disable our wrapper around the pow() function, since CoStLy still throws an exception for invocations like, e.g., pow(([-0.5, -0.5],[0, 0]), ([-1, -1],[0, 0])). Do you plan to address also this case? Thanks a lot 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 diff -rcp2N CoStLy/cimath.cpp newCoStLy/cimath.cpp *** CoStLy/cimath.cpp Fri Dec 6 09:50:16 2002 --- newCoStLy/cimath.cpp Sat Dec 7 19:45:57 2002 *************** typedef Interval interval; *** 53,56 **** --- 53,60 ---- typedef double flnumber; + #ifdef FILIB_NAMESPACES + using namespace filib; + #endif + #else //C-XSC Version *************** cinterval coth(const cinterval& z) *** 1104,1108 **** interval arg(const cinterval& z) { - interval wert(0.0); interval hxl(0.0),hyl(0.0); interval hxu(0.0),hyu(0.0); --- 1108,1111 ---- *************** ln *** 1318,1323 **** else { ! re_max = max( srez, -irez ); ! im_max = max( simz, -iimz ); re_ln = interval( inf( log( abs(interval(re_min),interval(im_min)) ) ), --- 1321,1326 ---- else { ! re_max = std::max( srez, -irez ); ! im_max = std::max( simz, -iimz ); re_ln = interval( inf( log( abs(interval(re_min),interval(im_min)) ) ), *************** cinterval sqrt(const cinterval& z) *** 1465,1469 **** //right half plane (no intersection) rwert = sqrt(hxl); //Minimum ! rwert = interval(inf(rwert),max(sup(re_sqrt(hxu,hyu)),sup(re_sqrt(hxu,hyl)))); iwert = im_sqrt(hxl,hyl); //Minimum iwert = interval(inf(iwert),sup(im_sqrt(hxl,hyu))); //with maximum --- 1468,1472 ---- //right half plane (no intersection) rwert = sqrt(hxl); //Minimum ! rwert = interval(inf(rwert),std::max(sup(re_sqrt(hxu,hyu)),sup(re_sqrt(hxu,hyl)))); iwert = im_sqrt(hxl,hyl); //Minimum iwert = interval(inf(iwert),sup(im_sqrt(hxl,hyu))); //with maximum *************** cinterval sqrt(const cinterval& z) *** 1508,1512 **** } ! list<cinterval> sqrt_all(const cinterval& z) { #ifdef FILIB_VERSION --- 1511,1515 ---- } ! std::list<cinterval> sqrt_all(const cinterval& z) { #ifdef FILIB_VERSION *************** list<cinterval> sqrt_all(const cinterval *** 1535,1539 **** //right half plane (no intersection) rwert = re_sqrt(hxl,ZERO_INTERVAL()); //Minimum ! rwert = interval(inf(rwert),max(sup(re_sqrt(hxu,hyu)),sup(re_sqrt(hxu,hyl)))); iwert = im_sqrt(hxl,hyl); //Minimum iwert = interval(inf(iwert),sup(im_sqrt(hxl,hyu))); //with maximum --- 1538,1542 ---- //right half plane (no intersection) rwert = re_sqrt(hxl,ZERO_INTERVAL()); //Minimum ! rwert = interval(inf(rwert),std::max(sup(re_sqrt(hxu,hyu)),sup(re_sqrt(hxu,hyl)))); iwert = im_sqrt(hxl,hyl); //Minimum iwert = interval(inf(iwert),sup(im_sqrt(hxl,hyu))); //with maximum *************** list<cinterval> sqrt_all(const cinterval *** 1573,1577 **** cinterval w(rwert,iwert); ! list<cinterval> res; res.push_back( w); res.push_back(-w); --- 1576,1580 ---- cinterval w(rwert,iwert); ! std::list<cinterval> res; res.push_back( w); res.push_back(-w); *************** interval root(const interval& i,unsigned *** 1630,1636 **** } ! list<cinterval> root_all(const cinterval& z,unsigned int n) { ! list<cinterval> res; if( n == 0 ) --- 1633,1639 ---- } ! std::list<cinterval> root_all(const cinterval& z,unsigned int n) { ! std::list<cinterval> res; if( n == 0 ) *************** list<cinterval> root_all(const cinterval *** 1656,1660 **** interval root_abs_z(root(abs(z),n)), arg_z(arg(z)); ! for(int i = 0; i < n; i++) { interval arg_z_plus_2_i_pi_div_n( (arg_z+2*i*PI())/ --- 1659,1663 ---- interval root_abs_z(root(abs(z),n)), arg_z(arg(z)); ! for(unsigned int i = 0; i < n; i++) { interval arg_z_plus_2_i_pi_div_n( (arg_z+2*i*PI())/ *************** interval im_arcsin(const flnumber& x,con *** 1771,1780 **** hilf5(0.0), hilf6(0.0), - nenner(0.0), - zaehler(0.0), ix( abs(interval(x)) ), iy( y ); interval t(0.0),r(0.0); ! flnumber xc( abs(x) ); //Imaginary part --- 1774,1781 ---- hilf5(0.0), hilf6(0.0), ix( abs(interval(x)) ), iy( y ); interval t(0.0),r(0.0); ! flnumber xc( fabs(x) ); //Imaginary part *************** interval real_asin(const cinterval& c) *** 2022,2026 **** interval imag_asin(const cinterval& c) { ! flnumber xl(0.0),xu(0.0),yl(0.0),yu(0.0),maxx(0.0),maxy(0.0),max(0.0); flnumber null( 0.0 ),eins( 1.0 ); bool im_spiegel( false ); --- 2023,2027 ---- interval imag_asin(const cinterval& c) { ! flnumber xl(0.0),xu(0.0),yl(0.0),yu(0.0),maxx(0.0),maxy(0.0); flnumber null( 0.0 ),eins( 1.0 ); bool im_spiegel( false ); *************** void harctan(const interval& x,const int *** 2455,2459 **** { re_arct = interval( inf(re_fun(inf(x),sup(y))),sup(re_fun(sup(x),sup(y))) ); ! im_arct = interval( inf(im_fun(max(sup(x),-inf(x)),inf(y))),sup(im_fun(0.0,sup(y))) ); } else //Intersection in argument --- 2456,2460 ---- { re_arct = interval( inf(re_fun(inf(x),sup(y))),sup(re_fun(sup(x),sup(y))) ); ! im_arct = interval( inf(im_fun(std::max(sup(x),-inf(x)),inf(y))),sup(im_fun(0.0,sup(y))) ); } else //Intersection in argument *************** void h_arctan(const interval& x,const in *** 2495,2499 **** { #ifdef FILIB_VERSION - cinterval z( ZERO_INTERVAL(), ZERO_INTERVAL() ); interval imarct(0.0),rearct(0.0); --- 2496,2499 ---- *************** cinterval pow(const cinterval& z,const i *** 2753,2757 **** } ! list<cinterval> pow_all(const cinterval& z, const interval& n) { interval abs_z(abs(z)); --- 2753,2757 ---- } ! std::list<cinterval> pow_all(const cinterval& z, const interval& n) { interval abs_z(abs(z)); *************** list<cinterval> pow_all(const cinterval& *** 2788,2792 **** ; ! list<cinterval> res; #ifdef FILIB_VERSION --- 2788,2792 ---- ; ! std::list<cinterval> res; #ifdef FILIB_VERSION *************** list<cinterval> pow_all(const cinterval& *** 2822,2826 **** < 0 ) //return entire complex plane { ! list<cinterval> res; #ifdef FILIB_VERSION --- 2822,2826 ---- < 0 ) //return entire complex plane { ! std::list<cinterval> res; #ifdef FILIB_VERSION *************** list<cinterval> pow_all(const cinterval& *** 2841,2845 **** ) { ! list<cinterval> res; res.push_back( cinterval( interval( 0.0, 1.0 ), --- 2841,2845 ---- ) { ! std::list<cinterval> res; res.push_back( cinterval( interval( 0.0, 1.0 ), *************** list<cinterval> pow_all(const cinterval& *** 2865,2869 **** ; ! list<cinterval> res; res.push_back( cinterval( interval( 0.0, d_2 ), --- 2865,2869 ---- ; ! std::list<cinterval> res; res.push_back( cinterval( interval( 0.0, d_2 ), *************** cinterval pow(const cinterval& z,const c *** 2890,2899 **** } ! list<cinterval> pow_all(const cinterval& z, const cinterval& n) { if( Im(n) == ZERO_INTERVAL() ) return pow_all( z, Re(n) ); else { // nonreal eponent: return entire complex plane ! list<cinterval> res; #ifdef FILIB_VERSION --- 2890,2899 ---- } ! std::list<cinterval> pow_all(const cinterval& z, const cinterval& n) { if( Im(n) == ZERO_INTERVAL() ) return pow_all( z, Re(n) ); else { // nonreal eponent: return entire complex plane ! std::list<cinterval> res; #ifdef FILIB_VERSION diff -rcp2N CoStLy/cimath.h newCoStLy/cimath.h *** CoStLy/cimath.h Fri Dec 6 16:04:33 2002 --- newCoStLy/cimath.h Sat Dec 7 19:45:57 2002 *************** cinterval sqrt(const cinterval&); *** 65,70 **** cinterval root(const cinterval&,unsigned int); ! list<cinterval> sqrt_all(const cinterval&); ! list<cinterval> root_all(const cinterval&,unsigned int); cinterval asin(const cinterval&); --- 65,70 ---- cinterval root(const cinterval&,unsigned int); ! std::list<cinterval> sqrt_all(const cinterval&); ! std::list<cinterval> root_all(const cinterval&,unsigned int); cinterval asin(const cinterval&); *************** cinterval pow(const cinterval&,const int *** 82,87 **** cinterval pow(const cinterval&,const cinterval&); ! list<cinterval> pow_all(const cinterval&,const interval&); ! list<cinterval> pow_all(const cinterval&,const cinterval&); #endif --- 82,87 ---- cinterval pow(const cinterval&,const cinterval&); ! std::list<cinterval> pow_all(const cinterval&,const interval&); ! std::list<cinterval> pow_all(const cinterval&,const cinterval&); #endif diff -rcp2N CoStLy/cinterval.cpp newCoStLy/cinterval.cpp *** CoStLy/cinterval.cpp Fri Dec 6 16:04:52 2002 --- newCoStLy/cinterval.cpp Sat Dec 7 19:45:57 2002 *************** *** 26,29 **** --- 26,32 ---- #include "cinterval.h" #include "error.h" + #ifdef FILIB_NAMESPACES + using filib::Double; + #endif cinterval operator - (const cinterval& z) *************** double MinMax(bool minimum, const double *** 166,170 **** // Fallunterscheidung zur Min-,Max- Best: ! double x, ay0 = abs(y0); if( ( sign(b)==sign(y0) )== minimum ) { --- 169,173 ---- // Fallunterscheidung zur Min-,Max- Best: ! double x, ay0 = fabs(y0); if( ( sign(b)==sign(y0) )== minimum ) { *************** Interval abs(const cinterval& z) *** 462,465 **** --- 465,469 ---- } + #ifdef HAS_Complex Complex mid (const cinterval& z) { *************** Complex diam(const cinterval& z) *** 471,474 **** --- 475,479 ---- return Complex( diam(z.re()), diam(z.im()) ); } + #endif //Binary operators *************** cinterval operator / (const cinterval& z *** 590,594 **** // output ! ostream& operator << (ostream& os, const cinterval& z) { os << "(" << z.re() << "," << z.im() << ")"; --- 595,599 ---- // output ! std::ostream& operator << (std::ostream& os, const cinterval& z) { os << "(" << z.re() << "," << z.im() << ")"; diff -rcp2N CoStLy/cinterval.h newCoStLy/cinterval.h *** CoStLy/cinterval.h Fri Dec 6 16:05:10 2002 --- newCoStLy/cinterval.h Sat Dec 7 19:45:57 2002 *************** *** 29,33 **** --- 29,39 ---- #include "Interval.h" //filib++ Header: Macro Version #include <iostream> + #ifdef HAS_Complex #include <Complex.h> + #endif + + #ifdef FILIB_NAMESPACES + typedef filib::Interval Interval; + #endif class cinterval *************** public: *** 45,49 **** --- 51,57 ---- explicit cinterval(const Interval& x) : real_part(x), imag_part(0.0) {} explicit cinterval(const double& d) : real_part(d), imag_part(0.0) {} + #ifdef HAS_Complex explicit cinterval(const Complex& c) : real_part(c.real()), imag_part(c.imag()) {} + #endif cinterval& operator = (const cinterval& z) *************** inline const Interval& Im(const cinterva *** 93,98 **** --- 101,108 ---- Interval abs (const cinterval&); + #ifdef HAS_Complex Complex mid (const cinterval&); Complex diam(const cinterval&); + #endif cinterval operator - (const cinterval&); *************** cinterval operator / (const cinterval&, *** 135,140 **** // output ! ostream& operator << (ostream&, const cinterval&); ! string& operator << (string&, const cinterval&); // --- 145,150 ---- // output ! std::ostream& operator << (std::ostream&, const cinterval&); ! std::string& operator << (std::string&, const cinterval&); // diff -rcp2N CoStLy/error.h newCoStLy/error.h *** CoStLy/error.h Fri Dec 6 16:05:25 2002 --- newCoStLy/error.h Sat Dec 7 19:45:57 2002 *************** *** 29,45 **** #include <stdexcept> ! class syntax_error : public logic_error { public: ! syntax_error(const string& what_arg) : logic_error(what_arg) {} }; ! class division_by_zero : public logic_error { public: ! division_by_zero() : logic_error("") {} }; ! class function_not_defined : public invalid_argument { public: --- 29,45 ---- #include <stdexcept> ! class syntax_error : public std::logic_error { public: ! syntax_error(const std::string& what_arg) : std::logic_error(what_arg) {} }; ! class division_by_zero : public std::logic_error { public: ! division_by_zero() : std::logic_error("") {} }; ! class function_not_defined : public std::invalid_argument { public: *************** class function_not_defined : public inva *** 47,51 **** }; ! class wrong_dimensions : public length_error { public: --- 47,51 ---- }; ! class wrong_dimensions : public std::length_error { public:
Dear Roberto,
thanks for the beta. In order to be able to use it we had to make the corrections you will find in the attached patch. They fall into 5 categories:
Thank you very much. We have updated our code accordingly.
1) namespace errors whereby some std:: qualifications were missing; 2) namespace errors when compiling (as we have to) with FILIB_NAMESPACES defined; 3) incorrect uses of abs() (an int -> int function) instead of fabs() (double -> double); 5) unused variables (causing annoying warnings).
Thank you for the information and for the elimnation of the unused variables. Our compiler (gcc 2.95) did not show us any errors. Which compiler are you using? It seems to be better than ours.
4) undefined type `Complex';
This error is puzzling. The Complex data type is contained in the C++ standard (see File Complex.h). It is simply the double version of the complex data type: typedef complex<double> Complex;
We are now using the revised beta version. However, we could not disable our wrapper around the pow() function, since CoStLy still throws an exception for invocations like, e.g., pow(([-0.5, -0.5],[0, 0]), ([-1, -1],[0, 0])).
In CoStLy 0.3, we have several power functions (see documentation) which perform differently. There is a special function for integer exponents, namely power(Z,n) (Z: cinterval, n: integer) Here, the logarithm is not used. power( ([-0.5, -0.5],[0, 0]), -1 ) is well defined. There are also new root functions that compute all roots, even for intervals that intersect the negative real axis. Do you have a priori information on the type of exponents that occur in your computations? For example, if you knew that only integer exponents occurred then you could replace all calls of pow(Z,P) (Z,P: cintervals) by calls of power(Z,n) (Z: cinterval; n:integer). Or, for real or complex interval exponents, you could use pow_all instead of pow. pow_all computes an enclosure of all solutions of Z^P. pow_all never throws an exception, but sometimes a list of intervals (which cover the solution set) and sometimes the whole complex plane is returned. The pow function is implemented as pow(Z,P) = e^(P*ln(Z)) with the principal branch of the logarithm. For our applications, we need this function to be analytic on its domain. That's why Z must not intersect the negative real axis. Best regards, Ingo Eble and Markus Neher
Markus Neher wrote:
1) namespace errors whereby some std:: qualifications were missing; 2) namespace errors when compiling (as we have to) with FILIB_NAMESPACES defined;
3) incorrect uses of abs() (an int -> int function) instead of fabs() (double -> double); 5) unused variables (causing annoying warnings).
Thank you for the information and for the elimnation of the unused variables. Our compiler (gcc 2.95) did not show us any errors. Which compiler are you using? It seems to be better than ours.
Hi Markus and Ingo, I am using gcc 3.2.1. Other members of the team use other versions, but all above version 3.0 (which is the first one implementing enough standard C++ for our purposes).
4) undefined type `Complex';
This error is puzzling. The Complex data type is contained in the C++ standard (see File Complex.h). It is simply the double version of the complex data type: typedef complex<double> Complex;
The header `Complex.h' is not part of the C++ standard. And, as far as I can tell, it is not part of CoStLy 0.3 beta. Perhaps a packaging problem? Or else, do you mean that `Complex.h' is part of the GCC 2.95 distribution? If so it was a non-standard feature that has been removed. Maybe the solution is simply to 1) #include <complex> in the appropriate CoStLy sources; and 2) use std::complex<double> instead of Complex.
We are now using the revised beta version. However, we could not disable our wrapper around the pow() function, since CoStLy still throws an exception for invocations like, e.g., pow(([-0.5, -0.5],[0, 0]), ([-1, -1],[0, 0])).
In CoStLy 0.3, we have several power functions (see documentation) which perform differently. There is a special function for integer exponents, namely
power(Z,n) (Z: cinterval, n: integer)
Here, the logarithm is not used. power( ([-0.5, -0.5],[0, 0]), -1 ) is well defined. There are also new root functions that compute all roots, even for intervals that intersect the negative real axis.
Do you have a priori information on the type of exponents that occur in your computations? For example, if you knew that only integer exponents occurred then you could replace all calls of pow(Z,P) (Z,P: cintervals) by calls of power(Z,n) (Z: cinterval; n:integer).
No, we do not have this knowledge a priori. Of course we can test the type of the exponent at run-time and then invoke the appropriate function.
Or, for real or complex interval exponents, you could use pow_all instead of pow. pow_all computes an enclosure of all solutions of Z^P. pow_all never throws an exception, but sometimes a list of intervals (which cover the solution set) and sometimes the whole complex plane is returned.
The pow function is implemented as pow(Z,P) = e^(P*ln(Z)) with the principal branch of the logarithm. For our applications, we need this function to be analytic on its domain. That's why Z must not intersect the negative real axis.
Thanks a lot. Please let us know when it is time to beta-test a new version of CoStLy. Cheers 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
Dear Roberto,
do you mean that `Complex.h' is part of the GCC 2.95 distribution? If so it was a non-standard feature that has been removed.
Exactly.
Maybe the solution is simply to
1) #include <complex> in the appropriate CoStLy sources; and 2) use std::complex<double> instead of Complex.
That has been implemented in cinterval.h. You can use the two attached patches to make the changes.
Do you have a priori information on the type of exponents that occur in your computations? For example, if you knew that only integer exponents occurred then you could replace all calls of pow(Z,P) (Z,P: cintervals) by calls of power(Z,n) (Z: cinterval; n:integer).
No, we do not have this knowledge a priori. Of course we can test the type of the exponent at run-time and then invoke the appropriate function.
At the moment, we do not plan to change the pow function. Regards, Ingo and Markus
Markus Neher wrote:
That has been implemented in cinterval.h. You can use the two attached patches to make the changes.
Thanks! I have applied them without problems.
Do you have a priori information on the type of exponents that occur in your computations? For example, if you knew that only integer exponents occurred then you could replace all calls of pow(Z,P) (Z,P: cintervals) by calls of power(Z,n) (Z: cinterval; n:integer).
No, we do not have this knowledge a priori. Of course we can test the type of the exponent at run-time and then invoke the appropriate function.
At the moment, we do not plan to change the pow function.
Well, OK. Thanks again 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
participants (2)
-
Markus Neher -
Roberto Bagnara