Re: pow(double, int)

[Petter Urkedal <petter at matfys dot lth dot se>]

I must have been tired yesterday, forgot to send the
patch for src/Makefile.am.  Here it is again, with
some GNUification of `cmath.cc' and a testcase.  Hope
it's right now.

-petter.


2000-04-05  Petter Urkedal  <petter@matfys.lth.se>

        * src/cmath.cc: New file...
        (pow(*, int)): Define functions...
        * src/Makefile.am (sources): ...register file.
        * src/Makefile.in: Regenerate.
        * bits/std_cmath.h (pow): ...declare functions here, and remove
        old definitions.
	* testsuite/26_numerics/cmath_pow_int.cc: Test int versions
	of pow against FP versions.
Index: bits/std_cmath.h
===================================================================
RCS file: /cvs/libstdc++/libstdc++/bits/std_cmath.h,v
retrieving revision 1.18
diff -c -r1.18 std_cmath.h
*** std_cmath.h	2000/03/21 06:24:04	1.18
--- std_cmath.h	2000/04/05 08:18:06
***************
*** 52,59 ****
      inline double
      abs(double __x) { return ::fabs(__x); }
  
!     inline double
!     pow(double __x, int __i) { return ::pow(__x, static_cast<double>(__i)); }
  
      // float
      inline float
--- 52,59 ----
      inline double
      abs(double __x) { return ::fabs(__x); }
  
!     double
!     pow(double, int);
  
      // float
      inline float
***************
*** 78,86 ****
      pow(float __x, float __y) 
      { return ::pow(static_cast<double>(__x), static_cast<double>(__y)); }
  
!     inline float
!     pow(float __x, int __i) 
!     { return ::pow(static_cast<double>(__x), static_cast<double>(__i)); }
  
  #if _GLIBCPP_HAVE_ABSF
      inline float
--- 78,85 ----
      pow(float __x, float __y) 
      { return ::pow(static_cast<double>(__x), static_cast<double>(__y)); }
  
!     float
!     pow(float, int);
  
  #if _GLIBCPP_HAVE_ABSF
      inline float
***************
*** 239,244 ****
--- 238,246 ----
  
  
      // XXX long double
+ 
+     long double
+     pow(long double, int);
  
  } // std
  
Index: src/Makefile.am
===================================================================
RCS file: /cvs/libstdc++/libstdc++/src/Makefile.am,v
retrieving revision 1.89
diff -c -r1.89 Makefile.am
*** Makefile.am	2000/03/27 19:53:00	1.89
--- Makefile.am	2000/04/05 08:18:06
***************
*** 158,164 ****
  
  sources = \
  	limitsMEMBERS.cc c++io.cc \
! 	complex.cc complexf.cc complexl.cc complex_io.cc \
  	stdexcept.cc ios.cc stdstreams.cc strstream.cc \
  	locale.cc localename.cc  \
  	locale-inst.cc stl-inst.cc misc-inst.cc valarray-inst.cc \
--- 158,164 ----
  
  sources = \
  	limitsMEMBERS.cc c++io.cc \
! 	cmath.cc complex.cc complexf.cc complexl.cc complex_io.cc \
  	stdexcept.cc ios.cc stdstreams.cc strstream.cc \
  	locale.cc localename.cc  \
  	locale-inst.cc stl-inst.cc misc-inst.cc valarray-inst.cc \

// -*- C++ -*- C math library.

// Copyright (C) 2000 Free Software Foundation, Inc.
//
// This file is part of the GNU ISO C++ Library.  This library 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, or (at your option)
// any later version.

// This library 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 library; see the file COPYING.  If not, write to the Free
// Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
// USA.

// As a special exception, you may use this file as part of a free software
// library without restriction.  Specifically, if other files instantiate
// templates or use macros or inline functions from this file, or you compile
// this file and link it with other files to produce an executable, this
// file does not by itself cause the resulting executable to be covered by
// the GNU General Public License.  This exception does not however
// invalidate any other reasons why the executable file might be covered by
// the GNU General Public License.

//
// ISO C++ 14882: 26.5  C library
// Code for signatures not found in the C library
//

#include <bits/std_cmath.h>


namespace std
{

  namespace {
    template<typename T>
    inline T
    pow_helper(T x, unsigned int y)
    {
      T z = y&1? x : 1;
      while (y >>= 1)
        {
          x *= x;
          if (y & 1) z *= x;
        }
      return z;
    }
  }

  float
  pow(float x, int y)
  {
    if (y < 0)
      return 1.0f/pow_helper(x, -y);
    else
      return pow_helper(x, y);
  }

  double
  pow(double x, int y)
  {
    if (y < 0)
      return 1.0/pow_helper(x, -y);
    else
      return pow_helper(x, y);
  }

  long double
  pow(long double x, int y)
  {
    if (y < 0)
      return 1.0l/pow_helper(x, -y);
    else
      return pow_helper(x, y);
  }

}

#include <limits>
#include <iostream>
#include <iomanip>
#include <cmath>

template<typename T>
bool
flteq(T x, T y)
{
  const T eps = std::numeric_limits<T>::epsilon();
  if (std::fabs(x) < eps && std::fabs(y) < eps)
    return true;
  else
    return (std::fabs(x-y) < eps*(std::fabs(x) + std::fabs(y)));
}

template<typename T>
int testit()
{
  const int dig10 = std::numeric_limits<T>::digits10;
  int ret = 0;
  for (int i = 0; i < 1000; ++i)
    {
      T x = 10.0*(drand48()-.5);
      int y = (int)(std::numeric_limits<T>::max_exponent10*(drand48()-.5));
      T z1, z2;
      z1 = std::pow(x, y);
      z2 = std::pow(x, T(y));

      // Note inverted logic due to NaN and Inf
      if (!flteq(z1, z2))
	{
	  std::cout
	    << "Equality failed for pow(" << x << ',' << y << "):\n"
	    << std::setprecision(dig10) << z1 << ' '
	    << std::setprecision(dig10) << z2 << std::endl;
	  ret = 1;
	}
    }
  if (!(std::pow(T(0), 0) == std::pow(T(0), T(0))))
    {
      std::cout
	<< "Versions don't agree at the singularity:\n"
	<< "  pow(0.0, 0) = "
	<< std::pow(T(0), 0) << '\n'
	<< "  pow(0.0, 0.0) =  "
	<< std::pow(T(0), T(0)) << '\n'
	<< std::endl;
    }
  return ret;
}

int main()
{
  int ret = 0;
  std::cout << "Testing float versions.\n";
  ret += testit<float>();
  // XXX uncomment when provided in std
//   std::cout << "Testing double versions.\n";
//   ret += testit<double>();
//   std::cout << "Testing long double versions.\n";
//   ret += testit<long double>();
  return ret;
}