CompositeNumericalTypes.h

Go to the documentation of this file.

#ifndef SimTK_SIMMATRIX_COMPOSITE_NUMERICAL_TYPES_H_
#define SimTK_SIMMATRIX_COMPOSITE_NUMERICAL_TYPES_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/internal/common.h"
    
namespace SimTK {

// These are CNT "depths". 0 means the corresponding CNT is a scalar,
// 1 means it is a composite with scalar elements, 2 means a composite
// with composite elements, and 3 means a composite with depth-2
// composite elements. Beyond that the user will have to diambiguate
// operations by using named routines rather than operators.
enum  {
    SCALAR_DEPTH              = 0,
    SCALAR_COMPOSITE_DEPTH    = 1,
    COMPOSITE_COMPOSITE_DEPTH = 2,
    COMPOSITE_3_DEPTH         = 3,
    MAX_RESOLVED_DEPTH        = COMPOSITE_3_DEPTH
};

template <class K> class CNT : private K {
public:
    typedef K                        T;
    typedef typename K::TNeg         TNeg;
    typedef typename K::TWithoutNegator TWithoutNegator;
    typedef typename K::TReal        TReal;
    typedef typename K::TImag        TImag;
    typedef typename K::TComplex     TComplex;
    typedef typename K::THerm        THerm;
    typedef typename K::TPosTrans    TPosTrans;
    typedef typename K::TSqHermT     TSqHermT;
    typedef typename K::TSqTHerm     TSqTHerm;
    typedef typename K::TElement     TElement;
    typedef typename K::TRow         TRow;          // type of a row or column
    typedef typename K::TCol         TCol;

    // These are the results of calculations and should be packed regardless
    // of the spacing of this CNT.
    typedef typename K::TSqrt        TSqrt;         // also turns unit^2 to unit
    typedef typename K::TAbs         TAbs;
    typedef typename K::TStandard    TStandard;     // packed, StdNumbers
    typedef typename K::TInvert      TInvert;       // also turns units into 1/units
    typedef typename K::TNormalize   TNormalize;    // TODO: what effect on units?

    typedef typename K::Scalar       Scalar;        // quantity< units, <unitlessScalar> >
    typedef typename K::ULessScalar  ULessScalar;   // <number> or negator<number>
    typedef typename K::Number       Number;        // <real>, <complex> or <conjugate>
    typedef typename K::StdNumber    StdNumber;     // <real>, <complex>
    typedef typename K::Precision    Precision;     // float, double

    typedef typename K::ScalarNormSq ScalarNormSq;  // type of conjugate square of underlying scalar or
                                                    //   numeric value (squares the units too)

    template <class P> struct Result {
        typedef typename K::template Result<P>::Mul Mul;
        typedef typename K::template Result<P>::Dvd Dvd;
        typedef typename K::template Result<P>::Add Add;
        typedef typename K::template Result<P>::Sub Sub;
    };

    // Shape-preserving element substitution
    template <class P> struct Substitute {
        typedef typename K::template Substitute<P>::Type Type;
    };

    enum {
        NRows               = K::NRows,
        NCols               = K::NCols,
        RowSpacing          = K::RowSpacing,
        ColSpacing          = K::ColSpacing,
        NPackedElements     = K::NPackedElements,
        NActualElements     = K::NActualElements,
        NActualScalars      = K::NActualScalars,
        ImagOffset          = K::ImagOffset,
        RealStrideFactor    = K::RealStrideFactor,
        ArgDepth            = K::ArgDepth,
        IsScalar            = K::IsScalar,          // scalar with units, real, complex, conjugate, negator
        IsULessScalar       = K::IsULessScalar,     // real, complex, conjugate, negator
        IsNumber            = K::IsNumber,          // real, complex, conjugate
        IsStdNumber         = K::IsStdNumber,       // real, complex
        IsPrecision         = K::IsPrecision,       // real (float, double)
        SignInterpretation  = K::SignInterpretation // 1 normally, -1 if elements are negated
    };

    static const Scalar* getData(const T& t) { return t.getData(); }
    static       Scalar* updData(T& t)       { return t.updData(); }

    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(); }

    // We expect to be able to negate and transpose (hermitian or
    // positional) with just type casting; no need for help from class
    // K except to tell us the appropriate types.
    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 K& t)
      { return reinterpret_cast<const THerm&>(t); }
    static       THerm& transpose(K& t)
      { return reinterpret_cast<THerm&>(t); }

    static const TPosTrans& positionalTranspose(const K& t)
      { return reinterpret_cast<const TPosTrans&>(t); }
    static       TPosTrans& positionalTranspose(K& t)
      { return reinterpret_cast<TPosTrans&>(t); }

    // If the underlying scalars of this CNT are negator<N> for some numeric type N,
    // this method removes the negator<>, effectively negating the entire CNT. You
    // can still deal with the sign correctly by using the above enum SignInterpretation
    // which will be -1 in that case, 1 if there was no negator<> to remove. Note:
    // I'm not talking about TWithoutNegator::SignInterpretation -- that one is guaranteed
    // to be 1! T::SignInterpretation is the one you want.
    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 K& t) {return t.scalarNormSqr();}

    static TSqrt      sqrt(const K& t)          {return t.sqrt();}
    static TAbs       abs(const K& t)           {return t.abs();}
    static TStandard  standardize(const K& t)   {return t.standardize();}
    static TNormalize normalize(const K& t)     {return t.normalize();}
    static TInvert    invert(const K& t)        {return t.invert();}

    static K getInfinity() {return K::getInfinity();}
    static K getNaN()      {return K::getNaN();}

    static bool isNaN(const K& t) {return t.isNaN();}
    static bool isInf(const K& t) {return t.isInf();}
    static bool isFinite(const K& t) {return t.isFinite();}

    template <class K2> static bool 
    isNumericallyEqual(const K& t1, const K2& t2) 
    {   return t1.isNumericallyEqual(t2);}
    template <class K2> static bool 
    isNumericallyEqual(const K& t1, const K2& t2, double tol)
    {   return t1.isNumericallyEqual(t2,tol);}
    static double getDefaultTolerance() {return K::getDefaultTolerance();}

};

} // namespace SimTK

#endif // SimTK_SIMMATRIX_COMPOSITE_NUMERICAL_TYPES_H_