Simbody  3.5
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 
165  ScalarNormSq scalarNormSqr() const {return NTraits<N>::scalarNormSqr(v);}
166  TSqrt sqrt() const {return NTraits<N>::sqrt(N(v));}
167  TAbs abs() const {return NTraits<N>::abs(v);}
168  TStandard standardize() const {return -NTraits<N>::standardize(v);}
169  TNormalize normalize() const {return -NTraits<N>::normalize(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  negator(const long double& t) {v = -N((typename NTraits<N>::Precision)t);}
226 
227  // Some of these may not compile if instantiated -- you can't cast a complex
228  // to a float, for example.
229  template <class P> negator(const std::complex<P>& t) {v = -N(t);}
230  template <class P> negator(const conjugate<P>& t) {v = -N(t);}
231 
232  // This can be used to negate a value of type N at zero cost. It is typically
233  // used for recasting temporary expressions to apply a final negation. Note that
234  // this is *not* the same as constructing a negator<N> from an N, which actually
235  // peforms a floating point negation.
236  static const negator<N>& recast(const N& val)
237  { return reinterpret_cast<const negator<N>&>(val); }
238 
239  const N& operator-() const { return v; }
240  N& operator-() { return v; } // an lvalue!
241  N operator+() const { return N(-v); } // performs the actual negation (expensive)
242 
243  operator N() const { return N(-v); } // implicit conversion to N (expensive)
244 
245  template <class P> negator& operator =(const P& t) { v = -t; return *this; }
246  template <class P> negator& operator+=(const P& t) { v -= t; return *this; } //swap sign
247  template <class P> negator& operator-=(const P& t) { v += t; return *this; }
248  template <class P> negator& operator*=(const P& t) { v *= t; return *this; } //don't swap!
249  template <class P> negator& operator/=(const P& t) { v /= t; return *this; }
250 
251  // If we know we've got a negator as an argument, get rid of its negation
252  // and change signs as necessary. We're guaranteed to get rid of at least
253  // one negator<> this way. Nothing to gain for multiplication or division,
254  // though.
255  template <class NN> negator& operator =(const negator<NN>& t)
256  { v = -t; return *this; }
257  template <class NN> negator& operator+=(const negator<NN>& t)
258  { v += -t; return *this; } //swap sign
259  template <class NN> negator& operator-=(const negator<NN>& t)
260  { v -= -t; return *this; }
261 
262 private:
263  N v;
264 
265 template <class N2> friend class negator;
266 };
267 
268 // isNaN() for real, complex, and conjugate numbers is provided in
269 // NTraits. Here we add isNaN() for negated scalar types.
270 
272 
273 inline bool isNaN(const negator<float>& x) {return isNaN(-x);}
274 inline bool isNaN(const negator<double>& x) {return isNaN(-x);}
275 inline bool isNaN(const negator<long double>& x) {return isNaN(-x);}
276 template <class P> inline bool
277 isNaN(const negator< std::complex<P> >& x) {return isNaN(-x);}
278 template <class P> inline bool
279 isNaN(const negator< conjugate<P> >& x) {return isNaN(-x);}
281 
282 // isFinite() for real, complex, and conjugate numbers is provided in
283 // NTraits. Here we add isFinite() for negated scalar types.
284 
286 
287 inline bool isFinite(const negator<float>& x) {return isFinite(-x);}
288 inline bool isFinite(const negator<double>& x) {return isFinite(-x);}
289 inline bool isFinite(const negator<long double>& x) {return isFinite(-x);}
290 template <class P> inline bool
291 isFinite(const negator< std::complex<P> >& x) {return isFinite(-x);}
292 template <class P> inline bool
293 isFinite(const negator< conjugate<P> >& x) {return isFinite(-x);}
295 
296 // isInf(x) for real, complex, and conjugate numbers is provided in
297 // NTraits. Here we add isInf() for negated scalar types.
298 
300 
301 inline bool isInf(const negator<float>& x) {return isInf(-x);}
302 inline bool isInf(const negator<double>& x) {return isInf(-x);}
303 inline bool isInf(const negator<long double>& x) {return isInf(-x);}
304 template <class P> inline bool
305 isInf(const negator< std::complex<P> >& x) {return isInf(-x);}
306 template <class P> inline bool
307 isInf(const negator< conjugate<P> >& x) {return isInf(-x);}
309 
310 // The member functions call the global ones just defined.
311 template <class N> inline bool
312 negator<N>::isFinite() const {return SimTK::isFinite(*this);}
313 template <class N> inline bool
314 negator<N>::isNaN() const {return SimTK::isNaN(*this);}
315 template <class N> inline bool
316 negator<N>::isInf() const {return SimTK::isInf(*this);}
317 
318 // Handle all binary numerical operators involving a negator<A> and a B, or negator<A>
319 // and negator<B>, obtaining results by stripping away the negator<>s and fiddling
320 // with signs appropriately.
321 // Careful: don't remove both negators in one step because Result isn't specialized
322 // for negators so it might not predict the same result type as you actually get.
323 // Be patient and let it strip one negator at a time -- in Release mode that won't
324 // add any code since all this stuff drops away.
325 //
326 // To appreciate the beauty of these operators, remember that -x is free when x
327 // is a negator.
328 
329 template <class DEST, class SRC> static inline const DEST&
330 negRecast(const SRC& s) { return reinterpret_cast<const DEST&>(s); }
331 
332  // Binary '+' with a negator as one or both arguments //
333 template <class A, class B> inline typename negator<A>::template Result<B>::Add
334 operator+(const negator<A>& l, const B& r)
335  {return negRecast<typename negator<A>::template Result<B>::Add>(r-(-l));}
336 template <class A, class B> inline typename CNT<A>::template Result<negator<B> >::Add
337 operator+(const A& l, const negator<B>& r)
338  {return negRecast<typename CNT<A>::template Result<negator<B> >::Add>(l-(-r));}
339 // One step at a time here (see above).
340 template <class A, class B> inline typename negator<A>::template Result<negator<B> >::Add
341 operator+(const negator<A>& l, const negator<B>& r)
342  {return negRecast<typename negator<A>::template Result<negator<B> >::Add>(r-(-l)); }
343 
344  // Binary '-' with a negator as one or both arguments //
345 template <class A, class B> inline typename negator<A>::template Result<B>::Sub
346 operator-(const negator<A>& l, const B& r)
347  {return negRecast<typename negator<A>::template Result<B>::Sub>(r+(-l));}
348 template <class A, class B> inline typename CNT<A>::template Result<negator<B> >::Sub
349 operator-(const A& l, const negator<B>& r)
350  {return negRecast<typename CNT<A>::template Result<negator<B> >::Sub>(l+(-r));}
351 // One step at a time here (see above).
352 template <class A, class B> inline typename negator<A>::template Result<negator<B> >::Sub
353 operator-(const negator<A>& l, const negator<B>& r)
354  {return negRecast<typename negator<A>::template Result<negator<B> >::Sub>(r+(-l));}
355 
356 // Binary '*' with a negator as one or both arguments
357 template <class A, class B> inline typename negator<A>::template Result<B>::Mul
358 operator*(const negator<A>& l, const B& r)
359  {return negRecast<typename negator<A>::template Result<B>::Mul>((-l)*r);}
360 template <class A, class B> inline typename CNT<A>::template Result<negator<B> >::Mul
361 operator*(const A& l, const negator<B>& r)
362  {return negRecast<typename CNT<A>::template Result<negator<B> >::Mul>(l*(-r));}
363 // One step at a time here (see above).
364 template <class A, class B> inline typename negator<A>::template Result<negator<B> >::Mul
365 operator*(const negator<A>& l, const negator<B>& r)
366  {return negRecast<typename negator<A>::template Result<negator<B> >::Mul>((-l)*r);}
367 
368 // Binary '/' with a negator as one or both arguments
369 template <class A, class B> inline typename negator<A>::template Result<B>::Dvd
370 operator/(const negator<A>& l, const B& r)
371  {return negRecast<typename negator<A>::template Result<B>::Dvd>((-l)/r);}
372 template <class A, class B> inline typename CNT<A>::template Result<negator<B> >::Dvd
373 operator/(const A& l, const negator<B>& r)
374  {return negRecast<typename CNT<A>::template Result<negator<B> >::Dvd>(l/(-r));}
375 // One step at a time here (see above).
376 template <class A, class B> inline typename negator<A>::template Result<negator<B> >::Dvd
377 operator/(const negator<A>& l, const negator<B>& r)
378  {return negRecast<typename negator<A>::template Result<negator<B> >::Dvd>((-l)/r);}
379 
380 // Binary '==' with a negator as one or both arguments
381 template <class A, class B> inline bool
382 operator==(const negator<A>& l, const B& r) { return (A)l == r; }
383 template <class A, class B> inline bool
384 operator==(const A& l, const negator<B>& r) { return l == (B)r; }
385 template <class A, class B> inline bool
386 operator==(const negator<A>& l, const negator<B>& r) { return (-l) == (-r); } // cheap!
387 
388 // The lazy man's '!=' operator.
389 template <class A, class B> inline bool
390 operator!=(const negator<A>& l, const B& r) { return !(l==r); }
391 template <class A, class B> inline bool
392 operator!=(const A& l, const negator<B>& r) { return !(l==r); }
393 template <class A, class B> inline bool
394 operator!=(const negator<A>& l, const negator<B>& r) { return !(l==r); }
395 
396 // And a final touch of elegance allowing smooth interoperability with iostream.
397 template <class NUM, class CHAR, class TRAITS> inline std::basic_istream<CHAR,TRAITS>&
398 operator>>(std::basic_istream<CHAR,TRAITS>& is, negator<NUM>& nn) {
399  NUM z; is >> z; nn=z;
400  return is;
401 }
402 template <class NUM, class CHAR, class TRAITS> inline std::basic_ostream<CHAR,TRAITS>&
403 operator<<(std::basic_ostream<CHAR,TRAITS>& os, const negator<NUM>& nn) {
404  return os << NUM(nn);
405 }
406 
407 } // namespace SimTK
408 
409 #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:218
NUMBER TWithoutNegator
Definition: negator.h:85
negator & operator/=(const P &t)
Definition: negator.h:249
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
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, long double.
Definition: String.h:45
negator< N > * updData()
Definition: negator.h:158
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
negator(const long double &t)
Definition: negator.h:225
const TImag & imag() const
Definition: negator.h:162
NUMBER TNeg
Definition: negator.h:84
std::basic_istream< CHAR, TRAITS > & operator>>(std::basic_istream< CHAR, TRAITS > &is, conjugate< R > &c)
Definition: conjugate.h:800
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
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:774
ScalarNormSq scalarNormSqr() const
Definition: negator.h:165
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:229
NTraits< N >::ScalarNormSq ScalarNormSq
Definition: negator.h:109
bool isFinite(const negator< float > &x)
Definition: negator.h:287
Definition: CompositeNumericalTypes.h:116
CNT< NImag >::TNeg TImag
Definition: negator.h:87
bool isInf() const
Returns true if the negated value contains an Inf or -Inf and does not contain a NaN.
Definition: negator.h:316
NTraits< N >::TSqHermT TSqHermT
Definition: negator.h:91
TNormalize normalize() const
Definition: negator.h:169
Matrix_< E > operator*(const MatrixBase< E > &l, const typename CNT< E >::StdNumber &r)
Definition: BigMatrix.h:605
static const DEST & negRecast(const SRC &s)
Definition: negator.h:330
CNT< NHerm >::TNeg THerm
Definition: negator.h:89
const negator< N > * getData() const
Definition: negator.h:157
negator & operator+=(const P &t)
Definition: negator.h:246
bool isNaN() const
Returns true if the negated value contains a NaN.
Definition: negator.h:314
const TReal & real() const
Definition: negator.h:160
negator< N > TPosTrans
Definition: negator.h:90
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
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
const float & real(const conjugate< float > &c)
Definition: conjugate.h:771
negator & operator*=(const P &t)
Definition: negator.h:248
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
RowVectorBase< typename CNT< ELEM >::TAbs > abs(const RowVectorBase< ELEM > &v)
Definition: VectorMath.h:120
bool isNaN(const negator< float > &x)
Definition: negator.h:273
bool operator!=(const conjugate< R > &a, const float &b)
Definition: conjugate.h:859
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 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
bool isInf(const negator< float > &x)
Definition: negator.h:301
NTraits< N >::TStandard TNormalize
Definition: negator.h:101
const N & operator-() const
Definition: negator.h:239
static const negator< N > & recast(const N &val)
Definition: negator.h:236
K::TNeg TNeg
Definition: CompositeNumericalTypes.h:139
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:257
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:772
CNT< NComplex >::TNeg TComplex
Definition: negator.h:88
N operator+() const
Definition: negator.h:241
TStandard standardize() const
Definition: negator.h:168
bool isFinite() const
Returns true if the negated value is finite (i.e., not NaN or Inf).
Definition: negator.h:312
TInvert invert() const
Definition: negator.h:175
negator(const double &t)
Definition: negator.h:224
CNT< PSub >::TNeg Sub
Definition: negator.h:131
TAbs abs() const
Definition: negator.h:167
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
static double getDefaultTolerance()
Definition: negator.h:188
negator & operator-=(const P &t)
Definition: negator.h:247
TSqrt sqrt() const
Definition: negator.h:166
negator< N > T
Definition: negator.h:83
N & operator-()
Definition: negator.h:240
NTraits< N >::Precision Precision
Definition: negator.h:108
Definition: negator.h:64
negator & operator-=(const negator< NN > &t)
Definition: negator.h:259
negator(const conjugate< P > &t)
Definition: negator.h:230
bool isNumericallyEqual(const float &a, const float &b, double tol=RTraits< float >::getDefaultTolerance())
Compare two floats for approximate equality.
Definition: NTraits.h:313
NTraits< N >::TSqrt TSqrt
Definition: negator.h:97
negator & operator=(const negator &n)
Definition: negator.h:217
P Type
Definition: negator.h:136