op_eps_meat.hpp

Go to the documentation of this file.
00001 // Copyright (C) 2009 NICTA
00002 // Copyright (C) 2009 Dimitrios Bouzas
00003 // 
00004 // Authors:
00005 // - Conrad Sanderson (conradsand at ieee dot org)
00006 // - Dimitrios Bouzas (dimitris dot mpouzas at gmail dot com)
00007 // 
00008 // This file is part of the Armadillo C++ library.
00009 // It is provided without any warranty of fitness
00010 // for any purpose. You can redistribute this file
00011 // and/or modify it under the terms of the GNU
00012 // Lesser General Public License (LGPL) as published
00013 // by the Free Software Foundation, either version 3
00014 // of the License or (at your option) any later version.
00015 // (see http://www.opensource.org/licenses for more info)
00016 
00017 
00018 
00019 //! \addtogroup op_eps
00020 //! @{
00021 
00022 
00023 
00024 template<typename eT>
00025 inline
00026 eT
00027 op_eps::direct_eps(const eT x)
00028   {
00029   //arma_extra_debug_sigprint();
00030   
00031   // acording to IEEE Standard for Floating-Point Arithmetic (IEEE 754)
00032   // the mantissa length for double is 53 bits = std::numeric_limits<double>::digits
00033   // the mantissa length for float  is 24 bits = std::numeric_limits<float >::digits
00034   
00035   //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)) );
00036   
00037   const eT radix_eT     = eT(std::numeric_limits<eT>::radix);
00038   const eT digits_m1_eT = eT(std::numeric_limits<eT>::digits - 1);
00039   
00040   return std::pow( radix_eT, (std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) );
00041   }
00042 
00043 
00044 
00045 template<typename T>
00046 inline
00047 T
00048 op_eps::direct_eps(const std::complex<T>& x)
00049   {
00050   //arma_extra_debug_sigprint();
00051   
00052   //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)) );
00053   
00054   const T radix_T     = T(std::numeric_limits<T>::radix);
00055   const T digits_m1_T = T(std::numeric_limits<T>::digits - 1);
00056   
00057   return std::pow( radix_T, (std::floor(std::log10(std::abs(x))/std::log10(radix_T)) - digits_m1_T) );
00058   }
00059 
00060 
00061 
00062 template<typename eT>
00063 inline
00064 void
00065 op_eps::direct_eps(Mat<eT>& out, const Mat<eT>& A)
00066   {
00067   arma_extra_debug_sigprint();
00068   
00069   out.copy_size(A);
00070   
00071         eT* out_mem = out.memptr();
00072   const eT* A_mem   = A.memptr();
00073   const u32 n_elem  = A.n_elem;
00074   
00075   for(u32 i=0; i<n_elem; ++i)
00076     {
00077     out_mem[i] = op_eps::direct_eps( A_mem[i] );
00078     }
00079   
00080   }
00081 
00082 
00083 
00084 template<typename T>
00085 inline
00086 void
00087 op_eps::direct_eps(Mat<T>& out, const Mat< std::complex<T> >& A)
00088   {
00089   arma_extra_debug_sigprint();
00090   
00091   typedef typename std::complex<T> eT;
00092   
00093   out.copy_size(A);
00094   
00095          T* out_mem = out.memptr();
00096   const eT* A_mem   = A.memptr();
00097   const u32 n_elem  = A.n_elem;
00098   
00099   for(u32 i=0; i<n_elem; ++i)
00100     {
00101     out_mem[i] = op_eps::direct_eps( A_mem[i] );
00102     }
00103   
00104   }
00105 
00106 
00107 
00108 template<typename T1>
00109 inline
00110 void
00111 op_eps::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_eps>& in)
00112   {
00113   arma_extra_debug_sigprint();
00114   
00115   typedef typename T1::elem_type  eT;
00116   
00117   // in this case it doesn't matter if there is aliasing
00118   // (i.e. &out = &in.m), hence we use unwrap rather than unwrap_check
00119   
00120   const unwrap<T1>   tmp(in.m);
00121   const Mat<eT>& A = tmp.M;
00122   
00123   op_eps::direct_eps(out, A);
00124   }
00125 
00126 
00127 
00128 //! @}