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   // TODO: implement a more efficient method for randn()
00050   //
00051   // possible option:
00052   // Marsaglia and Tsang Ziggurat technique to transform from a uniform to a normal distribution.
00053   // G. Marsaglia, W.W. Tsang.
00054   // "Ziggurat method for generating random variables",
00055   // J. Statistical Software, vol 5, 2000.
00056   // http://www.jstatsoft.org/v05/i08/
00057   
00058   inline
00059   operator eT () const
00060     {
00061     const u32 N  = 12;  // N must be >= 12 and an even number
00062     const u32 N2 = N/2;
00063     
00064     eT acc = eT(0);
00065     
00066     for(u32 i=0; i<N2; ++i)
00067       {
00068       const eT tmp1 = eT(std::rand()) / eT(RAND_MAX);
00069       const eT tmp2 = eT(std::rand()) / eT(RAND_MAX);
00070       acc += tmp1+tmp2;
00071       }
00072     
00073     return acc - eT(N2);
00074     }
00075   
00076   };
00077 
00078 
00079 
00080 template<typename T>
00081 struct eop_aux_randn< std::complex<T> >
00082   {
00083   arma_inline
00084   operator std::complex<T> () const
00085     {
00086     return std::complex<T>( T(eop_aux_randn<T>()), T(eop_aux_randn<T>()) );
00087     }
00088 
00089   };
00090 
00091 
00092 
00093 class eop_aux
00094   {
00095   public:
00096   
00097   #if defined(ARMA_USE_BOOST)
00098     #define arma_boost_wrap(trig_fn, val) ( (boost::math::trig_fn)(val) )
00099   #else
00100     #define arma_boost_wrap(trig_fn, val) ( arma_stop( #trig_fn "(): need Boost libraries" ), val )
00101   #endif
00102   
00103   template<typename eT> arma_inline static typename arma_integral_only<eT>::result acos (const eT& x) { return std::acos(double(x)); }
00104   template<typename eT> arma_inline static typename arma_integral_only<eT>::result asin (const eT& x) { return std::asin(double(x)); }
00105   template<typename eT> arma_inline static typename arma_integral_only<eT>::result atan (const eT& x) { return std::atan(double(x)); }
00106   
00107   template<typename eT> arma_inline static typename arma_float_only<eT>::result acos (const eT& x) { return std::acos(x); }
00108   template<typename eT> arma_inline static typename arma_float_only<eT>::result asin (const eT& x) { return std::asin(x); }
00109   template<typename eT> arma_inline static typename arma_float_only<eT>::result atan (const eT& x) { return std::atan(x); }
00110   
00111   template<typename  T> arma_inline static std::complex<T> acos (const std::complex<T>& x) { return arma_boost_wrap(acos,  x); }
00112   template<typename  T> arma_inline static std::complex<T> asin (const std::complex<T>& x) { return arma_boost_wrap(asin,  x); }
00113   template<typename  T> arma_inline static std::complex<T> atan (const std::complex<T>& x) { return arma_boost_wrap(atan,  x); }
00114 
00115   template<typename eT> arma_inline static typename arma_integral_only<eT>::result acosh (const eT& x) { return arma_boost_wrap(acosh, double(x)); }
00116   template<typename eT> arma_inline static typename arma_integral_only<eT>::result asinh (const eT& x) { return arma_boost_wrap(asinh, double(x)); }
00117   template<typename eT> arma_inline static typename arma_integral_only<eT>::result atanh (const eT& x) { return arma_boost_wrap(atanh, double(x)); }
00118 
00119   template<typename eT> arma_inline static typename arma_float_only<eT>::result acosh (const eT& x) { return arma_boost_wrap(acosh, x); }
00120   template<typename eT> arma_inline static typename arma_float_only<eT>::result asinh (const eT& x) { return arma_boost_wrap(asinh, x); }
00121   template<typename eT> arma_inline static typename arma_float_only<eT>::result atanh (const eT& x) { return arma_boost_wrap(atanh, x); }
00122   
00123   template<typename  T> arma_inline static std::complex<T> acosh (const std::complex<T>& x) { return arma_boost_wrap(acosh, x); }
00124   template<typename  T> arma_inline static std::complex<T> asinh (const std::complex<T>& x) { return arma_boost_wrap(asinh, x); }
00125   template<typename  T> arma_inline static std::complex<T> atanh (const std::complex<T>& x) { return arma_boost_wrap(atanh, x); }
00126   
00127   #undef arma_boost_wrap
00128   
00129   template<typename eT> arma_inline static eT              conj(const eT&              x) { return x;            }
00130   template<typename  T> arma_inline static std::complex<T> conj(const std::complex<T>& x) { return std::conj(x); }
00131   
00132   
00133   template<typename T1, typename T2> arma_inline static
00134   T1    pow(const T1    base, const T2 exponent) { return std::pow(base, exponent); }
00135   
00136   template<typename T2> arma_inline static
00137   char  pow(const char  base, const T2 exponent) { typedef char  out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00138   
00139   template<typename T2> arma_inline static
00140   short pow(const short base, const T2 exponent) { typedef short out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00141   
00142   template<typename T2> arma_inline static
00143   int   pow(const int   base, const T2 exponent) { typedef int   out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00144   
00145   template<typename T2> arma_inline static
00146   long  pow(const long  base, const T2 exponent) { typedef long  out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00147   
00148   
00149   template<typename T2> arma_inline static
00150   unsigned char  pow(const unsigned char  base, const T2 exponent) { typedef unsigned char  out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00151   
00152   template<typename T2> arma_inline static
00153   unsigned short pow(const unsigned short base, const T2 exponent) { typedef unsigned short out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00154   
00155   template<typename T2> arma_inline static
00156   unsigned int   pow(const unsigned int   base, const T2 exponent) { typedef unsigned int   out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00157   
00158   template<typename T2> arma_inline static
00159   unsigned long  pow(const unsigned long  base, const T2 exponent) { typedef unsigned long  out_t; return out_t( std::pow(double(base), double(exponent) ) ); }
00160 
00161   
00162 
00163   template<typename T1> arma_inline static
00164   T1 pow_int(const    T1    base, const int exponent) { return std::pow(base, exponent); }
00165   
00166   arma_inline static
00167   char  pow_int(const char  base, const int exponent) { typedef char  out_t; return out_t( std::pow(double(base), exponent) ); }
00168   
00169   arma_inline static
00170   short pow_int(const short base, const int exponent) { typedef short out_t; return out_t( std::pow(double(base), exponent) ); }
00171   
00172   arma_inline static
00173   int   pow_int(const int   base, const int exponent) { typedef int   out_t; return out_t( std::pow(double(base), exponent) ); }
00174   
00175   arma_inline static
00176   long  pow_int(const long  base, const int exponent) { typedef long  out_t; return out_t( std::pow(double(base), exponent) ); }
00177   
00178   
00179   arma_inline static
00180   unsigned char  pow_int(const unsigned char  base, const int exponent) { typedef unsigned char  out_t; return out_t( std::pow(double(base), exponent) ); }
00181 
00182   arma_inline static
00183   unsigned short pow_int(const unsigned short base, const int exponent) { typedef unsigned short out_t; return out_t( std::pow(double(base), exponent) ); }
00184   
00185   arma_inline static
00186   unsigned int   pow_int(const unsigned int   base, const int exponent) { typedef unsigned int   out_t; return out_t( std::pow(double(base), exponent) ); }
00187   
00188   arma_inline static
00189   unsigned long  pow_int(const unsigned long  base, const int exponent) { typedef unsigned long  out_t; return out_t( std::pow(double(base), exponent) ); }
00190   
00191   
00192   
00193   
00194   template<typename eT>
00195   inline
00196   static
00197   eT
00198   trunc_exp(const eT x)
00199     {
00200     if(std::numeric_limits<eT>::is_iec559 && (x >= Math<eT>::log_max() ))
00201       {
00202       return std::numeric_limits<eT>::max();
00203       }
00204     else
00205       {
00206       return std::exp(x);
00207       }
00208     }
00209   
00210   
00211   
00212   template<typename T>
00213   inline
00214   static
00215   std::complex<T>
00216   trunc_exp(const std::complex<T>& x)
00217     {
00218     return std::exp(x);
00219     }
00220   
00221   
00222   
00223   template<typename eT>
00224   inline 
00225   static
00226   eT
00227   trunc_log(const eT x)
00228     {
00229     if(std::numeric_limits<eT>::is_iec559)
00230       {
00231       if(x == std::numeric_limits<eT>::infinity())
00232         {
00233         return Math<eT>::log_max();
00234         }
00235       else
00236       if(x <= eT(0))
00237         {
00238         return Math<eT>::log_min();
00239         }
00240       else
00241         {
00242         return std::log(x);
00243         }
00244       }
00245     else
00246       {
00247       return std::log(x);
00248       }
00249     }
00250   
00251   
00252   
00253   template<typename T>
00254   inline 
00255   static
00256   std::complex<T>
00257   trunc_log(const std::complex<T>& x)
00258     {
00259     return std::log(x);
00260     }
00261   
00262   
00263   
00264   template<typename eT>
00265   arma_inline
00266   static
00267   typename arma_integral_only<eT>::result
00268   direct_eps(const eT& x)
00269     {
00270     return eT(0);
00271     }
00272   
00273   
00274   
00275   template<typename eT>
00276   inline
00277   static
00278   typename arma_float_only<eT>::result
00279   direct_eps(const eT& x)
00280     {
00281     //arma_extra_debug_sigprint();
00282     
00283     // acording to IEEE Standard for Floating-Point Arithmetic (IEEE 754)
00284     // the mantissa length for double is 53 bits = std::numeric_limits<double>::digits
00285     // the mantissa length for float  is 24 bits = std::numeric_limits<float >::digits
00286     
00287     //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)) );
00288     
00289     const eT radix_eT     = eT(std::numeric_limits<eT>::radix);
00290     const eT digits_m1_eT = eT(std::numeric_limits<eT>::digits - 1);
00291     
00292     // return std::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) );
00293     return eop_aux::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) );
00294     }
00295   
00296   
00297   
00298   template<typename T>
00299   inline
00300   static
00301   typename arma_float_only<T>::result
00302   direct_eps(const std::complex<T>& x)
00303     {
00304     //arma_extra_debug_sigprint();
00305     
00306     //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)) );
00307     
00308     const T radix_T     = T(std::numeric_limits<T>::radix);
00309     const T digits_m1_T = T(std::numeric_limits<T>::digits - 1);
00310     
00311     return std::pow( radix_T, T(std::floor(std::log10(std::abs(x))/std::log10(radix_T)) - digits_m1_T) );
00312     }
00313   
00314   
00315   
00316   //! work around a bug in GCC 4.4
00317   template<typename eT> arma_inline static
00318   typename arma_unsigned_integral_only<eT>::result arma_abs(const eT& x)              { return x;           }
00319   
00320   template<typename eT> arma_inline static
00321   typename arma_signed_integral_only<eT>::result   arma_abs(const eT& x)              { return std::abs(x); }
00322   
00323   template<typename eT> arma_inline static
00324   typename arma_float_only<eT>::result             arma_abs(const eT& x)              { return std::abs(x); }
00325   
00326   template<typename T> arma_inline static
00327   typename arma_float_only<T>::result              arma_abs(const std::complex<T>& x) { return std::abs(x); }
00328   
00329   
00330   
00331   template<typename eT, typename eop_type>
00332   arma_inline
00333   static
00334   eT
00335   generate()
00336     {
00337          if(is_same_type<eop_type, eop_rand          >::value == true) { return eT(eop_aux_rand<eT>());  }
00338     else if(is_same_type<eop_type, eop_randn         >::value == true) { return eT(eop_aux_randn<eT>()); }
00339     else if(is_same_type<eop_type, eop_zeros         >::value == true) { return eT(0);                   }
00340     else if(is_same_type<eop_type, eop_ones_full     >::value == true) { return eT(1);                   }
00341     else if(is_same_type<eop_type, eop_cube_rand     >::value == true) { return eT(eop_aux_rand<eT>());  }
00342     else if(is_same_type<eop_type, eop_cube_randn    >::value == true) { return eT(eop_aux_randn<eT>()); }
00343     else if(is_same_type<eop_type, eop_cube_zeros    >::value == true) { return eT(0);                   }
00344     else if(is_same_type<eop_type, eop_cube_ones_full>::value == true) { return eT(1);                   }
00345     else                                                               { return eT(0);                   }
00346     }
00347   
00348   };
00349 
00350 
00351 
00352 //! @}
00353