Simbody  3.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
NTraits.h
Go to the documentation of this file.
1 #ifndef SimTK_SIMMATRIX_NTRAITS_H_
2 #define SimTK_SIMMATRIX_NTRAITS_H_
3 
4 /* -------------------------------------------------------------------------- *
5  * Simbody(tm): SimTKcommon *
6  * -------------------------------------------------------------------------- *
7  * This is part of the SimTK biosimulation toolkit originating from *
8  * Simbios, the NIH National Center for Physics-Based Simulation of *
9  * Biological Structures at Stanford, funded under the NIH Roadmap for *
10  * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. *
11  * *
12  * Portions copyright (c) 2005-12 Stanford University and the Authors. *
13  * Authors: Michael Sherman *
14  * Contributors: *
15  * *
16  * Licensed under the Apache License, Version 2.0 (the "License"); you may *
17  * not use this file except in compliance with the License. You may obtain a *
18  * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
19  * *
20  * Unless required by applicable law or agreed to in writing, software *
21  * distributed under the License is distributed on an "AS IS" BASIS, *
22  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
23  * See the License for the specific language governing permissions and *
24  * limitations under the License. *
25  * -------------------------------------------------------------------------- */
26 
50 #include "SimTKcommon/Constants.h"
51 #include "SimTKcommon/internal/CompositeNumericalTypes.h"
52 
53 #include <cstddef>
54 #include <cassert>
55 #include <complex>
56 #include <iostream>
57 #include <limits>
58 
59 using std::complex;
60 
61 namespace SimTK {
62 
63 // This is the 3rd type of number, conjugate. It is like complex except that
64 // the represented value is the conjugate of the value represented by a complex
65 // number containing the same bit pattern. That is, complex and conjugate both
66 // contain two real numbers, re and im, with complex(re,im) meaning
67 // (re + im*i) while conjugate(re,im) means (re - im*i). It is guaranteed that
68 // our conjugate type and complex type have identical sizes and representations.
69 // Together, these defininitions and guarantees permit conjugation
70 // to be done by reinterpretation rather than be computation.
71 template <class R> class conjugate; // Only defined for float, double, long double
72 
73 // Specializations of this class provide information about Composite Numerical
74 // Types in the style of std::numeric_limits<T>. It is specialized for the
75 // numeric types but can be invoked on any composite numerical type as well.
76 template <class T> class CNT;
77 
78 // NTraits provides information like CNT, but for numeric types only.
79 template <class N> class NTraits; // allowed only for N=<number>
80 template <class R> class NTraits< complex<R> >;
81 template <class R> class NTraits< conjugate<R> >;
82 template <> class NTraits<float>;
83 template <> class NTraits<double>;
84 template <> class NTraits<long double>;
85 
86 // This is an adaptor for numeric types which negates the apparent values. A
87 // negator<N> has exactly the same internal representation as a numeric value N,
88 // but it is to be interpreted has having the negative of the value it would
89 // have if interpreted as an N. This permits negation to be done by
90 // reinterpretation rather than compuation. A full set of arithmetic operators
91 // are provided involving negator<N>'s and N's. Sometimes we can save an op or
92 // two this way. For example negator<N>*negator<N> can be performed as an N*N
93 // since the negations cancel, and we saved two floating point negations.
94 template <class N> class negator; // Only defined for numbers
95 
96 // This is here so we can provide references to 0's when needed, e.g.
97 // when returning the imaginary part of a real number. These are local to
98 // the compilation unit, so the returned addresses will differ in different
99 // files. There are enough zeroes here for the widest number,
100 // complex<long double> (or conjugate<long double>).
101 static const complex<long double> zeroes(0);
102 
110 template <class R1, class R2> struct Widest {/* Only defined for built-ins. */};
111 template <> struct Widest<float,float> {typedef float Type; typedef float Precision;};
112 template <> struct Widest<float,double> {typedef double Type; typedef double Precision;};
113 template <> struct Widest<float,long double> {typedef long double Type; typedef long double Precision;};
114 template <> struct Widest<double,float> {typedef double Type; typedef double Precision;};
115 template <> struct Widest<double,double> {typedef double Type; typedef double Precision;};
116 template <> struct Widest<double,long double> {typedef long double Type; typedef long double Precision;};
117 template <> struct Widest<long double,float> {typedef long double Type; typedef long double Precision;};
118 template <> struct Widest<long double,double> {typedef long double Type; typedef long double Precision;};
119 template <> struct Widest<long double,long double> {typedef long double Type; typedef long double Precision;};
120 template <class R1, class R2> struct Widest< complex<R1>,complex<R2> > {
121  typedef complex< typename Widest<R1,R2>::Type > Type;
122  typedef typename Widest<R1,R2>::Precision Precision;
123 };
124 template <class R1, class R2> struct Widest< complex<R1>,R2 > {
125  typedef complex< typename Widest<R1,R2>::Type > Type;
126  typedef typename Widest<R1,R2>::Precision Precision;
127 };
128 template <class R1, class R2> struct Widest< R1,complex<R2> > {
129  typedef complex< typename Widest<R1,R2>::Type > Type;
130  typedef typename Widest<R1,R2>::Precision Precision;
131 };
132 
142 template <class R1, class R2> struct Narrowest {/* Only defined for built-ins. */};
143 template <> struct Narrowest<float,float> {typedef float Type; typedef float Precision;};
144 template <> struct Narrowest<float,double> {typedef float Type; typedef float Precision;};
145 template <> struct Narrowest<float,long double> {typedef float Type; typedef float Precision;};
146 template <> struct Narrowest<double,float> {typedef float Type; typedef float Precision;};
147 template <> struct Narrowest<double,double> {typedef double Type; typedef double Precision;};
148 template <> struct Narrowest<double,long double> {typedef double Type; typedef double Precision;};
149 template <> struct Narrowest<long double,float> {typedef float Type; typedef float Precision;};
150 template <> struct Narrowest<long double,double> {typedef double Type; typedef double Precision;};
151 template <> struct Narrowest<long double,long double> {typedef long double Type; typedef long double Precision;};
152 template <class R1, class R2> struct Narrowest< complex<R1>,complex<R2> > {
153  typedef complex< typename Narrowest<R1,R2>::Type > Type;
154  typedef typename Narrowest<R1,R2>::Precision Precision;
155 };
156 template <class R1, class R2> struct Narrowest< complex<R1>,R2 > {
157  typedef complex< typename Narrowest<R1,R2>::Type > Type;
158  typedef typename Narrowest<R1,R2>::Precision Precision;
159 };
160 template <class R1, class R2> struct Narrowest< R1,complex<R2> > {
161  typedef complex< typename Narrowest<R1,R2>::Type > Type;
162  typedef typename Narrowest<R1,R2>::Precision Precision;
163 };
164 
166 template <class R> class RTraits {/* Only defined for real types */};
167 template <> class RTraits<float> {
168 public:
170  static const float& getEps() {static const float c=std::numeric_limits<float>::epsilon(); return c;}
172  static const float& getSignificant() {static const float c=std::pow(getEps(), 0.875f); return c;}
174  static double getDefaultTolerance() {return (double)getSignificant();}
175 };
176 template <> class RTraits<double> {
177 public:
178  static const double& getEps() {static const double c=std::numeric_limits<double>::epsilon(); return c;}
179  static const double& getSignificant() {static const double c=std::pow(getEps(), 0.875); return c;}
180  static double getDefaultTolerance() {return getSignificant();}
181 };
182 template <> class RTraits<long double> {
183 public:
184  static const long double& getEps() {static const long double c=std::numeric_limits<long double>::epsilon(); return c;}
185  static const long double& getSignificant() {static const long double c=std::pow(getEps(), 0.875L); return c;}
186  static double getDefaultTolerance() {return (double)getSignificant();}
187 };
188 
205 // See negator.h for isNaN() applied to negated scalars.
207 inline bool isNaN(const float& x) {return std::isnan(x);}
208 inline bool isNaN(const double& x) {return std::isnan(x);}
209 inline bool isNaN(const long double& x) {return std::isnan(x);}
210 
211 template <class P> inline bool
212 isNaN(const std::complex<P>& x)
213 { return isNaN(x.real()) || isNaN(x.imag());}
214 
215 template <class P> inline bool
216 isNaN(const conjugate<P>& x)
217 { return isNaN(x.real()) || isNaN(x.negImag());}
219 
232 // See negator.h for isFinite() applied to negated scalars.
234 inline bool isFinite(const float& x) {return std::isfinite(x);}
235 inline bool isFinite(const double& x) {return std::isfinite(x);}
236 inline bool isFinite(const long double& x) {return std::isfinite(x);}
237 
238 template <class P> inline bool
239 isFinite(const std::complex<P>& x)
240 { return isFinite(x.real()) && isFinite(x.imag());}
241 
242 template <class P> inline bool
243 isFinite(const conjugate<P>& x)
244 { return isFinite(x.real()) && isFinite(x.negImag());}
246 
261 // See negator.h for isInf() applied to negated scalars.
263 inline bool isInf(const float& x) {return std::isinf(x);}
264 inline bool isInf(const double& x) {return std::isinf(x);}
265 inline bool isInf(const long double& x) {return std::isinf(x);}
266 
267 template <class P> inline bool
268 isInf(const std::complex<P>& x) {
269  return (isInf(x.real()) && !isNaN(x.imag()))
270  || (isInf(x.imag()) && !isNaN(x.real()));
271 }
272 
273 template <class P> inline bool
274 isInf(const conjugate<P>& x) {
275  return (isInf(x.real()) && !isNaN(x.negImag()))
276  || (isInf(x.negImag()) && !isNaN(x.real()));
277 }
279 
319 
320 inline bool isNumericallyEqual(const float& a, const float& b,
321  double tol = RTraits<float>::getDefaultTolerance())
322 { if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false;
323  const float scale = std::max(std::max(std::abs(a),std::abs(b)), 1.f);
324  return std::abs(a-b) <= scale*(float)tol; }
326 inline bool isNumericallyEqual(const double& a, const double& b,
327  double tol = RTraits<double>::getDefaultTolerance())
328 { if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false;
329  const double scale = std::max(std::max(std::abs(a),std::abs(b)), 1.);
330  return std::abs(a-b) <= scale*tol; }
332 inline bool isNumericallyEqual(const long double& a, const long double& b,
333  double tol = RTraits<long double>::getDefaultTolerance())
334 { if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false;
335  const long double scale = std::max(std::max(std::abs(a),std::abs(b)), 1.L);
336  return std::abs(a-b) <= scale*(long double)tol; }
337 
339 inline bool isNumericallyEqual(const float& a, const double& b,
340  double tol = RTraits<float>::getDefaultTolerance())
341 { return isNumericallyEqual((double)a, b, tol); }
343 inline bool isNumericallyEqual(const double& a, const float& b,
344  double tol = RTraits<float>::getDefaultTolerance())
345 { return isNumericallyEqual(a, (double)b, tol); }
347 inline bool isNumericallyEqual(const float& a, const long double& b,
348  double tol = RTraits<float>::getDefaultTolerance())
349 { return isNumericallyEqual((long double)a, b, tol); }
351 inline bool isNumericallyEqual(const long double& a, const float& b,
352  double tol = RTraits<float>::getDefaultTolerance())
353 { return isNumericallyEqual(a, (long double)b, tol); }
355 inline bool isNumericallyEqual(const double& a, const long double& b,
356  double tol = RTraits<double>::getDefaultTolerance())
357 { return isNumericallyEqual((long double)a, b, tol); }
359 inline bool isNumericallyEqual(const long double& a, const double& b,
360  double tol = RTraits<double>::getDefaultTolerance())
361 { return isNumericallyEqual(a, (long double)b, tol); }
362 
364 inline bool isNumericallyEqual(const float& a, int b,
365  double tol = RTraits<float>::getDefaultTolerance())
366 { return isNumericallyEqual(a, (double)b, tol); }
368 inline bool isNumericallyEqual(int a, const float& b,
369  double tol = RTraits<float>::getDefaultTolerance())
370 { return isNumericallyEqual((double)a, b, tol); }
372 inline bool isNumericallyEqual(const double& a, int b,
373  double tol = RTraits<double>::getDefaultTolerance())
374 { return isNumericallyEqual(a, (double)b, tol); }
376 inline bool isNumericallyEqual(int a, const double& b,
377  double tol = RTraits<double>::getDefaultTolerance())
378 { return isNumericallyEqual((double)a, b, tol); }
380 inline bool isNumericallyEqual(const long double& a, int b,
381  double tol = RTraits<long double>::getDefaultTolerance())
382 { return isNumericallyEqual(a, (long double)b, tol); }
384 inline bool isNumericallyEqual(int a, const long double& b,
385  double tol = RTraits<long double>::getDefaultTolerance())
386 { return isNumericallyEqual((long double)a, b, tol); }
387 
391 template <class P, class Q>
392 inline bool isNumericallyEqual
393  ( const std::complex<P>& a, const std::complex<Q>& b,
394  double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance())
395 { return isNumericallyEqual(a.real(),b.real(),tol)
396  && isNumericallyEqual(a.imag(),b.imag(),tol); }
397 
401 template <class P, class Q>
402 inline bool isNumericallyEqual
403  ( const conjugate<P>& a, const conjugate<Q>& b,
404  double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance())
405 { return isNumericallyEqual(a.real(),b.real(),tol)
406  && isNumericallyEqual(a.imag(),b.imag(),tol); }
407 
411 template <class P, class Q>
412 inline bool isNumericallyEqual
413  ( const std::complex<P>& a, const conjugate<Q>& b,
414  double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance())
415 { return isNumericallyEqual(a.real(),b.real(),tol)
416  && isNumericallyEqual(a.imag(),b.imag(),tol); }
417 
421 template <class P, class Q>
422 inline bool isNumericallyEqual
423  ( const conjugate<P>& a, const std::complex<Q>& b,
424  double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance())
425 { return isNumericallyEqual(a.real(),b.real(),tol)
426  && isNumericallyEqual(a.imag(),b.imag(),tol); }
427 
429 template <class P> inline bool
430 isNumericallyEqual(const std::complex<P>& a, const float& b,
431  double tol = RTraits<float>::getDefaultTolerance())
432 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.f,tol); }
434 template <class P> inline bool
435 isNumericallyEqual(const float& a, const std::complex<P>& b,
436  double tol = RTraits<float>::getDefaultTolerance())
437 { return isNumericallyEqual(b,a,tol); }
439 template <class P> inline bool
440 isNumericallyEqual(const std::complex<P>& a, const double& b,
441  double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance())
442 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.,tol); }
444 template <class P> inline bool
445 isNumericallyEqual(const double& a, const std::complex<P>& b,
446  double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance())
447 { return isNumericallyEqual(b,a,tol); }
449 template <class P> inline bool
450 isNumericallyEqual(const std::complex<P>& a, const long double& b,
451  double tol = RTraits<P>::getDefaultTolerance())
452 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.L,tol); }
454 template <class P> inline bool
455 isNumericallyEqual(const long double& a, const std::complex<P>& b,
456  double tol = RTraits<P>::getDefaultTolerance())
457 { return isNumericallyEqual(b,a,tol); }
459 template <class P> inline bool
460 isNumericallyEqual(const std::complex<P>& a, int b,
461  double tol = RTraits<P>::getDefaultTolerance())
462 { typedef typename Widest<P,double>::Precision W; return isNumericallyEqual(a,(W)b,tol); }
464 template <class P> inline bool
465 isNumericallyEqual(int a, const std::complex<P>& b,
466  double tol = RTraits<P>::getDefaultTolerance())
467 { return isNumericallyEqual(b,a,tol); }
468 
470 template <class P> inline bool
471 isNumericallyEqual(const conjugate<P>& a, const float& b,
472  double tol = RTraits<float>::getDefaultTolerance())
473 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.f,tol); }
475 template <class P> inline bool
476 isNumericallyEqual(const float& a, const conjugate<P>& b,
477  double tol = RTraits<float>::getDefaultTolerance())
478 { return isNumericallyEqual(b,a,tol); }
480 template <class P> inline bool
481 isNumericallyEqual(const conjugate<P>& a, const double& b,
482  double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance())
483 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.,tol); }
485 template <class P> inline bool
486 isNumericallyEqual(const double& a, const conjugate<P>& b,
487  double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance())
488 { return isNumericallyEqual(b,a,tol); }
490 template <class P> inline bool
491 isNumericallyEqual(const conjugate<P>& a, const long double& b,
492  double tol = RTraits<P>::getDefaultTolerance())
493 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.L,tol); }
495 template <class P> inline bool
496 isNumericallyEqual(const long double& a, const conjugate<P>& b,
497  double tol = RTraits<P>::getDefaultTolerance())
498 { return isNumericallyEqual(b,a,tol); }
500 template <class P> inline bool
501 isNumericallyEqual(const conjugate<P>& a, int b,
502  double tol = RTraits<P>::getDefaultTolerance())
503 { typedef typename Widest<P,double>::Precision W; return isNumericallyEqual(a,(W)b,tol); }
505 template <class P> inline bool
506 isNumericallyEqual(int a, const conjugate<P>& b,
507  double tol = RTraits<P>::getDefaultTolerance())
508 { return isNumericallyEqual(b,a,tol); }
509 
511 
512 
513 template <class N> class NTraits {
514  // only the specializations below are allowed
515 };
516 
519 template <class R> class NTraits< complex<R> > {
520  typedef complex<R> C;
521 public:
522  typedef C T;
523  typedef negator<C> TNeg; // type of this after *cast* to its negative
524  typedef C TWithoutNegator; // type of this ignoring negator (there isn't one!)
525 
526  typedef R TReal;
527  typedef R TImag;
528  typedef C TComplex;
529  typedef conjugate<R> THerm; // this is a recast
530  typedef C TPosTrans;
531  typedef R TSqHermT; // ~C*C
532  typedef R TSqTHerm; // C*~C (same)
533  typedef C TElement;
534  typedef C TRow;
535  typedef C TCol;
536 
537  typedef C TSqrt;
538  typedef R TAbs;
539  typedef C TStandard; // complex is a standard type
540  typedef C TInvert; // this is a calculation, so use standard number
541  typedef C TNormalize;
542 
543  typedef C Scalar;
544  typedef C ULessScalar;
545  typedef C Number;
546  typedef C StdNumber;
547  typedef R Precision;
548  typedef R ScalarNormSq;
549 
550  // For complex scalar C, op result types are:
551  // Typeof(C*P) = Typeof(P*C)
552  // Typeof(C/P) = Typeof(inv(P)*C)
553  // Typeof(C+P) = Typeof(P+C)
554  // typeof(C-P) = Typeof(P::TNeg + C)
555  // These must be specialized for P=real, complex, conjugate.
556  template <class P> struct Result {
557  typedef typename CNT<P>::template Result<C>::Mul Mul;
558  typedef typename CNT< typename CNT<P>::THerm >::template Result<C>::Mul Dvd;
559  typedef typename CNT<P>::template Result<C>::Add Add;
560  typedef typename CNT< typename CNT<P>::TNeg >::template Result<C>::Add Sub;
561  };
562 
563  // Shape-preserving element substitution (easy for scalars!)
564  template <class P> struct Substitute {
565  typedef P Type;
566  };
567 
568  enum {
569  NRows = 1,
570  NCols = 1,
571  RowSpacing = 1,
572  ColSpacing = 1,
573  NPackedElements = 1, // not two!
574  NActualElements = 1,
575  NActualScalars = 1,
576  ImagOffset = 1,
577  RealStrideFactor = 2, // double stride when casting to real or imaginary
578  ArgDepth = SCALAR_DEPTH,
579  IsScalar = 1,
580  IsULessScalar = 1,
581  IsNumber = 1,
582  IsStdNumber = 1,
583  IsPrecision = 0,
584  SignInterpretation = 1 // after cast to Number, what is the sign?
585  };
586  static const T* getData(const T& t) { return &t; }
587  static T* updData(T& t) { return &t; }
588  static const R& real(const T& t) { return (reinterpret_cast<const R*>(&t))[0]; }
589  static R& real(T& t) { return (reinterpret_cast<R*>(&t))[0]; }
590  static const R& imag(const T& t) { return (reinterpret_cast<const R*>(&t))[1]; }
591  static R& imag(T& t) { return (reinterpret_cast<R*>(&t))[1]; }
592 
593  static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);}
594  static TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);}
595 
596  static const THerm& transpose(const T& t) {return reinterpret_cast<const THerm&>(t);}
597  static THerm& transpose(T& t) {return reinterpret_cast<THerm&>(t);}
598 
599  static const TPosTrans& positionalTranspose(const T& t)
600  {return reinterpret_cast<const TPosTrans&>(t);}
601  static TPosTrans& positionalTranspose(T& t)
602  {return reinterpret_cast<TPosTrans&>(t);}
603 
604  static const TWithoutNegator& castAwayNegatorIfAny(const T& t)
605  {return reinterpret_cast<const TWithoutNegator&>(t);}
606  static TWithoutNegator& updCastAwayNegatorIfAny(T& t)
607  {return reinterpret_cast<TWithoutNegator&>(t);}
608 
609  static ScalarNormSq scalarNormSqr(const T& t)
610  { return t.real()*t.real() + t.imag()*t.imag(); }
611  static TSqrt sqrt(const T& t)
612  { return std::sqrt(t); }
613  static TAbs abs(const T& t)
614  { return std::abs(t); } // no, not just sqrt of scalarNormSqr()!
615  static const TStandard& standardize(const T& t) {return t;} // already standard
616  static TNormalize normalize(const T& t) {return t/abs(t);}
617  static TInvert invert(const T& t) {return TReal(1)/t;}
618 
619  static const T& getNaN() {
620  static const T c=T(NTraits<R>::getNaN(), NTraits<R>::getNaN());
621  return c;
622  }
623  static const T& getInfinity() {
624  static const T c=T(NTraits<R>::getInfinity(),NTraits<R>::getInfinity());
625  return c;
626  }
627 
628  static const T& getI() {
629  static const T c = T(0,1);
630  return c;
631  }
632 
633  static bool isFinite(const T& t) {return SimTK::isFinite(t);}
634  static bool isNaN(const T& t) {return SimTK::isNaN(t);}
635  static bool isInf(const T& t) {return SimTK::isInf(t);}
636 
637  static double getDefaultTolerance() {return RTraits<R>::getDefaultTolerance();}
638 
639  template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b)
640  { return SimTK::isNumericallyEqual(a,b); }
641  template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b, double tol)
642  { return SimTK::isNumericallyEqual(a,b,tol); }
643  template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b)
644  { return SimTK::isNumericallyEqual(a,b); }
645  template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b, double tol)
646  { return SimTK::isNumericallyEqual(a,b,tol); }
647 
648  static bool isNumericallyEqual(const T& a, const float& b) {return SimTK::isNumericallyEqual(a,b);}
649  static bool isNumericallyEqual(const T& a, const float& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
650  static bool isNumericallyEqual(const T& a, const double& b) {return SimTK::isNumericallyEqual(a,b);}
651  static bool isNumericallyEqual(const T& a, const double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
652  static bool isNumericallyEqual(const T& a, const long double& b) {return SimTK::isNumericallyEqual(a,b);}
653  static bool isNumericallyEqual(const T& a, const long double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
654  static bool isNumericallyEqual(const T& a, int b) {return SimTK::isNumericallyEqual(a,b);}
655  static bool isNumericallyEqual(const T& a, int b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
656 
657  // The rest are the same as the real equivalents, with zero imaginary part.
658  static const T& getZero() {static const T c(NTraits<R>::getZero()); return c;}
659  static const T& getOne() {static const T c(NTraits<R>::getOne()); return c;}
660  static const T& getMinusOne() {static const T c(NTraits<R>::getMinusOne()); return c;}
661  static const T& getTwo() {static const T c(NTraits<R>::getTwo()); return c;}
662  static const T& getThree() {static const T c(NTraits<R>::getThree()); return c;}
663  static const T& getOneHalf() {static const T c(NTraits<R>::getOneHalf()); return c;}
664  static const T& getOneThird() {static const T c(NTraits<R>::getOneThird()); return c;}
665  static const T& getOneFourth() {static const T c(NTraits<R>::getOneFourth()); return c;}
666  static const T& getOneFifth() {static const T c(NTraits<R>::getOneFifth()); return c;}
667  static const T& getOneSixth() {static const T c(NTraits<R>::getOneSixth()); return c;}
668  static const T& getOneSeventh() {static const T c(NTraits<R>::getOneSeventh()); return c;}
669  static const T& getOneEighth() {static const T c(NTraits<R>::getOneEighth()); return c;}
670  static const T& getOneNinth() {static const T c(NTraits<R>::getOneNinth()); return c;}
671  static const T& getPi() {static const T c(NTraits<R>::getPi()); return c;}
672  static const T& getOneOverPi() {static const T c(NTraits<R>::getOneOverPi()); return c;}
673  static const T& getE() {static const T c(NTraits<R>::getE()); return c;}
674  static const T& getLog2E() {static const T c(NTraits<R>::getLog2E()); return c;}
675  static const T& getLog10E() {static const T c(NTraits<R>::getLog10E()); return c;}
676  static const T& getSqrt2() {static const T c(NTraits<R>::getSqrt2()); return c;}
677  static const T& getOneOverSqrt2() {static const T c(NTraits<R>::getOneOverSqrt2()); return c;}
678  static const T& getSqrt3() {static const T c(NTraits<R>::getSqrt3()); return c;}
679  static const T& getOneOverSqrt3() {static const T c(NTraits<R>::getOneOverSqrt3()); return c;}
680  static const T& getCubeRoot2() {static const T c(NTraits<R>::getCubeRoot2()); return c;}
681  static const T& getCubeRoot3() {static const T c(NTraits<R>::getCubeRoot3()); return c;}
682  static const T& getLn2() {static const T c(NTraits<R>::getLn2()); return c;}
683  static const T& getLn10() {static const T c(NTraits<R>::getLn10()); return c;}
684 };
685 
686 
687 // Specialize NTraits<complex>::Results for <complex> OP <scalar>. Result type is
688 // always just the complex type of sufficient precision.
689 #define SimTK_BNTCMPLX_SPEC(T1,T2) \
690 template<> template<> struct NTraits< complex<T1> >::Result<T2> { \
691  typedef Widest< complex<T1>,T2 >::Type W; \
692  typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \
693 }; \
694 template<> template<> struct NTraits< complex<T1> >::Result< complex<T2> > { \
695  typedef Widest< complex<T1>,complex<T2> >::Type W; \
696  typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \
697 }; \
698 template<> template<> struct NTraits< complex<T1> >::Result< conjugate<T2> > { \
699  typedef Widest< complex<T1>,complex<T2> >::Type W; \
700  typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \
701 }
702 SimTK_BNTCMPLX_SPEC(float,float);SimTK_BNTCMPLX_SPEC(float,double);SimTK_BNTCMPLX_SPEC(float,long double);
703 SimTK_BNTCMPLX_SPEC(double,float);SimTK_BNTCMPLX_SPEC(double,double);SimTK_BNTCMPLX_SPEC(double,long double);
704 SimTK_BNTCMPLX_SPEC(long double,float);SimTK_BNTCMPLX_SPEC(long double,double);SimTK_BNTCMPLX_SPEC(long double,long double);
705 #undef SimTK_BNTCMPLX_SPEC
706 
707 
708 // conjugate -- should be instantiated only for float, double, long double.
709 template <class R> class NTraits< conjugate<R> > {
710  typedef complex<R> C;
711 public:
712  typedef conjugate<R> T;
713  typedef negator<T> TNeg; // type of this after *cast* to its negative
714  typedef conjugate<R> TWithoutNegator; // type of this ignoring negator (there isn't one!)
715  typedef R TReal;
716  typedef negator<R> TImag;
717  typedef conjugate<R> TComplex;
718  typedef complex<R> THerm; // conjugate evaporates
719  typedef conjugate<R> TPosTrans; // Positional transpose of scalar does nothing
720  typedef R TSqHermT; // C*C~
721  typedef R TSqTHerm; // ~C*C (same)
722  typedef conjugate<R> TElement;
723  typedef conjugate<R> TRow;
724  typedef conjugate<R> TCol;
725 
726  typedef complex<R> TSqrt;
727  typedef R TAbs;
728  typedef complex<R> TStandard;
729  typedef conjugate<R> TInvert;
730  typedef conjugate<R> TNormalize;
731 
732  typedef conjugate<R> Scalar;
733  typedef conjugate<R> ULessScalar;
734  typedef conjugate<R> Number;
735  typedef complex<R> StdNumber;
736  typedef R Precision;
737  typedef R ScalarNormSq;
738 
739  // Typeof( Conj<S>*P ) is Typeof(P*Conj<S>)
740  // Typeof( Conj<S>/P ) is Typeof(inv(P)*Conj<S>)
741  // Typeof( Conj<S>+P ) is Typeof(P+Conj<S>)
742  // Typeof( Conj<S>-P ) is Typeof(P::TNeg+Conj<S>)
743  // Must specialize for P=real or P=complex or P=conjugate
744  template <class P> struct Result {
745  typedef typename CNT<P>::template Result<T>::Mul Mul;
746  typedef typename CNT<typename CNT<P>::THerm>::template Result<T>::Mul Dvd;
747  typedef typename CNT<P>::template Result<T>::Add Add;
748  typedef typename CNT<typename CNT<P>::TNeg>::template Result<T>::Add Sub;
749  };
750 
751  // Shape-preserving element substitution (easy for scalars!)
752  template <class P> struct Substitute {
753  typedef P Type;
754  };
755 
756  enum {
757  NRows = 1,
758  NCols = 1,
759  RowSpacing = 1,
760  ColSpacing = 1,
761  NPackedElements = 1, // not two!
762  NActualElements = 1,
763  NActualScalars = 1,
764  ImagOffset = 1,
765  RealStrideFactor = 2, // double stride when casting to real or imaginary
766  ArgDepth = SCALAR_DEPTH,
767  IsScalar = 1,
768  IsULessScalar = 1,
769  IsNumber = 1,
770  IsStdNumber = 0,
771  IsPrecision = 0,
772  SignInterpretation = 1 // after cast to Number, what is the sign?
773  };
774 
775  static const T* getData(const T& t) { return &t; }
776  static T* updData(T& t) { return &t; }
777  static const TReal& real(const T& t) { return t.real(); }
778  static TReal& real(T& t) { return t.real(); }
779  static const TImag& imag(const T& t) { return t.imag(); }
780  static TImag& imag(T& t) { return t.imag(); }
781 
782  static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);}
783  static TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);}
784 
785  static const THerm& transpose(const T& t) {return t.conj();}
786  static THerm& transpose(T& t) {return t.conj();}
787 
788  static const TPosTrans& positionalTranspose(const T& t)
789  {return reinterpret_cast<const TPosTrans&>(t);}
790  static TPosTrans& positionalTranspose(T& t)
791  {return reinterpret_cast<TPosTrans&>(t);}
792 
793  static const TWithoutNegator& castAwayNegatorIfAny(const T& t)
794  {return reinterpret_cast<const TWithoutNegator&>(t);}
795  static TWithoutNegator& updCastAwayNegatorIfAny(T& t)
796  {return reinterpret_cast<TWithoutNegator&>(t);}
797 
798  static ScalarNormSq scalarNormSqr(const T& t)
799  { return t.real()*t.real() + t.negImag()*t.negImag(); }
800  static TSqrt sqrt(const T& t)
801  { return std::sqrt(C(t)); } // cast to complex (one negation)
802  static TAbs abs(const T& t)
803  { return std::abs(t.conj()); } // no, not just sqrt of scalarNormSqr()!
804  static TStandard standardize(const T& t)
805  { return TStandard(t); } // i.e., convert to complex
806  static TNormalize normalize(const T& t) {return TNormalize(t/abs(t));}
807 
808  // 1/conj(z) = conj(1/z), for complex z.
809  static TInvert invert(const T& t)
810  { const typename NTraits<THerm>::TInvert cmplx(NTraits<THerm>::invert(t.conj()));
811  return reinterpret_cast<const TInvert&>(cmplx); } // recast complex to conjugate it
812 
813  // We want a "conjugate NaN", NaN - NaN*i, meaning both reals should
814  // be positive NaN.
815  static const T& getNaN() {
816  static const T c=T(NTraits<R>::getNaN(),NTraits<R>::getNaN());
817  return c;
818  }
819  // We want a "conjugate infinity", Inf - Inf*i, meaning both stored reals
820  // are positive Inf.
821  static const T& getInfinity() {
822  static const T c=T(NTraits<R>::getInfinity(),NTraits<R>::getInfinity());
823  return c;
824  }
825  // But we want the constant i (=sqrt(-1)) to be the same however we represent it,
826  // so for conjugate i = 0 - (-1)i.
827  static const T& getI() {
828  static const T c = T(0,-1);
829  return c;
830  }
831 
832  static bool isFinite(const T& t) {return SimTK::isFinite(t);}
833  static bool isNaN(const T& t) {return SimTK::isNaN(t);}
834  static bool isInf(const T& t) {return SimTK::isInf(t);}
835 
836  static double getDefaultTolerance() {return RTraits<R>::getDefaultTolerance();}
837 
838  template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b)
839  { return SimTK::isNumericallyEqual(a,b); }
840  template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b, double tol)
841  { return SimTK::isNumericallyEqual(a,b,tol); }
842  template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b)
843  { return SimTK::isNumericallyEqual(a,b); }
844  template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b, double tol)
845  { return SimTK::isNumericallyEqual(a,b,tol); }
846 
847  static bool isNumericallyEqual(const T& a, const float& b) {return SimTK::isNumericallyEqual(a,b);}
848  static bool isNumericallyEqual(const T& a, const float& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
849  static bool isNumericallyEqual(const T& a, const double& b) {return SimTK::isNumericallyEqual(a,b);}
850  static bool isNumericallyEqual(const T& a, const double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
851  static bool isNumericallyEqual(const T& a, const long double& b) {return SimTK::isNumericallyEqual(a,b);}
852  static bool isNumericallyEqual(const T& a, const long double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
853  static bool isNumericallyEqual(const T& a, int b) {return SimTK::isNumericallyEqual(a,b);}
854  static bool isNumericallyEqual(const T& a, int b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
855 
856  // The rest are the same as the real equivalents, with zero imaginary part.
857  static const T& getZero() {static const T c(NTraits<R>::getZero()); return c;}
858  static const T& getOne() {static const T c(NTraits<R>::getOne()); return c;}
859  static const T& getMinusOne() {static const T c(NTraits<R>::getMinusOne()); return c;}
860  static const T& getTwo() {static const T c(NTraits<R>::getTwo()); return c;}
861  static const T& getThree() {static const T c(NTraits<R>::getThree()); return c;}
862  static const T& getOneHalf() {static const T c(NTraits<R>::getOneHalf()); return c;}
863  static const T& getOneThird() {static const T c(NTraits<R>::getOneThird()); return c;}
864  static const T& getOneFourth() {static const T c(NTraits<R>::getOneFourth()); return c;}
865  static const T& getOneFifth() {static const T c(NTraits<R>::getOneFifth()); return c;}
866  static const T& getOneSixth() {static const T c(NTraits<R>::getOneSixth()); return c;}
867  static const T& getOneSeventh() {static const T c(NTraits<R>::getOneSeventh()); return c;}
868  static const T& getOneEighth() {static const T c(NTraits<R>::getOneEighth()); return c;}
869  static const T& getOneNinth() {static const T c(NTraits<R>::getOneNinth()); return c;}
870  static const T& getPi() {static const T c(NTraits<R>::getPi()); return c;}
871  static const T& getOneOverPi() {static const T c(NTraits<R>::getOneOverPi()); return c;}
872  static const T& getE() {static const T c(NTraits<R>::getE()); return c;}
873  static const T& getLog2E() {static const T c(NTraits<R>::getLog2E()); return c;}
874  static const T& getLog10E() {static const T c(NTraits<R>::getLog10E()); return c;}
875  static const T& getSqrt2() {static const T c(NTraits<R>::getSqrt2()); return c;}
876  static const T& getOneOverSqrt2() {static const T c(NTraits<R>::getOneOverSqrt2()); return c;}
877  static const T& getSqrt3() {static const T c(NTraits<R>::getSqrt3()); return c;}
878  static const T& getOneOverSqrt3() {static const T c(NTraits<R>::getOneOverSqrt3()); return c;}
879  static const T& getCubeRoot2() {static const T c(NTraits<R>::getCubeRoot2()); return c;}
880  static const T& getCubeRoot3() {static const T c(NTraits<R>::getCubeRoot3()); return c;}
881  static const T& getLn2() {static const T c(NTraits<R>::getLn2()); return c;}
882  static const T& getLn10() {static const T c(NTraits<R>::getLn10()); return c;}
883 };
884 
885 // Any op involving conjugate & a real is best left as a conjugate. However,
886 // an op involving conjugate & a complex or conjugate can lose the conjugate at zero cost
887 // and return just a complex in some cases. Also, we prefer negator<complex> to conjugate.
888 //
889 // Conj op complex
890 // a-bi * r+si = (ar+bs) + (as-br)i (complex)
891 // a-bi / r+si = hairy and slow anyway; we'll convert to complex
892 // a-bi + r+si = (a+r) + (s-b)i (complex)
893 // a-bi - r+si = (a-r) - (b+s)i = -[(r-a)+(b+s)i] (negator<complex>)
894 //
895 // Conj op Conj
896 // a-bi * r-si = (ar-bs) - (as+br)i = -[(bs-ar)+(as+br)i] (negator<complex>)
897 // a-bi / r-si = hairy and slow anyway; we'll convert to complex
898 // a-bi + r-si = (a+r) - (b+s)i (conjugate)
899 // a-bi - r-si = (a-r) + (s-b)i (complex)
900 
901 #define SimTK_NTRAITS_CONJ_SPEC(T1,T2) \
902 template<> template<> struct NTraits< conjugate<T1> >::Result<T2> { \
903  typedef conjugate<Widest<T1,T2>::Type> W; \
904  typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \
905 }; \
906 template<> template<> struct NTraits< conjugate<T1> >::Result<complex<T2> >{\
907  typedef Widest<complex<T1>,complex<T2> >::Type W; \
908  typedef W Mul; typedef W Dvd; typedef W Add; typedef negator<W> Sub; \
909 }; \
910 template<> template<> struct NTraits< conjugate<T1> >::Result<conjugate<T2> >{\
911  typedef Widest<T1,T2>::Type W; typedef complex<W> WC; \
912  typedef negator<WC> Mul; typedef WC Dvd; typedef conjugate<W> Add; typedef WC Sub;\
913 }
914 SimTK_NTRAITS_CONJ_SPEC(float,float);SimTK_NTRAITS_CONJ_SPEC(float,double);
915 SimTK_NTRAITS_CONJ_SPEC(float,long double);
916 SimTK_NTRAITS_CONJ_SPEC(double,float);SimTK_NTRAITS_CONJ_SPEC(double,double);
917 SimTK_NTRAITS_CONJ_SPEC(double,long double);
918 SimTK_NTRAITS_CONJ_SPEC(long double,float);SimTK_NTRAITS_CONJ_SPEC(long double,double);
919 SimTK_NTRAITS_CONJ_SPEC(long double,long double);
920 #undef SimTK_NTRAITS_CONJ_SPEC
921 
922 
923 // Specializations for real numbers.
924 // For real scalar R, op result types are:
925 // Typeof(R*P) = Typeof(P*R)
926 // Typeof(R/P) = Typeof(inv(P)*R)
927 // Typeof(R+P) = Typeof(P+R)
928 // typeof(R-P) = Typeof(P::TNeg + R)
929 // These must be specialized for P=Real and P=Complex.
930 #define SimTK_DEFINE_REAL_NTRAITS(R) \
931 template <> class NTraits<R> { \
932 public: \
933  typedef R T; \
934  typedef negator<T> TNeg; \
935  typedef T TWithoutNegator; \
936  typedef T TReal; \
937  typedef T TImag; \
938  typedef complex<T> TComplex; \
939  typedef T THerm; \
940  typedef T TPosTrans; \
941  typedef T TSqHermT; \
942  typedef T TSqTHerm; \
943  typedef T TElement; \
944  typedef T TRow; \
945  typedef T TCol; \
946  typedef T TSqrt; \
947  typedef T TAbs; \
948  typedef T TStandard; \
949  typedef T TInvert; \
950  typedef T TNormalize; \
951  typedef T Scalar; \
952  typedef T ULessScalar; \
953  typedef T Number; \
954  typedef T StdNumber; \
955  typedef T Precision; \
956  typedef T ScalarNormSq; \
957  template <class P> struct Result { \
958  typedef typename CNT<P>::template Result<R>::Mul Mul; \
959  typedef typename CNT< typename CNT<P>::THerm >::template Result<R>::Mul Dvd; \
960  typedef typename CNT<P>::template Result<R>::Add Add; \
961  typedef typename CNT< typename CNT<P>::TNeg >::template Result<R>::Add Sub; \
962  }; \
963  template <class P> struct Substitute { \
964  typedef P Type; \
965  }; \
966  enum { \
967  NRows = 1, \
968  NCols = 1, \
969  RowSpacing = 1, \
970  ColSpacing = 1, \
971  NPackedElements = 1, \
972  NActualElements = 1, \
973  NActualScalars = 1, \
974  ImagOffset = 0, \
975  RealStrideFactor = 1, \
976  ArgDepth = SCALAR_DEPTH, \
977  IsScalar = 1, \
978  IsULessScalar = 1, \
979  IsNumber = 1, \
980  IsStdNumber = 1, \
981  IsPrecision = 1, \
982  SignInterpretation = 1 \
983  }; \
984  static const T* getData(const T& t) { return &t; } \
985  static T* updData(T& t) { return &t; } \
986  static const T& real(const T& t) { return t; } \
987  static T& real(T& t) { return t; } \
988  static const T& imag(const T&) { return reinterpret_cast<const T&>(zeroes); } \
989  static T& imag(T&) { assert(false); return *reinterpret_cast<T*>(0); } \
990  static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);} \
991  static TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);} \
992  static const THerm& transpose(const T& t) {return reinterpret_cast<const THerm&>(t);} \
993  static THerm& transpose(T& t) {return reinterpret_cast<THerm&>(t);} \
994  static const TPosTrans& positionalTranspose(const T& t) \
995  {return reinterpret_cast<const TPosTrans&>(t);} \
996  static TPosTrans& positionalTranspose(T& t) \
997  {return reinterpret_cast<TPosTrans&>(t);} \
998  static const TWithoutNegator& castAwayNegatorIfAny(const T& t) \
999  {return reinterpret_cast<const TWithoutNegator&>(t);} \
1000  static TWithoutNegator& updCastAwayNegatorIfAny(T& t) \
1001  {return reinterpret_cast<TWithoutNegator&>(t);} \
1002  static ScalarNormSq scalarNormSqr(const T& t) {return t*t;} \
1003  static TSqrt sqrt(const T& t) {return std::sqrt(t);} \
1004  static TAbs abs(const T& t) {return std::abs(t);} \
1005  static const TStandard& standardize(const T& t) {return t;} \
1006  static TNormalize normalize(const T& t) {return (t>0?T(1):(t<0?T(-1):getNaN()));} \
1007  static TInvert invert(const T& t) {return T(1)/t;} \
1008  /* properties of this floating point representation, with memory addresses */ \
1009  static const T& getEps() {return RTraits<T>::getEps();} \
1010  static const T& getSignificant() {return RTraits<T>::getSignificant();} \
1011  static const T& getNaN() {static const T c=std::numeric_limits<T>::quiet_NaN(); return c;} \
1012  static const T& getInfinity() {static const T c=std::numeric_limits<T>::infinity(); return c;} \
1013  static const T& getLeastPositive(){static const T c=std::numeric_limits<T>::min(); return c;} \
1014  static const T& getMostPositive() {static const T c=std::numeric_limits<T>::max(); return c;} \
1015  static const T& getLeastNegative(){static const T c=-std::numeric_limits<T>::min(); return c;} \
1016  static const T& getMostNegative() {static const T c=-std::numeric_limits<T>::max(); return c;} \
1017  static const T& getSqrtEps() {static const T c=std::sqrt(getEps()); return c;} \
1018  static const T& getTiny() {static const T c=std::pow(getEps(), (T)1.25L); return c;} \
1019  static bool isFinite(const T& t) {return SimTK::isFinite(t);} \
1020  static bool isNaN (const T& t) {return SimTK::isNaN(t);} \
1021  static bool isInf (const T& t) {return SimTK::isInf(t);} \
1022  /* Methods to use for approximate comparisons. Perform comparison in the wider of the two */ \
1023  /* precisions, using the default tolerance from the narrower of the two precisions. */ \
1024  static double getDefaultTolerance() {return RTraits<T>::getDefaultTolerance();} \
1025  static bool isNumericallyEqual(const T& t, const float& f) {return SimTK::isNumericallyEqual(t,f);} \
1026  static bool isNumericallyEqual(const T& t, const double& d) {return SimTK::isNumericallyEqual(t,d);} \
1027  static bool isNumericallyEqual(const T& t, const long double& l) {return SimTK::isNumericallyEqual(t,l);} \
1028  static bool isNumericallyEqual(const T& t, int i) {return SimTK::isNumericallyEqual(t,i);} \
1029  /* Here the tolerance is given so we don't have to figure it out. */ \
1030  static bool isNumericallyEqual(const T& t, const float& f, double tol){return SimTK::isNumericallyEqual(t,f,tol);} \
1031  static bool isNumericallyEqual(const T& t, const double& d, double tol){return SimTK::isNumericallyEqual(t,d,tol);} \
1032  static bool isNumericallyEqual(const T& t, const long double& l, double tol){return SimTK::isNumericallyEqual(t,l,tol);} \
1033  static bool isNumericallyEqual(const T& t, int i, double tol){return SimTK::isNumericallyEqual(t,i,tol);} \
1034  /* Carefully calculated constants with convenient memory addresses. */ \
1035  static const T& getZero() {static const T c=(T)(0); return c;} \
1036  static const T& getOne() {static const T c=(T)(1); return c;} \
1037  static const T& getMinusOne() {static const T c=(T)(-1); return c;} \
1038  static const T& getTwo() {static const T c=(T)(2); return c;} \
1039  static const T& getThree() {static const T c=(T)(3); return c;} \
1040  static const T& getOneHalf() {static const T c=(T)(0.5L); return c;} \
1041  static const T& getOneThird() {static const T c=(T)(1.L/3.L); return c;} \
1042  static const T& getOneFourth() {static const T c=(T)(0.25L); return c;} \
1043  static const T& getOneFifth() {static const T c=(T)(0.2L); return c;} \
1044  static const T& getOneSixth() {static const T c=(T)(1.L/6.L); return c;} \
1045  static const T& getOneSeventh() {static const T c=(T)(1.L/7.L); return c;} \
1046  static const T& getOneEighth() {static const T c=(T)(0.125L); return c;} \
1047  static const T& getOneNinth() {static const T c=(T)(1.L/9.L); return c;} \
1048  static const T& getPi() {static const T c=(T)(SimTK_PI); return c;} \
1049  static const T& getOneOverPi() {static const T c=(T)(1.L/SimTK_PI); return c;} \
1050  static const T& getE() {static const T c=(T)(SimTK_E); return c;} \
1051  static const T& getLog2E() {static const T c=(T)(SimTK_LOG2E); return c;} \
1052  static const T& getLog10E() {static const T c=(T)(SimTK_LOG10E); return c;} \
1053  static const T& getSqrt2() {static const T c=(T)(SimTK_SQRT2); return c;} \
1054  static const T& getOneOverSqrt2() {static const T c=(T)(1.L/SimTK_SQRT2); return c;} \
1055  static const T& getSqrt3() {static const T c=(T)(SimTK_SQRT3); return c;} \
1056  static const T& getOneOverSqrt3() {static const T c=(T)(1.L/SimTK_SQRT3); return c;} \
1057  static const T& getCubeRoot2() {static const T c=(T)(SimTK_CBRT2); return c;} \
1058  static const T& getCubeRoot3() {static const T c=(T)(SimTK_CBRT3); return c;} \
1059  static const T& getLn2() {static const T c=(T)(SimTK_LN2); return c;} \
1060  static const T& getLn10() {static const T c=(T)(SimTK_LN10); return c;} \
1061  /* integer digit counts useful for formatted input and output */ \
1062  static const int getNumDigits() {static const int c=(int)(std::log10(1/getEps()) -0.5); return c;} \
1063  static const int getLosslessNumDigits() {static const int c=(int)(std::log10(1/getTiny())+0.5); return c;} \
1064 }; \
1065 template<> struct NTraits<R>::Result<float> \
1066  {typedef Widest<R,float>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1067 template<> struct NTraits<R>::Result<double> \
1068  {typedef Widest<R,double>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1069 template<> struct NTraits<R>::Result<long double> \
1070  {typedef Widest<R,long double>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1071 template<> struct NTraits<R>::Result<complex<float> > \
1072  {typedef Widest<R,complex<float> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1073 template<> struct NTraits<R>::Result<complex<double> > \
1074  {typedef Widest<R,complex<double> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1075 template<> struct NTraits<R>::Result<complex<long double> > \
1076  {typedef Widest<R,complex<long double> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1077 template<> struct NTraits<R>::Result<conjugate<float> > \
1078  {typedef conjugate<Widest<R,float>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1079 template<> struct NTraits<R>::Result<conjugate<double> > \
1080  {typedef conjugate<Widest<R,double>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1081 template<> struct NTraits<R>::Result<conjugate<long double> > \
1082  {typedef conjugate<Widest<R,long double>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}
1083 SimTK_DEFINE_REAL_NTRAITS(float);
1084 SimTK_DEFINE_REAL_NTRAITS(double);
1085 SimTK_DEFINE_REAL_NTRAITS(long double);
1086 #undef SimTK_DEFINE_REAL_NTRAITS
1087 
1089 template <class R> class CNT< complex<R> > : public NTraits< complex<R> > { };
1090 template <class R> class CNT< conjugate<R> > : public NTraits< conjugate<R> > { };
1091 template <> class CNT<float> : public NTraits<float> { };
1092 template <> class CNT<double> : public NTraits<double> { };
1093 template <> class CNT<long double> : public NTraits<long double> { };
1094 
1095 
1096 } // namespace SimTK
1097 
1098 #endif //SimTK_SIMMATRIX_NTRAITS_H_