$search
00001 // Copyright (C) 2010-2011 NICTA (www.nicta.com.au) 00002 // Copyright (C) 2010-2011 Conrad Sanderson 00003 // 00004 // This file is part of the Armadillo C++ library. 00005 // It is provided without any warranty of fitness 00006 // for any purpose. You can redistribute this file 00007 // and/or modify it under the terms of the GNU 00008 // Lesser General Public License (LGPL) as published 00009 // by the Free Software Foundation, either version 3 00010 // of the License or (at your option) any later version. 00011 // (see http://www.opensource.org/licenses for more info) 00012 00013 00016 00017 00018 00019 template<typename eT> 00020 struct eop_aux_randu 00021 { 00022 arma_inline 00023 operator eT () 00024 { 00025 return eT(std::rand()) / eT(RAND_MAX); 00026 } 00027 }; 00028 00029 00030 00031 template<typename T> 00032 struct eop_aux_randu< std::complex<T> > 00033 { 00034 arma_inline 00035 operator std::complex<T> () 00036 { 00037 return std::complex<T>( T(eop_aux_randu<T>()), T(eop_aux_randu<T>()) ); 00038 } 00039 }; 00040 00041 00042 00043 template<typename eT> 00044 struct eop_aux_randn 00045 { 00046 // // rudimentary method, based on the central limit theorem 00047 // // http://en.wikipedia.org/wiki/Central_limit_theorem 00048 // inline 00049 // operator eT () const 00050 // { 00051 // const uword N = 12; // N must be >= 12 and an even number 00052 // const uword N2 = N/2; 00053 // 00054 // eT acc = eT(0); 00055 // 00056 // for(uword i=0; i<N2; ++i) 00057 // { 00058 // const eT tmp1 = eT(std::rand()) / eT(RAND_MAX); 00059 // const eT tmp2 = eT(std::rand()) / eT(RAND_MAX); 00060 // acc += tmp1+tmp2; 00061 // } 00062 // 00063 // return acc - eT(N2); 00064 // } 00065 00066 00067 // polar form of the Box-Muller transformation 00068 // http://en.wikipedia.org/wiki/Box-Muller_transformation 00069 // http://en.wikipedia.org/wiki/Marsaglia_polar_method 00070 inline 00071 operator eT () const 00072 { 00073 // make sure we are internally using at least floats 00074 typedef typename promote_type<eT,float>::result eTp; 00075 00076 eTp tmp1; 00077 eTp tmp2; 00078 eTp w; 00079 00080 do 00081 { 00082 tmp1 = eTp(2) * eTp(std::rand()) / eTp(RAND_MAX) - eTp(1); 00083 tmp2 = eTp(2) * eTp(std::rand()) / eTp(RAND_MAX) - eTp(1); 00084 w = tmp1*tmp1 + tmp2*tmp2; 00085 } 00086 while ( w >= eTp(1) ); 00087 00088 return eT( tmp1 * std::sqrt( (eTp(-2) * std::log(w)) / w) ); 00089 } 00090 00091 00092 // other methods: 00093 // http://en.wikipedia.org/wiki/Ziggurat_algorithm 00094 // 00095 // Marsaglia and Tsang Ziggurat technique to transform from a uniform to a normal distribution. 00096 // G. Marsaglia, W.W. Tsang. 00097 // "Ziggurat method for generating random variables", 00098 // J. Statistical Software, vol 5, 2000. 00099 // http://www.jstatsoft.org/v05/i08/ 00100 }; 00101 00102 00103 00104 template<typename T> 00105 struct eop_aux_randn< std::complex<T> > 00106 { 00107 arma_inline 00108 operator std::complex<T> () const 00109 { 00110 return std::complex<T>( T(eop_aux_randn<T>()), T(eop_aux_randn<T>()) ); 00111 } 00112 00113 }; 00114 00115 00118 00119 class eop_aux 00120 { 00121 public: 00122 00123 template<typename eT> arma_inline static typename arma_integral_only<eT>::result acos (const eT x) { return eT( std::acos(double(x)) ); } 00124 template<typename eT> arma_inline static typename arma_integral_only<eT>::result asin (const eT x) { return eT( std::asin(double(x)) ); } 00125 template<typename eT> arma_inline static typename arma_integral_only<eT>::result atan (const eT x) { return eT( std::atan(double(x)) ); } 00126 00127 template<typename eT> arma_inline static typename arma_float_only<eT>::result acos (const eT x) { return std::acos(x); } 00128 template<typename eT> arma_inline static typename arma_float_only<eT>::result asin (const eT x) { return std::asin(x); } 00129 template<typename eT> arma_inline static typename arma_float_only<eT>::result atan (const eT x) { return std::atan(x); } 00130 00131 template<typename eT> arma_inline static typename arma_cx_only<eT>::result acos (const eT x) { return arma_acos(x); } 00132 template<typename eT> arma_inline static typename arma_cx_only<eT>::result asin (const eT x) { return arma_asin(x); } 00133 template<typename eT> arma_inline static typename arma_cx_only<eT>::result atan (const eT x) { return arma_atan(x); } 00134 00135 template<typename eT> arma_inline static typename arma_integral_only<eT>::result acosh (const eT x) { return eT( arma_acosh(double(x)) ); } 00136 template<typename eT> arma_inline static typename arma_integral_only<eT>::result asinh (const eT x) { return eT( arma_asinh(double(x)) ); } 00137 template<typename eT> arma_inline static typename arma_integral_only<eT>::result atanh (const eT x) { return eT( arma_atanh(double(x)) ); } 00138 00139 template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result acosh (const eT x) { return arma_acosh(x); } 00140 template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result asinh (const eT x) { return arma_asinh(x); } 00141 template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result atanh (const eT x) { return arma_atanh(x); } 00142 00143 template<typename eT> arma_inline static typename arma_not_cx<eT>::result conj(const eT x) { return x; } 00144 template<typename T> arma_inline static std::complex<T> conj(const std::complex<T>& x) { return std::conj(x); } 00145 00146 template<typename eT> arma_inline static typename arma_integral_only<eT>::result sqrt (const eT x) { return eT( std::sqrt (double(x)) ); } 00147 template<typename eT> arma_inline static typename arma_integral_only<eT>::result log10 (const eT x) { return eT( std::log10(double(x)) ); } 00148 template<typename eT> arma_inline static typename arma_integral_only<eT>::result log (const eT x) { return eT( std::log (double(x)) ); } 00149 template<typename eT> arma_inline static typename arma_integral_only<eT>::result exp (const eT x) { return eT( std::exp (double(x)) ); } 00150 template<typename eT> arma_inline static typename arma_integral_only<eT>::result cos (const eT x) { return eT( std::cos (double(x)) ); } 00151 template<typename eT> arma_inline static typename arma_integral_only<eT>::result sin (const eT x) { return eT( std::sin (double(x)) ); } 00152 template<typename eT> arma_inline static typename arma_integral_only<eT>::result tan (const eT x) { return eT( std::tan (double(x)) ); } 00153 template<typename eT> arma_inline static typename arma_integral_only<eT>::result cosh (const eT x) { return eT( std::cosh (double(x)) ); } 00154 template<typename eT> arma_inline static typename arma_integral_only<eT>::result sinh (const eT x) { return eT( std::sinh (double(x)) ); } 00155 template<typename eT> arma_inline static typename arma_integral_only<eT>::result tanh (const eT x) { return eT( std::tanh (double(x)) ); } 00156 00157 template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result sqrt (const eT x) { return std::sqrt (x); } 00158 template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result log10 (const eT x) { return std::log10(x); } 00159 template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result log (const eT x) { return std::log (x); } 00160 template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result exp (const eT x) { return std::exp (x); } 00161 template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result cos (const eT x) { return std::cos (x); } 00162 template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result sin (const eT x) { return std::sin (x); } 00163 template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result tan (const eT x) { return std::tan (x); } 00164 template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result cosh (const eT x) { return std::cosh (x); } 00165 template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result sinh (const eT x) { return std::sinh (x); } 00166 template<typename eT> arma_inline static typename arma_float_or_cx_only<eT>::result tanh (const eT x) { return std::tanh (x); } 00167 00168 template<typename eT> arma_inline static typename arma_unsigned_integral_only<eT>::result neg (const eT x) { return x; } 00169 template<typename eT> arma_inline static typename arma_signed_only<eT>::result neg (const eT x) { return -x; } 00170 00171 template<typename eT> arma_inline static typename arma_integral_only<eT>::result floor(const eT x) { return x; } 00172 template<typename eT> arma_inline static typename arma_float_only<eT>::result floor(const eT x) { return std::floor(x); } 00173 template<typename eT> arma_inline static typename arma_cx_only<eT>::result floor(const eT& x) { return eT( std::floor(x.real()), std::floor(x.imag()) ); } 00174 00175 template<typename eT> arma_inline static typename arma_integral_only<eT>::result ceil(const eT x) { return x; } 00176 template<typename eT> arma_inline static typename arma_float_only<eT>::result ceil(const eT x) { return std::ceil(x); } 00177 template<typename eT> arma_inline static typename arma_cx_only<eT>::result ceil(const eT& x) { return eT( std::ceil(x.real()), std::ceil(x.imag()) ); } 00178 00179 template<typename eT> 00180 arma_inline 00181 static 00182 typename arma_integral_only<eT>::result 00183 log2 (const eT x) 00184 { 00185 return eT( std::log(double(x))/ double(0.69314718055994530942) ); 00186 } 00187 00188 00189 template<typename eT> 00190 arma_inline 00191 static 00192 typename arma_float_or_cx_only<eT>::result 00193 log2 (const eT x) 00194 { 00195 typedef typename get_pod_type<eT>::result T; 00196 return std::log(x) / T(0.69314718055994530942); 00197 } 00198 00199 00200 template<typename eT> 00201 arma_inline 00202 static 00203 typename arma_integral_only<eT>::result 00204 exp10 (const eT x) 00205 { 00206 return eT( std::pow(double(10), double(x)) ); 00207 } 00208 00209 00210 template<typename eT> 00211 arma_inline 00212 static 00213 typename 00214 arma_float_or_cx_only<eT>::result 00215 exp10 (const eT x) 00216 { 00217 typedef typename get_pod_type<eT>::result T; 00218 return std::pow( T(10), x); 00219 } 00220 00221 00222 template<typename eT> 00223 arma_inline 00224 static 00225 typename arma_integral_only<eT>::result 00226 exp2 (const eT x) 00227 { 00228 return eT( std::pow(double(2), double(x)) ); 00229 } 00230 00231 00232 template<typename eT> 00233 arma_inline 00234 static 00235 typename arma_float_or_cx_only<eT>::result 00236 exp2 (const eT x) 00237 { 00238 typedef typename get_pod_type<eT>::result T; 00239 return std::pow( T(2), x); 00240 } 00241 00242 00243 template<typename T1, typename T2> 00244 arma_inline 00245 static 00246 typename arma_float_or_cx_only<T1>::result 00247 pow(const T1 base, const T2 exponent) 00248 { 00249 return std::pow(base, exponent); 00250 } 00251 00252 00253 00254 template<typename T1, typename T2> 00255 arma_inline 00256 static 00257 typename arma_integral_only<T1>::result 00258 pow(const T1 base, const T2 exponent) 00259 { 00260 return T1( std::pow( double(base), double(exponent) ) ); 00261 } 00262 00263 00264 00265 template<typename eT> 00266 arma_inline 00267 static 00268 typename arma_integral_only<eT>::result 00269 direct_eps(const eT) 00270 { 00271 return eT(0); 00272 } 00273 00274 00275 00276 template<typename eT> 00277 inline 00278 static 00279 typename arma_float_only<eT>::result 00280 direct_eps(const eT x) 00281 { 00282 //arma_extra_debug_sigprint(); 00283 00284 // acording to IEEE Standard for Floating-Point Arithmetic (IEEE 754) 00285 // the mantissa length for double is 53 bits = std::numeric_limits<double>::digits 00286 // the mantissa length for float is 24 bits = std::numeric_limits<float >::digits 00287 00288 //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)) ); 00289 00290 const eT radix_eT = eT(std::numeric_limits<eT>::radix); 00291 const eT digits_m1_eT = eT(std::numeric_limits<eT>::digits - 1); 00292 00293 // return std::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) ); 00294 return eop_aux::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) ); 00295 } 00296 00297 00298 00299 template<typename T> 00300 inline 00301 static 00302 typename arma_float_only<T>::result 00303 direct_eps(const std::complex<T> x) 00304 { 00305 //arma_extra_debug_sigprint(); 00306 00307 //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)) ); 00308 00309 const T radix_T = T(std::numeric_limits<T>::radix); 00310 const T digits_m1_T = T(std::numeric_limits<T>::digits - 1); 00311 00312 return std::pow( radix_T, T(std::floor(std::log10(std::abs(x))/std::log10(radix_T)) - digits_m1_T) ); 00313 } 00314 00315 00316 00318 template<typename eT> arma_inline static 00319 typename arma_unsigned_integral_only<eT>::result arma_abs(const eT x) { return x; } 00320 00321 template<typename eT> arma_inline static 00322 typename arma_signed_integral_only<eT>::result arma_abs(const eT x) { return std::abs(x); } 00323 00324 template<typename eT> arma_inline static 00325 typename arma_float_only<eT>::result arma_abs(const eT x) { return std::abs(x); } 00326 00327 template<typename T> arma_inline static 00328 typename arma_float_only<T>::result arma_abs(const std::complex<T> x) { return std::abs(x); } 00329 00330 }; 00331 00332 00333 00335