Chris@49: // Copyright (C) 2010-2013 NICTA (www.nicta.com.au)
Chris@49: // Copyright (C) 2010-2013 Conrad Sanderson
Chris@49: // 
Chris@49: // This Source Code Form is subject to the terms of the Mozilla Public
Chris@49: // License, v. 2.0. If a copy of the MPL was not distributed with this
Chris@49: // file, You can obtain one at http://mozilla.org/MPL/2.0/.
Chris@49: 
Chris@49: 
Chris@49: //! \addtogroup eop_aux
Chris@49: //! @{
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: struct eop_aux_randu
Chris@49:   {
Chris@49:   arma_inline
Chris@49:   operator eT ()
Chris@49:     {
Chris@49:     // make sure we are internally using at least floats
Chris@49:     typedef typename promote_type<eT,float>::result eTp;
Chris@49:     
Chris@49:     return eT( eTp(std::rand()) * ( eTp(1) / eTp(RAND_MAX) ) );
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   inline
Chris@49:   static
Chris@49:   void
Chris@49:   fill(eT* mem, const uword N)
Chris@49:     {
Chris@49:     uword i,j;
Chris@49:     
Chris@49:     for(i=0, j=1; j < N; i+=2, j+=2)
Chris@49:       {
Chris@49:       const eT tmp_i = eT(eop_aux_randu<eT>());
Chris@49:       const eT tmp_j = eT(eop_aux_randu<eT>());
Chris@49:       
Chris@49:       mem[i] = tmp_i;
Chris@49:       mem[j] = tmp_j;
Chris@49:       }
Chris@49:     
Chris@49:     if(i < N)
Chris@49:       {
Chris@49:       mem[i] = eT(eop_aux_randu<eT>());
Chris@49:       }
Chris@49:     }
Chris@49:   };
Chris@49: 
Chris@49: 
Chris@49:   
Chris@49: template<typename T>
Chris@49: struct eop_aux_randu< std::complex<T> >
Chris@49:   {
Chris@49:   arma_inline
Chris@49:   operator std::complex<T> ()
Chris@49:     {
Chris@49:     return std::complex<T>( T(eop_aux_randu<T>()), T(eop_aux_randu<T>()) );
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   inline
Chris@49:   static
Chris@49:   void
Chris@49:   fill(std::complex<T>* mem, const uword N)
Chris@49:     {
Chris@49:     for(uword i=0; i < N; ++i)
Chris@49:       {
Chris@49:       mem[i] = std::complex<T>( eop_aux_randu< std::complex<T> >() );
Chris@49:       }
Chris@49:     }
Chris@49:   };
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: struct eop_aux_randn
Chris@49:   {
Chris@49:   // rudimentary method, based on the central limit theorem:
Chris@49:   // http://en.wikipedia.org/wiki/Central_limit_theorem
Chris@49:   
Chris@49:   // polar form of the Box-Muller transformation:
Chris@49:   // http://en.wikipedia.org/wiki/Box-Muller_transformation
Chris@49:   // http://en.wikipedia.org/wiki/Marsaglia_polar_method
Chris@49:   
Chris@49:   // other methods:
Chris@49:   // http://en.wikipedia.org/wiki/Ziggurat_algorithm
Chris@49:   //
Chris@49:   // Marsaglia and Tsang Ziggurat technique to transform from a uniform to a normal distribution.
Chris@49:   // G. Marsaglia, W.W. Tsang.
Chris@49:   // "Ziggurat method for generating random variables",
Chris@49:   // J. Statistical Software, vol 5, 2000.
Chris@49:   // http://www.jstatsoft.org/v05/i08/
Chris@49:   
Chris@49:   
Chris@49:   // currently using polar form of the Box-Muller transformation
Chris@49:   inline
Chris@49:   operator eT () const
Chris@49:     {
Chris@49:     // make sure we are internally using at least floats
Chris@49:     typedef typename promote_type<eT,float>::result eTp;
Chris@49:     
Chris@49:     eTp tmp1;
Chris@49:     eTp tmp2;
Chris@49:     eTp w;
Chris@49:     
Chris@49:     do
Chris@49:       {
Chris@49:       tmp1 = eTp(2) * eTp(std::rand()) * (eTp(1) / eTp(RAND_MAX)) - eTp(1);
Chris@49:       tmp2 = eTp(2) * eTp(std::rand()) * (eTp(1) / eTp(RAND_MAX)) - eTp(1);
Chris@49:       
Chris@49:       w = tmp1*tmp1 + tmp2*tmp2;
Chris@49:       }
Chris@49:     while ( w >= eTp(1) );
Chris@49:     
Chris@49:     return eT( tmp1 * std::sqrt( (eTp(-2) * std::log(w)) / w) );
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   
Chris@49:   inline
Chris@49:   static
Chris@49:   void
Chris@49:   generate(eT& out1, eT& out2)
Chris@49:     {
Chris@49:     // make sure we are internally using at least floats
Chris@49:     typedef typename promote_type<eT,float>::result eTp;
Chris@49:     
Chris@49:     eTp tmp1;
Chris@49:     eTp tmp2;
Chris@49:     eTp w;
Chris@49:     
Chris@49:     do
Chris@49:       {
Chris@49:       tmp1 = eTp(2) * eTp(std::rand()) * (eTp(1) / eTp(RAND_MAX)) - eTp(1);
Chris@49:       tmp2 = eTp(2) * eTp(std::rand()) * (eTp(1) / eTp(RAND_MAX)) - eTp(1);
Chris@49:       
Chris@49:       w = tmp1*tmp1 + tmp2*tmp2;
Chris@49:       }
Chris@49:     while ( w >= eTp(1) );
Chris@49:     
Chris@49:     const eTp k = std::sqrt( (eTp(-2) * std::log(w)) / w);
Chris@49:     
Chris@49:     out1 = eT(tmp1*k);
Chris@49:     out2 = eT(tmp2*k);
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   
Chris@49:   inline
Chris@49:   static
Chris@49:   void
Chris@49:   fill(eT* mem, const uword N)
Chris@49:     {
Chris@49:     uword i, j;
Chris@49:     
Chris@49:     for(i=0, j=1; j < N; i+=2, j+=2)
Chris@49:       {
Chris@49:       eop_aux_randn<eT>::generate( mem[i], mem[j] );
Chris@49:       }
Chris@49:     
Chris@49:     if(i < N)
Chris@49:       {
Chris@49:       mem[i] = eT(eop_aux_randn<eT>());
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   };
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename T>
Chris@49: struct eop_aux_randn< std::complex<T> >
Chris@49:   {
Chris@49:   inline
Chris@49:   operator std::complex<T> () const
Chris@49:     {
Chris@49:     T a, b;
Chris@49:     
Chris@49:     eop_aux_randn<T>::generate(a, b);
Chris@49:     
Chris@49:     return std::complex<T>(a, b);
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   inline
Chris@49:   static
Chris@49:   void
Chris@49:   fill(std::complex<T>* mem, const uword N)
Chris@49:     {
Chris@49:     for(uword i=0; i < N; ++i)
Chris@49:       {
Chris@49:       mem[i] = std::complex<T>( eop_aux_randn< std::complex<T> >() );
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   };
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! use of the SFINAE approach to work around compiler limitations
Chris@49: //! http://en.wikipedia.org/wiki/SFINAE
Chris@49: 
Chris@49: class eop_aux
Chris@49:   {
Chris@49:   public:
Chris@49:   
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result    acos  (const eT x) { return eT( std::acos(double(x)) ); }
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result    asin  (const eT x) { return eT( std::asin(double(x)) ); }
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result    atan  (const eT x) { return eT( std::atan(double(x)) ); }
Chris@49:   
Chris@49:   template<typename eT> arma_inline static typename arma_real_only<eT>::result        acos  (const eT x) { return std::acos(x); }
Chris@49:   template<typename eT> arma_inline static typename arma_real_only<eT>::result        asin  (const eT x) { return std::asin(x); }
Chris@49:   template<typename eT> arma_inline static typename arma_real_only<eT>::result        atan  (const eT x) { return std::atan(x); }
Chris@49:   
Chris@49:   template<typename eT> arma_inline static typename arma_cx_only<eT>::result          acos  (const eT x) { return arma_acos(x); }
Chris@49:   template<typename eT> arma_inline static typename arma_cx_only<eT>::result          asin  (const eT x) { return arma_asin(x); }
Chris@49:   template<typename eT> arma_inline static typename arma_cx_only<eT>::result          atan  (const eT x) { return arma_atan(x); }
Chris@49:   
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result    acosh (const eT x) { return eT( arma_acosh(double(x)) ); }
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result    asinh (const eT x) { return eT( arma_asinh(double(x)) ); }
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result    atanh (const eT x) { return eT( arma_atanh(double(x)) ); }
Chris@49:   
Chris@49:   template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result acosh (const eT x) { return arma_acosh(x); }
Chris@49:   template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result asinh (const eT x) { return arma_asinh(x); }
Chris@49:   template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result atanh (const eT x) { return arma_atanh(x); }
Chris@49:   
Chris@49:   template<typename eT> arma_inline static typename arma_not_cx<eT>::result conj(const eT               x) { return x;            }
Chris@49:   template<typename  T> arma_inline static          std::complex<T>         conj(const std::complex<T>& x) { return std::conj(x); }
Chris@49:   
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result sqrt  (const eT x) { return eT( std::sqrt (double(x)) ); }
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result log10 (const eT x) { return eT( std::log10(double(x)) ); }
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result log   (const eT x) { return eT( std::log  (double(x)) ); }
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result exp   (const eT x) { return eT( std::exp  (double(x)) ); }
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result cos   (const eT x) { return eT( std::cos  (double(x)) ); }
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result sin   (const eT x) { return eT( std::sin  (double(x)) ); }
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result tan   (const eT x) { return eT( std::tan  (double(x)) ); }
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result cosh  (const eT x) { return eT( std::cosh (double(x)) ); }
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result sinh  (const eT x) { return eT( std::sinh (double(x)) ); }
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result tanh  (const eT x) { return eT( std::tanh (double(x)) ); }
Chris@49:   
Chris@49:   template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result sqrt  (const eT x) { return std::sqrt (x); }
Chris@49:   template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result log10 (const eT x) { return std::log10(x); }
Chris@49:   template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result log   (const eT x) { return std::log  (x); }
Chris@49:   template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result exp   (const eT x) { return std::exp  (x); }
Chris@49:   template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result cos   (const eT x) { return std::cos  (x); }
Chris@49:   template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result sin   (const eT x) { return std::sin  (x); }
Chris@49:   template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result tan   (const eT x) { return std::tan  (x); }
Chris@49:   template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result cosh  (const eT x) { return std::cosh (x); }
Chris@49:   template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result sinh  (const eT x) { return std::sinh (x); }
Chris@49:   template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result tanh  (const eT x) { return std::tanh (x); }
Chris@49:   
Chris@49:   template<typename eT> arma_inline static typename arma_unsigned_integral_only<eT>::result neg (const eT x) { return  x; }
Chris@49:   template<typename eT> arma_inline static typename arma_signed_only<eT>::result            neg (const eT x) { return -x; }
Chris@49:   
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result floor(const eT  x) { return x;                                                }
Chris@49:   template<typename eT> arma_inline static typename arma_real_only<eT>::result     floor(const eT  x) { return std::floor(x);                                    }
Chris@49:   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()) ); }
Chris@49:   
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result  ceil(const eT  x) { return x;                                                }
Chris@49:   template<typename eT> arma_inline static typename arma_real_only<eT>::result      ceil(const eT  x) { return std::ceil(x);                                     }
Chris@49:   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()) );   }
Chris@49:   
Chris@49:   template<typename eT> arma_inline static typename arma_integral_only<eT>::result round(const eT  x) { return x;                                                        }
Chris@49:   template<typename eT> arma_inline static typename arma_real_only<eT>::result     round(const eT  x) { return (x >= eT(0)) ? std::floor(x+0.5) : std::ceil(x-0.5);      }
Chris@49:   template<typename eT> arma_inline static typename arma_cx_only<eT>::result       round(const eT& x) { return eT( eop_aux::round(x.real()), eop_aux::round(x.imag()) ); }
Chris@49:   
Chris@49:   template<typename eT>
Chris@49:   arma_inline
Chris@49:   static
Chris@49:   typename arma_integral_only<eT>::result
Chris@49:   log2 (const eT x)
Chris@49:     {
Chris@49:     return eT( std::log(double(x))/ double(0.69314718055994530942) );
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   template<typename eT>
Chris@49:   arma_inline
Chris@49:   static
Chris@49:   typename arma_real_or_cx_only<eT>::result
Chris@49:   log2 (const eT x)
Chris@49:     {
Chris@49:     typedef typename get_pod_type<eT>::result T;
Chris@49:     return std::log(x) / T(0.69314718055994530942);
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   template<typename eT>
Chris@49:   arma_inline
Chris@49:   static
Chris@49:   typename arma_integral_only<eT>::result
Chris@49:   exp10 (const eT x)
Chris@49:     {
Chris@49:     return eT( std::pow(double(10), double(x)) );
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   template<typename eT>
Chris@49:   arma_inline
Chris@49:   static
Chris@49:   typename
Chris@49:   arma_real_or_cx_only<eT>::result
Chris@49:   exp10 (const eT x)
Chris@49:     {
Chris@49:     typedef typename get_pod_type<eT>::result T;
Chris@49:     return std::pow( T(10), x);
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   template<typename eT>
Chris@49:   arma_inline
Chris@49:   static
Chris@49:   typename arma_integral_only<eT>::result
Chris@49:   exp2 (const eT x)
Chris@49:     {
Chris@49:     return eT( std::pow(double(2), double(x)) );
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   template<typename eT>
Chris@49:   arma_inline
Chris@49:   static
Chris@49:   typename arma_real_or_cx_only<eT>::result
Chris@49:   exp2 (const eT x)
Chris@49:     {
Chris@49:     typedef typename get_pod_type<eT>::result T;
Chris@49:     return std::pow( T(2), x);
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   template<typename T1, typename T2>
Chris@49:   arma_inline
Chris@49:   static
Chris@49:   typename arma_real_or_cx_only<T1>::result
Chris@49:   pow(const T1 base, const T2 exponent)
Chris@49:     {
Chris@49:     return std::pow(base, exponent);
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   
Chris@49:   template<typename T1, typename T2>
Chris@49:   arma_inline
Chris@49:   static
Chris@49:   typename arma_integral_only<T1>::result
Chris@49:   pow(const T1 base, const T2 exponent)
Chris@49:     {
Chris@49:     return T1( std::pow( double(base), double(exponent) ) );
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   
Chris@49:   template<typename eT>
Chris@49:   arma_inline
Chris@49:   static
Chris@49:   typename arma_integral_only<eT>::result
Chris@49:   direct_eps(const eT)
Chris@49:     {
Chris@49:     return eT(0);
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   
Chris@49:   template<typename eT>
Chris@49:   inline
Chris@49:   static
Chris@49:   typename arma_real_only<eT>::result
Chris@49:   direct_eps(const eT x)
Chris@49:     {
Chris@49:     //arma_extra_debug_sigprint();
Chris@49:     
Chris@49:     // acording to IEEE Standard for Floating-Point Arithmetic (IEEE 754)
Chris@49:     // the mantissa length for double is 53 bits = std::numeric_limits<double>::digits
Chris@49:     // the mantissa length for float  is 24 bits = std::numeric_limits<float >::digits
Chris@49:     
Chris@49:     //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)) );
Chris@49:     
Chris@49:     const eT radix_eT     = eT(std::numeric_limits<eT>::radix);
Chris@49:     const eT digits_m1_eT = eT(std::numeric_limits<eT>::digits - 1);
Chris@49:     
Chris@49:     // return std::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) );
Chris@49:     return eop_aux::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) );
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   
Chris@49:   template<typename T>
Chris@49:   inline
Chris@49:   static
Chris@49:   typename arma_real_only<T>::result
Chris@49:   direct_eps(const std::complex<T> x)
Chris@49:     {
Chris@49:     //arma_extra_debug_sigprint();
Chris@49:     
Chris@49:     //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)) );
Chris@49:     
Chris@49:     const T radix_T     = T(std::numeric_limits<T>::radix);
Chris@49:     const T digits_m1_T = T(std::numeric_limits<T>::digits - 1);
Chris@49:     
Chris@49:     return std::pow( radix_T, T(std::floor(std::log10(std::abs(x))/std::log10(radix_T)) - digits_m1_T) );
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   
Chris@49:   //! work around a bug in GCC 4.4
Chris@49:   template<typename eT> arma_inline static
Chris@49:   typename arma_unsigned_integral_only<eT>::result arma_abs(const eT x)              { return x;           }
Chris@49:   
Chris@49:   template<typename eT> arma_inline static
Chris@49:   typename arma_signed_integral_only<eT>::result   arma_abs(const eT x)              { return std::abs(x); }
Chris@49:   
Chris@49:   template<typename eT> arma_inline static
Chris@49:   typename arma_real_only<eT>::result              arma_abs(const eT x)              { return std::abs(x); }
Chris@49:   
Chris@49:   template<typename T> arma_inline static
Chris@49:   typename arma_real_only<T>::result               arma_abs(const std::complex<T> x) { return std::abs(x); }
Chris@49:   
Chris@49:   };
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! @}
Chris@49: