Teuchos_ScalarTraits.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//                    Teuchos: Common Tools Package
//                 Copyright (2004) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, 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
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
// @HEADER

// Kris
// 06.18.03 -- Minor formatting changes
//          -- Changed calls to LAPACK objects to use new <OType, SType> templates
// 07.08.03 -- Move into Teuchos package/namespace
// 07.11.03 -- Added ScalarTraits for ARPREC mp_real
// 07.14.03 -- Fixed int rand() function (was set up to return a floating-point style random number)
// 07.17.03 -- Added squareroot() function

// NOTE: Before adding specializations of ScalarTraits, make sure that they do not duplicate 
// specializations already present in PyTrilinos (see packages/PyTrilinos/src/Teuchos_Traits.i)

// NOTE: halfPrecision and doublePrecision are not currently implemented for ARPREC, GMP or the ordinal types (e.g., int, char)

#ifndef _TEUCHOS_SCALARTRAITS_HPP_
#define _TEUCHOS_SCALARTRAITS_HPP_

#include "Teuchos_ConfigDefs.hpp"

#ifdef HAVE_TEUCHOS_ARPREC
#include <arprec/mp_real.h>
#endif

#ifdef HAVE_TEUCHOS_QD
#include <qd/qd_real.h>
#include <qd/dd_real.h>
#endif

#ifdef HAVE_TEUCHOS_GNU_MP
#include <gmp.h>
#include <gmpxx.h>
#endif


#include "Teuchos_ScalarTraitsDecl.hpp"


namespace Teuchos {


#ifndef DOXYGEN_SHOULD_SKIP_THIS


void throwScalarTraitsNanInfError( const std::string &errMsg );


template<class Scalar>
bool generic_real_isnaninf(const Scalar &x)
{
  typedef std::numeric_limits<Scalar> NL;
  // IEEE says this should fail for NaN (not all compilers do not obey IEEE)
  const Scalar tol = 1.0; // Any (bounded) number should do!
  if (!(x <= tol) && !(x > tol)) return true;
  // Use fact that Inf*0 = NaN (again, all compilers do not obey IEEE)
  Scalar z = static_cast<Scalar>(0.0) * x;
  if (!(z <= tol) && !(z > tol)) return true;
  // As a last result use comparisons
  if (x == NL::infinity() || x == -NL::infinity()) return true;
  // We give up and assume the number is finite
  return false;
}


#define TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR( VALUE, MSG ) \
  if (isnaninf(VALUE)) { \
    std::ostringstream omsg; \
    omsg << MSG; \
    Teuchos::throwScalarTraitsNanInfError(omsg.str());  \
  }


template<>
struct ScalarTraits<char>
{
  typedef char magnitudeType;
  typedef char halfPrecision;
  typedef char doublePrecision;
  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
  static inline magnitudeType magnitude(char a) { return static_cast<char>(std::fabs(static_cast<double>(a))); }
  static inline char zero()  { return 0; }
  static inline char one()   { return 1; }
  static inline char conjugate(char x) { return x; }
  static inline char real(char x) { return x; }
  static inline char imag(char) { return 0; }
  static inline bool isnaninf(char ) { return false; }
  static inline void seedrandom(unsigned int s) { 
    std::srand(s); 
#ifdef __APPLE__
    // throw away first random number to address bug 3655
    // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
    random();
#endif
  }
  //static inline char random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
  static inline char random() { return std::rand(); } // RAB: This version should be used for an unsigned char, not char
  static inline std::string name() { return "char"; }
  static inline char squareroot(char x) { return (char) std::sqrt((double) x); }
  static inline char pow(char x, char y) { return (char) std::pow((double)x,(double)y); }
};


template<>
struct ScalarTraits<short int>
{
  typedef short int magnitudeType;
  typedef short int halfPrecision;
  typedef short int doublePrecision;
  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
  static inline magnitudeType magnitude(short int a) { return static_cast<short int>(std::fabs(static_cast<double>(a))); }
  static inline short int zero()  { return 0; }
  static inline short int one()   { return 1; }
  static inline short int conjugate(short int x) { return x; }
  static inline short int real(short int x) { return x; }
  static inline short int imag(short int) { return 0; }
  static inline bool isnaninf(short int) { return false; }
  static inline void seedrandom(unsigned int s) { 
    std::srand(s); 
#ifdef __APPLE__
    // throw away first random number to address bug 3655
    // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
    random();
#endif
  }
  //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
  static inline short int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
  static inline std::string name() { return "short int"; }
  static inline short int squareroot(short int x) { return (short int) std::sqrt((double) x); }
  static inline short int pow(short int x, short int y) { return (short int) std::pow((double)x,(double)y); }
};


template<>
struct ScalarTraits<int>
{
  typedef int magnitudeType;
  typedef int halfPrecision;
  typedef int doublePrecision;
  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
  static inline magnitudeType magnitude(int a) { return static_cast<int>(std::fabs(static_cast<double>(a))); }
  static inline int zero()  { return 0; }
  static inline int one()   { return 1; }
  static inline int conjugate(int x) { return x; }
  static inline int real(int x) { return x; }
  static inline int imag(int) { return 0; }
  static inline bool isnaninf(int) { return false; }
  static inline void seedrandom(unsigned int s) { 
    std::srand(s); 
#ifdef __APPLE__
    // throw away first random number to address bug 3655
    // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
    random();
#endif
  }
  //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
  static inline int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
  static inline std::string name() { return "int"; }
  static inline int squareroot(int x) { return (int) std::sqrt((double) x); }
  static inline int pow(int x, int y) { return (int) std::pow((double)x,(double)y); }
};


template<>
struct ScalarTraits<unsigned int>
{
  typedef unsigned int magnitudeType;
  typedef unsigned int halfPrecision;
  typedef unsigned int doublePrecision;
  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
  static inline magnitudeType magnitude(unsigned int a) { return static_cast<unsigned int>(std::fabs(static_cast<double>(a))); }
  static inline unsigned int zero()  { return 0; }
  static inline unsigned int one()   { return 1; }
  static inline unsigned int conjugate(unsigned int x) { return x; }
  static inline unsigned int real(unsigned int x) { return x; }
  static inline unsigned int imag(unsigned int) { return 0; }
  static inline void seedrandom(unsigned int s) { 
    std::srand(s); 
#ifdef __APPLE__
    // throw away first random number to address bug 3655
    // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
    random();
#endif
  }
  //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
  static inline unsigned int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
  static inline std::string name() { return "unsigned int"; }
  static inline unsigned int squareroot(unsigned int x) { return (unsigned int) std::sqrt((double) x); }
  static inline unsigned int pow(unsigned int x, unsigned int y) { return (unsigned int) std::pow((double)x,(double)y); }
};


template<>
struct ScalarTraits<long int>
{
  typedef long int magnitudeType;
  typedef long int halfPrecision;
  typedef long int doublePrecision;
  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
  static inline magnitudeType magnitude(long int a) { return static_cast<long int>(std::fabs(static_cast<double>(a))); }
  static inline long int zero()  { return 0; }
  static inline long int one()   { return 1; }
  static inline long int conjugate(long int x) { return x; }
  static inline long int real(long int x) { return x; }
  static inline long int imag(long int) { return 0; }
  static inline void seedrandom(unsigned int s) { 
    std::srand(s); 
#ifdef __APPLE__
    // throw away first random number to address bug 3655
    // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
    random();
#endif
  }
  //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
  static inline long int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
  static inline std::string name() { return "long int"; }
  static inline long int squareroot(long int x) { return (long int) std::sqrt((double) x); }
  static inline long int pow(long int x, long int y) { return (long int) std::pow((double)x,(double)y); }
};


template<>
struct ScalarTraits<long unsigned int>
{
  typedef long unsigned int magnitudeType;
  typedef long unsigned int halfPrecision;
  typedef long unsigned int doublePrecision;
  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
  static inline magnitudeType magnitude(long unsigned int a) { return static_cast<long unsigned int>(std::fabs(static_cast<double>(a))); }
  static inline long unsigned int zero()  { return 0; }
  static inline long unsigned int one()   { return 1; }
  static inline long unsigned int conjugate(long unsigned int x) { return x; }
  static inline long unsigned int real(long unsigned int x) { return x; }
  static inline long unsigned int imag(long unsigned int) { return 0; }
  static inline void seedrandom(unsigned int s) { 
    std::srand(s); 
#ifdef __APPLE__
    // throw away first random number to address bug 3655
    // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
    random();
#endif
  }
  //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
  static inline long unsigned int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
  static inline std::string name() { return "long unsigned int"; }
  static inline long unsigned int squareroot(long unsigned int x) { return (long unsigned int) std::sqrt((double) x); }
  static inline long unsigned int pow(long unsigned int x, long unsigned int y) { return (long unsigned int) std::pow((double)x,(double)y); }
};


#ifdef HAVE_TEUCHOS_LONG_LONG_INT
template<>
struct ScalarTraits<long long int>
{
  typedef long long int magnitudeType;
  typedef long long int halfPrecision;
  typedef long long int doublePrecision;
  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
  static inline magnitudeType magnitude(long long int a) { return static_cast<long long int>(std::fabs(static_cast<double>(a))); }
  static inline long long int zero()  { return 0; }
  static inline long long int one()   { return 1; }
  static inline long long int conjugate(long long int x) { return x; }
  static inline long long int real(long long int x) { return x; }
  static inline long long int imag(long long int) { return 0; }
  static inline void seedrandom(unsigned int s) { 
    std::srand(s); 
#ifdef __APPLE__
    // throw away first random number to address bug 3655
    // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
    random();
#endif
  }
  //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
  static inline long long int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
  static inline std::string name() { return "long long int"; }
  static inline long long int squareroot(long long int x) { return (long long int) std::sqrt((double) x); }
  static inline long long int pow(long long int x, long long int y) { return (long long int) std::pow((double)x,(double)y); }
};
#endif // HAVE_TEUCHOS_LONG_LONG_INT


#ifndef __sun
extern TEUCHOS_LIB_DLL_EXPORT const float flt_nan;
#endif

 
template<>
struct ScalarTraits<float>
{
  typedef float magnitudeType;
  typedef float halfPrecision; // should become IEEE754-2008 binary16 or fp16 later, perhaps specified at configure according to architectural support
  typedef double doublePrecision;
  static const bool isComplex = false;
  static const bool isOrdinal = false;
  static const bool isComparable = true;
  static const bool hasMachineParameters = true;
  static inline float eps()   {
    return std::numeric_limits<float>::epsilon();
  }
  static inline float sfmin() {
    return std::numeric_limits<float>::min();
  }
  static inline float base()  {
    return static_cast<float>(std::numeric_limits<float>::radix);
  }
  static inline float prec()  {
    return eps()*base();
  }
  static inline float t()     {
    return static_cast<float>(std::numeric_limits<float>::digits);
  }
  static inline float rnd()   {
    return ( std::numeric_limits<float>::round_style == std::round_to_nearest ? one() : zero() );
  }
  static inline float emin()  {
    return static_cast<float>(std::numeric_limits<float>::min_exponent);
  }
  static inline float rmin()  {
    return std::numeric_limits<float>::min();
  }
  static inline float emax()  {
    return static_cast<float>(std::numeric_limits<float>::max_exponent);
  }
  static inline float rmax()  {
    return std::numeric_limits<float>::max();
  }
  static inline magnitudeType magnitude(float a)
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
        a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
#endif      
      return std::fabs(a);
    }    
  static inline float zero()  { return(0.0f); }
  static inline float one()   { return(1.0f); }    
  static inline float conjugate(float x)   { return(x); }    
  static inline float real(float x) { return x; }
  static inline float imag(float) { return zero(); }
  static inline float nan() {
#ifdef __sun
    return 0.0f/std::sin(0.0f);
#else
    return flt_nan;
#endif
  }
  static inline bool isnaninf(float x) {
    return generic_real_isnaninf<float>(x);
  }
  static inline void seedrandom(unsigned int s) { 
    std::srand(s); 
#ifdef __APPLE__
    // throw away first random number to address bug 3655
    // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
    random();
#endif
  }
  static inline float random() { float rnd = (float) std::rand() / RAND_MAX; return (-1.0f + 2.0f * rnd); }
  static inline std::string name() { return "float"; }
  static inline float squareroot(float x)
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
        x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
#endif
      errno = 0;
      const float rtn = std::sqrt(x);
      if (errno)
        return nan();
      return rtn;
    }
  static inline float pow(float x, float y) { return std::pow(x,y); }
};


#ifndef __sun
extern TEUCHOS_LIB_DLL_EXPORT const double dbl_nan;
#endif

 
template<>
struct ScalarTraits<double>
{
  typedef double magnitudeType;
  typedef float halfPrecision;
  /* there are different options as to how to double "double"
     - QD's DD (if available)
     - ARPREC
     - GNU MP
     - a true hardware quad

     in the shortterm, this should be specified at configure time. I have inserted a configure-time option (--enable-teuchos-double-to-dd) 
     which uses QD's DD when available. This must be used alongside --enable-teuchos-qd.
   */
#if defined(HAVE_TEUCHOS_DOUBLE_TO_QD)
  typedef dd_real doublePrecision;
#elif defined(HAVE_TEUCHOS_DOUBLE_TO_ARPREC)
  typedef mp_real doublePrecision;
#else
  typedef double doublePrecision;     // don't double "double" in this case
#endif
  static const bool isComplex = false;
  static const bool isOrdinal = false;
  static const bool isComparable = true;
  static const bool hasMachineParameters = true;
  static inline double eps()   {
    return std::numeric_limits<double>::epsilon();
  }
  static inline double sfmin() {
    return std::numeric_limits<double>::min();
  }
  static inline double base()  {
    return std::numeric_limits<double>::radix;
  }
  static inline double prec()  {
    return eps()*base();
  }
  static inline double t()     {
    return std::numeric_limits<double>::digits;
  }
  static inline double rnd()   {
    return ( std::numeric_limits<double>::round_style == std::round_to_nearest ? double(1.0) : double(0.0) );
  }
  static inline double emin()  {
    return std::numeric_limits<double>::min_exponent;
  }
  static inline double rmin()  {
    return std::numeric_limits<double>::min();
  }
  static inline double emax()  {
    return std::numeric_limits<double>::max_exponent;
  }
  static inline double rmax()  {
    return std::numeric_limits<double>::max();
  }
  static inline magnitudeType magnitude(double a)
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
        a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
#endif      
      return std::fabs(a);
    }
  static inline double zero()  { return 0.0; }
  static inline double one()   { return 1.0; }
  static inline double conjugate(double x)   { return(x); }    
  static inline double real(double x) { return(x); }
  static inline double imag(double) { return(0); }
  static inline double nan() {
#ifdef __sun
    return 0.0/std::sin(0.0);
#else
    return dbl_nan;
#endif
  }
  static inline bool isnaninf(double x) {
    return generic_real_isnaninf<double>(x);
  }
  static inline void seedrandom(unsigned int s) { 
    std::srand(s); 
#ifdef __APPLE__
    // throw away first random number to address bug 3655
    // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
    random();
#endif
  }
  static inline double random() { double rnd = (double) std::rand() / RAND_MAX; return (double)(-1.0 + 2.0 * rnd); }
  static inline std::string name() { return "double"; }
  static inline double squareroot(double x)
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
        x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
#endif      
      errno = 0;
      const double rtn = std::sqrt(x);
      if (errno)
        return nan();
      return rtn;
    }
  static inline double pow(double x, double y) { return std::pow(x,y); }
};


#ifdef HAVE_TEUCHOS_QD


bool operator&&(const dd_real &a, const dd_real &b);
bool operator&&(const qd_real &a, const qd_real &b);


template<>
struct ScalarTraits<dd_real>
{
  typedef dd_real magnitudeType;
  typedef double halfPrecision;
  typedef qd_real doublePrecision;
  static const bool isComplex = false;
  static const bool isOrdinal = false;
  static const bool isComparable = true;
  static const bool hasMachineParameters = true;
  static inline dd_real eps()   { return std::numeric_limits<dd_real>::epsilon(); }
  static inline dd_real sfmin() { return std::numeric_limits<dd_real>::min(); }
  static inline dd_real base()  { return std::numeric_limits<dd_real>::radix; }
  static inline dd_real prec()  { return eps()*base(); }
  static inline dd_real t()     { return std::numeric_limits<dd_real>::digits; }
  static inline dd_real rnd()   { return ( std::numeric_limits<dd_real>::round_style == std::round_to_nearest ? dd_real(1.0) : dd_real(0.0) ); }
  static inline dd_real emin()  { return std::numeric_limits<dd_real>::min_exponent; }
  static inline dd_real rmin()  { return std::numeric_limits<dd_real>::min(); }
  static inline dd_real emax()  { return std::numeric_limits<dd_real>::max_exponent; }
  static inline dd_real rmax()  { return std::numeric_limits<dd_real>::max(); }
  static inline magnitudeType magnitude(dd_real a)
  {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
        a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
#endif      
      return abs(a);
  }
  static inline dd_real zero()  { return dd_real(0.0); }
  static inline dd_real one()   { return dd_real(1.0); }
  static inline dd_real conjugate(dd_real x)   { return(x); }    
  static inline dd_real real(dd_real x) { return x ; }
  static inline dd_real imag(dd_real) { return zero(); }
  static inline dd_real nan() { return dd_real::_nan; }
  static inline bool isnaninf(dd_real x) { return isnan(x) || isinf(x); }
  static inline void seedrandom(unsigned int s) {
    // ddrand() uses std::rand(), so the std::srand() is our seed
    std::srand(s); 
#ifdef __APPLE__
    // throw away first random number to address bug 3655
    // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
    random();
#endif
  }
  static inline dd_real random() { return ddrand(); }
  static inline std::string name() { return "dd_real"; }
  static inline dd_real squareroot(dd_real x)
  {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
        x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
#endif      
      return sqrt(x);
  }
  static inline dd_real pow(dd_real x, dd_real y) { return pow(x,y); }
};


template<>
struct ScalarTraits<qd_real>
{
  typedef qd_real magnitudeType;
  typedef dd_real halfPrecision;
  typedef qd_real doublePrecision;
  static const bool isComplex = false;
  static const bool isOrdinal = false;
  static const bool isComparable = true;
  static const bool hasMachineParameters = true;
  static inline qd_real eps()   { return std::numeric_limits<qd_real>::epsilon(); }
  static inline qd_real sfmin() { return std::numeric_limits<qd_real>::min(); }
  static inline qd_real base()  { return std::numeric_limits<qd_real>::radix; }
  static inline qd_real prec()  { return eps()*base(); }
  static inline qd_real t()     { return std::numeric_limits<qd_real>::digits; }
  static inline qd_real rnd()   { return ( std::numeric_limits<qd_real>::round_style == std::round_to_nearest ? qd_real(1.0) : qd_real(0.0) ); }
  static inline qd_real emin()  { return std::numeric_limits<qd_real>::min_exponent; }
  static inline qd_real rmin()  { return std::numeric_limits<qd_real>::min(); }
  static inline qd_real emax()  { return std::numeric_limits<qd_real>::max_exponent; }
  static inline qd_real rmax()  { return std::numeric_limits<qd_real>::max(); }
  static inline magnitudeType magnitude(qd_real a)
  {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
        a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
#endif      
      return abs(a);
  }
  static inline qd_real zero()  { return qd_real(0.0); }
  static inline qd_real one()   { return qd_real(1.0); }
  static inline qd_real conjugate(qd_real x)   { return(x); }    
  static inline qd_real real(qd_real x) { return x ; }
  static inline qd_real imag(qd_real) { return zero(); }
  static inline qd_real nan() { return qd_real::_nan; }
  static inline bool isnaninf(qd_real x) { return isnan(x) || isinf(x); }
  static inline void seedrandom(unsigned int s) {
    // qdrand() uses std::rand(), so the std::srand() is our seed
    std::srand(s); 
#ifdef __APPLE__
    // throw away first random number to address bug 3655
    // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
    random();
#endif
  }
  static inline qd_real random() { return qdrand(); }
  static inline std::string name() { return "qd_real"; }
  static inline qd_real squareroot(qd_real x)
  {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
        x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
#endif      
      return sqrt(x);
  }
  static inline qd_real pow(qd_real x, qd_real y) { return pow(x,y); }
};


#endif  // HAVE_TEUCHOS_QD


#ifdef HAVE_TEUCHOS_GNU_MP


extern gmp_randclass gmp_rng; 


/* Regarding halfPrecision, doublePrecision and mpf_class: 
   Because the precision of an mpf_class float is not determined by the data type, 
   there is no way to fill the typedefs for this object. 

   Instead, we could create new data classes (e.g., Teuchos::MPF128, Teuchos::MPF256) for 
   commonly used levels of precision, and fill out ScalarTraits for these. This would allow us
   to typedef the promotions and demotions in the appropriate way. These classes would serve to 
   wrap an mpf_class object, calling the constructor for the appropriate precision, exposing the 
   arithmetic routines but hiding the precision-altering routines.
   
   Alternatively (perhaps, preferably), would could create a single class templated on the precision (e.g., Teuchos::MPF<N>). 
   Then we have a single (partially-specialized) implementation of ScalarTraits. This class, as above, must expose all of the 
   operations expected of a scalar type; however, most of these can be trivially stolen from the gmpcxx.h header file

   CGB/RAB, 01/05/2009
*/
template<>
struct ScalarTraits<mpf_class>
{
  typedef mpf_class magnitudeType;
  typedef mpf_class halfPrecision;
  typedef mpf_class doublePrecision;
  static const bool isComplex = false;
  static const bool hasMachineParameters = false;
  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
  static magnitudeType magnitude(mpf_class a) { return std::abs(a); }
  static inline mpf_class zero() { mpf_class zero = 0.0; return zero; }
  static inline mpf_class one() { mpf_class one = 1.0; return one; }    
  static inline mpf_class conjugate(mpf_class x) { return x; }
  static inline mpf_class real(mpf_class x) { return(x); }
  static inline mpf_class imag(mpf_class x) { return(0); }
  static inline bool isnaninf(mpf_class x) { return false; } // mpf_class currently can't handle nan or inf!
  static inline void seedrandom(unsigned int s) { 
    unsigned long int seedVal = static_cast<unsigned long int>(s);
    gmp_rng.seed( seedVal );  
  }
  static inline mpf_class random() { 
    return gmp_rng.get_f(); 
  }
  static inline std::string name() { return "mpf_class"; }
  static inline mpf_class squareroot(mpf_class x) { return std::sqrt(x); }
  static inline mpf_class pow(mpf_class x, mpf_class y) { return pow(x,y); }
  // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
};

#endif  // HAVE_TEUCHOS_GNU_MP

#ifdef HAVE_TEUCHOS_ARPREC

/* See discussion above for mpf_class, regarding halfPrecision and doublePrecision. Something similar will need to be done
   for ARPREC. */
template<>
struct ScalarTraits<mp_real>
{
  typedef mp_real magnitudeType;
  typedef double halfPrecision;
  typedef mp_real doublePrecision;
  static const bool isComplex = false;
  static const bool isComparable = true;
  static const bool isOrdinal = false;
  static const bool hasMachineParameters = false;
  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
  static magnitudeType magnitude(mp_real a) { return abs(a); }
  static inline mp_real zero() { mp_real zero = 0.0; return zero; }
  static inline mp_real one() { mp_real one = 1.0; return one; }    
  static inline mp_real conjugate(mp_real x) { return x; }
  static inline mp_real real(mp_real x) { return(x); }
  static inline mp_real imag(mp_real x) { return zero(); }
  static inline bool isnaninf(mp_real x) { return false; } // ToDo: Change this?
  static inline void seedrandom(unsigned int s) { 
    long int seedVal = static_cast<long int>(s);
    srand48(seedVal);
  }
  static inline mp_real random() { return mp_rand(); }
  static inline std::string name() { return "mp_real"; }
  static inline mp_real squareroot(mp_real x) { return sqrt(x); }
  static inline mp_real pow(mp_real x, mp_real y) { return pow(x,y); }
  // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
};

  
#endif // HAVE_TEUCHOS_ARPREC

 
#ifdef HAVE_TEUCHOS_COMPLEX


// Partial specialization for std::complex numbers templated on real type T
template<class T> 
struct ScalarTraits<
  std::complex<T>
>
{
  typedef std::complex<T>  ComplexT;
  typedef std::complex<typename ScalarTraits<T>::halfPrecision> halfPrecision;
  typedef std::complex<typename ScalarTraits<T>::doublePrecision> doublePrecision;
  typedef typename ScalarTraits<T>::magnitudeType magnitudeType;
  static const bool isComplex = true;
  static const bool isOrdinal = ScalarTraits<T>::isOrdinal;
  static const bool isComparable = false;
  static const bool hasMachineParameters = true;
  static inline magnitudeType eps()          { return ScalarTraits<magnitudeType>::eps(); }
  static inline magnitudeType sfmin()        { return ScalarTraits<magnitudeType>::sfmin(); }
  static inline magnitudeType base()         { return ScalarTraits<magnitudeType>::base(); }
  static inline magnitudeType prec()         { return ScalarTraits<magnitudeType>::prec(); }
  static inline magnitudeType t()            { return ScalarTraits<magnitudeType>::t(); }
  static inline magnitudeType rnd()          { return ScalarTraits<magnitudeType>::rnd(); }
  static inline magnitudeType emin()         { return ScalarTraits<magnitudeType>::emin(); }
  static inline magnitudeType rmin()         { return ScalarTraits<magnitudeType>::rmin(); }
  static inline magnitudeType emax()         { return ScalarTraits<magnitudeType>::emax(); }
  static inline magnitudeType rmax()         { return ScalarTraits<magnitudeType>::rmax(); }
  static magnitudeType magnitude(ComplexT a)
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
        a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
#endif      
      return std::abs(a);
    }
  static inline ComplexT zero()              { return ComplexT(ScalarTraits<magnitudeType>::zero(),ScalarTraits<magnitudeType>::zero()); }
  static inline ComplexT one()               { return ComplexT(ScalarTraits<magnitudeType>::one(),ScalarTraits<magnitudeType>::zero()); }
  static inline ComplexT conjugate(ComplexT a){ return ComplexT(a.real(),-a.imag()); }
  static inline magnitudeType real(ComplexT a) { return a.real(); }
  static inline magnitudeType imag(ComplexT a) { return a.imag(); }
  static inline ComplexT nan()               { return ComplexT(ScalarTraits<magnitudeType>::nan(),ScalarTraits<magnitudeType>::nan()); }
  static inline bool isnaninf(ComplexT x)    { return ScalarTraits<magnitudeType>::isnaninf(x.real()) || ScalarTraits<magnitudeType>::isnaninf(x.imag()); }
  static inline void seedrandom(unsigned int s) { ScalarTraits<magnitudeType>::seedrandom(s); }
  static inline ComplexT random()
    {
      const T rnd1 = ScalarTraits<magnitudeType>::random();
      const T rnd2 = ScalarTraits<magnitudeType>::random();
      return ComplexT(rnd1,rnd2);
    }
  static inline std::string name() { return std::string("std::complex<")+std::string(ScalarTraits<magnitudeType>::name())+std::string(">"); }
  // This will only return one of the square roots of x, the other can be obtained by taking its conjugate
  static inline ComplexT squareroot(ComplexT x)
    {
#ifdef TEUCHOS_DEBUG
      TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
        x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
#endif
      typedef ScalarTraits<magnitudeType>  STMT;
      const T r  = x.real(), i = x.imag(), zero = STMT::zero(), two = 2.0;
      const T a  = STMT::squareroot((r*r)+(i*i));
      const T nr = STMT::squareroot((a+r)/two);
      const T ni = ( i == zero ? zero : STMT::squareroot((a-r)/two) );
      return ComplexT(nr,ni);
    }
    // 2010/03/19: rabartl: Above, I had to add the check for i == zero
    // to avoid a returned NaN on the Intel 10.1 compiler.  For some
    // reason, having these two squareroot calls in a row produce a NaN
    // in an optimized build with this compiler.  Amazingly, when I put
    // in print statements (i.e. std::cout << ...) the NaN went away and
    // the second squareroot((a-r)/two) returned zero correctly.  I
    // guess that calling the output routine flushed the registers or
    // something and restarted the squareroot rountine for this compiler
    // or something similar.  Actually, due to roundoff, it is possible that a-r
    // might be slightly less than zero (i.e. -1e-16) and that would cause
    // a possbile NaN return.  THe above if test is the right thing to do
    // I think and is very cheap.
  static inline ComplexT pow(ComplexT x, ComplexT y) { return pow(x,y); }
};


#endif //  HAVE_TEUCHOS_COMPLEX


#endif // DOXYGEN_SHOULD_SKIP_THIS


} // Teuchos namespace


#endif // _TEUCHOS_SCALARTRAITS_HPP_