Scalar.h

Go to the documentation of this file.
00001 #ifndef SimTK_SIMMATRIX_SCALAR_H_
00002 #define SimTK_SIMMATRIX_SCALAR_H_
00003 
00004 /* -------------------------------------------------------------------------- *
00005  *                      SimTK Core: SimTK Simmatrix(tm)                       *
00006  * -------------------------------------------------------------------------- *
00007  * This is part of the SimTK Core biosimulation toolkit originating from      *
00008  * Simbios, the NIH National Center for Physics-Based Simulation of           *
00009  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
00010  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
00011  *                                                                            *
00012  * Portions copyright (c) 2005-9 Stanford University and the Authors.         *
00013  * Authors: Michael Sherman                                                   *
00014  * Contributors:                                                              *
00015  *                                                                            *
00016  * Permission is hereby granted, free of charge, to any person obtaining a    *
00017  * copy of this software and associated documentation files (the "Software"), *
00018  * to deal in the Software without restriction, including without limitation  *
00019  * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
00020  * and/or sell copies of the Software, and to permit persons to whom the      *
00021  * Software is furnished to do so, subject to the following conditions:       *
00022  *                                                                            *
00023  * The above copyright notice and this permission notice shall be included in *
00024  * all copies or substantial portions of the Software.                        *
00025  *                                                                            *
00026  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
00027  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
00028  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
00029  * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
00030  * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
00031  * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
00032  * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
00033  * -------------------------------------------------------------------------- */
00034 
00041 #include "SimTKcommon/internal/common.h"
00042 #include "SimTKcommon/Constants.h"
00043 #include "SimTKcommon/internal/Exception.h"
00044 #include "SimTKcommon/internal/ExceptionMacros.h"
00045 #include "SimTKcommon/internal/String.h"
00046 
00047 #include "SimTKcommon/internal/conjugate.h"
00048 #include "SimTKcommon/internal/CompositeNumericalTypes.h"
00049 #include "SimTKcommon/internal/NTraits.h"
00050 #include "SimTKcommon/internal/negator.h"
00051 
00052 #include <complex>
00053 #include <cmath>
00054 #include <climits>
00055 
00056 namespace SimTK {
00057 
00059     // Handy default-precision definitions //
00061 
00062 typedef conjugate<Real> Conjugate;  // like Complex
00063 
00064     // Note that these constants have memory addresses, so you can
00065     // return references to them.
00066     // These are static variables rather than static members to avoid
00067     // problems with static initialization order.
00068 
00069     // Properties of the floaing point representation.
00070 
00071 static const Real& NaN               = NTraits<Real>::getNaN();      // "not a number"
00072 static const Real& Infinity          = NTraits<Real>::getInfinity();
00073 
00074 // Epsilon is the size of roundoff noise; it is the smallest positive number
00075 // such that 1+Eps != 1.
00076 static const Real& Eps               = NTraits<Real>::getEps();         // double ~1e-16, float ~1e-7
00077 static const Real& SqrtEps           = NTraits<Real>::getSqrtEps();     // eps^(1/2): double ~1e- 8, float ~3e-4
00078 static const Real& TinyReal          = NTraits<Real>::getTiny();        // eps^(5/4): double ~1e-20, float ~1e-9
00079 static const Real& SignificantReal   = NTraits<Real>::getSignificant(); // eps^(7/8): double ~1e-14, float ~1e-6
00080 
00081 static const Real& LeastPositiveReal = NTraits<Real>::getLeastPositive(); // double ~1e-308, float ~1e-38
00082 static const Real& MostPositiveReal  = NTraits<Real>::getMostPositive();  // double ~1e+308, float ~1e+38
00083 static const Real& LeastNegativeReal = NTraits<Real>::getLeastNegative();
00084 static const Real& MostNegativeReal  = NTraits<Real>::getMostNegative();
00085 
00086 // This is the number of decimal digits that can be reliably stored and
00087 // retrieved in the default Real precision (typically log10(1/eps)-1).
00088 static const int NumDigitsReal = NTraits<Real>::getNumDigits(); // double ~15, float ~6
00089 
00090 // This is the smallest number of decimal digits you should store in a file
00091 // if you want to be able to get exactly the same bit pattern back when you
00092 // read it in. Typically, this is about log10(1/tiny).
00093 static const int LosslessNumDigitsReal = NTraits<Real>::getLosslessNumDigits(); // double ~20, float ~9
00094 
00095     // Carefully calculated constants, with convenient memory addresses.
00096 
00097 static const Real& Zero         = NTraits<Real>::getZero();
00098 static const Real& One          = NTraits<Real>::getOne();
00099 static const Real& MinusOne     = NTraits<Real>::getMinusOne();
00100 static const Real& Two          = NTraits<Real>::getTwo();
00101 static const Real& Three        = NTraits<Real>::getThree();
00102 
00103 static const Real& OneHalf      = NTraits<Real>::getOneHalf();
00104 static const Real& OneThird     = NTraits<Real>::getOneThird();
00105 static const Real& OneFourth    = NTraits<Real>::getOneFourth();
00106 static const Real& OneFifth     = NTraits<Real>::getOneFifth();
00107 static const Real& OneSixth     = NTraits<Real>::getOneSixth();
00108 static const Real& OneSeventh   = NTraits<Real>::getOneSeventh();
00109 static const Real& OneEighth    = NTraits<Real>::getOneEighth();
00110 static const Real& OneNinth     = NTraits<Real>::getOneNinth();
00111 static const Real& Pi           = NTraits<Real>::getPi();
00112 static const Real& OneOverPi    = NTraits<Real>::getOneOverPi();
00113 static const Real& E            = NTraits<Real>::getE();
00114 static const Real& Log2E        = NTraits<Real>::getLog2E();
00115 static const Real& Log10E       = NTraits<Real>::getLog10E();
00116 static const Real& Sqrt2        = NTraits<Real>::getSqrt2();
00117 static const Real& OneOverSqrt2 = NTraits<Real>::getOneOverSqrt2();  // also sqrt(2)/2
00118 static const Real& Sqrt3        = NTraits<Real>::getSqrt3();
00119 static const Real& OneOverSqrt3 = NTraits<Real>::getOneOverSqrt3();
00120 static const Real& CubeRoot2    = NTraits<Real>::getCubeRoot2();
00121 static const Real& CubeRoot3    = NTraits<Real>::getCubeRoot3();
00122 static const Real& Ln2          = NTraits<Real>::getLn2();
00123 static const Real& Ln10         = NTraits<Real>::getLn10();
00124 
00125 // We only need one complex constant. For the rest just use
00126 // Complex(the Real constant), or if you need an address use
00127 // NTraits<Complex>::getPi(), etc.
00128 static const Complex& I = NTraits<Complex>::getI();
00129 
00130 
00132     // SOME SCALAR UTILITIES //
00134 
00135 
00149 // We use the trick that v & v-1 returns the value that is v with its
00150 // rightmost bit cleared (if it has a rightmost bit set).
00151 inline bool atMostOneBitIsSet(unsigned char v)      {return (v&v-1)==0;}
00152 inline bool atMostOneBitIsSet(unsigned short v)     {return (v&v-1)==0;}
00153 inline bool atMostOneBitIsSet(unsigned int v)       {return (v&v-1)==0;}
00154 inline bool atMostOneBitIsSet(unsigned long v)      {return (v&v-1)==0;}
00155 inline bool atMostOneBitIsSet(unsigned long long v) {return (v&v-1)==0;}
00156 inline bool atMostOneBitIsSet(signed char v)        {return (v&v-1)==0;}
00157 inline bool atMostOneBitIsSet(char v)               {return (v&v-1)==0;}
00158 inline bool atMostOneBitIsSet(short v)              {return (v&v-1)==0;}
00159 inline bool atMostOneBitIsSet(int v)                {return (v&v-1)==0;}
00160 inline bool atMostOneBitIsSet(long v)               {return (v&v-1)==0;}
00161 inline bool atMostOneBitIsSet(long long v)          {return (v&v-1)==0;}
00163 
00178 inline bool exactlyOneBitIsSet(unsigned char v)      {return v && atMostOneBitIsSet(v);}
00179 inline bool exactlyOneBitIsSet(unsigned short v)     {return v && atMostOneBitIsSet(v);}
00180 inline bool exactlyOneBitIsSet(unsigned int v)       {return v && atMostOneBitIsSet(v);}
00181 inline bool exactlyOneBitIsSet(unsigned long v)      {return v && atMostOneBitIsSet(v);}
00182 inline bool exactlyOneBitIsSet(unsigned long long v) {return v && atMostOneBitIsSet(v);}
00183 inline bool exactlyOneBitIsSet(signed char v)        {return v && atMostOneBitIsSet(v);}
00184 inline bool exactlyOneBitIsSet(char v)               {return v && atMostOneBitIsSet(v);}
00185 inline bool exactlyOneBitIsSet(short v)              {return v && atMostOneBitIsSet(v);}
00186 inline bool exactlyOneBitIsSet(int v)                {return v && atMostOneBitIsSet(v);}
00187 inline bool exactlyOneBitIsSet(long v)               {return v && atMostOneBitIsSet(v);}
00188 inline bool exactlyOneBitIsSet(long long v)          {return v && atMostOneBitIsSet(v);}
00190 
00217 inline bool signBit(unsigned char)      {return 0;}
00218 inline bool signBit(unsigned short)     {return 0;}
00219 inline bool signBit(unsigned int)       {return 0;}
00220 inline bool signBit(unsigned long)      {return 0;}
00221 inline bool signBit(unsigned long long) {return 0;}
00222 
00223 // Note that plain 'char' type is not overloaded -- see above.
00224 
00225 // We're assuming sizeof(char)==1, short==2, int==4, long long==8 which is safe
00226 // for all our anticipated platforms. But some 64 bit implementations have
00227 // sizeof(long)==4 while others have sizeof(long)==8 so we'll use the ANSI
00228 // standard value LONG_MIN which is a long value with just the high bit set.
00229 // (We're assuming two's complement integers here; I haven't seen anything
00230 // else in decades.)
00231 inline bool signBit(signed char i) {return ((unsigned char)i      & 0x80U) != 0;}
00232 inline bool signBit(short       i) {return ((unsigned short)i     & 0x8000U) != 0;}
00233 inline bool signBit(int         i) {return ((unsigned int)i       & 0x80000000U) != 0;}
00234 inline bool signBit(long long   i) {return ((unsigned long long)i & 0x8000000000000000ULL) != 0;}
00235 
00236 inline bool signBit(long        i) {return ((unsigned long)i
00237                                             & (unsigned long)LONG_MIN) != 0;}
00238 
00239 inline bool signBit(const float&  f) {return std::signbit(f);}
00240 inline bool signBit(const double& d) {return std::signbit(d);}
00241 inline bool signBit(const negator<float>&  nf) {return std::signbit(-nf);} // !!
00242 inline bool signBit(const negator<double>& nd) {return std::signbit(-nd);} // !!
00244 
00257 inline unsigned int sign(unsigned char      u) {return u==0 ? 0 : 1;}
00258 inline unsigned int sign(unsigned short     u) {return u==0 ? 0 : 1;}
00259 inline unsigned int sign(unsigned int       u) {return u==0 ? 0 : 1;}
00260 inline unsigned int sign(unsigned long      u) {return u==0 ? 0 : 1;}
00261 inline unsigned int sign(unsigned long long u) {return u==0 ? 0 : 1;}
00262 
00263 // Don't overload for plain "char" because it may be signed or unsigned
00264 // depending on the compiler.
00265 
00266 inline int sign(signed char i) {return i>0 ? 1 : (i<0 ? -1 : 0);}
00267 inline int sign(short       i) {return i>0 ? 1 : (i<0 ? -1 : 0);}
00268 inline int sign(int         i) {return i>0 ? 1 : (i<0 ? -1 : 0);}
00269 inline int sign(long        i) {return i>0 ? 1 : (i<0 ? -1 : 0);}
00270 inline int sign(long long   i) {return i>0 ? 1 : (i<0 ? -1 : 0);}
00271 
00272 inline int sign(const float&       x) {return x>0 ? 1 : (x<0 ? -1 : 0);}
00273 inline int sign(const double&      x) {return x>0 ? 1 : (x<0 ? -1 : 0);}
00274 inline int sign(const long double& x) {return x>0 ? 1 : (x<0 ? -1 : 0);}
00275 
00276 inline int sign(const negator<float>&       x) {return -sign(-x);} // -x is free
00277 inline int sign(const negator<double>&      x) {return -sign(-x);}
00278 inline int sign(const negator<long double>& x) {return -sign(-x);}
00280 
00297 inline unsigned char  square(unsigned char  u) {return u*u;}
00298 inline unsigned short square(unsigned short u) {return u*u;}
00299 inline unsigned int   square(unsigned int   u) {return u*u;}
00300 inline unsigned long  square(unsigned long  u) {return u*u;}
00301 
00302 inline char        square(char c) {return c*c;}
00303 
00304 inline signed char square(signed char i) {return i*i;}
00305 inline short       square(short       i) {return i*i;}
00306 inline int         square(int         i) {return i*i;}
00307 inline long        square(long        i) {return i*i;}
00308 
00309 inline float       square(const float&       x) {return x*x;}
00310 inline double      square(const double&      x) {return x*x;}
00311 inline long double square(const long double& x) {return x*x;}
00312 
00313 // Negation is free for negators, so we can square them and clean
00314 // them up at the same time at no extra cost.
00315 inline float       square(const negator<float>&       x) {return square(-x);}
00316 inline double      square(const negator<double>&      x) {return square(-x);}
00317 inline long double square(const negator<long double>& x) {return square(-x);}
00318 
00319 // It is safer to templatize using complex classes, and doesn't make
00320 // debugging any worse since complex is already templatized. 
00321 // 5 flops vs. 6 for general complex multiply.
00322 template <class P> inline 
00323 std::complex<P> square(const std::complex<P>& x) {
00324     const P re=x.real(), im=x.imag();
00325     return std::complex<P>(re*re-im*im, 2*re*im);
00326 }
00327 
00328 // We can square a conjugate and clean it up back to complex at
00329 // the same time at zero cost (or maybe 1 negation depending
00330 // on how the compiler handles the "-2" below).
00331 template <class P> inline 
00332 std::complex<P> square(const conjugate<P>& x) {
00333     const P re=x.real(), nim=x.negImag();
00334     return std::complex<P>(re*re-nim*nim, -2*re*nim);
00335 }
00336 
00337 template <class P> inline
00338 std::complex<P> square(const negator< std::complex<P> >& x) {
00339     return square(-x); // negation is free for negators
00340 }
00341 
00342 // Note return type here after squaring negator<conjugate>
00343 // is complex, not conjugate.
00344 template <class P> inline
00345 std::complex<P> square(const negator< conjugate<P> >& x) {
00346     return square(-x); // negation is free for negators
00347 }
00349 
00368 inline unsigned char  cube(unsigned char  u) {return u*u*u;}
00369 inline unsigned short cube(unsigned short u) {return u*u*u;}
00370 inline unsigned int   cube(unsigned int   u) {return u*u*u;}
00371 inline unsigned long  cube(unsigned long  u) {return u*u*u;}
00372 
00373 inline char        cube(char c) {return c*c*c;}
00374 
00375 inline signed char cube(signed char i) {return i*i*i;}
00376 inline short       cube(short       i) {return i*i*i;}
00377 inline int         cube(int         i) {return i*i*i;}
00378 inline long        cube(long        i) {return i*i*i;}
00379 
00380 inline float       cube(const float&       x) {return x*x*x;}
00381 inline double      cube(const double&      x) {return x*x*x;}
00382 inline long double cube(const long double& x) {return x*x*x;}
00383 
00384 // To keep this cheap we'll defer getting rid of the negator<> until
00385 // some other operation. We cube -x and then recast that to negator<>
00386 // on return, for a total cost of 2 flops.
00387 inline negator<float> cube(const negator<float>& x) {
00388     return negator<float>::recast(cube(-x));
00389 }
00390 inline negator<double> cube(const negator<double>& x) {
00391     return negator<double>::recast(cube(-x));
00392 }
00393 inline negator<long double> cube(const negator<long double>& x) {
00394     return negator<long double>::recast(cube(-x));
00395 }
00396 
00397 // Cubing a complex this way is cheaper than doing it by
00398 // multiplication. Cost here is 8 flops vs. 11 for a square
00399 // followed by a multiply.
00400 template <class P> inline
00401 std::complex<P> cube(const std::complex<P>& x) {
00402     const P re=x.real(), im=x.imag(), rr=re*re, ii=im*im;
00403     return std::complex<P>(re*(rr-3*ii), im*(3*rr-ii));
00404 }
00405 
00406 // Cubing a negated complex allows us to cube and eliminate the
00407 // negator at the same time for zero extra cost. Compare the 
00408 // expressions here to the normal cube above to see the free
00409 // sign changes in both parts. 8 flops here.
00410 template <class P> inline
00411 std::complex<P> cube(const negator< std::complex<P> >& x) {
00412     // -x is free for a negator
00413     const P nre=(-x).real(), nim=(-x).imag(), rr=nre*nre, ii=nim*nim;
00414     return std::complex<P>(nre*(3*ii-rr), nim*(ii-3*rr));
00415 }
00416 
00417 // Cubing a conjugate this way saves a lot over multiplying, and
00418 // also lets us convert the result to complex for free.
00419 template <class P> inline
00420 std::complex<P> cube(const conjugate<P>& x) {
00421     const P re=x.real(), nim=x.negImag(), rr=re*re, ii=nim*nim;
00422     return std::complex<P>(re*(rr-3*ii), nim*(ii-3*rr));
00423 }
00424 
00425 
00426 // Cubing a negated conjugate this way saves a lot over multiplying, and
00427 // also lets us convert the result to non-negated complex for free.
00428 template <class P> inline
00429 std::complex<P> cube(const negator< conjugate<P> >& x) {
00430     // -x is free for a negator
00431     const P nre=(-x).real(), im=(-x).negImag(), rr=nre*nre, ii=im*im;
00432     return std::complex<P>(nre*(3*ii-rr), im*(3*rr-ii));
00433 }
00435 
00436 } // namespace SimTK
00437 
00438 #endif //SimTK_SIMMATRIX_SCALAR_H_

Generated by  doxygen 1.6.2