[PURRS-devel] Re: CoStLy 0.3

7 Dec 2002

      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: