NTraits.h

Go to the documentation of this file.

#ifndef SimTK_SIMMATRIX_NTRAITS_H_
#define SimTK_SIMMATRIX_NTRAITS_H_
 
/* -------------------------------------------------------------------------- *
 *                       Simbody(tm): SimTKcommon                             *
 * -------------------------------------------------------------------------- *
 * This is part of the SimTK biosimulation toolkit originating from           *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody.  *
 *                                                                            *
 * Portions copyright (c) 2005-12 Stanford University and the Authors.        *
 * Authors: Michael Sherman                                                   *
 * Contributors:                                                              *
 *                                                                            *
 * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
 * not use this file except in compliance with the License. You may obtain a  *
 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
 *                                                                            *
 * Unless required by applicable law or agreed to in writing, software        *
 * distributed under the License is distributed on an "AS IS" BASIS,          *
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
 * See the License for the specific language governing permissions and        *
 * limitations under the License.                                             *
 * -------------------------------------------------------------------------- */
 
#include "SimTKcommon/Constants.h"
#include "SimTKcommon/internal/CompositeNumericalTypes.h"
 
#include <bit>
#include <cassert>
#include <cmath>
#include <complex>
#include <cstddef>
#include <cstdint>
#include <limits>
 
using std::complex;
 
namespace SimTK::detail {
// Consteval square root. Uses `std::sqrt` when `__cpp_lib_constexpr_cmath` is
// available (C++23). Otherwise uses a Newton-Raphson solver, which may be
// off by 1 ULP due to differences in rounding.
consteval double consteval_sqrt(double x) noexcept {
#if defined(__cpp_lib_constexpr_cmath) && __cpp_lib_constexpr_cmath >= 202202L
    return std::sqrt(x);  // In new C++es, just call `std::sqrt`.
#else
    static_assert(std::numeric_limits<double>::is_iec559);
 
    constexpr double max_double = std::numeric_limits<double>::max();
    constexpr double nan = std::numeric_limits<double>::quiet_NaN();
 
    // Trivial/edge cases.
    if (x != x)          return x;    // forward NaNs
    if (x <  0.0)        return nan;  // sqrt(negative_number) -> NaN
    if (x == 0.0)        return x;    // preserve -0 and +0
    if (x >  max_double) return x;    // +inf
 
    // Subnormal scaling.
    double scale = 1.0;
    if (x < std::numeric_limits<double>::min()) {
        x     *= 4503599627370496.0;  // 2^52
        scale  = 1.0 / 67108864.0;    // 2^-26
    }
 
    auto b = std::bit_cast<std::uint64_t>(x);
    b = (b + (1023ULL << 52)) >> 1;
    auto g = std::bit_cast<double>(b);
 
    // Newton-Raphson solver: might be off by 1ULP
    for (;;) {
        const auto next = (g + x / g) * 0.5;
        if (next >= g) break;
        g = next;
    }
 
    return g * scale;
#endif
}
}
 
namespace SimTK {
 
// This is the 3rd type of number, conjugate. It is like complex except that
// the represented value is the conjugate of the value represented by a complex
// number containing the same bit pattern. That is, complex and conjugate both
// contain two real numbers, re and im, with complex(re,im) meaning
// (re + im*i) while conjugate(re,im) means (re - im*i). It is guaranteed that
// our conjugate type and complex type have identical sizes and representations.
// Together, these defininitions and guarantees permit conjugation
// to be done by reinterpretation rather than be computation.
template <class R> class conjugate; // Only defined for float, double
 
// Specializations of this class provide information about Composite Numerical 
// Types in the style of std::numeric_limits<T>. It is specialized for the 
// numeric types but can be invoked on any composite numerical type as well.
template <class T> class CNT;
 
// NTraits provides information like CNT, but for numeric types only.
template <class N> class NTraits; // allowed only for N=<number>
template <class R> class NTraits< complex<R> >;
template <class R> class NTraits< conjugate<R> >;
template <> class NTraits<float>;
template <> class NTraits<double>;
 
// This is an adaptor for numeric types which negates the apparent values. A
// negator<N> has exactly the same internal representation as a numeric value N, 
// but it is to be interpreted has having the negative of the value it would 
// have if interpreted as an N. This permits negation to be done by 
// reinterpretation rather than computation. A full set of arithmetic operators
// are provided involving negator<N>'s and N's. Sometimes we can save an op or
// two this way. For example negator<N>*negator<N> can be performed as an N*N
// since the negations cancel, and we saved two floating point negations.
template <class N> class negator;      // Only defined for numbers
 
template <class R1, class R2> struct Widest {/* Only defined for built-ins. */};
template <> struct Widest<float,float>              {typedef float       Type;  typedef float       Precision;};
template <> struct Widest<float,double>             {typedef double      Type;  typedef double      Precision;};
template <> struct Widest<double,float>             {typedef double      Type;  typedef double      Precision;};
template <> struct Widest<double,double>            {typedef double      Type;  typedef double      Precision;};
template <class R1, class R2> struct Widest< complex<R1>,complex<R2> > { 
    typedef complex< typename Widest<R1,R2>::Type > Type; 
    typedef typename Widest<R1,R2>::Precision       Precision; 
};
template <class R1, class R2> struct Widest< complex<R1>,R2 > { 
    typedef complex< typename Widest<R1,R2>::Type > Type; 
    typedef typename Widest<R1,R2>::Precision       Precision;  
};
template <class R1, class R2> struct Widest< R1,complex<R2> > { 
    typedef complex< typename Widest<R1,R2>::Type > Type; 
    typedef typename Widest<R1,R2>::Precision       Precision; 
};
 
template <class R1, class R2> struct Narrowest {/* Only defined for built-ins. */};
template <> struct Narrowest<float,float>              {typedef float  Type; typedef float Precision;};
template <> struct Narrowest<float,double>             {typedef float  Type; typedef float Precision;};
template <> struct Narrowest<double,float>             {typedef float  Type; typedef float Precision;};
template <> struct Narrowest<double,double>            {typedef double Type; typedef double Precision;};
template <class R1, class R2> struct Narrowest< complex<R1>,complex<R2> > { 
    typedef complex< typename Narrowest<R1,R2>::Type >  Type; 
    typedef typename Narrowest<R1,R2>::Precision        Precision;
};
template <class R1, class R2> struct Narrowest< complex<R1>,R2 > { 
    typedef complex< typename Narrowest<R1,R2>::Type >  Type; 
    typedef typename Narrowest<R1,R2>::Precision        Precision; 
};
template <class R1, class R2> struct Narrowest< R1,complex<R2> > { 
    typedef complex< typename Narrowest<R1,R2>::Type >  Type; 
    typedef typename Narrowest<R1,R2>::Precision        Precision; 
};
 
template <class R> class RTraits {/* Only defined for real types */};
template <> class RTraits<float> {
public:
    static constexpr float getEps()               {return std::numeric_limits<float>::epsilon();}
    static constexpr float getSignificant() {
        // Three nested square roots compute x / x^(1/8) = x^(7/8) = x^0.875
        return getEps() / detail::consteval_sqrt(detail::consteval_sqrt(detail::consteval_sqrt(getEps())));
    }
    static constexpr double getDefaultTolerance() {return static_cast<double>(getSignificant());}
};
template <> class RTraits<double> {
public:
    static constexpr double getEps()              {return std::numeric_limits<double>::epsilon();}
    static constexpr double getSignificant() {
        // Three nested square roots compute x / x^(1/8) = x^(7/8) = x^0.875
        return getEps() / detail::consteval_sqrt(detail::consteval_sqrt(detail::consteval_sqrt(getEps())));
    }
    static constexpr double getDefaultTolerance() {return getSignificant();}
};
 
// See negator.h for isNaN() applied to negated scalars.
inline bool isNaN(const float& x)  {return std::isnan(x);}
inline bool isNaN(const double& x) {return std::isnan(x);}
 
template <class P> inline bool
isNaN(const std::complex<P>& x)
{   return isNaN(x.real()) || isNaN(x.imag());}
 
template <class P> inline bool
isNaN(const conjugate<P>& x)
{   return isNaN(x.real()) || isNaN(x.negImag());}
 
// See negator.h for isFinite() applied to negated scalars.
inline bool isFinite(const float&  x) {return std::isfinite(x);}
inline bool isFinite(const double& x) {return std::isfinite(x);}
 
template <class P> inline bool
isFinite(const std::complex<P>& x)
{   return isFinite(x.real()) && isFinite(x.imag());}
 
template <class P> inline bool
isFinite(const conjugate<P>& x)
{   return isFinite(x.real()) && isFinite(x.negImag());}
 
// See negator.h for isInf() applied to negated scalars.
inline bool isInf(const float&  x) {return std::isinf(x);}
inline bool isInf(const double& x) {return std::isinf(x);}
 
template <class P> inline bool
isInf(const std::complex<P>& x) {
    return (isInf(x.real()) && !isNaN(x.imag()))
        || (isInf(x.imag()) && !isNaN(x.real()));
}
 
template <class P> inline bool
isInf(const conjugate<P>& x) {
    return (isInf(x.real())    && !isNaN(x.negImag()))
        || (isInf(x.negImag()) && !isNaN(x.real()));
}
 
inline bool isNumericallyEqual(const float& a, const float& b, 
                               double tol = RTraits<float>::getDefaultTolerance())
{   if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false;
    const float scale = std::max(std::max(std::abs(a),std::abs(b)), 1.f);
    return std::abs(a-b) <= scale*(float)tol; }
inline bool isNumericallyEqual(const double& a, const double& b, 
                               double tol = RTraits<double>::getDefaultTolerance())
{   if (isNaN(a)) return isNaN(b); else if (isNaN(b)) return false;
    const double scale = std::max(std::max(std::abs(a),std::abs(b)), 1.);
    return std::abs(a-b) <= scale*tol; }
 
inline bool isNumericallyEqual(const float& a, const double& b, 
                               double tol = RTraits<float>::getDefaultTolerance())
{   return isNumericallyEqual((double)a, b, tol); }
inline bool isNumericallyEqual(const double& a, const float& b, 
                               double tol = RTraits<float>::getDefaultTolerance())
{   return isNumericallyEqual(a, (double)b, tol); }
 
inline bool isNumericallyEqual(const float& a, int b,
                               double tol = RTraits<float>::getDefaultTolerance())
{   return isNumericallyEqual(a, (double)b, tol); }
inline bool isNumericallyEqual(int a, const float& b,
                               double tol = RTraits<float>::getDefaultTolerance())
{   return isNumericallyEqual((double)a, b, tol); }
inline bool isNumericallyEqual(const double& a, int b,
                               double tol = RTraits<double>::getDefaultTolerance())
{   return isNumericallyEqual(a, (double)b, tol); }
inline bool isNumericallyEqual(int a, const double& b,
                               double tol = RTraits<double>::getDefaultTolerance())
{   return isNumericallyEqual((double)a, b, tol); }
 
template <class P, class Q>
inline bool isNumericallyEqual
  ( const std::complex<P>& a, const std::complex<Q>& b, 
    double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance())
{   return isNumericallyEqual(a.real(),b.real(),tol)
        && isNumericallyEqual(a.imag(),b.imag(),tol); }
 
template <class P, class Q>
inline bool isNumericallyEqual
  ( const conjugate<P>& a, const conjugate<Q>& b, 
    double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance())
{   return isNumericallyEqual(a.real(),b.real(),tol)
        && isNumericallyEqual(a.imag(),b.imag(),tol); }
 
template <class P, class Q>
inline bool isNumericallyEqual
  ( const std::complex<P>& a, const conjugate<Q>& b, 
    double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance())
{   return isNumericallyEqual(a.real(),b.real(),tol)
        && isNumericallyEqual(a.imag(),b.imag(),tol); }
 
template <class P, class Q>
inline bool isNumericallyEqual
  ( const conjugate<P>& a, const std::complex<Q>& b, 
    double tol = RTraits<typename Narrowest<P,Q>::Precision>::getDefaultTolerance())
{   return isNumericallyEqual(a.real(),b.real(),tol)
        && isNumericallyEqual(a.imag(),b.imag(),tol); }
 
template <class P> inline bool 
isNumericallyEqual(const std::complex<P>& a, const float& b, 
                   double tol = RTraits<float>::getDefaultTolerance())
{   return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.f,tol); }
template <class P> inline bool 
isNumericallyEqual(const float& a, const std::complex<P>& b,
                   double tol = RTraits<float>::getDefaultTolerance())
{   return isNumericallyEqual(b,a,tol); }
template <class P> inline bool 
isNumericallyEqual(const std::complex<P>& a, const double& b, 
                   double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance())
{   return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.,tol); }
template <class P> inline bool 
isNumericallyEqual(const double& a, const std::complex<P>& b,
                   double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance())
{   return isNumericallyEqual(b,a,tol); }
template <class P> inline bool 
isNumericallyEqual(const std::complex<P>& a, int b, 
                   double tol = RTraits<P>::getDefaultTolerance())
{   typedef typename Widest<P,double>::Precision W; return isNumericallyEqual(a,(W)b,tol); }
template <class P> inline bool 
isNumericallyEqual(int a, const std::complex<P>& b, 
                   double tol = RTraits<P>::getDefaultTolerance())
{   return isNumericallyEqual(b,a,tol); }
 
template <class P> inline bool 
isNumericallyEqual(const conjugate<P>& a, const float& b, 
                   double tol = RTraits<float>::getDefaultTolerance())
{   return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.f,tol); }
template <class P> inline bool 
isNumericallyEqual(const float& a, const conjugate<P>& b,
                   double tol = RTraits<float>::getDefaultTolerance())
{   return isNumericallyEqual(b,a,tol); }
template <class P> inline bool 
isNumericallyEqual(const conjugate<P>& a, const double& b, 
                   double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance())
{   return isNumericallyEqual(a.real(),b,tol) && isNumericallyEqual(a.imag(),0.,tol); }
template <class P> inline bool 
isNumericallyEqual(const double& a, const conjugate<P>& b,
                   double tol = RTraits<typename Narrowest<P,double>::Precision>::getDefaultTolerance())
{   return isNumericallyEqual(b,a,tol); }
template <class P> inline bool 
isNumericallyEqual(const conjugate<P>& a, int b, 
                   double tol = RTraits<P>::getDefaultTolerance())
{   typedef typename Widest<P,double>::Precision W; return isNumericallyEqual(a,(W)b,tol); }
template <class P> inline bool 
isNumericallyEqual(int a, const conjugate<P>& b, 
                   double tol = RTraits<P>::getDefaultTolerance())
{   return isNumericallyEqual(b,a,tol); }
 
 
 
template <class N> class NTraits { 
    // only the specializations below are allowed 
};
 
template <class R> class NTraits< complex<R> > {
    typedef complex<R>  C;
public:
    typedef C                T;
    typedef negator<C>       TNeg;            // type of this after *cast* to its negative
    typedef C                TWithoutNegator; // type of this ignoring negator (there isn't one!)
 
    typedef R                TReal;
    typedef R                TImag;
    typedef C                TComplex;    
    typedef conjugate<R>     THerm;     // this is a recast
    typedef C                TPosTrans;
    typedef R                TSqHermT;  // ~C*C
    typedef R                TSqTHerm;  // C*~C (same)
    typedef C                TElement;
    typedef C                TRow;
    typedef C                TCol;
 
    typedef C                TSqrt;
    typedef R                TAbs;
    typedef C                TStandard; // complex is a standard type
    typedef C                TInvert;   // this is a calculation, so use standard number
    typedef C                TNormalize;
 
    typedef C                Scalar;
    typedef C                ULessScalar;
    typedef C                Number;
    typedef C                StdNumber;
    typedef R                Precision;
    typedef R                ScalarNormSq;
 
    // For complex scalar C, op result types are:
    //   Typeof(C*P) = Typeof(P*C)
    //   Typeof(C/P) = Typeof(inv(P)*C)
    //   Typeof(C+P) = Typeof(P+C)
    //   typeof(C-P) = Typeof(P::TNeg + C)
    // These must be specialized for P=real, complex, conjugate.
    template <class P> struct Result { 
        typedef typename CNT<P>::template Result<C>::Mul Mul;
        typedef typename CNT< typename CNT<P>::THerm >::template Result<C>::Mul Dvd;
        typedef typename CNT<P>::template Result<C>::Add Add;
        typedef typename CNT< typename CNT<P>::TNeg >::template Result<C>::Add Sub;
    };
 
    // Shape-preserving element substitution (easy for scalars!)
    template <class P> struct Substitute {
        typedef P Type;
    };
 
    enum {
        NRows               = 1,
        NCols               = 1,
        RowSpacing          = 1,
        ColSpacing          = 1,
        NPackedElements     = 1,      // not two!
        NActualElements     = 1,
        NActualScalars      = 1,
        ImagOffset          = 1,
        RealStrideFactor    = 2,      // double stride when casting to real or imaginary
        ArgDepth            = SCALAR_DEPTH,
        IsScalar            = 1,
        IsULessScalar       = 1,
        IsNumber            = 1,
        IsStdNumber         = 1,
        IsPrecision         = 0,
        SignInterpretation  = 1       // after cast to Number, what is the sign?
    }; 
    static const T* getData(const T& t) { return &t; } 
    static T*       updData(T& t)       { return &t; }
    static const R& real(const T& t) { return (reinterpret_cast<const R*>(&t))[0]; }
    static R&       real(T& t)       { return (reinterpret_cast<R*>(&t))[0]; }
    static const R& imag(const T& t) { return (reinterpret_cast<const R*>(&t))[1]; }
    static R&       imag(T& t)       { return (reinterpret_cast<R*>(&t))[1]; }
 
    static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);}
    static       TNeg& negate(T& t)       {return reinterpret_cast<TNeg&>(t);}
 
    static const THerm& transpose(const T& t) {return reinterpret_cast<const THerm&>(t);}
    static       THerm& transpose(T& t)       {return reinterpret_cast<THerm&>(t);}
 
    static const TPosTrans& positionalTranspose(const T& t)
        {return reinterpret_cast<const TPosTrans&>(t);}
    static       TPosTrans& positionalTranspose(T& t)
        {return reinterpret_cast<TPosTrans&>(t);} 
 
    static const TWithoutNegator& castAwayNegatorIfAny(const T& t)
        {return reinterpret_cast<const TWithoutNegator&>(t);}
    static       TWithoutNegator& updCastAwayNegatorIfAny(T& t)
        {return reinterpret_cast<TWithoutNegator&>(t);}
 
    static ScalarNormSq scalarNormSqr(const T& t)
        { return t.real()*t.real() + t.imag()*t.imag(); }
    static TSqrt    sqrt(const T& t)
        { return std::sqrt(t); }
    static TAbs     abs(const T& t)
        { return std::abs(t); } // no, not just sqrt of scalarNormSqr()!
    static const TStandard& standardize(const T& t) {return t;} // already standard
    static TNormalize normalize(const T& t) {return t/abs(t);}
    static TInvert    invert(const T& t)    {return TReal(1)/t;}
 
    static constexpr T getNaN() {return T(NTraits<R>::getNaN(), NTraits<R>::getNaN());}
    static constexpr T getInfinity() {return T(NTraits<R>::getInfinity(),NTraits<R>::getInfinity());}
 
    static constexpr T getI() {return T(0,1);}
 
    static bool isFinite(const T& t) {return SimTK::isFinite(t);}
    static bool isNaN(const T& t) {return SimTK::isNaN(t);}
    static bool isInf(const T& t) {return SimTK::isInf(t);}
 
    static double getDefaultTolerance() {return RTraits<R>::getDefaultTolerance();}
 
    template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b)
    {   return SimTK::isNumericallyEqual(a,b); }
    template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b, double tol)
    {   return SimTK::isNumericallyEqual(a,b,tol); }
    template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b)
    {   return SimTK::isNumericallyEqual(a,b); }
    template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b, double tol)
    {   return SimTK::isNumericallyEqual(a,b,tol); }
 
    static bool isNumericallyEqual(const T& a, const float& b) {return SimTK::isNumericallyEqual(a,b);}
    static bool isNumericallyEqual(const T& a, const float& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
    static bool isNumericallyEqual(const T& a, const double& b) {return SimTK::isNumericallyEqual(a,b);}
    static bool isNumericallyEqual(const T& a, const double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
    static bool isNumericallyEqual(const T& a, int b) {return SimTK::isNumericallyEqual(a,b);}
    static bool isNumericallyEqual(const T& a, int b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
 
    // The rest are the same as the real equivalents, with zero imaginary part.              
    static constexpr T getZero()         {return T(NTraits<R>::getZero());         }
    static constexpr T getOne()          {return T(NTraits<R>::getOne());          }
    static constexpr T getMinusOne()     {return T(NTraits<R>::getMinusOne());     }
    static constexpr T getTwo()          {return T(NTraits<R>::getTwo());          }
    static constexpr T getThree()        {return T(NTraits<R>::getThree());        }
    static constexpr T getOneHalf()      {return T(NTraits<R>::getOneHalf());      }
    static constexpr T getOneThird()     {return T(NTraits<R>::getOneThird());     }
    static constexpr T getOneFourth()    {return T(NTraits<R>::getOneFourth());    }
    static constexpr T getOneFifth()     {return T(NTraits<R>::getOneFifth());     }
    static constexpr T getOneSixth()     {return T(NTraits<R>::getOneSixth());     }
    static constexpr T getOneSeventh()   {return T(NTraits<R>::getOneSeventh());   }
    static constexpr T getOneEighth()    {return T(NTraits<R>::getOneEighth());    }
    static constexpr T getOneNinth()     {return T(NTraits<R>::getOneNinth());     }
    static constexpr T getPi()           {return T(NTraits<R>::getPi());           }
    static constexpr T getOneOverPi()    {return T(NTraits<R>::getOneOverPi());    }
    static constexpr T getE()            {return T(NTraits<R>::getE());            }
    static constexpr T getLog2E()        {return T(NTraits<R>::getLog2E());        }
    static constexpr T getLog10E()       {return T(NTraits<R>::getLog10E());       }
    static constexpr T getSqrt2()        {return T(NTraits<R>::getSqrt2());        }
    static constexpr T getOneOverSqrt2() {return T(NTraits<R>::getOneOverSqrt2()); }
    static constexpr T getSqrt3()        {return T(NTraits<R>::getSqrt3());        }
    static constexpr T getOneOverSqrt3() {return T(NTraits<R>::getOneOverSqrt3()); }
    static constexpr T getCubeRoot2()    {return T(NTraits<R>::getCubeRoot2());    }
    static constexpr T getCubeRoot3()    {return T(NTraits<R>::getCubeRoot3());    }
    static constexpr T getLn2()          {return T(NTraits<R>::getLn2());          }
    static constexpr T getLn10()         {return T(NTraits<R>::getLn10());         }
};
 
 
// Specialize NTraits<complex>::Results for <complex> OP <scalar>. Result type is
// always just the complex type of sufficient precision.
#define SimTK_BNTCMPLX_SPEC(T1,T2)  \
template<> template<> struct NTraits< complex<T1> >::Result<T2> {      \
    typedef Widest< complex<T1>,T2 >::Type W;                      \
    typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub;         \
};                                                                      \
template<> template<> struct NTraits< complex<T1> >::Result< complex<T2> > {  \
    typedef Widest< complex<T1>,complex<T2> >::Type W;        \
    typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub;         \
};                                                                      \
template<> template<> struct NTraits< complex<T1> >::Result< conjugate<T2> > {  \
    typedef Widest< complex<T1>,complex<T2> >::Type W;        \
    typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub;         \
}
SimTK_BNTCMPLX_SPEC(float,float);SimTK_BNTCMPLX_SPEC(float,double);
SimTK_BNTCMPLX_SPEC(double,float);SimTK_BNTCMPLX_SPEC(double,double);
#undef SimTK_BNTCMPLX_SPEC
 
 
// conjugate -- should be instantiated only for float, double.
template <class R> class NTraits< conjugate<R> > {
    typedef complex<R>          C;
public:
    typedef conjugate<R>        T;
    typedef negator<T>          TNeg;            // type of this after *cast* to its negative
    typedef conjugate<R>        TWithoutNegator; // type of this ignoring negator (there isn't one!)
    typedef R                   TReal;
    typedef negator<R>          TImag;
    typedef conjugate<R>        TComplex;     
    typedef complex<R>          THerm;      // conjugate evaporates
    typedef conjugate<R>        TPosTrans;  // Positional transpose of scalar does nothing
    typedef R                   TSqHermT;   // C*C~
    typedef R                   TSqTHerm;   // ~C*C (same)
    typedef conjugate<R>        TElement;
    typedef conjugate<R>        TRow;
    typedef conjugate<R>        TCol;
 
    typedef complex<R>          TSqrt;
    typedef R                   TAbs;
    typedef complex<R>          TStandard;
    typedef conjugate<R>        TInvert;
    typedef conjugate<R>        TNormalize;
 
    typedef conjugate<R>        Scalar;
    typedef conjugate<R>        ULessScalar;
    typedef conjugate<R>        Number;
    typedef complex<R>          StdNumber;
    typedef R                   Precision;
    typedef R                   ScalarNormSq;
 
    // Typeof( Conj<S>*P ) is Typeof(P*Conj<S>)
    // Typeof( Conj<S>/P ) is Typeof(inv(P)*Conj<S>)
    // Typeof( Conj<S>+P ) is Typeof(P+Conj<S>)
    // Typeof( Conj<S>-P ) is Typeof(P::TNeg+Conj<S>)
    // Must specialize for P=real or P=complex or P=conjugate
    template <class P> struct Result {
        typedef typename CNT<P>::template Result<T>::Mul Mul;
        typedef typename CNT<typename CNT<P>::THerm>::template Result<T>::Mul Dvd;
        typedef typename CNT<P>::template Result<T>::Add Add;
        typedef typename CNT<typename CNT<P>::TNeg>::template Result<T>::Add Sub;
    };
 
    // Shape-preserving element substitution (easy for scalars!)
    template <class P> struct Substitute {
        typedef P Type;
    };
 
    enum {
        NRows               = 1,
        NCols               = 1,
        RowSpacing          = 1,
        ColSpacing          = 1,
        NPackedElements     = 1,      // not two!
        NActualElements     = 1,
        NActualScalars      = 1,
        ImagOffset          = 1,
        RealStrideFactor    = 2,      // double stride when casting to real or imaginary
        ArgDepth            = SCALAR_DEPTH,
        IsScalar            = 1,
        IsULessScalar       = 1,
        IsNumber            = 1,
        IsStdNumber         = 0,
        IsPrecision         = 0,
        SignInterpretation  = 1       // after cast to Number, what is the sign?
    }; 
 
    static const T*     getData(const T& t) { return &t; } 
    static T*           updData(T& t)       { return &t; }
    static const TReal& real(const T& t) { return t.real(); }
    static TReal&       real(T& t)       { return t.real(); }
    static const TImag& imag(const T& t) { return t.imag(); }
    static TImag&       imag(T& t)       { return t.imag(); }
 
    static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);}
    static       TNeg& negate(T& t)       {return reinterpret_cast<TNeg&>(t);}
 
    static const THerm& transpose(const T& t) {return t.conj();}
    static       THerm& transpose(T& t)       {return t.conj();}
 
    static const TPosTrans& positionalTranspose(const T& t)
        {return reinterpret_cast<const TPosTrans&>(t);}
    static       TPosTrans& positionalTranspose(T& t)
        {return reinterpret_cast<TPosTrans&>(t);} 
 
    static const TWithoutNegator& castAwayNegatorIfAny(const T& t)
        {return reinterpret_cast<const TWithoutNegator&>(t);}
    static       TWithoutNegator& updCastAwayNegatorIfAny(T& t)
        {return reinterpret_cast<TWithoutNegator&>(t);}
 
    static ScalarNormSq scalarNormSqr(const T& t)
        { return t.real()*t.real() + t.negImag()*t.negImag(); }
    static TSqrt    sqrt(const T& t)
        { return std::sqrt(C(t)); } // cast to complex (one negation)
    static TAbs     abs(const T& t)
        { return std::abs(t.conj()); }  // no, not just sqrt of scalarNormSqr()!
    static TStandard standardize(const T& t)
        { return TStandard(t); }        // i.e., convert to complex
    static TNormalize normalize(const T& t) {return TNormalize(t/abs(t));}
 
    // 1/conj(z) = conj(1/z), for complex z.
    static TInvert invert(const T& t)    
    {   const typename NTraits<THerm>::TInvert cmplx(NTraits<THerm>::invert(t.conj()));
        return reinterpret_cast<const TInvert&>(cmplx); } // recast complex to conjugate it
 
    // We want a "conjugate NaN", NaN - NaN*i, meaning both reals should
    // be positive NaN.
    static constexpr T getNaN() {return T(NTraits<R>::getNaN(),NTraits<R>::getNaN());}
    // We want a "conjugate infinity", Inf - Inf*i, meaning both stored reals
    // are positive Inf.
    static constexpr T getInfinity() {return T(NTraits<R>::getInfinity(),NTraits<R>::getInfinity());}
    // But we want the constant i (=sqrt(-1)) to be the same however we represent it,
    // so for conjugate i = 0 - (-1)i.
    static constexpr T getI() {return T(0,-1);}
 
    static bool isFinite(const T& t) {return SimTK::isFinite(t);}
    static bool isNaN(const T& t) {return SimTK::isNaN(t);}
    static bool isInf(const T& t) {return SimTK::isInf(t);}
 
    static double getDefaultTolerance() {return RTraits<R>::getDefaultTolerance();}
 
    template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b)
    {   return SimTK::isNumericallyEqual(a,b); }
    template <class R2> static bool isNumericallyEqual(const T& a, const conjugate<R2>& b, double tol)
    {   return SimTK::isNumericallyEqual(a,b,tol); }
    template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b)
    {   return SimTK::isNumericallyEqual(a,b); }
    template <class R2> static bool isNumericallyEqual(const T& a, const complex<R2>& b, double tol)
    {   return SimTK::isNumericallyEqual(a,b,tol); }
 
    static bool isNumericallyEqual(const T& a, const float& b) {return SimTK::isNumericallyEqual(a,b);}
    static bool isNumericallyEqual(const T& a, const float& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
    static bool isNumericallyEqual(const T& a, const double& b) {return SimTK::isNumericallyEqual(a,b);}
    static bool isNumericallyEqual(const T& a, const double& b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
    static bool isNumericallyEqual(const T& a, int b) {return SimTK::isNumericallyEqual(a,b);}
    static bool isNumericallyEqual(const T& a, int b, double tol) {return SimTK::isNumericallyEqual(a,b,tol);}
 
    // The rest are the same as the real equivalents, with zero imaginary part.              
    static constexpr T getZero()         {return T(NTraits<R>::getZero());         }
    static constexpr T getOne()          {return T(NTraits<R>::getOne());          }
    static constexpr T getMinusOne()     {return T(NTraits<R>::getMinusOne());     }
    static constexpr T getTwo()          {return T(NTraits<R>::getTwo());          }
    static constexpr T getThree()        {return T(NTraits<R>::getThree());        }
    static constexpr T getOneHalf()      {return T(NTraits<R>::getOneHalf());      }
    static constexpr T getOneThird()     {return T(NTraits<R>::getOneThird());     }
    static constexpr T getOneFourth()    {return T(NTraits<R>::getOneFourth());    }
    static constexpr T getOneFifth()     {return T(NTraits<R>::getOneFifth());     }
    static constexpr T getOneSixth()     {return T(NTraits<R>::getOneSixth());     }
    static constexpr T getOneSeventh()   {return T(NTraits<R>::getOneSeventh());   }
    static constexpr T getOneEighth()    {return T(NTraits<R>::getOneEighth());    }
    static constexpr T getOneNinth()     {return T(NTraits<R>::getOneNinth());     }
    static constexpr T getPi()           {return T(NTraits<R>::getPi());           }
    static constexpr T getOneOverPi()    {return T(NTraits<R>::getOneOverPi());    }
    static constexpr T getE()            {return T(NTraits<R>::getE());            }
    static constexpr T getLog2E()        {return T(NTraits<R>::getLog2E());        }
    static constexpr T getLog10E()       {return T(NTraits<R>::getLog10E());       }
    static constexpr T getSqrt2()        {return T(NTraits<R>::getSqrt2());        }
    static constexpr T getOneOverSqrt2() {return T(NTraits<R>::getOneOverSqrt2()); }
    static constexpr T getSqrt3()        {return T(NTraits<R>::getSqrt3());        }
    static constexpr T getOneOverSqrt3() {return T(NTraits<R>::getOneOverSqrt3()); }
    static constexpr T getCubeRoot2()    {return T(NTraits<R>::getCubeRoot2());    }
    static constexpr T getCubeRoot3()    {return T(NTraits<R>::getCubeRoot3());    }
    static constexpr T getLn2()          {return T(NTraits<R>::getLn2());          }
    static constexpr T getLn10()         {return T(NTraits<R>::getLn10());         }
};
 
// Any op involving conjugate & a real is best left as a conjugate. However,
// an op involving conjugate & a complex or conjugate can lose the conjugate at zero cost
// and return just a complex in some cases. Also, we prefer negator<complex> to conjugate.
//
// Conj op complex
//   a-bi * r+si = (ar+bs) + (as-br)i               (complex)
//   a-bi / r+si = hairy and slow anyway; we'll convert to complex
//   a-bi + r+si = (a+r) + (s-b)i                   (complex)
//   a-bi - r+si = (a-r) - (b+s)i = -[(r-a)+(b+s)i] (negator<complex>)
//
// Conj op Conj
//   a-bi * r-si = (ar-bs) - (as+br)i = -[(bs-ar)+(as+br)i] (negator<complex>)
//   a-bi / r-si = hairy and slow anyway; we'll convert to complex
//   a-bi + r-si = (a+r) - (b+s)i  (conjugate)
//   a-bi - r-si = (a-r) + (s-b)i  (complex)
 
#define SimTK_NTRAITS_CONJ_SPEC(T1,T2)                                      \
template<> template<> struct NTraits< conjugate<T1> >::Result<T2> {         \
  typedef conjugate<Widest<T1,T2>::Type> W;                                 \
  typedef W Mul; typedef W Dvd; typedef W Add; typedef W Sub;               \
};                                                                          \
template<> template<> struct NTraits< conjugate<T1> >::Result<complex<T2> >{\
  typedef Widest<complex<T1>,complex<T2> >::Type W;               \
  typedef W Mul; typedef W Dvd; typedef W Add; typedef negator<W> Sub;      \
};                                                                          \
template<> template<> struct NTraits< conjugate<T1> >::Result<conjugate<T2> >{\
    typedef Widest<T1,T2>::Type W; typedef complex<W> WC;              \
    typedef negator<WC> Mul; typedef WC Dvd; typedef conjugate<W> Add; typedef WC Sub;\
}
SimTK_NTRAITS_CONJ_SPEC(float,float);SimTK_NTRAITS_CONJ_SPEC(float,double);
SimTK_NTRAITS_CONJ_SPEC(double,float);SimTK_NTRAITS_CONJ_SPEC(double,double);
#undef SimTK_NTRAITS_CONJ_SPEC
 
// Specializations for real numbers.
// For real scalar R, op result types are:
//   Typeof(R*P) = Typeof(P*R)
//   Typeof(R/P) = Typeof(inv(P)*R)
//   Typeof(R+P) = Typeof(P+R)
//   typeof(R-P) = Typeof(P::TNeg + R)
// These must be specialized for P=Real and P=Complex.
#define SimTK_DEFINE_REAL_NTRAITS(R)            \
template <> class NTraits<R> {                  \
public:                                         \
    typedef R                T;                 \
    typedef negator<T>       TNeg;              \
    typedef T                TWithoutNegator;   \
    typedef T                TReal;             \
    typedef T                TImag;             \
    typedef complex<T>       TComplex;          \
    typedef T                THerm;             \
    typedef T                TPosTrans;         \
    typedef T                TSqHermT;          \
    typedef T                TSqTHerm;          \
    typedef T                TElement;          \
    typedef T                TRow;              \
    typedef T                TCol;              \
    typedef T                TSqrt;             \
    typedef T                TAbs;              \
    typedef T                TStandard;         \
    typedef T                TInvert;           \
    typedef T                TNormalize;        \
    typedef T                Scalar;            \
    typedef T                ULessScalar;       \
    typedef T                Number;            \
    typedef T                StdNumber;         \
    typedef T                Precision;         \
    typedef T                ScalarNormSq;      \
    template <class P> struct Result {          \
        typedef typename CNT<P>::template Result<R>::Mul Mul;   \
        typedef typename CNT< typename CNT<P>::THerm >::template Result<R>::Mul Dvd;    \
        typedef typename CNT<P>::template Result<R>::Add Add;   \
        typedef typename CNT< typename CNT<P>::TNeg >::template Result<R>::Add Sub;     \
    };                                          \
    template <class P> struct Substitute {      \
        typedef P Type;                         \
    };                                          \
    enum {                                      \
        NRows               = 1,                \
        NCols               = 1,                \
        RowSpacing          = 1,                \
        ColSpacing          = 1,                \
        NPackedElements     = 1,                \
        NActualElements     = 1,                \
        NActualScalars      = 1,                \
        ImagOffset          = 0,                \
        RealStrideFactor    = 1,                \
        ArgDepth            = SCALAR_DEPTH,     \
        IsScalar            = 1,                \
        IsULessScalar       = 1,                \
        IsNumber            = 1,                \
        IsStdNumber         = 1,                \
        IsPrecision         = 1,                \
        SignInterpretation  = 1                 \
    };                                          \
    static const T* getData(const T& t) { return &t; }  \
    static T*       updData(T& t)       { return &t; }  \
    static const T& real(const T& t) { return t; }      \
    static T&       real(T& t)       { return t; }      \
    static const T& imag(const T&)   { static const T v{};                return v; }     \
    static T&       imag(T&)         { static T v{};       assert(false); return v; }     \
    static const TNeg& negate(const T& t) {return reinterpret_cast<const TNeg&>(t);}      \
    static       TNeg& negate(T& t) {return reinterpret_cast<TNeg&>(t);}                  \
    static const THerm& transpose(const T& t) {return reinterpret_cast<const THerm&>(t);} \
    static       THerm& transpose(T& t) {return reinterpret_cast<THerm&>(t);}             \
    static const TPosTrans& positionalTranspose(const T& t)                               \
        {return reinterpret_cast<const TPosTrans&>(t);}                                   \
    static       TPosTrans& positionalTranspose(T& t)                                     \
        {return reinterpret_cast<TPosTrans&>(t);}                                         \
    static const TWithoutNegator& castAwayNegatorIfAny(const T& t)                        \
        {return reinterpret_cast<const TWithoutNegator&>(t);}                             \
    static       TWithoutNegator& updCastAwayNegatorIfAny(T& t)                           \
        {return reinterpret_cast<TWithoutNegator&>(t);}                                   \
    static ScalarNormSq scalarNormSqr(const T& t) {return t*t;}                           \
    static TSqrt        sqrt(const T& t) {return std::sqrt(t);}                           \
    static TAbs         abs(const T& t) {return std::abs(t);}                             \
    static const TStandard& standardize(const T& t) {return t;}                           \
    static TNormalize normalize(const T& t) {return (t>0?T(1):(t<0?T(-1):getNaN()));}     \
    static TInvert invert(const T& t) {return T(1)/t;}                                    \
    /* properties of this floating point representation */                                \
    static constexpr T getEps()         {return RTraits<T>::getEps();}                    \
    static constexpr T getSignificant() {return RTraits<T>::getSignificant();}            \
    static constexpr T getNaN()         {return std::numeric_limits<T>::quiet_NaN();}     \
    static constexpr T getInfinity()    {return std::numeric_limits<T>::infinity();}      \
    static constexpr T getLeastPositive()  {return std::numeric_limits<T>::min();}        \
    static constexpr T getMostPositive()   {return std::numeric_limits<T>::max();}        \
    static constexpr T getLeastNegative()  {return -std::numeric_limits<T>::min();}       \
    static constexpr T getMostNegative()   {return -std::numeric_limits<T>::max();}       \
    static constexpr T getSqrtEps()     {return detail::consteval_sqrt(getEps());}                                      \
    static constexpr T getTiny()        {constexpr auto s=detail::consteval_sqrt(getEps()); return getEps()*detail::consteval_sqrt(s);}   \
    static bool isFinite(const T& t) {return SimTK::isFinite(t);}                                               \
    static bool isNaN   (const T& t) {return SimTK::isNaN(t);}                                                  \
    static bool isInf   (const T& t) {return SimTK::isInf(t);}                                                  \
    /* Methods to use for approximate comparisons. Perform comparison in the wider of the two */                \
    /* precisions, using the default tolerance from the narrower of the two precisions.       */                \
    static double getDefaultTolerance() {return RTraits<T>::getDefaultTolerance();}                             \
    static bool isNumericallyEqual(const T& t, const float& f) {return SimTK::isNumericallyEqual(t,f);}         \
    static bool isNumericallyEqual(const T& t, const double& d) {return SimTK::isNumericallyEqual(t,d);}        \
    static bool isNumericallyEqual(const T& t, int i) {return SimTK::isNumericallyEqual(t,i);}                  \
    /* Here the tolerance is given so we don't have to figure it out. */                                                        \
    static bool isNumericallyEqual(const T& t, const float& f, double tol){return SimTK::isNumericallyEqual(t,f,tol);}          \
    static bool isNumericallyEqual(const T& t, const double& d, double tol){return SimTK::isNumericallyEqual(t,d,tol);}         \
    static bool isNumericallyEqual(const T& t, int i, double tol){return SimTK::isNumericallyEqual(t,i,tol);}                   \
    /* Carefully calculated constants. */                             \
    static constexpr T getZero()         {return T(0);}               \
    static constexpr T getOne()          {return T(1);}               \
    static constexpr T getMinusOne()     {return T(-1);}              \
    static constexpr T getTwo()          {return T(2);}               \
    static constexpr T getThree()        {return T(3);}               \
    static constexpr T getOneHalf()      {return T(0.5L);}            \
    static constexpr T getOneThird()     {return T(1.L/3.L);}         \
    static constexpr T getOneFourth()    {return T(0.25L);}           \
    static constexpr T getOneFifth()     {return T(0.2L);}            \
    static constexpr T getOneSixth()     {return T(1.L/6.L);}         \
    static constexpr T getOneSeventh()   {return T(1.L/7.L);}         \
    static constexpr T getOneEighth()    {return T(0.125L);}          \
    static constexpr T getOneNinth()     {return T(1.L/9.L);}         \
    static constexpr T getPi()           {return T(SimTK_PI);}        \
    static constexpr T getOneOverPi()    {return T(1.L/SimTK_PI);}    \
    static constexpr T getE()            {return T(SimTK_E);}         \
    static constexpr T getLog2E()        {return T(SimTK_LOG2E);}     \
    static constexpr T getLog10E()       {return T(SimTK_LOG10E);}    \
    static constexpr T getSqrt2()        {return T(SimTK_SQRT2);}     \
    static constexpr T getOneOverSqrt2() {return T(1.L/SimTK_SQRT2);} \
    static constexpr T getSqrt3()        {return T(SimTK_SQRT3);}     \
    static constexpr T getOneOverSqrt3() {return T(1.L/SimTK_SQRT3);} \
    static constexpr T getCubeRoot2()    {return T(SimTK_CBRT2);}     \
    static constexpr T getCubeRoot3()    {return T(SimTK_CBRT3);}     \
    static constexpr T getLn2()          {return T(SimTK_LN2);}       \
    static constexpr T getLn10()         {return T(SimTK_LN10);}      \
    /* integer digit counts useful for formatted input and output */                     \
    static constexpr int getNumDigits()         {return std::numeric_limits<T>::digits10;}               \
    static constexpr int getLosslessNumDigits() {return (int)(5.0L*(std::numeric_limits<T>::digits-1)/4.0L*(SimTK_LN2/SimTK_LN10)+0.5L);} \
}; \
template<> struct NTraits<R>::Result<float> \
  {typedef Widest<R,float>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;};    \
template<> struct NTraits<R>::Result<double> \
  {typedef Widest<R,double>::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;};    \
template<> struct NTraits<R>::Result<complex<float> > \
  {typedef Widest<R,complex<float> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
template<> struct NTraits<R>::Result<complex<double> > \
  {typedef Widest<R,complex<double> >::Type Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
template<> struct NTraits<R>::Result<conjugate<float> > \
  {typedef conjugate<Widest<R,float>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}; \
template<> struct NTraits<R>::Result<conjugate<double> > \
  {typedef conjugate<Widest<R,double>::Type> Mul;typedef Mul Dvd;typedef Mul Add;typedef Mul Sub;}
 
#if defined(__clang__)
#pragma clang diagnostic push
// The function `T& imag(T&)` generates a null-dereference warning.
#pragma clang diagnostic ignored "-Wnull-dereference"
#endif
SimTK_DEFINE_REAL_NTRAITS(float);
SimTK_DEFINE_REAL_NTRAITS(double);
#if defined(__clang__)
#pragma clang diagnostic pop
#endif
 
#undef SimTK_DEFINE_REAL_NTRAITS
 
template <class R> class CNT< complex<R> > : public NTraits< complex<R> > { };
template <class R> class CNT< conjugate<R> > : public NTraits< conjugate<R> > { };
template <> class CNT<float> : public NTraits<float> { };
template <> class CNT<double> : public NTraits<double> { };
 
 
} // namespace SimTK
 
#endif //SimTK_SIMMATRIX_NTRAITS_H_