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_randu
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_randu< std::complex<T> >
00036   {
00037   arma_inline
00038   operator std::complex<T> ()
00039     {
00040     return std::complex<T>( T(eop_aux_randu<T>()), T(eop_aux_randu<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 //! use of the SFINAE approach to work around compiler limitations
00120 //! http://en.wikipedia.org/wiki/SFINAE
00121 
00122 class eop_aux
00123   {
00124   public:
00125   
00126   template<typename eT> arma_inline static typename arma_integral_only<eT>::result    acos  (const eT& x) { return std::acos(double(x)); }
00127   template<typename eT> arma_inline static typename arma_integral_only<eT>::result    asin  (const eT& x) { return std::asin(double(x)); }
00128   template<typename eT> arma_inline static typename arma_integral_only<eT>::result    atan  (const eT& x) { return std::atan(double(x)); }
00129   
00130   template<typename eT> arma_inline static typename arma_float_only<eT>::result       acos  (const eT& x) { return std::acos(x); }
00131   template<typename eT> arma_inline static typename arma_float_only<eT>::result       asin  (const eT& x) { return std::asin(x); }
00132   template<typename eT> arma_inline static typename arma_float_only<eT>::result       atan  (const eT& x) { return std::atan(x); }
00133   
00134   template<typename eT> arma_inline static typename arma_cx_only<eT>::result          acos  (const eT& x) { return arma_acos(x); }
00135   template<typename eT> arma_inline static typename arma_cx_only<eT>::result          asin  (const eT& x) { return arma_asin(x); }
00136   template<typename eT> arma_inline static typename arma_cx_only<eT>::result          atan  (const eT& x) { return arma_atan(x); }
00137   
00138   template<typename eT> arma_inline static typename arma_integral_only<eT>::result    acosh (const eT& x) { return arma_acosh(double(x)); }
00139   template<typename eT> arma_inline static typename arma_integral_only<eT>::result    asinh (const eT& x) { return arma_asinh(double(x)); }
00140   template<typename eT> arma_inline static typename arma_integral_only<eT>::result    atanh (const eT& x) { return arma_atanh(double(x)); }
00141   
00142   template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result acosh (const eT& x) { return arma_acosh(x); }
00143   template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result asinh (const eT& x) { return arma_asinh(x); }
00144   template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result atanh (const eT& x) { return arma_atanh(x); }
00145   
00146   template<typename eT> arma_inline static eT              conj(const eT&              x) { return x;            }
00147   template<typename  T> arma_inline static std::complex<T> conj(const std::complex<T>& x) { return std::conj(x); }
00148   
00149   template<typename eT> arma_inline static typename arma_integral_only<eT>::result sqrt  (const eT& x) { return std::sqrt (double(x)); }
00150   template<typename eT> arma_inline static typename arma_integral_only<eT>::result log10 (const eT& x) { return std::log10(double(x)); }
00151   template<typename eT> arma_inline static typename arma_integral_only<eT>::result log   (const eT& x) { return std::log  (double(x)); }
00152   template<typename eT> arma_inline static typename arma_integral_only<eT>::result exp   (const eT& x) { return std::exp  (double(x)); }
00153   template<typename eT> arma_inline static typename arma_integral_only<eT>::result cos   (const eT& x) { return std::cos  (double(x)); }
00154   template<typename eT> arma_inline static typename arma_integral_only<eT>::result sin   (const eT& x) { return std::sin  (double(x)); }
00155   template<typename eT> arma_inline static typename arma_integral_only<eT>::result tan   (const eT& x) { return std::tan  (double(x)); }
00156   template<typename eT> arma_inline static typename arma_integral_only<eT>::result cosh  (const eT& x) { return std::cosh (double(x)); }
00157   template<typename eT> arma_inline static typename arma_integral_only<eT>::result sinh  (const eT& x) { return std::sinh (double(x)); }
00158   template<typename eT> arma_inline static typename arma_integral_only<eT>::result tanh  (const eT& x) { return std::tanh (double(x)); }
00159   
00160   template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result sqrt  (const eT& x) { return std::sqrt (x); }
00161   template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result log10 (const eT& x) { return std::log10(x); }
00162   template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result log   (const eT& x) { return std::log  (x); }
00163   template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result exp   (const eT& x) { return std::exp  (x); }
00164   template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result cos   (const eT& x) { return std::cos  (x); }
00165   template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result sin   (const eT& x) { return std::sin  (x); }
00166   template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result tan   (const eT& x) { return std::tan  (x); }
00167   template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result cosh  (const eT& x) { return std::cosh (x); }
00168   template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result sinh  (const eT& x) { return std::sinh (x); }
00169   template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result tanh  (const eT& x) { return std::tanh (x); }
00170   
00171   
00172   
00173   template<typename T1, typename T2>
00174   arma_inline
00175   static
00176   typename arma_float_or_cx_only<T1>::result
00177   pow(const T1 base, const T2 exponent)
00178     {
00179     return std::pow(base, exponent);
00180     }
00181   
00182   
00183   
00184   template<typename T1, typename T2>
00185   arma_inline
00186   static
00187   typename arma_integral_only<T1>::result
00188   pow(const T1 base, const T2 exponent)
00189     {
00190     return T1( std::pow( double(base), double(exponent) ) );
00191     }
00192   
00193   
00194   
00195   template<typename T1>
00196   arma_inline
00197   static
00198   typename arma_float_or_cx_only<T1>::result
00199   pow_int(const T1 base, const int exponent)
00200     { 
00201     return std::pow(base, exponent);
00202     }
00203   
00204   
00205   
00206   template<typename T1>
00207   arma_inline
00208   static
00209   typename arma_integral_only<T1>::result
00210   pow_int(const T1 base, const int exponent)
00211     { 
00212     return T1( std::pow( double(base), exponent) );
00213     }
00214   
00215   
00216   
00217   template<typename eT>
00218   arma_inline
00219   static
00220   typename arma_integral_only<eT>::result
00221   direct_eps(const eT& x)
00222     {
00223     return eT(0);
00224     }
00225   
00226   
00227   
00228   template<typename eT>
00229   inline
00230   static
00231   typename arma_float_only<eT>::result
00232   direct_eps(const eT& x)
00233     {
00234     //arma_extra_debug_sigprint();
00235     
00236     // acording to IEEE Standard for Floating-Point Arithmetic (IEEE 754)
00237     // the mantissa length for double is 53 bits = std::numeric_limits<double>::digits
00238     // the mantissa length for float  is 24 bits = std::numeric_limits<float >::digits
00239     
00240     //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)) );
00241     
00242     const eT radix_eT     = eT(std::numeric_limits<eT>::radix);
00243     const eT digits_m1_eT = eT(std::numeric_limits<eT>::digits - 1);
00244     
00245     // return std::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) );
00246     return eop_aux::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) );
00247     }
00248   
00249   
00250   
00251   template<typename T>
00252   inline
00253   static
00254   typename arma_float_only<T>::result
00255   direct_eps(const std::complex<T>& x)
00256     {
00257     //arma_extra_debug_sigprint();
00258     
00259     //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)) );
00260     
00261     const T radix_T     = T(std::numeric_limits<T>::radix);
00262     const T digits_m1_T = T(std::numeric_limits<T>::digits - 1);
00263     
00264     return std::pow( radix_T, T(std::floor(std::log10(std::abs(x))/std::log10(radix_T)) - digits_m1_T) );
00265     }
00266   
00267   
00268   
00269   //! work around a bug in GCC 4.4
00270   template<typename eT> arma_inline static
00271   typename arma_unsigned_integral_only<eT>::result arma_abs(const eT& x)              { return x;           }
00272   
00273   template<typename eT> arma_inline static
00274   typename arma_signed_integral_only<eT>::result   arma_abs(const eT& x)              { return std::abs(x); }
00275   
00276   template<typename eT> arma_inline static
00277   typename arma_float_only<eT>::result             arma_abs(const eT& x)              { return std::abs(x); }
00278   
00279   template<typename T> arma_inline static
00280   typename arma_float_only<T>::result              arma_abs(const std::complex<T>& x) { return std::abs(x); }
00281   
00282   
00283   
00284   template<typename eT, typename eop_type>
00285   arma_inline
00286   static
00287   eT
00288   generate()
00289     {
00290          if(is_same_type<eop_type, eop_randu         >::value == true) { return eT(eop_aux_randu<eT>()); }
00291     else if(is_same_type<eop_type, eop_randn         >::value == true) { return eT(eop_aux_randn<eT>()); }
00292     else if(is_same_type<eop_type, eop_zeros         >::value == true) { return eT(0);                   }
00293     else if(is_same_type<eop_type, eop_ones_full     >::value == true) { return eT(1);                   }
00294     else if(is_same_type<eop_type, eop_cube_randu    >::value == true) { return eT(eop_aux_randu<eT>()); }
00295     else if(is_same_type<eop_type, eop_cube_randn    >::value == true) { return eT(eop_aux_randn<eT>()); }
00296     else if(is_same_type<eop_type, eop_cube_zeros    >::value == true) { return eT(0);                   }
00297     else if(is_same_type<eop_type, eop_cube_ones_full>::value == true) { return eT(1);                   }
00298     else                                                               { return eT(0);                   }
00299     }
00300   
00301   };
00302 
00303 
00304 
00305 //! @}
00306