Simbody  3.7
negator.h
Go to the documentation of this file.
1 #ifndef SimTK_SIMMATRIX_NEGATOR_H_
2 #define SimTK_SIMMATRIX_NEGATOR_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 
57 #include <iostream>
58 
59 namespace SimTK {
60 
61 // Specializations of this class provide information about Composite Numerical Types
62 // (i.e. composite types) in the style of std::numeric_limits<T>.
63 template <class T> class CNT;
64 template <class N> class NTraits; // Same, but only defined for numeric types.
65 template <class N> class negator; // negator is only defined for numbers.
66 
67 
74 template <class NUMBER>
76  typedef NUMBER N;
77  typedef typename NTraits<N>::TReal NReal;
78  typedef typename NTraits<N>::TImag NImag;
79  typedef typename NTraits<N>::TComplex NComplex;
80  typedef typename NTraits<N>::THerm NHerm;
81  typedef typename NTraits<N>::TInvert NInvert;
82 public:
83  typedef negator<N> T;
84  typedef NUMBER TNeg; // negator evaporates
85  typedef NUMBER TWithoutNegator; // "
86  typedef typename CNT<NReal>::TNeg TReal;
87  typedef typename CNT<NImag>::TNeg TImag;
88  typedef typename CNT<NComplex>::TNeg TComplex;
89  typedef typename CNT<NHerm>::TNeg THerm;
91  typedef typename NTraits<N>::TSqHermT TSqHermT;
92  typedef typename NTraits<N>::TSqTHerm TSqTHerm;
94  typedef negator<N> TRow;
95  typedef negator<N> TCol;
96 
97  typedef typename NTraits<N>::TSqrt TSqrt;
98  typedef typename NTraits<N>::TAbs TAbs;
99  typedef typename NTraits<N>::TStandard TStandard;
100  typedef typename CNT<NInvert>::TNeg TInvert;
101  typedef typename NTraits<N>::TStandard TNormalize; // neg<conj> -> complex
102 
103 
106  typedef NUMBER Number;
110 
111  // negator may be used in combination with any composite numerical type, not just
112  // numbers. Hence we must use CNT<P> here rather than NTraits<P> (they are
113  // the same when P turns out to be a numeric type).
114  // Typeof( Neg<N>*P ) is Typeof(P*N)::TNeg
115  // Typeof( Neg<N>/P ) is Typeof(N/P)::TNeg
116  // Typeof( Neg<N>+P ) is Typeof(P-N)
117  // Typeof( Neg<N>-P ) is Typeof(P+N)::TNeg
118  // No need to specialize because we strip off the "negator" here.
119  template <class P> struct Result {
120  private:
121  // These are the type of the calculations we actually perform.
122  typedef typename CNT<P>::template Result<N>::Mul PMul;
123  typedef typename NTraits<N>::template Result<P>::Dvd PDvd;
124  typedef typename CNT<P>::template Result<N>::Sub PAdd;
125  typedef typename CNT<P>::template Result<N>::Add PSub;
126  public:
127  // These are the types to which we must recast the results.
128  typedef typename CNT<PMul>::TNeg Mul;
129  typedef typename CNT<PDvd>::TNeg Dvd;
130  typedef PAdd Add;
131  typedef typename CNT<PSub>::TNeg Sub;
132  };
133 
134  // Shape-preserving element substitution (easy for scalars!)
135  template <class P> struct Substitute {
136  typedef P Type;
137  };
138 
139  enum {
140  NRows = 1, // Negators are always scalars
141  NCols = 1,
142  RowSpacing = 1,
143  ColSpacing = 1,
144  NPackedElements = 1,
145  NActualElements = 1,
146  NActualScalars = 1,
148  RealStrideFactor = NTraits<N>::RealStrideFactor,
149  ArgDepth = SCALAR_DEPTH,
150  IsScalar = 1,
151  IsULessScalar = 1,
152  IsNumber = 0,
153  IsStdNumber = 0,
154  IsPrecision = 0,
155  SignInterpretation = -1 // if you cast away the negator, don't forget this!
156  };
157  const negator<N>* getData() const {return this;}
158  negator<N>* updData() {return this;}
159 
160  const TReal& real() const {return reinterpret_cast<const TReal&>(NTraits<N>::real(v));}
161  TReal& real() {return reinterpret_cast< TReal&>(NTraits<N>::real(v));}
162  const TImag& imag() const {return reinterpret_cast<const TImag&>(NTraits<N>::imag(v));}
163  TImag& imag() {return reinterpret_cast< TImag&>(NTraits<N>::imag(v));}
164 
166  TSqrt sqrt() const {return NTraits<N>::sqrt(N(v));}
167  TAbs abs() const {return NTraits<N>::abs(v);}
170 
171  // Inverse (1/x) of a non-negated type N will return a non-negated type, so we
172  // can cast it to a negated type here to save a flop. The return type might
173  // not be N (a negated conjugate comes back as a complex), so there may be
174  // a flop done in the final conversion to TInvert.
175  TInvert invert() const {return recast(NTraits<N>::invert(v));}
176 
177  static negator<N> getNaN() {return recast(NTraits<N>::getNaN());}
178  static negator<N> getInfinity() {return recast(NTraits<N>::getInfinity());}
179 
181  inline bool isFinite() const;
183  inline bool isNaN() const;
186  inline bool isInf() const;
187 
189 
192  template <class T2> bool isNumericallyEqual(const T2& t2) const
193  { return CNT<T2>::isNumericallyEqual(t2, -v); } // perform negation
194 
198  template <class N2> bool isNumericallyEqual(const negator<N2>& t2) const
199  { return SimTK::isNumericallyEqual(v, t2.v); }
200 
202  template <class T2> bool isNumericallyEqual(const T2& t2, double tol) const
203  { return CNT<T2>::isNumericallyEqual(t2, -v, tol); } // perform negation
204 
207  template <class N2> bool isNumericallyEqual(const negator<N2>& t2, double tol) const
208  { return SimTK::isNumericallyEqual(v, t2.v, tol); }
209 
210 
212  #ifndef NDEBUG
213  v = NTraits<N>::getNaN();
214  #endif
215  }
216  negator(const negator& n) : v(n.v) { }
217  negator& operator=(const negator& n) { v=n.v; return *this; }
218 
219  // These are implicit conversions from numeric type NN to negator<N>. The
220  // value must be unchanged, so we must negate. Note that NN==N is a
221  // certainty for one of these cases.
222  negator(int t) {v = -N((typename NTraits<N>::Precision)t);}
223  negator(const float& t) {v = -N((typename NTraits<N>::Precision)t);}
224  negator(const double& t) {v = -N((typename NTraits<N>::Precision)t);}
225 
226  // Some of these may not compile if instantiated -- you can't cast a complex
227  // to a float, for example.
228  template <class P> negator(const std::complex<P>& t) {v = -N(t);}
229  template <class P> negator(const conjugate<P>& t) {v = -N(t);}
230 
231  // This can be used to negate a value of type N at zero cost. It is typically
232  // used for recasting temporary expressions to apply a final negation. Note that
233  // this is *not* the same as constructing a negator<N> from an N, which actually
234  // peforms a floating point negation.
235  static const negator<N>& recast(const N& val)
236  { return reinterpret_cast<const negator<N>&>(val); }
237 
238  const N& operator-() const { return v; }
239  N& operator-() { return v; } // an lvalue!
240  N operator+() const { return N(-v); } // performs the actual negation (expensive)
241 
242  operator N() const { return N(-v); } // implicit conversion to N (expensive)
243 
244  template <class P> negator& operator =(const P& t) { v = -t; return *this; }
245  template <class P> negator& operator+=(const P& t) { v -= t; return *this; } //swap sign
246  template <class P> negator& operator-=(const P& t) { v += t; return *this; }
247  template <class P> negator& operator*=(const P& t) { v *= t; return *this; } //don't swap!
248  template <class P> negator& operator/=(const P& t) { v /= t; return *this; }
249 
250  // If we know we've got a negator as an argument, get rid of its negation
251  // and change signs as necessary. We're guaranteed to get rid of at least
252  // one negator<> this way. Nothing to gain for multiplication or division,
253  // though.
254  template <class NN> negator& operator =(const negator<NN>& t)
255  { v = -t; return *this; }
256  template <class NN> negator& operator+=(const negator<NN>& t)
257  { v += -t; return *this; } //swap sign
258  template <class NN> negator& operator-=(const negator<NN>& t)
259  { v -= -t; return *this; }
260 
261 private:
262  N v;
263 
264 template <class N2> friend class negator;
265 };
266 
267 // isNaN() for real, complex, and conjugate numbers is provided in
268 // NTraits. Here we add isNaN() for negated scalar types.
269 
271 
272 inline bool isNaN(const negator<float>& x) {return isNaN(-x);}
273 inline bool isNaN(const negator<double>& x) {return isNaN(-x);}
274 template <class P> inline bool
275 isNaN(const negator< std::complex<P> >& x) {return isNaN(-x);}
276 template <class P> inline bool
277 isNaN(const negator< conjugate<P> >& x) {return isNaN(-x);}
279 
280 // isFinite() for real, complex, and conjugate numbers is provided in
281 // NTraits. Here we add isFinite() for negated scalar types.
282 
284 
285 inline bool isFinite(const negator<float>& x) {return isFinite(-x);}
286 inline bool isFinite(const negator<double>& x) {return isFinite(-x);}
287 template <class P> inline bool
288 isFinite(const negator< std::complex<P> >& x) {return isFinite(-x);}
289 template <class P> inline bool
290 isFinite(const negator< conjugate<P> >& x) {return isFinite(-x);}
292 
293 // isInf(x) for real, complex, and conjugate numbers is provided in
294 // NTraits. Here we add isInf() for negated scalar types.
295 
297 
298 inline bool isInf(const negator<float>& x) {return isInf(-x);}
299 inline bool isInf(const negator<double>& x) {return isInf(-x);}
300 template <class P> inline bool
301 isInf(const negator< std::complex<P> >& x) {return isInf(-x);}
302 template <class P> inline bool
303 isInf(const negator< conjugate<P> >& x) {return isInf(-x);}
305 
306 // The member functions call the global ones just defined.
307 template <class N> inline bool
308 negator<N>::isFinite() const {return SimTK::isFinite(*this);}
309 template <class N> inline bool
310 negator<N>::isNaN() const {return SimTK::isNaN(*this);}
311 template <class N> inline bool
312 negator<N>::isInf() const {return SimTK::isInf(*this);}
313 
314 // Handle all binary numerical operators involving a negator<A> and a B, or negator<A>
315 // and negator<B>, obtaining results by stripping away the negator<>s and fiddling
316 // with signs appropriately.
317 // Careful: don't remove both negators in one step because Result isn't specialized
318 // for negators so it might not predict the same result type as you actually get.
319 // Be patient and let it strip one negator at a time -- in Release mode that won't
320 // add any code since all this stuff drops away.
321 //
322 // To appreciate the beauty of these operators, remember that -x is free when x
323 // is a negator.
324 
325 template <class DEST, class SRC> static inline const DEST&
326 negRecast(const SRC& s) { return reinterpret_cast<const DEST&>(s); }
327 
328  // Binary '+' with a negator as one or both arguments //
329 template <class A, class B> inline typename negator<A>::template Result<B>::Add
330 operator+(const negator<A>& l, const B& r)
331  {return negRecast<typename negator<A>::template Result<B>::Add>(r-(-l));}
332 template <class A, class B> inline typename CNT<A>::template Result<negator<B> >::Add
333 operator+(const A& l, const negator<B>& r)
334  {return negRecast<typename CNT<A>::template Result<negator<B> >::Add>(l-(-r));}
335 // One step at a time here (see above).
336 template <class A, class B> inline typename negator<A>::template Result<negator<B> >::Add
337 operator+(const negator<A>& l, const negator<B>& r)
338  {return negRecast<typename negator<A>::template Result<negator<B> >::Add>(r-(-l)); }
339 
340  // Binary '-' with a negator as one or both arguments //
341 template <class A, class B> inline typename negator<A>::template Result<B>::Sub
342 operator-(const negator<A>& l, const B& r)
343  {return negRecast<typename negator<A>::template Result<B>::Sub>(r+(-l));}
344 template <class A, class B> inline typename CNT<A>::template Result<negator<B> >::Sub
345 operator-(const A& l, const negator<B>& r)
346  {return negRecast<typename CNT<A>::template Result<negator<B> >::Sub>(l+(-r));}
347 // One step at a time here (see above).
348 template <class A, class B> inline typename negator<A>::template Result<negator<B> >::Sub
349 operator-(const negator<A>& l, const negator<B>& r)
350  {return negRecast<typename negator<A>::template Result<negator<B> >::Sub>(r+(-l));}
351 
352 // Binary '*' with a negator as one or both arguments
353 template <class A, class B> inline typename negator<A>::template Result<B>::Mul
354 operator*(const negator<A>& l, const B& r)
355  {return negRecast<typename negator<A>::template Result<B>::Mul>((-l)*r);}
356 template <class A, class B> inline typename CNT<A>::template Result<negator<B> >::Mul
357 operator*(const A& l, const negator<B>& r)
358  {return negRecast<typename CNT<A>::template Result<negator<B> >::Mul>(l*(-r));}
359 // One step at a time here (see above).
360 template <class A, class B> inline typename negator<A>::template Result<negator<B> >::Mul
361 operator*(const negator<A>& l, const negator<B>& r)
362  {return negRecast<typename negator<A>::template Result<negator<B> >::Mul>((-l)*r);}
363 
364 // Binary '/' with a negator as one or both arguments
365 template <class A, class B> inline typename negator<A>::template Result<B>::Dvd
366 operator/(const negator<A>& l, const B& r)
367  {return negRecast<typename negator<A>::template Result<B>::Dvd>((-l)/r);}
368 template <class A, class B> inline typename CNT<A>::template Result<negator<B> >::Dvd
369 operator/(const A& l, const negator<B>& r)
370  {return negRecast<typename CNT<A>::template Result<negator<B> >::Dvd>(l/(-r));}
371 // One step at a time here (see above).
372 template <class A, class B> inline typename negator<A>::template Result<negator<B> >::Dvd
373 operator/(const negator<A>& l, const negator<B>& r)
374  {return negRecast<typename negator<A>::template Result<negator<B> >::Dvd>((-l)/r);}
375 
376 // Binary '==' with a negator as one or both arguments
377 template <class A, class B> inline bool
378 operator==(const negator<A>& l, const B& r) { return (A)l == r; }
379 template <class A, class B> inline bool
380 operator==(const A& l, const negator<B>& r) { return l == (B)r; }
381 template <class A, class B> inline bool
382 operator==(const negator<A>& l, const negator<B>& r) { return (-l) == (-r); } // cheap!
383 
384 // The lazy man's '!=' operator.
385 template <class A, class B> inline bool
386 operator!=(const negator<A>& l, const B& r) { return !(l==r); }
387 template <class A, class B> inline bool
388 operator!=(const A& l, const negator<B>& r) { return !(l==r); }
389 template <class A, class B> inline bool
390 operator!=(const negator<A>& l, const negator<B>& r) { return !(l==r); }
391 
392 // And a final touch of elegance allowing smooth interoperability with iostream.
393 template <class NUM, class CHAR, class TRAITS> inline std::basic_istream<CHAR,TRAITS>&
394 operator>>(std::basic_istream<CHAR,TRAITS>& is, negator<NUM>& nn) {
395  NUM z; is >> z; nn=z;
396  return is;
397 }
398 template <class NUM, class CHAR, class TRAITS> inline std::basic_ostream<CHAR,TRAITS>&
399 operator<<(std::basic_ostream<CHAR,TRAITS>& os, const negator<NUM>& nn) {
400  return os << NUM(nn);
401 }
402 
403 } // namespace SimTK
404 
405 #endif //SimTK_SIMMATRIX_NEGATOR_H_
Matrix_< E > operator/(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition: BigMatrix.h:613
TImag & imag()
Definition: negator.h:163
#define SimTK_SimTKCOMMON_EXPORT
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:224
bool isNumericallyEqual(const T2 &t2, double tol) const
This is the generic case (see above) but with an explicitly-provided tolerance.
Definition: negator.h:202
NUMBER TWithoutNegator
Definition: negator.h:85
negator & operator/=(const P &t)
Definition: negator.h:248
negator< N > TElement
Definition: negator.h:93
PAdd Add
Definition: negator.h:130
negator< N > TCol
Definition: negator.h:95
CNT< PDvd >::TNeg Dvd
Definition: negator.h:129
TAbs abs() const
Definition: negator.h:167
static negator< N > getInfinity()
Definition: negator.h:178
negator< N > ULessScalar
Definition: negator.h:105
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition: Assembler.h:37
NTraits< N >::StdNumber StdNumber
Definition: negator.h:107
SimTK::conjugate<R> should be instantiated only for float, double.
Definition: String.h:45
negator< N > * updData()
Definition: negator.h:158
bool isInf() const
Returns true if the negated value contains an Inf or -Inf and does not contain a NaN.
Definition: negator.h:312
Matrix_< typename CNT< E1 >::template Result< E2 >::Sub > operator-(const MatrixBase< E1 > &l, const MatrixBase< E2 > &r)
Definition: BigMatrix.h:584
negator(int t)
Definition: negator.h:222
TReal & real()
Definition: negator.h:161
NUMBER TNeg
Definition: negator.h:84
std::basic_istream< CHAR, TRAITS > & operator>>(std::basic_istream< CHAR, TRAITS > &is, conjugate< R > &c)
Definition: conjugate.h:505
negator< N > TRow
Definition: negator.h:94
negator<N>, where N is a number type (real, complex, conjugate), is represented in memory identically...
Definition: String.h:44
bool isNumericallyEqual(const negator< N2 > &t2) const
In this partial specialization we know that both types have negators so we can just compare the under...
Definition: negator.h:198
NUMBER Number
Definition: negator.h:106
static negator< N > getNaN()
Definition: negator.h:177
bool operator==(const PhiMatrix &p1, const PhiMatrix &p2)
Definition: SpatialAlgebra.h:791
static bool isNumericallyEqual(const K &t1, const K2 &t2)
CNTs are expected to support an "==" operator for exact, bitwise equality.
Definition: CompositeNumericalTypes.h:264
negator(const std::complex< P > &t)
Definition: negator.h:228
NTraits< N >::ScalarNormSq ScalarNormSq
Definition: negator.h:109
bool isFinite(const negator< float > &x)
Definition: negator.h:285
Definition: CompositeNumericalTypes.h:116
CNT< NImag >::TNeg TImag
Definition: negator.h:87
const TReal & real() const
Definition: negator.h:160
NTraits< N >::TSqHermT TSqHermT
Definition: negator.h:91
N operator+() const
Definition: negator.h:240
Matrix_< E > operator*(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition: BigMatrix.h:605
bool isNumericallyEqual(const negator< N2 > &t2, double tol) const
This is the partially specialized case again (see above) but with an explicitly-provided tolerance...
Definition: negator.h:207
static const DEST & negRecast(const SRC &s)
Definition: negator.h:326
CNT< NHerm >::TNeg THerm
Definition: negator.h:89
negator & operator+=(const P &t)
Definition: negator.h:245
TStandard standardize() const
Definition: negator.h:168
const negator< N > * getData() const
Definition: negator.h:157
negator< N > TPosTrans
Definition: negator.h:90
NTraits< N >::TAbs TAbs
Definition: negator.h:98
negator()
Definition: negator.h:211
Definition: negator.h:135
CNT< NReal >::TNeg TReal
Definition: negator.h:86
TInvert invert() const
Definition: negator.h:175
const float & real(const conjugate< float > &c)
Definition: conjugate.h:482
negator & operator*=(const P &t)
Definition: negator.h:247
RowVectorBase< typename CNT< ELEM >::TAbs > abs(const RowVectorBase< ELEM > &v)
Definition: VectorMath.h:120
bool isNaN(const negator< float > &x)
Definition: negator.h:272
negator(const negator &n)
Definition: negator.h:216
negator< N > Scalar
Definition: negator.h:104
Definition: negator.h:119
Specialized information about Composite Numerical Types which allows us to define appropriate templat...
Definition: CompositeNumericalTypes.h:136
bool isInf(const negator< float > &x)
Definition: negator.h:298
NTraits< N >::TStandard TNormalize
Definition: negator.h:101
static const negator< N > & recast(const N &val)
Definition: negator.h:235
K::TNeg TNeg
Definition: CompositeNumericalTypes.h:139
TNormalize normalize() const
Definition: negator.h:169
CNT< NInvert >::TNeg TInvert
Definition: negator.h:100
Matrix_< typename CNT< E1 >::template Result< E2 >::Add > operator+(const MatrixBase< E1 > &l, const MatrixBase< E2 > &r)
Definition: BigMatrix.h:568
CNT< PMul >::TNeg Mul
Definition: negator.h:128
NTraits< N >::TStandard TStandard
Definition: negator.h:99
negator & operator+=(const negator< NN > &t)
Definition: negator.h:256
NTraits< N >::TSqTHerm TSqTHerm
Definition: negator.h:92
negator(const float &t)
Definition: negator.h:223
const negator< float > & imag(const conjugate< float > &c)
Definition: conjugate.h:483
bool isFinite() const
Returns true if the negated value is finite (i.e., not NaN or Inf).
Definition: negator.h:308
bool isNaN() const
Returns true if the negated value contains a NaN.
Definition: negator.h:310
CNT< NComplex >::TNeg TComplex
Definition: negator.h:88
bool operator!=(const L &left, const R &right)
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:641
bool isNumericallyEqual(const T2 &t2) const
In the generic case we&#39;ll perform the negation here to get a number, and then delegate to the other t...
Definition: negator.h:192
const N & operator-() const
Definition: negator.h:238
const TImag & imag() const
Definition: negator.h:162
TSqrt sqrt() const
Definition: negator.h:166
negator(const double &t)
Definition: negator.h:224
CNT< PSub >::TNeg Sub
Definition: negator.h:131
static double getDefaultTolerance()
Definition: negator.h:188
negator & operator-=(const P &t)
Definition: negator.h:246
negator< N > T
Definition: negator.h:83
N & operator-()
Definition: negator.h:239
NTraits< N >::Precision Precision
Definition: negator.h:108
Definition: negator.h:64
negator & operator-=(const negator< NN > &t)
Definition: negator.h:258
negator(const conjugate< P > &t)
Definition: negator.h:229
bool isNumericallyEqual(const float &a, const float &b, double tol=RTraits< float >::getDefaultTolerance())
Compare two floats for approximate equality.
Definition: NTraits.h:293
ScalarNormSq scalarNormSqr() const
Definition: negator.h:165
NTraits< N >::TSqrt TSqrt
Definition: negator.h:97
negator & operator=(const negator &n)
Definition: negator.h:217
P Type
Definition: negator.h:136