eop_aux.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2010 NICTA and the authors listed below
00002 // http://nicta.com.au
00003 // 
00004 // Authors:
00005 // - Conrad Sanderson (conradsand at ieee dot org)
00006 // 
00007 // This file is part of the Armadillo C++ library.
00008 // It is provided without any warranty of fitness
00009 // for any purpose. You can redistribute this file
00010 // and/or modify it under the terms of the GNU
00011 // Lesser General Public License (LGPL) as published
00012 // by the Free Software Foundation, either version 3
00013 // of the License or (at your option) any later version.
00014 // (see http://www.opensource.org/licenses for more info)
00015 
00016 
00017 //! \addtogroup eop_aux
00018 //! @{
00019 
00020 
00021 
00022 template<typename eT>
00023 struct eop_aux_rand
00024   {
00025   arma_inline
00026   operator eT ()
00027     {
00028     return eT(std::rand()) / eT(RAND_MAX);
00029     }
00030   };
00031 
00032 
00033   
00034 template<typename T>
00035 struct eop_aux_rand< std::complex<T> >
00036   {
00037   arma_inline
00038   operator std::complex<T> ()
00039     {
00040     return std::complex<T>( T(eop_aux_rand<T>()), T(eop_aux_rand<T>()) );
00041     }
00042   };
00043 
00044 
00045 
00046 template<typename eT>
00047 struct eop_aux_randn
00048   {
00049   // // rudimentary method, based on the central limit theorem
00050   // // http://en.wikipedia.org/wiki/Central_limit_theorem
00051   // inline
00052   // operator eT () const
00053   //   {
00054   //   const u32 N  = 12;  // N must be >= 12 and an even number
00055   //   const u32 N2 = N/2;
00056   //   
00057   //   eT acc = eT(0);
00058   //   
00059   //   for(u32 i=0; i<N2; ++i)
00060   //     {
00061   //     const eT tmp1 = eT(std::rand()) / eT(RAND_MAX);
00062   //     const eT tmp2 = eT(std::rand()) / eT(RAND_MAX);
00063   //     acc += tmp1+tmp2;
00064   //     }
00065   //   
00066   //   return acc - eT(N2);
00067   //   }
00068   
00069   
00070   // polar form of the Box-Muller transformation
00071   // http://en.wikipedia.org/wiki/Box-Muller_transformation
00072   // http://en.wikipedia.org/wiki/Marsaglia_polar_method
00073   inline
00074   operator eT () const
00075     {
00076     // make sure we are internally using at least floats
00077     typedef typename promote_type<eT,float>::result eTp;
00078     
00079     eTp tmp1;
00080     eTp tmp2;
00081     eTp w;
00082     
00083     do
00084       {
00085       tmp1 = eTp(2) * eTp(std::rand()) / eTp(RAND_MAX) - eTp(1);
00086       tmp2 = eTp(2) * eTp(std::rand()) / eTp(RAND_MAX) - eTp(1);
00087       w = tmp1*tmp1 + tmp2*tmp2;
00088       }
00089     while ( w >= eTp(1) );
00090     
00091     return eT( tmp1 * std::sqrt( (eTp(-2) * std::log(w)) / w) );
00092     }
00093   
00094   
00095   // other methods:
00096   // http://en.wikipedia.org/wiki/Ziggurat_algorithm
00097   //
00098   // Marsaglia and Tsang Ziggurat technique to transform from a uniform to a normal distribution.
00099   // G. Marsaglia, W.W. Tsang.
00100   // "Ziggurat method for generating random variables",
00101   // J. Statistical Software, vol 5, 2000.
00102   // http://www.jstatsoft.org/v05/i08/
00103   };
00104 
00105 
00106 
00107 template<typename T>
00108 struct eop_aux_randn< std::complex<T> >
00109   {
00110   arma_inline
00111   operator std::complex<T> () const
00112     {
00113     return std::complex<T>( T(eop_aux_randn<T>()), T(eop_aux_randn<T>()) );
00114     }
00115 
00116   };
00117 
00118 
00119 
00120 class eop_aux
00121   {
00122   public:
00123   
00124   #if defined(ARMA_USE_BOOST)
00125     #define arma_boost_wrap(trig_fn, val) ( (boost::math::trig_fn)(val) )
00126   #else
00127     #define arma_boost_wrap(trig_fn, val) ( arma_stop( #trig_fn "(): need Boost libraries" ), val )
00128   #endif
00129   
00130   template<typename eT> arma_inline static typename arma_integral_only<eT>::result acos (const eT& x) { return std::acos(double(x)); }
00131   template<typename eT> arma_inline static typename arma_integral_only<eT>::result asin (const eT& x) { return std::asin(double(x)); }
00132   template<typename eT> arma_inline static typename arma_integral_only<eT>::result atan (const eT& x) { return std::atan(double(x)); }
00133   
00134   template<typename eT> arma_inline static typename arma_float_only<eT>::result acos (const eT& x) { return std::acos(x); }
00135   template<typename eT> arma_inline static typename arma_float_only<eT>::result asin (const eT& x) { return std::asin(x); }
00136   template<typename eT> arma_inline static typename arma_float_only<eT>::result atan (const eT& x) { return std::atan(x); }
00137   
00138   template<typename  T> arma_inline static std::complex<T> acos (const std::complex<T>& x) { return arma_boost_wrap(acos,  x); }
00139   template<typename  T> arma_inline static std::complex<T> asin (const std::complex<T>& x) { return arma_boost_wrap(asin,  x); }
00140   template<typename  T> arma_inline static std::complex<T> atan (const std::complex<T>& x) { return arma_boost_wrap(atan,  x); }
00141 
00142   template<typename eT> arma_inline static typename arma_integral_only<eT>::result acosh (const eT& x) { return arma_boost_wrap(acosh, double(x)); }
00143   template<typename eT> arma_inline static typename arma_integral_only<eT>::result asinh (const eT& x) { return arma_boost_wrap(asinh, double(x)); }
00144   template<typename eT> arma_inline static typename arma_integral_only<eT>::result atanh (const eT& x) { return arma_boost_wrap(atanh, double(x)); }
00145 
00146   template<typename eT> arma_inline static typename arma_float_only<eT>::result acosh (const eT& x) { return arma_boost_wrap(acosh, x); }
00147   template<typename eT> arma_inline static typename arma_float_only<eT>::result asinh (const eT& x) { return arma_boost_wrap(asinh, x); }
00148   template<typename eT> arma_inline static typename arma_float_only<eT>::result atanh (const eT& x) { return arma_boost_wrap(atanh, x); }
00149   
00150   template<typename  T> arma_inline static std::complex<T> acosh (const std::complex<T>& x) { return arma_boost_wrap(acosh, x); }
00151   template<typename  T> arma_inline static std::complex<T> asinh (const std::complex<T>& x) { return arma_boost_wrap(asinh, x); }
00152   template<typename  T> arma_inline static std::complex<T> atanh (const std::complex<T>& x) { return arma_boost_wrap(atanh, x); }
00153   
00154   #undef arma_boost_wrap
00155   
00156   template<typename eT> arma_inline static eT              conj(const eT&              x) { return x;            }
00157   template<typename  T> arma_inline static std::complex<T> conj(const std::complex<T>& x) { return std::conj(x); }
00158   
00159   
00160   template<typename T1, typename T2> arma_inline static
00161   T1    pow(const T1    base, const T2 exponent) { return std::pow(base, exponent); }
00162   
00163   template<typename T2> arma_inline static
00164   char  pow(const char  base, const T2 exponent) { typedef char  out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00165   
00166   template<typename T2> arma_inline static
00167   short pow(const short base, const T2 exponent) { typedef short out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00168   
00169   template<typename T2> arma_inline static
00170   int   pow(const int   base, const T2 exponent) { typedef int   out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00171   
00172   template<typename T2> arma_inline static
00173   long  pow(const long  base, const T2 exponent) { typedef long  out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00174   
00175   
00176   template<typename T2> arma_inline static
00177   unsigned char  pow(const unsigned char  base, const T2 exponent) { typedef unsigned char  out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00178   
00179   template<typename T2> arma_inline static
00180   unsigned short pow(const unsigned short base, const T2 exponent) { typedef unsigned short out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00181   
00182   template<typename T2> arma_inline static
00183   unsigned int   pow(const unsigned int   base, const T2 exponent) { typedef unsigned int   out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00184   
00185   template<typename T2> arma_inline static
00186   unsigned long  pow(const unsigned long  base, const T2 exponent) { typedef unsigned long  out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00187 
00188   
00189 
00190   template<typename T1> arma_inline static
00191   T1 pow_int(const    T1    base, const int exponent) { return std::pow(base, exponent); }
00192   
00193   arma_inline static
00194   char  pow_int(const char  base, const int exponent) { typedef char  out_t; return out_t( std::pow(double(base), exponent) ); }
00195   
00196   arma_inline static
00197   short pow_int(const short base, const int exponent) { typedef short out_t; return out_t( std::pow(double(base), exponent) ); }
00198   
00199   arma_inline static
00200   int   pow_int(const int   base, const int exponent) { typedef int   out_t; return out_t( std::pow(double(base), exponent) ); }
00201   
00202   arma_inline static
00203   long  pow_int(const long  base, const int exponent) { typedef long  out_t; return out_t( std::pow(double(base), exponent) ); }
00204   
00205   
00206   arma_inline static
00207   unsigned char  pow_int(const unsigned char  base, const int exponent) { typedef unsigned char  out_t; return out_t( std::pow(double(base), exponent) ); }
00208 
00209   arma_inline static
00210   unsigned short pow_int(const unsigned short base, const int exponent) { typedef unsigned short out_t; return out_t( std::pow(double(base), exponent) ); }
00211   
00212   arma_inline static
00213   unsigned int   pow_int(const unsigned int   base, const int exponent) { typedef unsigned int   out_t; return out_t( std::pow(double(base), exponent) ); }
00214   
00215   arma_inline static
00216   unsigned long  pow_int(const unsigned long  base, const int exponent) { typedef unsigned long  out_t; return out_t( std::pow(double(base), exponent) ); }
00217   
00218   
00219   
00220   
00221   template<typename eT>
00222   inline
00223   static
00224   eT
00225   trunc_exp(const eT x)
00226     {
00227     if(std::numeric_limits<eT>::is_iec559 && (x >= Math<eT>::log_max() ))
00228       {
00229       return std::numeric_limits<eT>::max();
00230       }
00231     else
00232       {
00233       return std::exp(x);
00234       }
00235     }
00236   
00237   
00238   
00239   template<typename T>
00240   inline
00241   static
00242   std::complex<T>
00243   trunc_exp(const std::complex<T>& x)
00244     {
00245     return std::exp(x);
00246     }
00247   
00248   
00249   
00250   template<typename eT>
00251   inline 
00252   static
00253   eT
00254   trunc_log(const eT x)
00255     {
00256     if(std::numeric_limits<eT>::is_iec559)
00257       {
00258       if(x == std::numeric_limits<eT>::infinity())
00259         {
00260         return Math<eT>::log_max();
00261         }
00262       else
00263       if(x <= eT(0))
00264         {
00265         return Math<eT>::log_min();
00266         }
00267       else
00268         {
00269         return std::log(x);
00270         }
00271       }
00272     else
00273       {
00274       return std::log(x);
00275       }
00276     }
00277   
00278   
00279   
00280   template<typename T>
00281   inline 
00282   static
00283   std::complex<T>
00284   trunc_log(const std::complex<T>& x)
00285     {
00286     return std::log(x);
00287     }
00288   
00289   
00290   
00291   template<typename eT>
00292   arma_inline
00293   static
00294   typename arma_integral_only<eT>::result
00295   direct_eps(const eT& x)
00296     {
00297     return eT(0);
00298     }
00299   
00300   
00301   
00302   template<typename eT>
00303   inline
00304   static
00305   typename arma_float_only<eT>::result
00306   direct_eps(const eT& x)
00307     {
00308     //arma_extra_debug_sigprint();
00309     
00310     // acording to IEEE Standard for Floating-Point Arithmetic (IEEE 754)
00311     // the mantissa length for double is 53 bits = std::numeric_limits<double>::digits
00312     // the mantissa length for float  is 24 bits = std::numeric_limits<float >::digits
00313     
00314     //return std::pow( std::numeric_limits<eT>::radix, (std::floor(std::log10(std::abs(x))/std::log10(std::numeric_limits<eT>::radix))-(std::numeric_limits<eT>::digits-1)) );
00315     
00316     const eT radix_eT     = eT(std::numeric_limits<eT>::radix);
00317     const eT digits_m1_eT = eT(std::numeric_limits<eT>::digits - 1);
00318     
00319     // return std::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) );
00320     return eop_aux::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) );
00321     }
00322   
00323   
00324   
00325   template<typename T>
00326   inline
00327   static
00328   typename arma_float_only<T>::result
00329   direct_eps(const std::complex<T>& x)
00330     {
00331     //arma_extra_debug_sigprint();
00332     
00333     //return std::pow( std::numeric_limits<T>::radix, (std::floor(std::log10(std::abs(x))/std::log10(std::numeric_limits<T>::radix))-(std::numeric_limits<T>::digits-1)) );
00334     
00335     const T radix_T     = T(std::numeric_limits<T>::radix);
00336     const T digits_m1_T = T(std::numeric_limits<T>::digits - 1);
00337     
00338     return std::pow( radix_T, T(std::floor(std::log10(std::abs(x))/std::log10(radix_T)) - digits_m1_T) );
00339     }
00340   
00341   
00342   
00343   //! work around a bug in GCC 4.4
00344   template<typename eT> arma_inline static
00345   typename arma_unsigned_integral_only<eT>::result arma_abs(const eT& x)              { return x;           }
00346   
00347   template<typename eT> arma_inline static
00348   typename arma_signed_integral_only<eT>::result   arma_abs(const eT& x)              { return std::abs(x); }
00349   
00350   template<typename eT> arma_inline static
00351   typename arma_float_only<eT>::result             arma_abs(const eT& x)              { return std::abs(x); }
00352   
00353   template<typename T> arma_inline static
00354   typename arma_float_only<T>::result              arma_abs(const std::complex<T>& x) { return std::abs(x); }
00355   
00356   
00357   
00358   template<typename eT, typename eop_type>
00359   arma_inline
00360   static
00361   eT
00362   generate()
00363     {
00364          if(is_same_type<eop_type, eop_rand          >::value == true) { return eT(eop_aux_rand<eT>());  }
00365     else if(is_same_type<eop_type, eop_randn         >::value == true) { return eT(eop_aux_randn<eT>()); }
00366     else if(is_same_type<eop_type, eop_zeros         >::value == true) { return eT(0);                   }
00367     else if(is_same_type<eop_type, eop_ones_full     >::value == true) { return eT(1);                   }
00368     else if(is_same_type<eop_type, eop_cube_rand     >::value == true) { return eT(eop_aux_rand<eT>());  }
00369     else if(is_same_type<eop_type, eop_cube_randn    >::value == true) { return eT(eop_aux_randn<eT>()); }
00370     else if(is_same_type<eop_type, eop_cube_zeros    >::value == true) { return eT(0);                   }
00371     else if(is_same_type<eop_type, eop_cube_ones_full>::value == true) { return eT(1);                   }
00372     else                                                               { return eT(0);                   }
00373     }
00374   
00375   };
00376 
00377 
00378 
00379 //! @}
00380