Teuchos - Trilinos Tools Package  Version of the Day
Teuchos_ScalarTraits.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Teuchos: Common Tools Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 // NOTE: Before adding specializations of ScalarTraits, make sure that they do
43 // not duplicate specializations already present in PyTrilinos (see
44 // packages/PyTrilinos/src/Teuchos_Traits.i)
45 
46 // NOTE: halfPrecision and doublePrecision are not currently implemented for
47 // ARPREC, GMP or the ordinal types (e.g., int, char)
48 
49 #ifndef _TEUCHOS_SCALARTRAITS_HPP_
50 #define _TEUCHOS_SCALARTRAITS_HPP_
51 
56 #include "Teuchos_ConfigDefs.hpp"
57 
58 #ifdef HAVE_TEUCHOS_ARPREC
59 #include <arprec/mp_real.h>
60 #endif
61 
62 #ifdef HAVE_TEUCHOSCORE_QUADMATH
63 #include <quadmath.h>
64 
65 // Teuchos_ConfigDefs.hpp includes <iostream>, which includes
66 // <ostream>. If this ever changes, include <ostream> here.
67 
68 namespace std {
69 
79 ostream&
80 operator<< (std::ostream& out, const __float128& x);
81 
91 istream&
92 operator>> (std::istream& in, __float128& x);
93 
94 } // namespace std
95 
96 #endif // HAVE_TEUCHOSCORE_QUADMATH
97 
98 #ifdef HAVE_TEUCHOS_QD
99 #include <qd/qd_real.h>
100 #include <qd/dd_real.h>
101 #endif
102 
103 #ifdef HAVE_TEUCHOS_GNU_MP
104 #include <gmp.h>
105 #include <gmpxx.h>
106 #endif
107 
108 
109 #include "Teuchos_ScalarTraitsDecl.hpp"
110 
111 
112 namespace Teuchos {
113 
114 
115 #ifndef DOXYGEN_SHOULD_SKIP_THIS
116 
117 
118 TEUCHOSCORE_LIB_DLL_EXPORT
119 void throwScalarTraitsNanInfError( const std::string &errMsg );
120 
121 
122 template<class Scalar>
123 bool generic_real_isnaninf(const Scalar &x)
124 {
125 #ifdef HAVE_TEUCHOSCORE_CXX11
126  if (std::isnan(x)) return true;
127  if (std::isinf(x)) return true;
128  return false;
129 #else
130  typedef std::numeric_limits<Scalar> STD_NL;
131  // IEEE says this should fail for NaN (not all compilers do not obey IEEE)
132  const Scalar tol = 1.0; // Any (bounded) number should do!
133  if (!(x <= tol) && !(x > tol)) return true;
134  // Use fact that Inf*0 = NaN (again, all compilers do not obey IEEE)
135  Scalar z = static_cast<Scalar>(0.0) * x;
136  if (!(z <= tol) && !(z > tol)) return true;
137  // As a last result use comparisons
138  if (x == STD_NL::infinity() || x == -STD_NL::infinity()) return true;
139  // We give up and assume the number is finite
140  return false;
141 #endif
142 }
143 
144 
145 #define TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR( VALUE, MSG ) \
146  if (isnaninf(VALUE)) { \
147  std::ostringstream omsg; \
148  omsg << MSG; \
149  Teuchos::throwScalarTraitsNanInfError(omsg.str()); \
150  }
151 
152 
153 template<>
154 struct ScalarTraits<char>
155 {
156  typedef char magnitudeType;
157  typedef char halfPrecision;
158  typedef char doublePrecision;
159  static const bool isComplex = false;
160  static const bool isOrdinal = true;
161  static const bool isComparable = true;
162  static const bool hasMachineParameters = false;
163  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
164  static inline magnitudeType magnitude(char a) { return static_cast<char>(std::fabs(static_cast<double>(a))); }
165  static inline char zero() { return 0; }
166  static inline char one() { return 1; }
167  static inline char conjugate(char x) { return x; }
168  static inline char real(char x) { return x; }
169  static inline char imag(char) { return 0; }
170  static inline bool isnaninf(char ) { return false; }
171  static inline void seedrandom(unsigned int s) {
172  std::srand(s);
173 #ifdef __APPLE__
174  // throw away first random number to address bug 3655
175  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
176  random();
177 #endif
178  }
179  //static inline char random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
180  static inline char random() { return std::rand(); } // RAB: This version should be used for an unsigned char, not char
181  static inline std::string name() { return "char"; }
182  static inline char squareroot(char x) { return (char) std::sqrt((double) x); }
183  static inline char pow(char x, char y) { return (char) std::pow((double)x,(double)y); }
184  static inline char log(char x) { return static_cast<char> (std::log (static_cast<double> (x))); }
185  static inline char log10(char x) { return static_cast<char> (std::log10 (static_cast<double> (x))); }
186 };
187 
188 
189 template<>
190 struct ScalarTraits<short int>
191 {
192  typedef short int magnitudeType;
193  typedef short int halfPrecision;
194  typedef short int doublePrecision;
195  static const bool isComplex = false;
196  static const bool isOrdinal = true;
197  static const bool isComparable = true;
198  static const bool hasMachineParameters = false;
199  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
200  static inline magnitudeType magnitude(short int a) { return static_cast<short int>(std::fabs(static_cast<double>(a))); }
201  static inline short int zero() { return 0; }
202  static inline short int one() { return 1; }
203  static inline short int conjugate(short int x) { return x; }
204  static inline short int real(short int x) { return x; }
205  static inline short int imag(short int) { return 0; }
206  static inline bool isnaninf(short int) { return false; }
207  static inline void seedrandom(unsigned int s) {
208  std::srand(s);
209 #ifdef __APPLE__
210  // throw away first random number to address bug 3655
211  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
212  random();
213 #endif
214  }
215  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
216  static inline short int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
217  static inline std::string name() { return "short int"; }
218  static inline short int squareroot(short int x) { return (short int) std::sqrt((double) x); }
219  static inline short int pow(short int x, short int y) { return (short int) std::pow((double)x,(double)y); }
220  static inline short int log(short int x) { return static_cast<short int> (std::log (static_cast<double> (x))); }
221  static inline short int log10(short int x) { return static_cast<short int> (std::log10 (static_cast<double> (x))); }
222 };
223 
224 template<>
225 struct ScalarTraits<unsigned short int>
226 {
227  typedef unsigned short int magnitudeType;
228  typedef unsigned short int halfPrecision;
229  typedef unsigned short int doublePrecision;
230  static const bool isComplex = false;
231  static const bool isOrdinal = true;
232  static const bool isComparable = true;
233  static const bool hasMachineParameters = false;
234  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
235  static inline magnitudeType magnitude(unsigned short int a) { return static_cast<unsigned short int>(std::fabs(static_cast<double>(a))); }
236  static inline unsigned short int zero() { return 0; }
237  static inline unsigned short int one() { return 1; }
238  static inline unsigned short int conjugate(unsigned short int x) { return x; }
239  static inline unsigned short int real(unsigned short int x) { return x; }
240  static inline unsigned short int imag(unsigned short int) { return 0; }
241  static inline bool isnaninf(unsigned short int) { return false; }
242  static inline void seedrandom(unsigned int s) {
243  std::srand(s);
244 #ifdef __APPLE__
245  // throw away first random number to address bug 3655
246  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
247  random();
248 #endif
249  }
250  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
251  static inline unsigned short int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
252  static inline std::string name() { return "unsigned short int"; }
253  static inline unsigned short int squareroot(unsigned short int x) { return (unsigned short int) std::sqrt((double) x); }
254  static inline unsigned short int pow(unsigned short int x, unsigned short int y) { return (unsigned short int) std::pow((double)x,(double)y); }
255  static inline unsigned short int log(unsigned short int x) { return static_cast<unsigned short int> (std::log (static_cast<double> (x))); }
256  static inline unsigned short int log10(unsigned short int x) { return static_cast<unsigned short int> (std::log10 (static_cast<double> (x))); }
257 };
258 
259 
260 template<>
261 struct ScalarTraits<int>
262 {
263  typedef int magnitudeType;
264  typedef int halfPrecision;
265  typedef int doublePrecision;
266  static const bool isComplex = false;
267  static const bool isOrdinal = true;
268  static const bool isComparable = true;
269  static const bool hasMachineParameters = false;
270  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
271  static inline magnitudeType magnitude(int a) { return static_cast<int>(std::fabs(static_cast<double>(a))); }
272  static inline int zero() { return 0; }
273  static inline int one() { return 1; }
274  static inline int conjugate(int x) { return x; }
275  static inline int real(int x) { return x; }
276  static inline int imag(int) { return 0; }
277  static inline bool isnaninf(int) { return false; }
278  static inline void seedrandom(unsigned int s) {
279  std::srand(s);
280 #ifdef __APPLE__
281  // throw away first random number to address bug 3655
282  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
283  random();
284 #endif
285  }
286  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
287  static inline int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
288  static inline std::string name() { return "int"; }
289  static inline int squareroot(int x) { return (int) std::sqrt((double) x); }
290  static inline int pow(int x, int y) { return (int) std::pow((double)x,(double)y); }
291  static inline int log(int x) { return static_cast<int> (std::log (static_cast<double> (x))); }
292  static inline int log10(int x) { return static_cast<int> (std::log10 (static_cast<double> (x))); }
293 };
294 
295 
296 template<>
297 struct ScalarTraits<unsigned int>
298 {
299  typedef unsigned int magnitudeType;
300  typedef unsigned int halfPrecision;
301  typedef unsigned int doublePrecision;
302  static const bool isComplex = false;
303  static const bool isOrdinal = true;
304  static const bool isComparable = true;
305  static const bool hasMachineParameters = false;
306  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
307  static inline magnitudeType magnitude(unsigned int a) { return static_cast<unsigned int>(std::fabs(static_cast<double>(a))); }
308  static inline unsigned int zero() { return 0; }
309  static inline unsigned int one() { return 1; }
310  static inline unsigned int conjugate(unsigned int x) { return x; }
311  static inline unsigned int real(unsigned int x) { return x; }
312  static inline unsigned int imag(unsigned int) { return 0; }
313  static inline bool isnaninf(unsigned int) { return false; }
314  static inline void seedrandom(unsigned int s) {
315  std::srand(s);
316 #ifdef __APPLE__
317  // throw away first random number to address bug 3655
318  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
319  random();
320 #endif
321  }
322  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
323  static inline unsigned int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
324  static inline std::string name() { return "unsigned int"; }
325  static inline unsigned int squareroot(unsigned int x) { return (unsigned int) std::sqrt((double) x); }
326  static inline unsigned int pow(unsigned int x, unsigned int y) { return (unsigned int) std::pow((double)x,(double)y); }
327  static inline unsigned int log(unsigned int x) { return static_cast<unsigned int> (std::log (static_cast<double> (x))); }
328  static inline unsigned int log10(unsigned int x) { return static_cast<unsigned int> (std::log10 (static_cast<double> (x))); }
329 };
330 
331 
332 template<>
333 struct ScalarTraits<long int>
334 {
335  typedef long int magnitudeType;
336  typedef long int halfPrecision;
337  typedef long int doublePrecision;
338  static const bool isComplex = false;
339  static const bool isOrdinal = true;
340  static const bool isComparable = true;
341  static const bool hasMachineParameters = false;
342  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
343  static inline magnitudeType magnitude(long int a) { return static_cast<long int>(std::fabs(static_cast<double>(a))); }
344  static inline long int zero() { return 0; }
345  static inline long int one() { return 1; }
346  static inline long int conjugate(long int x) { return x; }
347  static inline long int real(long int x) { return x; }
348  static inline long int imag(long int) { return 0; }
349  static inline bool isnaninf(long int) { return false; }
350  static inline void seedrandom(unsigned int s) {
351  std::srand(s);
352 #ifdef __APPLE__
353  // throw away first random number to address bug 3655
354  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
355  random();
356 #endif
357  }
358  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
359  static inline long int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
360  static inline std::string name() { return "long int"; }
361  static inline long int squareroot(long int x) { return (long int) std::sqrt((double) x); }
362  static inline long int pow(long int x, long int y) { return (long int) std::pow((double)x,(double)y); }
363  // Note: Depending on the number of bits in long int, the cast from
364  // long int to double may not be exact.
365  static inline long int log(long int x) { return static_cast<long int> (std::log (static_cast<double> (x))); }
366  static inline long int log10(long int x) { return static_cast<long int> (std::log10 (static_cast<double> (x))); }
367 };
368 
369 
370 template<>
371 struct ScalarTraits<long unsigned int>
372 {
373  typedef long unsigned int magnitudeType;
374  typedef long unsigned int halfPrecision;
375  typedef long unsigned int doublePrecision;
376  static const bool isComplex = false;
377  static const bool isOrdinal = true;
378  static const bool isComparable = true;
379  static const bool hasMachineParameters = false;
380  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
381  static inline magnitudeType magnitude(long unsigned int a) { return static_cast<long unsigned int>(std::fabs(static_cast<double>(a))); }
382  static inline long unsigned int zero() { return 0; }
383  static inline long unsigned int one() { return 1; }
384  static inline long unsigned int conjugate(long unsigned int x) { return x; }
385  static inline long unsigned int real(long unsigned int x) { return x; }
386  static inline long unsigned int imag(long unsigned int) { return 0; }
387  static inline bool isnaninf(long unsigned int) { return false; }
388  static inline void seedrandom(unsigned int s) {
389  std::srand(s);
390 #ifdef __APPLE__
391  // throw away first random number to address bug 3655
392  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
393  random();
394 #endif
395  }
396  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
397  static inline long unsigned int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
398  static inline std::string name() { return "long unsigned int"; }
399  static inline long unsigned int squareroot(long unsigned int x) { return (long unsigned int) std::sqrt((double) x); }
400  static inline long unsigned int pow(long unsigned int x, long unsigned int y) { return (long unsigned int) std::pow((double)x,(double)y); }
401  // Note: Depending on the number of bits in long unsigned int, the
402  // cast from long unsigned int to double may not be exact.
403  static inline long unsigned int log(long unsigned int x) { return static_cast<long unsigned int> (std::log (static_cast<double> (x))); }
404  static inline long unsigned int log10(long unsigned int x) { return static_cast<long unsigned int> (std::log10 (static_cast<double> (x))); }
405 };
406 
407 
408 #ifdef HAVE_TEUCHOS_LONG_LONG_INT
409 template<>
410 struct ScalarTraits<long long int>
411 {
412  typedef long long int magnitudeType;
413  typedef long long int halfPrecision;
414  typedef long long int doublePrecision;
415  static const bool isComplex = false;
416  static const bool isOrdinal = true;
417  static const bool isComparable = true;
418  static const bool hasMachineParameters = false;
419  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
420  static inline magnitudeType magnitude(long long int a) { return static_cast<long long int>(std::fabs(static_cast<double>(a))); }
421  static inline long long int zero() { return 0; }
422  static inline long long int one() { return 1; }
423  static inline long long int conjugate(long long int x) { return x; }
424  static inline long long int real(long long int x) { return x; }
425  static inline long long int imag(long long int) { return 0; }
426  static inline bool isnaninf(long long int) { return false; }
427  static inline void seedrandom(unsigned int s) {
428  std::srand(s);
429 #ifdef __APPLE__
430  // throw away first random number to address bug 3655
431  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
432  random();
433 #endif
434  }
435  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
436  static inline long long int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
437  static inline std::string name() { return "long long int"; }
438  static inline long long int squareroot(long long int x) { return (long long int) std::sqrt((double) x); }
439  static inline long long int pow(long long int x, long long int y) { return (long long int) std::pow((double)x,(double)y); }
440  // Note: Depending on the number of bits in long long int, the cast
441  // from long long int to double may not be exact.
442  static inline long long int log(long long int x) { return static_cast<long long int> (std::log (static_cast<double> (x))); }
443  static inline long long int log10(long long int x) { return static_cast<long long int> (std::log10 (static_cast<double> (x))); }
444 };
445 
446 template<>
447 struct ScalarTraits<unsigned long long int>
448 {
449  typedef unsigned long long int magnitudeType;
450  typedef unsigned long long int halfPrecision;
451  typedef unsigned long long int doublePrecision;
452  static const bool isComplex = false;
453  static const bool isOrdinal = true;
454  static const bool isComparable = true;
455  static const bool hasMachineParameters = false;
456  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
457  static inline magnitudeType magnitude(unsigned long long int a) { return static_cast<unsigned long long int>(std::fabs(static_cast<double>(a))); }
458  static inline unsigned long long int zero() { return 0; }
459  static inline unsigned long long int one() { return 1; }
460  static inline unsigned long long int conjugate(unsigned long long int x) { return x; }
461  static inline unsigned long long int real(unsigned long long int x) { return x; }
462  static inline unsigned long long int imag(unsigned long long int) { return 0; }
463  static inline bool isnaninf(unsigned long long int) { return false; }
464  static inline void seedrandom(unsigned int s) {
465  std::srand(s);
466 #ifdef __APPLE__
467  // throw away first random number to address bug 3655
468  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
469  random();
470 #endif
471  }
472  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
473  static inline unsigned long long int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
474  static inline std::string name() { return "unsigned long long int"; }
475  static inline unsigned long long int squareroot(unsigned long long int x) { return (unsigned long long int) std::sqrt((double) x); }
476  static inline unsigned long long int pow(unsigned long long int x, unsigned long long int y) { return (unsigned long long int) std::pow((double)x,(double)y); }
477  // Note: Depending on the number of bits in unsigned long long int,
478  // the cast from unsigned long long int to double may not be exact.
479  static inline unsigned long long int log(unsigned long long int x) { return static_cast<unsigned long long int> (std::log (static_cast<double> (x))); }
480  static inline unsigned long long int log10(unsigned long long int x) { return static_cast<unsigned long long int> (std::log10 (static_cast<double> (x))); }
481 };
482 #endif // HAVE_TEUCHOS_LONG_LONG_INT
483 
484 
485 #ifdef HAVE_TEUCHOS___INT64
486 
487 template<>
488 struct ScalarTraits<__int64>
489 {
490  typedef __int64 magnitudeType;
491  typedef __int64 halfPrecision;
492  typedef __int64 doublePrecision;
493  static const bool isComplex = false;
494  static const bool isOrdinal = true;
495  static const bool isComparable = true;
496  static const bool hasMachineParameters = false;
497  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
498  static inline magnitudeType magnitude(__int64 a) { return static_cast<__int64>(std::fabs(static_cast<double>(a))); }
499  static inline __int64 zero() { return 0; }
500  static inline __int64 one() { return 1; }
501  static inline __int64 conjugate(__int64 x) { return x; }
502  static inline __int64 real(__int64 x) { return x; }
503  static inline __int64 imag(__int64) { return 0; }
504  static inline void seedrandom(unsigned int s) {
505  std::srand(s);
506 #ifdef __APPLE__
507  // throw away first random number to address bug 3655
508  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
509  random();
510 #endif
511  }
512  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
513  static inline __int64 random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
514  static inline std::string name() { return "__int64"; }
515  static inline __int64 squareroot(__int64 x) { return (__int64) std::sqrt((double) x); }
516  static inline __int64 pow(__int64 x, __int64 y) { return (__int64) std::pow((double)x,(double)y); }
517  // Note: Depending on the number of bits in __int64, the cast
518  // from __int64 to double may not be exact.
519  static inline __int64 log(__int64 x) { return static_cast<__int64> (std::log (static_cast<double> (x))); }
520  static inline __int64 log10(__int64 x) { return static_cast<__int64> (std::log10 (static_cast<double> (x))); }
521 };
522 
523 template<>
524 struct ScalarTraits<unsigned __int64>
525 {
526  typedef unsigned __int64 magnitudeType;
527  typedef unsigned __int64 halfPrecision;
528  typedef unsigned __int64 doublePrecision;
529  static const bool isComplex = false;
530  static const bool isOrdinal = true;
531  static const bool isComparable = true;
532  static const bool hasMachineParameters = false;
533  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
534  static inline magnitudeType magnitude(unsigned __int64 a) { return static_cast<unsigned __int64>(std::fabs(static_cast<double>(a))); }
535  static inline unsigned __int64 zero() { return 0; }
536  static inline unsigned __int64 one() { return 1; }
537  static inline unsigned __int64 conjugate(unsigned __int64 x) { return x; }
538  static inline unsigned __int64 real(unsigned __int64 x) { return x; }
539  static inline unsigned __int64 imag(unsigned __int64) { return 0; }
540  static inline void seedrandom(unsigned int s) {
541  std::srand(s);
542 #ifdef __APPLE__
543  // throw away first random number to address bug 3655
544  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
545  random();
546 #endif
547  }
548  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
549  static inline unsigned __int64 random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
550  static inline std::string name() { return "unsigned __int64"; }
551  static inline unsigned __int64 squareroot(unsigned __int64 x) { return (unsigned __int64) std::sqrt((double) x); }
552  static inline unsigned __int64 pow(unsigned __int64 x, unsigned __int64 y) { return (unsigned __int64) std::pow((double)x,(double)y); }
553  // Note: Depending on the number of bits in unsigned __int64,
554  // the cast from unsigned __int64 to double may not be exact.
555  static inline unsigned __int64 log(unsigned __int64 x) { return static_cast<unsigned __int64> (std::log (static_cast<double> (x))); }
556  static inline unsigned __int64 log10(unsigned __int64 x) { return static_cast<unsigned __int64> (std::log10 (static_cast<double> (x))); }
557 };
558 
559 #endif // HAVE_TEUCHOS___INT64
560 
561 
562 #ifndef __sun
563 extern TEUCHOSCORE_LIB_DLL_EXPORT const float flt_nan;
564 #endif
565 
566 
567 template<>
568 struct ScalarTraits<float>
569 {
570  typedef float magnitudeType;
571  typedef float halfPrecision; // should become IEEE754-2008 binary16 or fp16 later, perhaps specified at configure according to architectural support
572  typedef double doublePrecision;
573  static const bool isComplex = false;
574  static const bool isOrdinal = false;
575  static const bool isComparable = true;
576  static const bool hasMachineParameters = true;
577  static inline float eps() {
578  return std::numeric_limits<float>::epsilon();
579  }
580  static inline float sfmin() {
581  return std::numeric_limits<float>::min();
582  }
583  static inline float base() {
584  return static_cast<float>(std::numeric_limits<float>::radix);
585  }
586  static inline float prec() {
587  return eps()*base();
588  }
589  static inline float t() {
590  return static_cast<float>(std::numeric_limits<float>::digits);
591  }
592  static inline float rnd() {
593  return ( std::numeric_limits<float>::round_style == std::round_to_nearest ? one() : zero() );
594  }
595  static inline float emin() {
596  return static_cast<float>(std::numeric_limits<float>::min_exponent);
597  }
598  static inline float rmin() {
599  return std::numeric_limits<float>::min();
600  }
601  static inline float emax() {
602  return static_cast<float>(std::numeric_limits<float>::max_exponent);
603  }
604  static inline float rmax() {
605  return std::numeric_limits<float>::max();
606  }
607  static inline magnitudeType magnitude(float a)
608  {
609 #ifdef TEUCHOS_DEBUG
610  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
611  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
612 #endif
613  return std::fabs(a);
614  }
615  static inline float zero() { return(0.0f); }
616  static inline float one() { return(1.0f); }
617  static inline float conjugate(float x) { return(x); }
618  static inline float real(float x) { return x; }
619  static inline float imag(float) { return zero(); }
620  static inline float nan() {
621 #ifdef __sun
622  return 0.0f/std::sin(0.0f);
623 #else
624  return flt_nan;
625 #endif
626  }
627  static inline bool isnaninf(float x) {
628  return generic_real_isnaninf<float>(x);
629  }
630  static inline void seedrandom(unsigned int s) {
631  std::srand(s);
632 #ifdef __APPLE__
633  // throw away first random number to address bug 3655
634  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
635  random();
636 #endif
637  }
638  static inline float random() { float rnd = (float) std::rand() / RAND_MAX; return (-1.0f + 2.0f * rnd); }
639  static inline std::string name() { return "float"; }
640  static inline float squareroot(float x)
641  {
642 #ifdef TEUCHOS_DEBUG
643  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
644  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
645 #endif
646  errno = 0;
647  const float rtn = std::sqrt(x);
648  if (errno)
649  return nan();
650  return rtn;
651  }
652  static inline float pow(float x, float y) { return std::pow(x,y); }
653  static inline float log(float x) { return std::log(x); }
654  static inline float log10(float x) { return std::log10(x); }
655 };
656 
657 
658 #ifndef __sun
659 extern TEUCHOSCORE_LIB_DLL_EXPORT const double dbl_nan;
660 #endif
661 
662 
663 template<>
664 struct ScalarTraits<double>
665 {
666  typedef double magnitudeType;
667  typedef float halfPrecision;
668  /* there are different options as to how to double "double"
669  - QD's DD (if available)
670  - ARPREC
671  - GNU MP
672  - a true hardware quad
673 
674  in the shortterm, this should be specified at configure time. I have inserted a configure-time option (--enable-teuchos-double-to-dd)
675  which uses QD's DD when available. This must be used alongside --enable-teuchos-qd.
676  */
677 #if defined(HAVE_TEUCHOS_DOUBLE_TO_QD)
678  typedef dd_real doublePrecision;
679 #elif defined(HAVE_TEUCHOS_DOUBLE_TO_ARPREC)
680  typedef mp_real doublePrecision;
681 #else
682  typedef double doublePrecision; // don't double "double" in this case
683 #endif
684  static const bool isComplex = false;
685  static const bool isOrdinal = false;
686  static const bool isComparable = true;
687  static const bool hasMachineParameters = true;
688  static inline double eps() {
689  return std::numeric_limits<double>::epsilon();
690  }
691  static inline double sfmin() {
692  return std::numeric_limits<double>::min();
693  }
694  static inline double base() {
695  return std::numeric_limits<double>::radix;
696  }
697  static inline double prec() {
698  return eps()*base();
699  }
700  static inline double t() {
701  return std::numeric_limits<double>::digits;
702  }
703  static inline double rnd() {
704  return ( std::numeric_limits<double>::round_style == std::round_to_nearest ? double(1.0) : double(0.0) );
705  }
706  static inline double emin() {
707  return std::numeric_limits<double>::min_exponent;
708  }
709  static inline double rmin() {
710  return std::numeric_limits<double>::min();
711  }
712  static inline double emax() {
713  return std::numeric_limits<double>::max_exponent;
714  }
715  static inline double rmax() {
716  return std::numeric_limits<double>::max();
717  }
718  static inline magnitudeType magnitude(double a)
719  {
720 #ifdef TEUCHOS_DEBUG
721  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
722  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
723 #endif
724  return std::fabs(a);
725  }
726  static inline double zero() { return 0.0; }
727  static inline double one() { return 1.0; }
728  static inline double conjugate(double x) { return(x); }
729  static inline double real(double x) { return(x); }
730  static inline double imag(double) { return(0); }
731  static inline double nan() {
732 #ifdef __sun
733  return 0.0/std::sin(0.0);
734 #else
735  return dbl_nan;
736 #endif
737  }
738  static inline bool isnaninf(double x) {
739  return generic_real_isnaninf<double>(x);
740  }
741  static inline void seedrandom(unsigned int s) {
742  std::srand(s);
743 #ifdef __APPLE__
744  // throw away first random number to address bug 3655
745  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
746  random();
747 #endif
748  }
749  static inline double random() { double rnd = (double) std::rand() / RAND_MAX; return (double)(-1.0 + 2.0 * rnd); }
750  static inline std::string name() { return "double"; }
751  static inline double squareroot(double x)
752  {
753 #ifdef TEUCHOS_DEBUG
754  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
755  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
756 #endif
757  errno = 0;
758  const double rtn = std::sqrt(x);
759  if (errno)
760  return nan();
761  return rtn;
762  }
763  static inline double pow(double x, double y) { return std::pow(x,y); }
764  static inline double log(double x) { return std::log(x); }
765  static inline double log10(double x) { return std::log10(x); }
766 };
767 
768 
769 #ifdef HAVE_TEUCHOSCORE_QUADMATH
770 
771 template<>
772 struct ScalarTraits<__float128> {
773  typedef __float128 magnitudeType;
774  // Unfortunately, we can't rely on a standard __float256 type. It
775  // might make sense to use qd_real, but mixing quadmath and QD might
776  // cause unforeseen issues.
777  typedef __float128 doublePrecision;
778  typedef double halfPrecision;
779 
780  static const bool isComplex = false;
781  static const bool isOrdinal = false;
782  static const bool isComparable = true;
783  static const bool hasMachineParameters = true;
784 
785  static __float128 eps () {
786  return FLT128_EPSILON;
787  }
788  static __float128 sfmin () {
789  return FLT128_MIN; // ???
790  }
791  static __float128 base () {
792  return 2;
793  }
794  static __float128 prec () {
795  return eps () * base ();
796  }
797  static __float128 t () {
798  return FLT128_MANT_DIG;
799  }
800  static __float128 rnd () {
801  return 1.0;
802  }
803  static __float128 emin () {
804  return FLT128_MIN_EXP;
805  }
806  static __float128 rmin () {
807  return FLT128_MIN; // ??? // should be base^(emin-1)
808  }
809  static __float128 emax () {
810  return FLT128_MAX_EXP;
811  }
812  static __float128 rmax () {
813  return FLT128_MAX; // ??? // should be (base^emax)*(1-eps)
814  }
815  static magnitudeType magnitude (const __float128& x) {
816  return fabsq (x);
817  }
818  static __float128 zero () {
819  return 0.0;
820  }
821  static __float128 one () {
822  return 1.0;
823  }
824  static __float128 conjugate (const __float128& x) {
825  return x;
826  }
827  static __float128 real (const __float128& x) {
828  return x;
829  }
830  static __float128 imag (const __float128& /* x */) {
831  return 0.0;
832  }
833  static __float128 nan () {
834  return strtoflt128 ("NAN()", NULL); // ???
835  }
836  static bool isnaninf (const __float128& x) {
837  return isinfq (x) || isnanq (x);
838  }
839  static inline void seedrandom (unsigned int s) {
840  std::srand (s);
841 #ifdef __APPLE__
842  // throw away first random number to address bug 3655
843  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
844  random ();
845 #endif
846  }
847  static __float128 random () {
848  // Half the smallest normalized double, is the scaling factor of
849  // the lower-order term in the double-double representation.
850  const __float128 scalingFactor =
851  static_cast<__float128> (std::numeric_limits<double>::min ()) /
852  static_cast<__float128> (2.0);
853  const __float128 higherOrderTerm =
854  static_cast<__float128> (ScalarTraits<double>::random ());
855  const __float128 lowerOrderTerm =
856  static_cast<__float128> (ScalarTraits<double>::random ()) *
857  scalingFactor;
858  return higherOrderTerm + lowerOrderTerm;
859  }
860  static std::string name () {
861  return "__float128";
862  }
863  static __float128 squareroot (const __float128& x) {
864  return sqrtq (x);
865  }
866  static __float128 pow (const __float128& x, const __float128& y) {
867  return powq (x, y);
868  }
869  static __float128 log (const __float128& x) {
870  return logq (x);
871  }
872  static __float128 log10 (const __float128& x) {
873  return log10q (x);
874  }
875 };
876 #endif // HAVE_TEUCHOSCORE_QUADMATH
877 
878 
879 
880 #ifdef HAVE_TEUCHOS_QD
881 
882 bool operator&&(const dd_real &a, const dd_real &b);
883 bool operator&&(const qd_real &a, const qd_real &b);
884 
885 template<>
886 struct ScalarTraits<dd_real>
887 {
888  typedef dd_real magnitudeType;
889  typedef double halfPrecision;
890  typedef qd_real doublePrecision;
891  static const bool isComplex = false;
892  static const bool isOrdinal = false;
893  static const bool isComparable = true;
894  static const bool hasMachineParameters = true;
895  static inline dd_real eps() { return std::numeric_limits<dd_real>::epsilon(); }
896  static inline dd_real sfmin() { return std::numeric_limits<dd_real>::min(); }
897  static inline dd_real base() { return std::numeric_limits<dd_real>::radix; }
898  static inline dd_real prec() { return eps()*base(); }
899  static inline dd_real t() { return std::numeric_limits<dd_real>::digits; }
900  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) ); }
901  static inline dd_real emin() { return std::numeric_limits<dd_real>::min_exponent; }
902  static inline dd_real rmin() { return std::numeric_limits<dd_real>::min(); }
903  static inline dd_real emax() { return std::numeric_limits<dd_real>::max_exponent; }
904  static inline dd_real rmax() { return std::numeric_limits<dd_real>::max(); }
905  static inline magnitudeType magnitude(dd_real a)
906  {
907 #ifdef TEUCHOS_DEBUG
908  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
909  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
910 #endif
911  return ::abs(a);
912  }
913  static inline dd_real zero() { return dd_real(0.0); }
914  static inline dd_real one() { return dd_real(1.0); }
915  static inline dd_real conjugate(dd_real x) { return(x); }
916  static inline dd_real real(dd_real x) { return x ; }
917  static inline dd_real imag(dd_real) { return zero(); }
918  static inline dd_real nan() { return dd_real::_nan; }
919  static inline bool isnaninf(dd_real x) { return isnan(x) || isinf(x); }
920  static inline void seedrandom(unsigned int s) {
921  // ddrand() uses std::rand(), so the std::srand() is our seed
922  std::srand(s);
923 #ifdef __APPLE__
924  // throw away first random number to address bug 3655
925  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
926  random();
927 #endif
928  }
929  static inline dd_real random() { return ddrand(); }
930  static inline std::string name() { return "dd_real"; }
931  static inline dd_real squareroot(dd_real x)
932  {
933 #ifdef TEUCHOS_DEBUG
934  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
935  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
936 #endif
937  return ::sqrt(x);
938  }
939  static inline dd_real pow(dd_real x, dd_real y) { return ::pow(x,y); }
940  // dd_real puts its transcendental functions in the global namespace.
941  static inline dd_real log(dd_real x) { return ::log(x); }
942  static inline dd_real log10(dd_real x) { return ::log10(x); }
943 };
944 
945 
946 template<>
947 struct ScalarTraits<qd_real>
948 {
949  typedef qd_real magnitudeType;
950  typedef dd_real halfPrecision;
951  typedef qd_real doublePrecision;
952  static const bool isComplex = false;
953  static const bool isOrdinal = false;
954  static const bool isComparable = true;
955  static const bool hasMachineParameters = true;
956  static inline qd_real eps() { return std::numeric_limits<qd_real>::epsilon(); }
957  static inline qd_real sfmin() { return std::numeric_limits<qd_real>::min(); }
958  static inline qd_real base() { return std::numeric_limits<qd_real>::radix; }
959  static inline qd_real prec() { return eps()*base(); }
960  static inline qd_real t() { return std::numeric_limits<qd_real>::digits; }
961  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) ); }
962  static inline qd_real emin() { return std::numeric_limits<qd_real>::min_exponent; }
963  static inline qd_real rmin() { return std::numeric_limits<qd_real>::min(); }
964  static inline qd_real emax() { return std::numeric_limits<qd_real>::max_exponent; }
965  static inline qd_real rmax() { return std::numeric_limits<qd_real>::max(); }
966  static inline magnitudeType magnitude(qd_real a)
967  {
968 #ifdef TEUCHOS_DEBUG
969  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
970  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
971 #endif
972  return ::abs(a);
973  }
974  static inline qd_real zero() { return qd_real(0.0); }
975  static inline qd_real one() { return qd_real(1.0); }
976  static inline qd_real conjugate(qd_real x) { return(x); }
977  static inline qd_real real(qd_real x) { return x ; }
978  static inline qd_real imag(qd_real) { return zero(); }
979  static inline qd_real nan() { return qd_real::_nan; }
980  static inline bool isnaninf(qd_real x) { return isnan(x) || isinf(x); }
981  static inline void seedrandom(unsigned int s) {
982  // qdrand() uses std::rand(), so the std::srand() is our seed
983  std::srand(s);
984 #ifdef __APPLE__
985  // throw away first random number to address bug 3655
986  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
987  random();
988 #endif
989  }
990  static inline qd_real random() { return qdrand(); }
991  static inline std::string name() { return "qd_real"; }
992  static inline qd_real squareroot(qd_real x)
993  {
994 #ifdef TEUCHOS_DEBUG
995  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
996  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
997 #endif
998  return ::sqrt(x);
999  }
1000  static inline qd_real pow(qd_real x, qd_real y) { return ::pow(x,y); }
1001  // qd_real puts its transcendental functions in the global namespace.
1002  static inline qd_real log(qd_real x) { return ::log(x); }
1003  static inline qd_real log10(qd_real x) { return ::log10(x); }
1004 };
1005 
1006 
1007 #endif // HAVE_TEUCHOS_QD
1008 
1009 
1010 #ifdef HAVE_TEUCHOS_GNU_MP
1011 
1012 
1013 extern gmp_randclass gmp_rng;
1014 
1015 
1016 /* Regarding halfPrecision, doublePrecision and mpf_class:
1017  Because the precision of an mpf_class float is not determined by the data type,
1018  there is no way to fill the typedefs for this object.
1019 
1020  Instead, we could create new data classes (e.g., Teuchos::MPF128, Teuchos::MPF256) for
1021  commonly used levels of precision, and fill out ScalarTraits for these. This would allow us
1022  to typedef the promotions and demotions in the appropriate way. These classes would serve to
1023  wrap an mpf_class object, calling the constructor for the appropriate precision, exposing the
1024  arithmetic routines but hiding the precision-altering routines.
1025 
1026  Alternatively (perhaps, preferably), would could create a single class templated on the precision (e.g., Teuchos::MPF<N>).
1027  Then we have a single (partially-specialized) implementation of ScalarTraits. This class, as above, must expose all of the
1028  operations expected of a scalar type; however, most of these can be trivially stolen from the gmpcxx.h header file
1029 
1030  CGB/RAB, 01/05/2009
1031 */
1032 template<>
1033 struct ScalarTraits<mpf_class>
1034 {
1035  typedef mpf_class magnitudeType;
1036  typedef mpf_class halfPrecision;
1037  typedef mpf_class doublePrecision;
1038  static const bool isComplex = false;
1039  static const bool hasMachineParameters = false;
1040  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
1041  static magnitudeType magnitude(mpf_class a) { return std::abs(a); }
1042  static inline mpf_class zero() { mpf_class zero = 0.0; return zero; }
1043  static inline mpf_class one() { mpf_class one = 1.0; return one; }
1044  static inline mpf_class conjugate(mpf_class x) { return x; }
1045  static inline mpf_class real(mpf_class x) { return(x); }
1046  static inline mpf_class imag(mpf_class x) { return(0); }
1047  static inline bool isnaninf(mpf_class x) { return false; } // mpf_class currently can't handle nan or inf!
1048  static inline void seedrandom(unsigned int s) {
1049  unsigned long int seedVal = static_cast<unsigned long int>(s);
1050  gmp_rng.seed( seedVal );
1051  }
1052  static inline mpf_class random() {
1053  return gmp_rng.get_f();
1054  }
1055  static inline std::string name() { return "mpf_class"; }
1056  static inline mpf_class squareroot(mpf_class x) { return std::sqrt(x); }
1057  static inline mpf_class pow(mpf_class x, mpf_class y) { return pow(x,y); }
1058  // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
1059 };
1060 
1061 #endif // HAVE_TEUCHOS_GNU_MP
1062 
1063 #ifdef HAVE_TEUCHOS_ARPREC
1064 
1065 /* See discussion above for mpf_class, regarding halfPrecision and doublePrecision. Something similar will need to be done
1066  for ARPREC. */
1067 template<>
1068 struct ScalarTraits<mp_real>
1069 {
1070  typedef mp_real magnitudeType;
1071  typedef double halfPrecision;
1072  typedef mp_real doublePrecision;
1073  static const bool isComplex = false;
1074  static const bool isComparable = true;
1075  static const bool isOrdinal = false;
1076  static const bool hasMachineParameters = false;
1077  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
1078  static magnitudeType magnitude(mp_real a) { return abs(a); }
1079  static inline mp_real zero() { mp_real zero = 0.0; return zero; }
1080  static inline mp_real one() { mp_real one = 1.0; return one; }
1081  static inline mp_real conjugate(mp_real x) { return x; }
1082  static inline mp_real real(mp_real x) { return(x); }
1083  static inline mp_real imag(mp_real x) { return zero(); }
1084  static inline bool isnaninf(mp_real x) { return false; } // ToDo: Change this?
1085  static inline void seedrandom(unsigned int s) {
1086  long int seedVal = static_cast<long int>(s);
1087  srand48(seedVal);
1088  }
1089  static inline mp_real random() { return mp_rand(); }
1090  static inline std::string name() { return "mp_real"; }
1091  static inline mp_real squareroot(mp_real x) { return sqrt(x); }
1092  static inline mp_real pow(mp_real x, mp_real y) { return pow(x,y); }
1093  // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
1094 };
1095 
1096 
1097 #endif // HAVE_TEUCHOS_ARPREC
1098 
1099 
1100 #ifdef HAVE_TEUCHOS_COMPLEX
1101 
1102 
1103 // Partial specialization for std::complex numbers templated on real type T
1104 template<class T>
1105 struct ScalarTraits<
1106  std::complex<T>
1107 >
1108 {
1109  typedef std::complex<T> ComplexT;
1110  typedef std::complex<typename ScalarTraits<T>::halfPrecision> halfPrecision;
1111  typedef std::complex<typename ScalarTraits<T>::doublePrecision> doublePrecision;
1112  typedef typename ScalarTraits<T>::magnitudeType magnitudeType;
1113  static const bool isComplex = true;
1114  static const bool isOrdinal = ScalarTraits<T>::isOrdinal;
1115  static const bool isComparable = false;
1116  static const bool hasMachineParameters = true;
1117  static inline magnitudeType eps() { return ScalarTraits<magnitudeType>::eps(); }
1118  static inline magnitudeType sfmin() { return ScalarTraits<magnitudeType>::sfmin(); }
1119  static inline magnitudeType base() { return ScalarTraits<magnitudeType>::base(); }
1120  static inline magnitudeType prec() { return ScalarTraits<magnitudeType>::prec(); }
1121  static inline magnitudeType t() { return ScalarTraits<magnitudeType>::t(); }
1122  static inline magnitudeType rnd() { return ScalarTraits<magnitudeType>::rnd(); }
1123  static inline magnitudeType emin() { return ScalarTraits<magnitudeType>::emin(); }
1124  static inline magnitudeType rmin() { return ScalarTraits<magnitudeType>::rmin(); }
1125  static inline magnitudeType emax() { return ScalarTraits<magnitudeType>::emax(); }
1126  static inline magnitudeType rmax() { return ScalarTraits<magnitudeType>::rmax(); }
1127  static magnitudeType magnitude(ComplexT a)
1128  {
1129 #ifdef TEUCHOS_DEBUG
1130  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
1131  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
1132 #endif
1133  return std::abs(a);
1134  }
1135  static inline ComplexT zero() { return ComplexT(ScalarTraits<magnitudeType>::zero(),ScalarTraits<magnitudeType>::zero()); }
1136  static inline ComplexT one() { return ComplexT(ScalarTraits<magnitudeType>::one(),ScalarTraits<magnitudeType>::zero()); }
1137  static inline ComplexT conjugate(ComplexT a){ return ComplexT(a.real(),-a.imag()); }
1138  static inline magnitudeType real(ComplexT a) { return a.real(); }
1139  static inline magnitudeType imag(ComplexT a) { return a.imag(); }
1140  static inline ComplexT nan() { return ComplexT(ScalarTraits<magnitudeType>::nan(),ScalarTraits<magnitudeType>::nan()); }
1141  static inline bool isnaninf(ComplexT x) { return ScalarTraits<magnitudeType>::isnaninf(x.real()) || ScalarTraits<magnitudeType>::isnaninf(x.imag()); }
1142  static inline void seedrandom(unsigned int s) { ScalarTraits<magnitudeType>::seedrandom(s); }
1143  static inline ComplexT random()
1144  {
1145  const T rnd1 = ScalarTraits<magnitudeType>::random();
1146  const T rnd2 = ScalarTraits<magnitudeType>::random();
1147  return ComplexT(rnd1,rnd2);
1148  }
1149  static inline std::string name() { return std::string("std::complex<")+std::string(ScalarTraits<magnitudeType>::name())+std::string(">"); }
1150  // This will only return one of the square roots of x, the other can be obtained by taking its conjugate
1151  static inline ComplexT squareroot(ComplexT x)
1152  {
1153 #ifdef TEUCHOS_DEBUG
1154  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
1155  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
1156 #endif
1157  typedef ScalarTraits<magnitudeType> STMT;
1158  const T r = x.real(), i = x.imag(), zero = STMT::zero(), two = 2.0;
1159  const T a = STMT::squareroot((r*r)+(i*i));
1160  const T nr = STMT::squareroot((a+r)/two);
1161  const T ni = ( i == zero ? zero : STMT::squareroot((a-r)/two) );
1162  return ComplexT(nr,ni);
1163  }
1164  // 2010/03/19: rabartl: Above, I had to add the check for i == zero
1165  // to avoid a returned NaN on the Intel 10.1 compiler. For some
1166  // reason, having these two squareroot calls in a row produce a NaN
1167  // in an optimized build with this compiler. Amazingly, when I put
1168  // in print statements (i.e. std::cout << ...) the NaN went away and
1169  // the second squareroot((a-r)/two) returned zero correctly. I
1170  // guess that calling the output routine flushed the registers or
1171  // something and restarted the squareroot rountine for this compiler
1172  // or something similar. Actually, due to roundoff, it is possible that a-r
1173  // might be slightly less than zero (i.e. -1e-16) and that would cause
1174  // a possbile NaN return. THe above if test is the right thing to do
1175  // I think and is very cheap.
1176  static inline ComplexT pow(ComplexT x, ComplexT y) { return pow(x,y); }
1177 };
1178 
1179 #endif // HAVE_TEUCHOS_COMPLEX
1180 #endif // DOXYGEN_SHOULD_SKIP_THIS
1181 
1182 } // Teuchos namespace
1183 
1184 #endif // _TEUCHOS_SCALARTRAITS_HPP_
Teuchos_ConfigDefs.hpp
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos_ScalarTraitsDecl.hpp
Declaration and default implementation for basic traits for the scalar field type.
Teuchos
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
Generated by   doxygen 1.8.13