Simbody  3.5
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"
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 
103 template <class R1, class R2> struct Widest {/* Only defined for built-ins. */};
104 template <> struct Widest<float,float> {typedef float Type; typedef float Precision;};
105 template <> struct Widest<float,double> {typedef double Type; typedef double Precision;};
106 template <> struct Widest<float,long double> {typedef long double Type; typedef long double Precision;};
107 template <> struct Widest<double,float> {typedef double Type; typedef double Precision;};
108 template <> struct Widest<double,double> {typedef double Type; typedef double Precision;};
109 template <> struct Widest<double,long double> {typedef long double Type; typedef long double Precision;};
110 template <> struct Widest<long double,float> {typedef long double Type; typedef long double Precision;};
111 template <> struct Widest<long double,double> {typedef long double Type; typedef long double Precision;};
112 template <> struct Widest<long double,long double> {typedef long double Type; typedef long double Precision;};
113 template <class R1, class R2> struct Widest< complex<R1>,complex<R2> > {
114  typedef complex< typename Widest<R1,R2>::Type > Type;
116 };
117 template <class R1, class R2> struct Widest< complex<R1>,R2 > {
118  typedef complex< typename Widest<R1,R2>::Type > Type;
120 };
121 template <class R1, class R2> struct Widest< R1,complex<R2> > {
122  typedef complex< typename Widest<R1,R2>::Type > Type;
124 };
125 
135 template <class R1, class R2> struct Narrowest {/* Only defined for built-ins. */};
136 template <> struct Narrowest<float,float> {typedef float Type; typedef float Precision;};
137 template <> struct Narrowest<float,double> {typedef float Type; typedef float Precision;};
138 template <> struct Narrowest<float,long double> {typedef float Type; typedef float Precision;};
139 template <> struct Narrowest<double,float> {typedef float Type; typedef float Precision;};
140 template <> struct Narrowest<double,double> {typedef double Type; typedef double Precision;};
141 template <> struct Narrowest<double,long double> {typedef double Type; typedef double Precision;};
142 template <> struct Narrowest<long double,float> {typedef float Type; typedef float Precision;};
143 template <> struct Narrowest<long double,double> {typedef double Type; typedef double Precision;};
144 template <> struct Narrowest<long double,long double> {typedef long double Type; typedef long double Precision;};
145 template <class R1, class R2> struct Narrowest< complex<R1>,complex<R2> > {
146  typedef complex< typename Narrowest<R1,R2>::Type > Type;
148 };
149 template <class R1, class R2> struct Narrowest< complex<R1>,R2 > {
150  typedef complex< typename Narrowest<R1,R2>::Type > Type;
152 };
153 template <class R1, class R2> struct Narrowest< R1,complex<R2> > {
154  typedef complex< typename Narrowest<R1,R2>::Type > Type;
156 };
157 
159 template <class R> class RTraits {/* Only defined for real types */};
160 template <> class RTraits<float> {
161 public:
163  static const float& getEps() {static const float c=std::numeric_limits<float>::epsilon(); return c;}
165  static const float& getSignificant() {static const float c=std::pow(getEps(), 0.875f); return c;}
167  static double getDefaultTolerance() {return (double)getSignificant();}
168 };
169 template <> class RTraits<double> {
170 public:
171  static const double& getEps() {static const double c=std::numeric_limits<double>::epsilon(); return c;}
172  static const double& getSignificant() {static const double c=std::pow(getEps(), 0.875); return c;}
173  static double getDefaultTolerance() {return getSignificant();}
174 };
175 template <> class RTraits<long double> {
176 public:
177  static const long double& getEps() {static const long double c=std::numeric_limits<long double>::epsilon(); return c;}
178  static const long double& getSignificant() {static const long double c=std::pow(getEps(), 0.875L); return c;}
179  static double getDefaultTolerance() {return (double)getSignificant();}
180 };
181 
198 // See negator.h for isNaN() applied to negated scalars.
200 inline bool isNaN(const float& x) {return std::isnan(x);}
201 inline bool isNaN(const double& x) {return std::isnan(x);}
202 inline bool isNaN(const long double& x) {return std::isnan(x);}
203 
204 template <class P> inline bool
205 isNaN(const std::complex<P>& x)
206 { return isNaN(x.real()) || isNaN(x.imag());}
207 
208 template <class P> inline bool
210 { return isNaN(x.real()) || isNaN(x.negImag());}
212 
225 // See negator.h for isFinite() applied to negated scalars.
227 inline bool isFinite(const float& x) {return std::isfinite(x);}
228 inline bool isFinite(const double& x) {return std::isfinite(x);}
229 inline bool isFinite(const long double& x) {return std::isfinite(x);}
230 
231 template <class P> inline bool
232 isFinite(const std::complex<P>& x)
233 { return isFinite(x.real()) && isFinite(x.imag());}
234 
235 template <class P> inline bool
237 { return isFinite(x.real()) && isFinite(x.negImag());}
239 
254 // See negator.h for isInf() applied to negated scalars.
256 inline bool isInf(const float& x) {return std::isinf(x);}
257 inline bool isInf(const double& x) {return std::isinf(x);}
258 inline bool isInf(const long double& x) {return std::isinf(x);}
259 
260 template <class P> inline bool
261 isInf(const std::complex<P>& x) {
262  return (isInf(x.real()) && !isNaN(x.imag()))
263  || (isInf(x.imag()) && !isNaN(x.real()));
264 }
265 
266 template <class P> inline bool
267 isInf(const conjugate<P>& x) {
268  return (isInf(x.real()) && !isNaN(x.negImag()))
269  || (isInf(x.negImag()) && !isNaN(x.real()));
270 }
272 
312 inline bool isNumericallyEqual(const float& a, const float& b,
315 { if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false;
316  const float scale = std::max(std::max(std::abs(a),std::abs(b)), 1.f);
317  return std::abs(a-b) <= scale*(float)tol; }
319 inline bool isNumericallyEqual(const double& a, const double& b,
321 { if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false;
322  const double scale = std::max(std::max(std::abs(a),std::abs(b)), 1.);
323  return std::abs(a-b) <= scale*tol; }
325 inline bool isNumericallyEqual(const long double& a, const long double& b,
327 { if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false;
328  const long double scale = std::max(std::max(std::abs(a),std::abs(b)), 1.L);
329  return std::abs(a-b) <= scale*(long double)tol; }
330 
332 inline bool isNumericallyEqual(const float& a, const double& b,
334 { return isNumericallyEqual((double)a, b, tol); }
336 inline bool isNumericallyEqual(const double& a, const float& b,
338 { return isNumericallyEqual(a, (double)b, tol); }
340 inline bool isNumericallyEqual(const float& a, const long double& b,
342 { return isNumericallyEqual((long double)a, b, tol); }
344 inline bool isNumericallyEqual(const long double& a, const float& b,
346 { return isNumericallyEqual(a, (long double)b, tol); }
348 inline bool isNumericallyEqual(const double& a, const long double& b,
350 { return isNumericallyEqual((long double)a, b, tol); }
352 inline bool isNumericallyEqual(const long double& a, const double& b,
354 { return isNumericallyEqual(a, (long double)b, tol); }
355 
357 inline bool isNumericallyEqual(const float& a, int b,
359 { return isNumericallyEqual(a, (double)b, tol); }
361 inline bool isNumericallyEqual(int a, const float& b,
363 { return isNumericallyEqual((double)a, b, tol); }
365 inline bool isNumericallyEqual(const double& a, int b,
367 { return isNumericallyEqual(a, (double)b, tol); }
369 inline bool isNumericallyEqual(int a, const double& b,
371 { return isNumericallyEqual((double)a, b, tol); }
373 inline bool isNumericallyEqual(const long double& a, int b,
375 { return isNumericallyEqual(a, (long double)b, tol); }
377 inline bool isNumericallyEqual(int a, const long double& b,
379 { return isNumericallyEqual((long double)a, b, tol); }
380 
384 template <class P, class Q>
385 inline bool isNumericallyEqual
386  ( const std::complex<P>& a, const std::complex<Q>& b,
387  double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance())
388 { return isNumericallyEqual(a.real(),b.real(),tol)
389  && isNumericallyEqual(a.imag(),b.imag(),tol); }
390 
394 template <class P, class Q>
395 inline bool isNumericallyEqual
396  ( const conjugate<P>& a, const conjugate<Q>& b,
397  double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance())
398 { return isNumericallyEqual(a.real(),b.real(),tol)
399  && isNumericallyEqual(a.imag(),b.imag(),tol); }
400 
404 template <class P, class Q>
405 inline bool isNumericallyEqual
406  ( const std::complex<P>& a, const conjugate<Q>& b,
407  double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance())
408 { return isNumericallyEqual(a.real(),b.real(),tol)
409  && isNumericallyEqual(a.imag(),b.imag(),tol); }
410 
414 template <class P, class Q>
415 inline bool isNumericallyEqual
416  ( const conjugate<P>& a, const std::complex<Q>& b,
417  double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance())
418 { return isNumericallyEqual(a.real(),b.real(),tol)
419  && isNumericallyEqual(a.imag(),b.imag(),tol); }
420 
422 template <class P> inline bool
423 isNumericallyEqual(const std::complex<P>& a, const float& b,
425 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.f,tol); }
427 template <class P> inline bool
428 isNumericallyEqual(const float& a, const std::complex<P>& b,
430 { return isNumericallyEqual(b,a,tol); }
432 template <class P> inline bool
433 isNumericallyEqual(const std::complex<P>& a, const double& b,
434  double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance())
435 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.,tol); }
437 template <class P> inline bool
438 isNumericallyEqual(const double& a, const std::complex<P>& b,
439  double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance())
440 { return isNumericallyEqual(b,a,tol); }
442 template <class P> inline bool
443 isNumericallyEqual(const std::complex<P>& a, const long double& b,
444  double tol = RTraits<P>::getDefaultTolerance())
445 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.L,tol); }
447 template <class P> inline bool
448 isNumericallyEqual(const long double& a, const std::complex<P>& b,
449  double tol = RTraits<P>::getDefaultTolerance())
450 { return isNumericallyEqual(b,a,tol); }
452 template <class P> inline bool
453 isNumericallyEqual(const std::complex<P>& a, int b,
454  double tol = RTraits<P>::getDefaultTolerance())
455 { typedef typename Widest<P,double>::Precision W; return isNumericallyEqual(a,(W)b,tol); }
457 template <class P> inline bool
458 isNumericallyEqual(int a, const std::complex<P>& b,
459  double tol = RTraits<P>::getDefaultTolerance())
460 { return isNumericallyEqual(b,a,tol); }
461 
463 template <class P> inline bool
464 isNumericallyEqual(const conjugate<P>& a, const float& b,
466 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.f,tol); }
468 template <class P> inline bool
469 isNumericallyEqual(const float& a, const conjugate<P>& b,
471 { return isNumericallyEqual(b,a,tol); }
473 template <class P> inline bool
474 isNumericallyEqual(const conjugate<P>& a, const double& b,
475  double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance())
476 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.,tol); }
478 template <class P> inline bool
479 isNumericallyEqual(const double& a, const conjugate<P>& b,
480  double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance())
481 { return isNumericallyEqual(b,a,tol); }
483 template <class P> inline bool
484 isNumericallyEqual(const conjugate<P>& a, const long double& b,
485  double tol = RTraits<P>::getDefaultTolerance())
486 { return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.L,tol); }
488 template <class P> inline bool
489 isNumericallyEqual(const long double& a, const conjugate<P>& b,
490  double tol = RTraits<P>::getDefaultTolerance())
491 { return isNumericallyEqual(b,a,tol); }
493 template <class P> inline bool
495  double tol = RTraits<P>::getDefaultTolerance())
496 { typedef typename Widest<P,double>::Precision W; return isNumericallyEqual(a,(W)b,tol); }
498 template <class P> inline bool
500  double tol = RTraits<P>::getDefaultTolerance())
501 { return isNumericallyEqual(b,a,tol); }
502 
504 
505 
506 template <class N> class NTraits {
507  // only the specializations below are allowed
508 };
509 
512 template <class R> class NTraits< complex<R> > {
513  typedef complex<R> C;
514 public:
515  typedef C T;
516  typedef negator<C> TNeg; // type of this after *cast* to its negative
517  typedef C TWithoutNegator; // type of this ignoring negator (there isn't one!)
518 
519  typedef R TReal;
520  typedef R TImag;
521  typedef C TComplex;
522  typedef conjugate<R> THerm; // this is a recast
523  typedef C TPosTrans;
524  typedef R TSqHermT; // ~C*C
525  typedef R TSqTHerm; // C*~C (same)
526  typedef C TElement;
527  typedef C TRow;
528  typedef C TCol;
529 
530  typedef C TSqrt;
531  typedef R TAbs;
532  typedef C TStandard; // complex is a standard type
533  typedef C TInvert; // this is a calculation, so use standard number
534  typedef C TNormalize;
535 
536  typedef C Scalar;
537  typedef C ULessScalar;
538  typedef C Number;
539  typedef C StdNumber;
540  typedef R Precision;
541  typedef R ScalarNormSq;
542 
543  // For complex scalar C, op result types are:
544  // Typeof(C*P) = Typeof(P*C)
545  // Typeof(C/P) = Typeof(inv(P)*C)
546  // Typeof(C+P) = Typeof(P+C)
547  // typeof(C-P) = Typeof(P::TNeg + C)
548  // These must be specialized for P=real, complex, conjugate.
549  template <class P> struct Result {
551  typedef typename CNT< typename CNT<P>::THerm >::template Result<C>::Mul Dvd;
553  typedef typename CNT< typename CNT<P>::TNeg >::template Result<C>::Add Sub;
554  };
555 
556  // Shape-preserving element substitution (easy for scalars!)
557  template <class P> struct Substitute {
558  typedef P Type;
559  };
560 
561  enum {
562  NRows = 1,
563  NCols = 1,
564  RowSpacing = 1,
565  ColSpacing = 1,
566  NPackedElements = 1, // not two!
567  NActualElements = 1,
568  NActualScalars = 1,
569  ImagOffset = 1,
570  RealStrideFactor = 2, // double stride when casting to real or imaginary
571  ArgDepth = SCALAR_DEPTH,
572  IsScalar = 1,
573  IsULessScalar = 1,
574  IsNumber = 1,
575  IsStdNumber = 1,
576  IsPrecision = 0,
577  SignInterpretation = 1 // after cast to Number, what is the sign?
578  };
579  static const T* getData(const T& t) { return &t; }
580  static T* updData(T& t) { return &t; }
581  static const R& real(const T& t) { return (reinterpret_cast<const R*>(&t))[0]; }
582  static R& real(T& t) { return (reinterpret_cast<R*>(&t))[0]; }
583  static const R& imag(const T& t) { return (reinterpret_cast<const R*>(&t))[1]; }
584  static R& imag(T& t) { return (reinterpret_cast<R*>(&t))[1]; }
585 
586  static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);}
587  static TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);}
588 
589  static const THerm& transpose(const T& t) {return reinterpret_cast<const THerm&>(t);}
590  static THerm& transpose(T& t) {return reinterpret_cast<THerm&>(t);}
591 
592  static const TPosTrans& positionalTranspose(const T& t)
593  {return reinterpret_cast<const TPosTrans&>(t);}
594  static TPosTrans& positionalTranspose(T& t)
595  {return reinterpret_cast<TPosTrans&>(t);}
596 
597  static const TWithoutNegator& castAwayNegatorIfAny(const T& t)
598  {return reinterpret_cast<const TWithoutNegator&>(t);}
599  static TWithoutNegator& updCastAwayNegatorIfAny(T& t)
600  {return reinterpret_cast<TWithoutNegator&>(t);}
601 
602  static ScalarNormSq scalarNormSqr(const T& t)
603  { return t.real()*t.real() + t.imag()*t.imag(); }
604  static TSqrt sqrt(const T& t)
605  { return std::sqrt(t); }
606  static TAbs abs(const T& t)
607  { return std::abs(t); } // no, not just sqrt of scalarNormSqr()!
608  static const TStandard& standardize(const T& t) {return t;} // already standard
609  static TNormalize normalize(const T& t) {return t/abs(t);}
610  static TInvert invert(const T& t) {return TReal(1)/t;}
611 
612  static const T& getNaN() {
613  static const T c=T(NTraits<R>::getNaN(), NTraits<R>::getNaN());
614  return c;
615  }
616  static const T& getInfinity() {
617  static const T c=T(NTraits<R>::getInfinity(),NTraits<R>::getInfinity());
618  return c;
619  }
620 
621  static const T& getI() {
622  static const T c = T(0,1);
623  return c;
624  }
625 
626  static bool isFinite(const T& t) {return SimTK::isFinite(t);}
627  static bool isNaN(const T& t) {return SimTK::isNaN(t);}
628  static bool isInf(const T& t) {return SimTK::isInf(t);}
629 
631 
632  template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b)
633  { return SimTK::isNumericallyEqual(a,b); }
634  template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b, double tol)
635  { return SimTK::isNumericallyEqual(a,b,tol); }
636  template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b)
637  { return SimTK::isNumericallyEqual(a,b); }
638  template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b, double tol)
639  { return SimTK::isNumericallyEqual(a,b,tol); }
640 
641  static bool isNumericallyEqual(const T& a, const float& b) {return SimTK::isNumericallyEqual(a,b);}
642  static bool isNumericallyEqual(const T& a, const float& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
643  static bool isNumericallyEqual(const T& a, const double& b) {return SimTK::isNumericallyEqual(a,b);}
644  static bool isNumericallyEqual(const T& a, const double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
645  static bool isNumericallyEqual(const T& a, const long double& b) {return SimTK::isNumericallyEqual(a,b);}
646  static bool isNumericallyEqual(const T& a, const long double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
647  static bool isNumericallyEqual(const T& a, int b) {return SimTK::isNumericallyEqual(a,b);}
648  static bool isNumericallyEqual(const T& a, int b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
649 
650  // The rest are the same as the real equivalents, with zero imaginary part.
651  static const T& getZero() {static const T c(NTraits<R>::getZero()); return c;}
652  static const T& getOne() {static const T c(NTraits<R>::getOne()); return c;}
653  static const T& getMinusOne() {static const T c(NTraits<R>::getMinusOne()); return c;}
654  static const T& getTwo() {static const T c(NTraits<R>::getTwo()); return c;}
655  static const T& getThree() {static const T c(NTraits<R>::getThree()); return c;}
656  static const T& getOneHalf() {static const T c(NTraits<R>::getOneHalf()); return c;}
657  static const T& getOneThird() {static const T c(NTraits<R>::getOneThird()); return c;}
658  static const T& getOneFourth() {static const T c(NTraits<R>::getOneFourth()); return c;}
659  static const T& getOneFifth() {static const T c(NTraits<R>::getOneFifth()); return c;}
660  static const T& getOneSixth() {static const T c(NTraits<R>::getOneSixth()); return c;}
661  static const T& getOneSeventh() {static const T c(NTraits<R>::getOneSeventh()); return c;}
662  static const T& getOneEighth() {static const T c(NTraits<R>::getOneEighth()); return c;}
663  static const T& getOneNinth() {static const T c(NTraits<R>::getOneNinth()); return c;}
664  static const T& getPi() {static const T c(NTraits<R>::getPi()); return c;}
665  static const T& getOneOverPi() {static const T c(NTraits<R>::getOneOverPi()); return c;}
666  static const T& getE() {static const T c(NTraits<R>::getE()); return c;}
667  static const T& getLog2E() {static const T c(NTraits<R>::getLog2E()); return c;}
668  static const T& getLog10E() {static const T c(NTraits<R>::getLog10E()); return c;}
669  static const T& getSqrt2() {static const T c(NTraits<R>::getSqrt2()); return c;}
670  static const T& getOneOverSqrt2() {static const T c(NTraits<R>::getOneOverSqrt2()); return c;}
671  static const T& getSqrt3() {static const T c(NTraits<R>::getSqrt3()); return c;}
672  static const T& getOneOverSqrt3() {static const T c(NTraits<R>::getOneOverSqrt3()); return c;}
673  static const T& getCubeRoot2() {static const T c(NTraits<R>::getCubeRoot2()); return c;}
674  static const T& getCubeRoot3() {static const T c(NTraits<R>::getCubeRoot3()); return c;}
675  static const T& getLn2() {static const T c(NTraits<R>::getLn2()); return c;}
676  static const T& getLn10() {static const T c(NTraits<R>::getLn10()); return c;}
677 };
678 
679 
680 // Specialize NTraits<complex>::Results for <complex> OP <scalar>. Result type is
681 // always just the complex type of sufficient precision.
682 #define SimTK_BNTCMPLX_SPEC(T1,T2) \
683 template<> template<> struct NTraits< complex<T1> >::Result<T2> { \
684  typedef Widest< complex<T1>,T2 >::Type W; \
685  typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \
686 }; \
687 template<> template<> struct NTraits< complex<T1> >::Result< complex<T2> > { \
688  typedef Widest< complex<T1>,complex<T2> >::Type W; \
689  typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \
690 }; \
691 template<> template<> struct NTraits< complex<T1> >::Result< conjugate<T2> > { \
692  typedef Widest< complex<T1>,complex<T2> >::Type W; \
693  typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \
694 }
695 SimTK_BNTCMPLX_SPEC(float,float);SimTK_BNTCMPLX_SPEC(float,double);SimTK_BNTCMPLX_SPEC(float,long double);
696 SimTK_BNTCMPLX_SPEC(double,float);SimTK_BNTCMPLX_SPEC(double,double);SimTK_BNTCMPLX_SPEC(double,long double);
697 SimTK_BNTCMPLX_SPEC(long double,float);SimTK_BNTCMPLX_SPEC(long double,double);SimTK_BNTCMPLX_SPEC(long double,long double);
698 #undef SimTK_BNTCMPLX_SPEC
699 
700 
701 // conjugate -- should be instantiated only for float, double, long double.
702 template <class R> class NTraits< conjugate<R> > {
703  typedef complex<R> C;
704 public:
705  typedef conjugate<R> T;
706  typedef negator<T> TNeg; // type of this after *cast* to its negative
707  typedef conjugate<R> TWithoutNegator; // type of this ignoring negator (there isn't one!)
708  typedef R TReal;
709  typedef negator<R> TImag;
711  typedef complex<R> THerm; // conjugate evaporates
712  typedef conjugate<R> TPosTrans; // Positional transpose of scalar does nothing
713  typedef R TSqHermT; // C*C~
714  typedef R TSqTHerm; // ~C*C (same)
718 
719  typedef complex<R> TSqrt;
720  typedef R TAbs;
721  typedef complex<R> TStandard;
724 
728  typedef complex<R> StdNumber;
729  typedef R Precision;
730  typedef R ScalarNormSq;
731 
732  // Typeof( Conj<S>*P ) is Typeof(P*Conj<S>)
733  // Typeof( Conj<S>/P ) is Typeof(inv(P)*Conj<S>)
734  // Typeof( Conj<S>+P ) is Typeof(P+Conj<S>)
735  // Typeof( Conj<S>-P ) is Typeof(P::TNeg+Conj<S>)
736  // Must specialize for P=real or P=complex or P=conjugate
737  template <class P> struct Result {
739  typedef typename CNT<typename CNT<P>::THerm>::template Result<T>::Mul Dvd;
741  typedef typename CNT<typename CNT<P>::TNeg>::template Result<T>::Add Sub;
742  };
743 
744  // Shape-preserving element substitution (easy for scalars!)
745  template <class P> struct Substitute {
746  typedef P Type;
747  };
748 
749  enum {
750  NRows = 1,
751  NCols = 1,
752  RowSpacing = 1,
753  ColSpacing = 1,
754  NPackedElements = 1, // not two!
755  NActualElements = 1,
756  NActualScalars = 1,
757  ImagOffset = 1,
758  RealStrideFactor = 2, // double stride when casting to real or imaginary
759  ArgDepth = SCALAR_DEPTH,
760  IsScalar = 1,
761  IsULessScalar = 1,
762  IsNumber = 1,
763  IsStdNumber = 0,
764  IsPrecision = 0,
765  SignInterpretation = 1 // after cast to Number, what is the sign?
766  };
767 
768  static const T* getData(const T& t) { return &t; }
769  static T* updData(T& t) { return &t; }
770  static const TReal& real(const T& t) { return t.real(); }
771  static TReal& real(T& t) { return t.real(); }
772  static const TImag& imag(const T& t) { return t.imag(); }
773  static TImag& imag(T& t) { return t.imag(); }
774 
775  static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);}
776  static TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);}
777 
778  static const THerm& transpose(const T& t) {return t.conj();}
779  static THerm& transpose(T& t) {return t.conj();}
780 
781  static const TPosTrans& positionalTranspose(const T& t)
782  {return reinterpret_cast<const TPosTrans&>(t);}
783  static TPosTrans& positionalTranspose(T& t)
784  {return reinterpret_cast<TPosTrans&>(t);}
785 
786  static const TWithoutNegator& castAwayNegatorIfAny(const T& t)
787  {return reinterpret_cast<const TWithoutNegator&>(t);}
788  static TWithoutNegator& updCastAwayNegatorIfAny(T& t)
789  {return reinterpret_cast<TWithoutNegator&>(t);}
790 
791  static ScalarNormSq scalarNormSqr(const T& t)
792  { return t.real()*t.real() + t.negImag()*t.negImag(); }
793  static TSqrt sqrt(const T& t)
794  { return std::sqrt(C(t)); } // cast to complex (one negation)
795  static TAbs abs(const T& t)
796  { return std::abs(t.conj()); } // no, not just sqrt of scalarNormSqr()!
797  static TStandard standardize(const T& t)
798  { return TStandard(t); } // i.e., convert to complex
799  static TNormalize normalize(const T& t) {return TNormalize(t/abs(t));}
800 
801  // 1/conj(z) = conj(1/z), for complex z.
802  static TInvert invert(const T& t)
803  { const typename NTraits<THerm>::TInvert cmplx(NTraits<THerm>::invert(t.conj()));
804  return reinterpret_cast<const TInvert&>(cmplx); } // recast complex to conjugate it
805 
806  // We want a "conjugate NaN", NaN - NaN*i, meaning both reals should
807  // be positive NaN.
808  static const T& getNaN() {
809  static const T c=T(NTraits<R>::getNaN(),NTraits<R>::getNaN());
810  return c;
811  }
812  // We want a "conjugate infinity", Inf - Inf*i, meaning both stored reals
813  // are positive Inf.
814  static const T& getInfinity() {
815  static const T c=T(NTraits<R>::getInfinity(),NTraits<R>::getInfinity());
816  return c;
817  }
818  // But we want the constant i (=sqrt(-1)) to be the same however we represent it,
819  // so for conjugate i = 0 - (-1)i.
820  static const T& getI() {
821  static const T c = T(0,-1);
822  return c;
823  }
824 
825  static bool isFinite(const T& t) {return SimTK::isFinite(t);}
826  static bool isNaN(const T& t) {return SimTK::isNaN(t);}
827  static bool isInf(const T& t) {return SimTK::isInf(t);}
828 
830 
831  template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b)
832  { return SimTK::isNumericallyEqual(a,b); }
833  template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b, double tol)
834  { return SimTK::isNumericallyEqual(a,b,tol); }
835  template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b)
836  { return SimTK::isNumericallyEqual(a,b); }
837  template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b, double tol)
838  { return SimTK::isNumericallyEqual(a,b,tol); }
839 
840  static bool isNumericallyEqual(const T& a, const float& b) {return SimTK::isNumericallyEqual(a,b);}
841  static bool isNumericallyEqual(const T& a, const float& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
842  static bool isNumericallyEqual(const T& a, const double& b) {return SimTK::isNumericallyEqual(a,b);}
843  static bool isNumericallyEqual(const T& a, const double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
844  static bool isNumericallyEqual(const T& a, const long double& b) {return SimTK::isNumericallyEqual(a,b);}
845  static bool isNumericallyEqual(const T& a, const long double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
846  static bool isNumericallyEqual(const T& a, int b) {return SimTK::isNumericallyEqual(a,b);}
847  static bool isNumericallyEqual(const T& a, int b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
848 
849  // The rest are the same as the real equivalents, with zero imaginary part.
850  static const T& getZero() {static const T c(NTraits<R>::getZero()); return c;}
851  static const T& getOne() {static const T c(NTraits<R>::getOne()); return c;}
852  static const T& getMinusOne() {static const T c(NTraits<R>::getMinusOne()); return c;}
853  static const T& getTwo() {static const T c(NTraits<R>::getTwo()); return c;}
854  static const T& getThree() {static const T c(NTraits<R>::getThree()); return c;}
855  static const T& getOneHalf() {static const T c(NTraits<R>::getOneHalf()); return c;}
856  static const T& getOneThird() {static const T c(NTraits<R>::getOneThird()); return c;}
857  static const T& getOneFourth() {static const T c(NTraits<R>::getOneFourth()); return c;}
858  static const T& getOneFifth() {static const T c(NTraits<R>::getOneFifth()); return c;}
859  static const T& getOneSixth() {static const T c(NTraits<R>::getOneSixth()); return c;}
860  static const T& getOneSeventh() {static const T c(NTraits<R>::getOneSeventh()); return c;}
861  static const T& getOneEighth() {static const T c(NTraits<R>::getOneEighth()); return c;}
862  static const T& getOneNinth() {static const T c(NTraits<R>::getOneNinth()); return c;}
863  static const T& getPi() {static const T c(NTraits<R>::getPi()); return c;}
864  static const T& getOneOverPi() {static const T c(NTraits<R>::getOneOverPi()); return c;}
865  static const T& getE() {static const T c(NTraits<R>::getE()); return c;}
866  static const T& getLog2E() {static const T c(NTraits<R>::getLog2E()); return c;}
867  static const T& getLog10E() {static const T c(NTraits<R>::getLog10E()); return c;}
868  static const T& getSqrt2() {static const T c(NTraits<R>::getSqrt2()); return c;}
869  static const T& getOneOverSqrt2() {static const T c(NTraits<R>::getOneOverSqrt2()); return c;}
870  static const T& getSqrt3() {static const T c(NTraits<R>::getSqrt3()); return c;}
871  static const T& getOneOverSqrt3() {static const T c(NTraits<R>::getOneOverSqrt3()); return c;}
872  static const T& getCubeRoot2() {static const T c(NTraits<R>::getCubeRoot2()); return c;}
873  static const T& getCubeRoot3() {static const T c(NTraits<R>::getCubeRoot3()); return c;}
874  static const T& getLn2() {static const T c(NTraits<R>::getLn2()); return c;}
875  static const T& getLn10() {static const T c(NTraits<R>::getLn10()); return c;}
876 };
877 
878 // Any op involving conjugate & a real is best left as a conjugate. However,
879 // an op involving conjugate & a complex or conjugate can lose the conjugate at zero cost
880 // and return just a complex in some cases. Also, we prefer negator<complex> to conjugate.
881 //
882 // Conj op complex
883 // a-bi * r+si = (ar+bs) + (as-br)i (complex)
884 // a-bi / r+si = hairy and slow anyway; we'll convert to complex
885 // a-bi + r+si = (a+r) + (s-b)i (complex)
886 // a-bi - r+si = (a-r) - (b+s)i = -[(r-a)+(b+s)i] (negator<complex>)
887 //
888 // Conj op Conj
889 // a-bi * r-si = (ar-bs) - (as+br)i = -[(bs-ar)+(as+br)i] (negator<complex>)
890 // a-bi / r-si = hairy and slow anyway; we'll convert to complex
891 // a-bi + r-si = (a+r) - (b+s)i (conjugate)
892 // a-bi - r-si = (a-r) + (s-b)i (complex)
893 
894 #define SimTK_NTRAITS_CONJ_SPEC(T1,T2) \
895 template<> template<> struct NTraits< conjugate<T1> >::Result<T2> { \
896  typedef conjugate<Widest<T1,T2>::Type> W; \
897  typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub; \
898 }; \
899 template<> template<> struct NTraits< conjugate<T1> >::Result<complex<T2> >{\
900  typedef Widest<complex<T1>,complex<T2> >::Type W; \
901  typedef W Mul; typedef W Dvd; typedef W Add; typedef negator<W> Sub; \
902 }; \
903 template<> template<> struct NTraits< conjugate<T1> >::Result<conjugate<T2> >{\
904  typedef Widest<T1,T2>::Type W; typedef complex<W> WC; \
905  typedef negator<WC> Mul; typedef WC Dvd; typedef conjugate<W> Add; typedef WC Sub;\
906 }
907 SimTK_NTRAITS_CONJ_SPEC(float,float);SimTK_NTRAITS_CONJ_SPEC(float,double);
908 SimTK_NTRAITS_CONJ_SPEC(float,long double);
909 SimTK_NTRAITS_CONJ_SPEC(double,float);SimTK_NTRAITS_CONJ_SPEC(double,double);
910 SimTK_NTRAITS_CONJ_SPEC(double,long double);
911 SimTK_NTRAITS_CONJ_SPEC(long double,float);SimTK_NTRAITS_CONJ_SPEC(long double,double);
912 SimTK_NTRAITS_CONJ_SPEC(long double,long double);
913 #undef SimTK_NTRAITS_CONJ_SPEC
914 
915 
916 // Specializations for real numbers.
917 // For real scalar R, op result types are:
918 // Typeof(R*P) = Typeof(P*R)
919 // Typeof(R/P) = Typeof(inv(P)*R)
920 // Typeof(R+P) = Typeof(P+R)
921 // typeof(R-P) = Typeof(P::TNeg + R)
922 // These must be specialized for P=Real and P=Complex.
923 #define SimTK_DEFINE_REAL_NTRAITS(R) \
924 template <> class NTraits<R> { \
925 public: \
926  typedef R T; \
927  typedef negator<T> TNeg; \
928  typedef T TWithoutNegator; \
929  typedef T TReal; \
930  typedef T TImag; \
931  typedef complex<T> TComplex; \
932  typedef T THerm; \
933  typedef T TPosTrans; \
934  typedef T TSqHermT; \
935  typedef T TSqTHerm; \
936  typedef T TElement; \
937  typedef T TRow; \
938  typedef T TCol; \
939  typedef T TSqrt; \
940  typedef T TAbs; \
941  typedef T TStandard; \
942  typedef T TInvert; \
943  typedef T TNormalize; \
944  typedef T Scalar; \
945  typedef T ULessScalar; \
946  typedef T Number; \
947  typedef T StdNumber; \
948  typedef T Precision; \
949  typedef T ScalarNormSq; \
950  template <class P> struct Result { \
951  typedef typename CNT<P>::template Result<R>::Mul Mul; \
952  typedef typename CNT< typename CNT<P>::THerm >::template Result<R>::Mul Dvd; \
953  typedef typename CNT<P>::template Result<R>::Add Add; \
954  typedef typename CNT< typename CNT<P>::TNeg >::template Result<R>::Add Sub; \
955  }; \
956  template <class P> struct Substitute { \
957  typedef P Type; \
958  }; \
959  enum { \
960  NRows = 1, \
961  NCols = 1, \
962  RowSpacing = 1, \
963  ColSpacing = 1, \
964  NPackedElements = 1, \
965  NActualElements = 1, \
966  NActualScalars = 1, \
967  ImagOffset = 0, \
968  RealStrideFactor = 1, \
969  ArgDepth = SCALAR_DEPTH, \
970  IsScalar = 1, \
971  IsULessScalar = 1, \
972  IsNumber = 1, \
973  IsStdNumber = 1, \
974  IsPrecision = 1, \
975  SignInterpretation = 1 \
976  }; \
977  static const T* getData(const T& t) { return &t; } \
978  static T* updData(T& t) { return &t; } \
979  static const T& real(const T& t) { return t; } \
980  static T& real(T& t) { return t; } \
981  static const T& imag(const T&) { return getZero(); } \
982  static T& imag(T&) { assert(false); return *reinterpret_cast<T*>(0); } \
983  static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);} \
984  static TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);} \
985  static const THerm& transpose(const T& t) {return reinterpret_cast<const THerm&>(t);} \
986  static THerm& transpose(T& t) {return reinterpret_cast<THerm&>(t);} \
987  static const TPosTrans& positionalTranspose(const T& t) \
988  {return reinterpret_cast<const TPosTrans&>(t);} \
989  static TPosTrans& positionalTranspose(T& t) \
990  {return reinterpret_cast<TPosTrans&>(t);} \
991  static const TWithoutNegator& castAwayNegatorIfAny(const T& t) \
992  {return reinterpret_cast<const TWithoutNegator&>(t);} \
993  static TWithoutNegator& updCastAwayNegatorIfAny(T& t) \
994  {return reinterpret_cast<TWithoutNegator&>(t);} \
995  static ScalarNormSq scalarNormSqr(const T& t) {return t*t;} \
996  static TSqrt sqrt(const T& t) {return std::sqrt(t);} \
997  static TAbs abs(const T& t) {return std::abs(t);} \
998  static const TStandard& standardize(const T& t) {return t;} \
999  static TNormalize normalize(const T& t) {return (t>0?T(1):(t<0?T(-1):getNaN()));} \
1000  static TInvert invert(const T& t) {return T(1)/t;} \
1001  /* properties of this floating point representation, with memory addresses */ \
1002  static const T& getEps() {return RTraits<T>::getEps();} \
1003  static const T& getSignificant() {return RTraits<T>::getSignificant();} \
1004  static const T& getNaN() {static const T c=std::numeric_limits<T>::quiet_NaN(); return c;} \
1005  static const T& getInfinity() {static const T c=std::numeric_limits<T>::infinity(); return c;} \
1006  static const T& getLeastPositive(){static const T c=std::numeric_limits<T>::min(); return c;} \
1007  static const T& getMostPositive() {static const T c=std::numeric_limits<T>::max(); return c;} \
1008  static const T& getLeastNegative(){static const T c=-std::numeric_limits<T>::min(); return c;} \
1009  static const T& getMostNegative() {static const T c=-std::numeric_limits<T>::max(); return c;} \
1010  static const T& getSqrtEps() {static const T c=std::sqrt(getEps()); return c;} \
1011  static const T& getTiny() {static const T c=std::pow(getEps(), (T)1.25L); return c;} \
1012  static bool isFinite(const T& t) {return SimTK::isFinite(t);} \
1013  static bool isNaN (const T& t) {return SimTK::isNaN(t);} \
1014  static bool isInf (const T& t) {return SimTK::isInf(t);} \
1015  /* Methods to use for approximate comparisons. Perform comparison in the wider of the two */ \
1016  /* precisions, using the default tolerance from the narrower of the two precisions. */ \
1017  static double getDefaultTolerance() {return RTraits<T>::getDefaultTolerance();} \
1018  static bool isNumericallyEqual(const T& t, const float& f) {return SimTK::isNumericallyEqual(t,f);} \
1019  static bool isNumericallyEqual(const T& t, const double& d) {return SimTK::isNumericallyEqual(t,d);} \
1020  static bool isNumericallyEqual(const T& t, const long double& l) {return SimTK::isNumericallyEqual(t,l);} \
1021  static bool isNumericallyEqual(const T& t, int i) {return SimTK::isNumericallyEqual(t,i);} \
1022  /* Here the tolerance is given so we don't have to figure it out. */ \
1023  static bool isNumericallyEqual(const T& t, const float& f, double tol){return SimTK::isNumericallyEqual(t,f,tol);} \
1024  static bool isNumericallyEqual(const T& t, const double& d, double tol){return SimTK::isNumericallyEqual(t,d,tol);} \
1025  static bool isNumericallyEqual(const T& t, const long double& l, double tol){return SimTK::isNumericallyEqual(t,l,tol);} \
1026  static bool isNumericallyEqual(const T& t, int i, double tol){return SimTK::isNumericallyEqual(t,i,tol);} \
1027  /* Carefully calculated constants with convenient memory addresses. */ \
1028  static const T& getZero() {static const T c=(T)(0); return c;} \
1029  static const T& getOne() {static const T c=(T)(1); return c;} \
1030  static const T& getMinusOne() {static const T c=(T)(-1); return c;} \
1031  static const T& getTwo() {static const T c=(T)(2); return c;} \
1032  static const T& getThree() {static const T c=(T)(3); return c;} \
1033  static const T& getOneHalf() {static const T c=(T)(0.5L); return c;} \
1034  static const T& getOneThird() {static const T c=(T)(1.L/3.L); return c;} \
1035  static const T& getOneFourth() {static const T c=(T)(0.25L); return c;} \
1036  static const T& getOneFifth() {static const T c=(T)(0.2L); return c;} \
1037  static const T& getOneSixth() {static const T c=(T)(1.L/6.L); return c;} \
1038  static const T& getOneSeventh() {static const T c=(T)(1.L/7.L); return c;} \
1039  static const T& getOneEighth() {static const T c=(T)(0.125L); return c;} \
1040  static const T& getOneNinth() {static const T c=(T)(1.L/9.L); return c;} \
1041  static const T& getPi() {static const T c=(T)(SimTK_PI); return c;} \
1042  static const T& getOneOverPi() {static const T c=(T)(1.L/SimTK_PI); return c;} \
1043  static const T& getE() {static const T c=(T)(SimTK_E); return c;} \
1044  static const T& getLog2E() {static const T c=(T)(SimTK_LOG2E); return c;} \
1045  static const T& getLog10E() {static const T c=(T)(SimTK_LOG10E); return c;} \
1046  static const T& getSqrt2() {static const T c=(T)(SimTK_SQRT2); return c;} \
1047  static const T& getOneOverSqrt2() {static const T c=(T)(1.L/SimTK_SQRT2); return c;} \
1048  static const T& getSqrt3() {static const T c=(T)(SimTK_SQRT3); return c;} \
1049  static const T& getOneOverSqrt3() {static const T c=(T)(1.L/SimTK_SQRT3); return c;} \
1050  static const T& getCubeRoot2() {static const T c=(T)(SimTK_CBRT2); return c;} \
1051  static const T& getCubeRoot3() {static const T c=(T)(SimTK_CBRT3); return c;} \
1052  static const T& getLn2() {static const T c=(T)(SimTK_LN2); return c;} \
1053  static const T& getLn10() {static const T c=(T)(SimTK_LN10); return c;} \
1054  /* integer digit counts useful for formatted input and output */ \
1055  static int getNumDigits() {static const int c=(int)(std::log10(1/getEps()) -0.5); return c;} \
1056  static int getLosslessNumDigits() {static const int c=(int)(std::log10(1/getTiny())+0.5); return c;} \
1057 }; \
1058 template<> struct NTraits<R>::Result<float> \
1059  {typedef Widest<R,float>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1060 template<> struct NTraits<R>::Result<double> \
1061  {typedef Widest<R,double>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1062 template<> struct NTraits<R>::Result<long double> \
1063  {typedef Widest<R,long double>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1064 template<> struct NTraits<R>::Result<complex<float> > \
1065  {typedef Widest<R,complex<float> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1066 template<> struct NTraits<R>::Result<complex<double> > \
1067  {typedef Widest<R,complex<double> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1068 template<> struct NTraits<R>::Result<complex<long double> > \
1069  {typedef Widest<R,complex<long double> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1070 template<> struct NTraits<R>::Result<conjugate<float> > \
1071  {typedef conjugate<Widest<R,float>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1072 template<> struct NTraits<R>::Result<conjugate<double> > \
1073  {typedef conjugate<Widest<R,double>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
1074 template<> struct NTraits<R>::Result<conjugate<long double> > \
1075  {typedef conjugate<Widest<R,long double>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}
1078 SimTK_DEFINE_REAL_NTRAITS(long double);
1079 #undef SimTK_DEFINE_REAL_NTRAITS
1080 
1082 template <class R> class CNT< complex<R> > : public NTraits< complex<R> > { };
1083 template <class R> class CNT< conjugate<R> > : public NTraits< conjugate<R> > { };
1084 template <> class CNT<float> : public NTraits<float> { };
1085 template <> class CNT<double> : public NTraits<double> { };
1086 template <> class CNT<long double> : public NTraits<long double> { };
1087 
1088 
1089 } // namespace SimTK
1090 
1091 #endif //SimTK_SIMMATRIX_NTRAITS_H_
static TNormalize normalize(const T &t)
Definition: NTraits.h:609
static TImag & imag(T &t)
Definition: NTraits.h:773
static bool isNumericallyEqual(const T &a, int b)
Definition: NTraits.h:647
static const T & getOneFourth()
Definition: NTraits.h:857
CNT< P >::template Result< C >::Mul Mul
Definition: NTraits.h:550
C ULessScalar
Definition: NTraits.h:537
C TStandard
Definition: NTraits.h:532
static const T & getInfinity()
Definition: NTraits.h:616
R TReal
Definition: NTraits.h:708
static const T & getOneEighth()
Definition: NTraits.h:662
double Type
Definition: NTraits.h:143
static bool isNumericallyEqual(const T &a, const double &b)
Definition: NTraits.h:643
static const T & getOneOverSqrt3()
Definition: NTraits.h:672
C TInvert
Definition: NTraits.h:533
double Type
Definition: NTraits.h:140
static const T & getLn10()
Definition: NTraits.h:676
static const T & getOne()
Definition: NTraits.h:851
static const T & getMinusOne()
Definition: NTraits.h:653
static const T & getSqrt3()
Definition: NTraits.h:671
static const T & getOneSixth()
Definition: NTraits.h:660
conjugate< R > TCol
Definition: NTraits.h:717
complex< typename Widest< R1, R2 >::Type > Type
Definition: NTraits.h:122
static const T & getLog10E()
Definition: NTraits.h:867
conjugate< R > TComplex
Definition: NTraits.h:710
static bool isNumericallyEqual(const T &a, const double &b, double tol)
Definition: NTraits.h:644
RTraits is a helper class for NTraits.
Definition: NTraits.h:159
SimTK_DEFINE_REAL_NTRAITS(float)
long double Type
Definition: NTraits.h:109
float Precision
Definition: NTraits.h:142
static const T & getOneSeventh()
Definition: NTraits.h:860
R TSqTHerm
Definition: NTraits.h:714
static const T & getThree()
Definition: NTraits.h:854
static bool isNumericallyEqual(const T &a, const long double &b, double tol)
Definition: NTraits.h:845
static TPosTrans & positionalTranspose(T &t)
Definition: NTraits.h:783
complex< R > TSqrt
Definition: NTraits.h:719
R ScalarNormSq
Definition: NTraits.h:730
static const T & getCubeRoot3()
Definition: NTraits.h:674
long double Type
Definition: NTraits.h:144
R TReal
Definition: NTraits.h:519
CNT< typename CNT< P >::TNeg >::template Result< C >::Add Sub
Definition: NTraits.h:553
static const T & getZero()
Definition: NTraits.h:651
static bool isFinite(const T &t)
Definition: NTraits.h:626
static const T & getOneThird()
Definition: NTraits.h:856
C TRow
Definition: NTraits.h:527
static const THerm & transpose(const T &t)
Definition: NTraits.h:778
static bool isNumericallyEqual(const T &a, const long double &b, double tol)
Definition: NTraits.h:646
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition: Assembler.h:37
SimTK::conjugate<R> should be instantiated only for float, double, long double.
Definition: String.h:45
static const THerm & transpose(const T &t)
Definition: NTraits.h:589
long double Type
Definition: NTraits.h:106
static bool isNumericallyEqual(const T &a, const float &b)
Definition: NTraits.h:641
static const T & getSqrt2()
Definition: NTraits.h:868
double Type
Definition: NTraits.h:105
static bool isNumericallyEqual(const T &a, const complex< R2 > &b, double tol)
Definition: NTraits.h:837
static const T & getOneSixth()
Definition: NTraits.h:859
float Type
Definition: NTraits.h:142
conjugate< R > TInvert
Definition: NTraits.h:722
static const T & getOneNinth()
Definition: NTraits.h:663
static bool isNaN(const T &t)
Definition: NTraits.h:826
static const double & getSignificant()
Definition: NTraits.h:172
complex< R > THerm
Definition: NTraits.h:711
static const T & getOneFifth()
Definition: NTraits.h:659
static const TPosTrans & positionalTranspose(const T &t)
Definition: NTraits.h:592
static const T & getThree()
Definition: NTraits.h:655
static bool isNaN(const T &t)
Definition: NTraits.h:627
static const T & getSqrt2()
Definition: NTraits.h:669
C TNormalize
Definition: NTraits.h:534
conjugate< R > ULessScalar
Definition: NTraits.h:726
static TInvert invert(const T &t)
Definition: NTraits.h:802
static bool isNumericallyEqual(const T &a, int b, double tol)
Definition: NTraits.h:847
static const T & getOneOverSqrt2()
Definition: NTraits.h:869
static const T & getOneNinth()
Definition: NTraits.h:862
static const T & getMinusOne()
Definition: NTraits.h:852
static const T & getOneOverSqrt2()
Definition: NTraits.h:670
static const T & getLog2E()
Definition: NTraits.h:667
CNT< P >::template Result< T >::Mul Mul
Definition: NTraits.h:738
float Precision
Definition: NTraits.h:138
float Type
Definition: NTraits.h:137
static const TNeg & negate(const T &t)
Definition: NTraits.h:775
static const TPosTrans & positionalTranspose(const T &t)
Definition: NTraits.h:781
static const float & getEps()
Attainable accuracy at this precision.
Definition: NTraits.h:163
static const T & getE()
Definition: NTraits.h:865
conjugate< R > T
Definition: NTraits.h:705
C TWithoutNegator
Definition: NTraits.h:517
negator<N>, where N is a number type (real, complex, conjugate), is represented in memory identically...
Definition: String.h:44
C T
Definition: NTraits.h:515
conjugate< R > Number
Definition: NTraits.h:727
Widest< R1, R2 >::Precision Precision
Definition: NTraits.h:123
static const T & getCubeRoot2()
Definition: NTraits.h:673
Widest< R1, R2 >::Precision Precision
Definition: NTraits.h:115
static const T & getNaN()
Definition: NTraits.h:808
CNT< P >::template Result< T >::Add Add
Definition: NTraits.h:740
static const TNeg & negate(const T &t)
Definition: NTraits.h:586
static bool isNumericallyEqual(const T &a, const conjugate< R2 > &b)
Definition: NTraits.h:636
static THerm & transpose(T &t)
Definition: NTraits.h:590
static const T & getLn2()
Definition: NTraits.h:675
static const T & getCubeRoot2()
Definition: NTraits.h:872
C Scalar
Definition: NTraits.h:536
negator< R > TImag
Definition: NTraits.h:709
static TSqrt sqrt(const T &t)
Definition: NTraits.h:793
static const T & getPi()
Definition: NTraits.h:664
static TInvert invert(const T &t)
Definition: NTraits.h:610
static const T & getPi()
Definition: NTraits.h:863
static TReal & real(T &t)
Definition: NTraits.h:771
float Precision
Definition: NTraits.h:136
SimTK_BNTCMPLX_SPEC(float, float)
float Type
Definition: NTraits.h:136
R TImag
Definition: NTraits.h:520
C TComplex
Definition: NTraits.h:521
double Type
Definition: NTraits.h:141
static TAbs abs(const T &t)
Definition: NTraits.h:606
static TNeg & negate(T &t)
Definition: NTraits.h:587
complex< typename Widest< R1, R2 >::Type > Type
Definition: NTraits.h:118
static TSqrt sqrt(const T &t)
Definition: NTraits.h:604
R Precision
Definition: NTraits.h:729
CNT< typename CNT< P >::TNeg >::template Result< T >::Add Sub
Definition: NTraits.h:741
static bool isNumericallyEqual(const T &a, const conjugate< R2 > &b, double tol)
Definition: NTraits.h:833
bool isFinite(const negator< float > &x)
Definition: negator.h:287
R TAbs
Definition: NTraits.h:720
static R & real(T &t)
Definition: NTraits.h:582
conjugate< R > TNormalize
Definition: NTraits.h:723
Definition: CompositeNumericalTypes.h:116
static const R & imag(const T &t)
Definition: NTraits.h:583
static bool isFinite(const T &t)
Definition: NTraits.h:825
Narrowest< R1, R2 >::Precision Precision
Definition: NTraits.h:151
conjugate< R > TWithoutNegator
Definition: NTraits.h:707
CNT< typename CNT< P >::THerm >::template Result< T >::Mul Dvd
Definition: NTraits.h:739
conjugate< R > TPosTrans
Definition: NTraits.h:712
static bool isNumericallyEqual(const T &a, int b)
Definition: NTraits.h:846
static TAbs abs(const T &t)
Definition: NTraits.h:795
float Precision
Definition: NTraits.h:137
double Type
Definition: NTraits.h:107
static TNormalize normalize(const T &t)
Definition: NTraits.h:799
static const T & getTwo()
Definition: NTraits.h:853
static const double & getEps()
Definition: NTraits.h:171
C TCol
Definition: NTraits.h:528
static const TWithoutNegator & castAwayNegatorIfAny(const T &t)
Definition: NTraits.h:786
static TStandard standardize(const T &t)
Definition: NTraits.h:797
static TWithoutNegator & updCastAwayNegatorIfAny(T &t)
Definition: NTraits.h:599
static const T & getTwo()
Definition: NTraits.h:654
float Type
Definition: NTraits.h:139
R TAbs
Definition: NTraits.h:531
This class is specialized for all 36 combinations of standard types (that is, real and complex types ...
Definition: NTraits.h:103
static TWithoutNegator & updCastAwayNegatorIfAny(T &t)
Definition: NTraits.h:788
static bool isNumericallyEqual(const T &a, const conjugate< R2 > &b)
Definition: NTraits.h:831
static double getDefaultTolerance()
Definition: NTraits.h:179
static double getDefaultTolerance()
Definition: NTraits.h:173
High precision mathematical and physical constants.
The purpose of the CNT<T> class is to hide the differences between built-in numerical types and compo...
static ScalarNormSq scalarNormSqr(const T &t)
Definition: NTraits.h:602
static const T & getOneThird()
Definition: NTraits.h:657
static bool isNumericallyEqual(const T &a, const double &b, double tol)
Definition: NTraits.h:843
complex< typename Widest< R1, R2 >::Type > Type
Definition: NTraits.h:114
double Precision
Definition: NTraits.h:140
conjugate< R > TRow
Definition: NTraits.h:716
static const T & getOneSeventh()
Definition: NTraits.h:661
long double Precision
Definition: NTraits.h:110
long double Type
Definition: NTraits.h:111
static bool isNumericallyEqual(const T &a, const float &b, double tol)
Definition: NTraits.h:841
ELEM max(const VectorBase< ELEM > &v)
Definition: VectorMath.h:251
long double Precision
Definition: NTraits.h:106
R ScalarNormSq
Definition: NTraits.h:541
static const TReal & real(const T &t)
Definition: NTraits.h:770
conjugate< R > THerm
Definition: NTraits.h:522
static TNeg & negate(T &t)
Definition: NTraits.h:776
static const T & getOne()
Definition: NTraits.h:652
static double getDefaultTolerance()
Definition: NTraits.h:630
static bool isNumericallyEqual(const T &a, int b, double tol)
Definition: NTraits.h:648
double Precision
Definition: NTraits.h:107
static const TImag & imag(const T &t)
Definition: NTraits.h:772
static const T & getInfinity()
Definition: NTraits.h:814
This class is specialized for all 36 combinations of standard types (that is, real and complex types ...
Definition: NTraits.h:135
static const T & getOneHalf()
Definition: NTraits.h:855
Widest< R1, R2 >::Precision Precision
Definition: NTraits.h:119
static const T & getE()
Definition: NTraits.h:666
static const float & getSignificant()
What multiple of attainable accuracy do we consider significant?
Definition: NTraits.h:165
complex< typename Narrowest< R1, R2 >::Type > Type
Definition: NTraits.h:150
long double Precision
Definition: NTraits.h:112
static const long double & getSignificant()
Definition: NTraits.h:178
static R & imag(T &t)
Definition: NTraits.h:584
RowVectorBase< typename CNT< ELEM >::TAbs > abs(const RowVectorBase< ELEM > &v)
Definition: VectorMath.h:120
bool isNaN(const negator< float > &x)
Definition: negator.h:273
Narrowest< R1, R2 >::Precision Precision
Definition: NTraits.h:147
static const T & getOneFourth()
Definition: NTraits.h:658
static T * updData(T &t)
Definition: NTraits.h:769
Specialized information about Composite Numerical Types which allows us to define appropriate templat...
Definition: CompositeNumericalTypes.h:136
static const T & getNaN()
Definition: NTraits.h:612
static const TStandard & standardize(const T &t)
Definition: NTraits.h:608
conjugate< R > TElement
Definition: NTraits.h:715
SimTK_NTRAITS_CONJ_SPEC(float, float)
static bool isInf(const T &t)
Definition: NTraits.h:628
static const T & getCubeRoot3()
Definition: NTraits.h:873
static double getDefaultTolerance()
Definition: NTraits.h:829
bool isInf(const negator< float > &x)
Definition: negator.h:301
double Type
Definition: NTraits.h:108
static ScalarNormSq scalarNormSqr(const T &t)
Definition: NTraits.h:791
static bool isNumericallyEqual(const T &a, const complex< R2 > &b)
Definition: NTraits.h:835
static const T & getSqrt3()
Definition: NTraits.h:870
double Precision
Definition: NTraits.h:143
static const T & getLog10E()
Definition: NTraits.h:668
static const T & getOneEighth()
Definition: NTraits.h:861
static const T & getLog2E()
Definition: NTraits.h:866
Narrowest< R1, R2 >::Precision Precision
Definition: NTraits.h:155
static T * updData(T &t)
Definition: NTraits.h:580
long double Precision
Definition: NTraits.h:111
C StdNumber
Definition: NTraits.h:539
long double Type
Definition: NTraits.h:112
negator< T > TNeg
Definition: NTraits.h:706
complex< typename Narrowest< R1, R2 >::Type > Type
Definition: NTraits.h:154
CNT< typename CNT< P >::THerm >::template Result< C >::Mul Dvd
Definition: NTraits.h:551
static bool isNumericallyEqual(const T &a, const float &b, double tol)
Definition: NTraits.h:642
R Precision
Definition: NTraits.h:540
static bool isInf(const T &t)
Definition: NTraits.h:827
static const T & getOneOverPi()
Definition: NTraits.h:864
double Precision
Definition: NTraits.h:105
static const T & getOneOverPi()
Definition: NTraits.h:665
conjugate< R > Scalar
Definition: NTraits.h:725
long double Type
Definition: NTraits.h:110
C TPosTrans
Definition: NTraits.h:523
static const T * getData(const T &t)
Definition: NTraits.h:579
long double Precision
Definition: NTraits.h:109
C TSqrt
Definition: NTraits.h:530
static const T & getI()
Definition: NTraits.h:820
static const T & getOneHalf()
Definition: NTraits.h:656
R TSqTHerm
Definition: NTraits.h:525
R TSqHermT
Definition: NTraits.h:713
static bool isNumericallyEqual(const T &a, const complex< R2 > &b)
Definition: NTraits.h:632
static bool isNumericallyEqual(const T &a, const conjugate< R2 > &b, double tol)
Definition: NTraits.h:638
C TElement
Definition: NTraits.h:526
static const T & getLn10()
Definition: NTraits.h:875
double Precision
Definition: NTraits.h:108
static double getDefaultTolerance()
The default numerical error tolerance is always given in double precision.
Definition: NTraits.h:167
static const T * getData(const T &t)
Definition: NTraits.h:768
complex< typename Narrowest< R1, R2 >::Type > Type
Definition: NTraits.h:146
complex< R > TStandard
Definition: NTraits.h:721
static const T & getZero()
Definition: NTraits.h:850
float Type
Definition: NTraits.h:138
static const long double & getEps()
Definition: NTraits.h:177
static bool isNumericallyEqual(const T &a, const double &b)
Definition: NTraits.h:842
static const TWithoutNegator & castAwayNegatorIfAny(const T &t)
Definition: NTraits.h:597
double Precision
Definition: NTraits.h:141
negator< C > TNeg
Definition: NTraits.h:516
float Precision
Definition: NTraits.h:104
static bool isNumericallyEqual(const T &a, const complex< R2 > &b, double tol)
Definition: NTraits.h:634
static TPosTrans & positionalTranspose(T &t)
Definition: NTraits.h:594
static const T & getOneOverSqrt3()
Definition: NTraits.h:871
R TSqHermT
Definition: NTraits.h:524
static const R & real(const T &t)
Definition: NTraits.h:581
CNT< P >::template Result< C >::Add Add
Definition: NTraits.h:552
long double Precision
Definition: NTraits.h:144
float Type
Definition: NTraits.h:104
complex< R > StdNumber
Definition: NTraits.h:728
Definition: negator.h:64
static const T & getLn2()
Definition: NTraits.h:874
static bool isNumericallyEqual(const T &a, const long double &b)
Definition: NTraits.h:844
static const T & getOneFifth()
Definition: NTraits.h:858
static const T & getI()
Definition: NTraits.h:621
float Precision
Definition: NTraits.h:139
bool isNumericallyEqual(const float &a, const float &b, double tol=RTraits< float >::getDefaultTolerance())
Compare two floats for approximate equality.
Definition: NTraits.h:313
C Number
Definition: NTraits.h:538
static bool isNumericallyEqual(const T &a, const long double &b)
Definition: NTraits.h:645
static bool isNumericallyEqual(const T &a, const float &b)
Definition: NTraits.h:840
static THerm & transpose(T &t)
Definition: NTraits.h:779