From fcd771fda7cace3816e00148a83dd05abfee7b1d Mon Sep 17 00:00:00 2001 From: conrad Date: Wed, 20 Aug 2025 00:44:22 +1000 Subject: [PATCH] sort out the eps() mess --- include/armadillo | 1 - include/armadillo_bits/eop_aux.hpp | 49 ---------- include/armadillo_bits/eop_core_bones.hpp | 1 - include/armadillo_bits/eop_core_meat.hpp | 3 - include/armadillo_bits/fn_elem.hpp | 29 ++++++ include/armadillo_bits/fn_eps.hpp | 106 ---------------------- include/armadillo_bits/op_misc_bones.hpp | 14 +++ include/armadillo_bits/op_misc_meat.hpp | 85 +++++++++++++++++ 8 files changed, 128 insertions(+), 160 deletions(-) delete mode 100644 include/armadillo_bits/fn_eps.hpp diff --git a/include/armadillo b/include/armadillo index 6a9d608d..d7f2f010 100644 --- a/include/armadillo +++ b/include/armadillo @@ -549,7 +549,6 @@ namespace arma #include "armadillo_bits/fn_circshift.hpp" #include "armadillo_bits/fn_shuffle.hpp" #include "armadillo_bits/fn_prod.hpp" - #include "armadillo_bits/fn_eps.hpp" #include "armadillo_bits/fn_pinv.hpp" #include "armadillo_bits/fn_rank.hpp" #include "armadillo_bits/fn_kron.hpp" diff --git a/include/armadillo_bits/eop_aux.hpp b/include/armadillo_bits/eop_aux.hpp index e1c3c9f5..f79b8e49 100644 --- a/include/armadillo_bits/eop_aux.hpp +++ b/include/armadillo_bits/eop_aux.hpp @@ -136,55 +136,6 @@ struct eop_aux template arma_inline static typename arma_integral_only::result pow (const T1 base, const T2 exponent) { return T1( std::pow( double(base), double(exponent) ) ); } template arma_inline static typename arma_real_or_cx_only::result pow (const T1 base, const T2 exponent) { return T1( std::pow( base, exponent ) ); } - - - template - arma_inline - static - typename arma_integral_only::result - direct_eps(const eT) - { - return eT(0); - } - - - template - inline - static - typename arma_real_only::result - direct_eps(const eT x) - { - //arma_debug_sigprint(); - - // according to IEEE Standard for Floating-Point Arithmetic (IEEE 754) - // the mantissa length for double is 53 bits = std::numeric_limits::digits - // the mantissa length for float is 24 bits = std::numeric_limits::digits - - //return std::pow( std::numeric_limits::radix, (std::floor(std::log10(std::abs(x))/std::log10(std::numeric_limits::radix))-(std::numeric_limits::digits-1)) ); - - const eT radix_eT = eT(std::numeric_limits::radix); - const eT digits_m1_eT = eT(std::numeric_limits::digits - 1); - - // return std::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) ); - return eop_aux::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) ); - } - - - template - inline - static - typename arma_real_only::result - direct_eps(const std::complex& x) - { - //arma_debug_sigprint(); - - //return std::pow( std::numeric_limits::radix, (std::floor(std::log10(std::abs(x))/std::log10(std::numeric_limits::radix))-(std::numeric_limits::digits-1)) ); - - const T radix_T = T(std::numeric_limits::radix); - const T digits_m1_T = T(std::numeric_limits::digits - 1); - - return std::pow( radix_T, T(std::floor(std::log10(std::abs(x))/std::log10(radix_T)) - digits_m1_T) ); - } }; diff --git a/include/armadillo_bits/eop_core_bones.hpp b/include/armadillo_bits/eop_core_bones.hpp index de8fe75c..c87504f1 100644 --- a/include/armadillo_bits/eop_core_bones.hpp +++ b/include/armadillo_bits/eop_core_bones.hpp @@ -87,7 +87,6 @@ struct eop_acosh : public eop_core , public e struct eop_asinh : public eop_core , public eop_use_mp_true {}; struct eop_atanh : public eop_core , public eop_use_mp_true {}; struct eop_sinc : public eop_core , public eop_use_mp_true {}; -struct eop_eps : public eop_core , public eop_use_mp_true {}; struct eop_abs : public eop_core , public eop_use_mp_false {}; struct eop_arg : public eop_core , public eop_use_mp_false {}; struct eop_conj : public eop_core , public eop_use_mp_false {}; diff --git a/include/armadillo_bits/eop_core_meat.hpp b/include/armadillo_bits/eop_core_meat.hpp index 499bb575..1dab15e1 100644 --- a/include/armadillo_bits/eop_core_meat.hpp +++ b/include/armadillo_bits/eop_core_meat.hpp @@ -1104,9 +1104,6 @@ eop_core::process(const eT val, const eT ) { return eop_ template<> template arma_inline eT eop_core::process(const eT val, const eT ) { return arma_sinc(val); } -template<> template arma_inline eT -eop_core::process(const eT val, const eT ) { return eop_aux::direct_eps(val); } - template<> template arma_inline eT eop_core::process(const eT val, const eT ) { return eop_aux::arma_abs(val); } diff --git a/include/armadillo_bits/fn_elem.hpp b/include/armadillo_bits/fn_elem.hpp index 2befe7b9..dd1c3f5b 100644 --- a/include/armadillo_bits/fn_elem.hpp +++ b/include/armadillo_bits/fn_elem.hpp @@ -1190,6 +1190,35 @@ tgamma(const BaseCube& A) +// +// eps + +template +arma_warn_unused +inline +typename enable_if2< (is_arma_type::value && is_real_or_cx::value), const mtOp >::result +eps(const T1& X) + { + arma_debug_sigprint(); + + return mtOp(X); + } + + + +template +arma_warn_unused +inline +typename arma_real_or_cx_only::result +eps(const eT& x) + { + arma_debug_sigprint(); + + return op_eps::direct_eps(x); + } + + + // the functions below are currently unused; reserved for potential future use template void exp_approx(const T1&) { arma_stop_logic_error("unimplemented"); } diff --git a/include/armadillo_bits/fn_eps.hpp b/include/armadillo_bits/fn_eps.hpp deleted file mode 100644 index cb71fdfd..00000000 --- a/include/armadillo_bits/fn_eps.hpp +++ /dev/null @@ -1,106 +0,0 @@ -// SPDX-License-Identifier: Apache-2.0 -// -// Copyright 2008-2016 Conrad Sanderson (https://conradsanderson.id.au) -// Copyright 2008-2016 National ICT Australia (NICTA) -// -// Licensed under the Apache License, Version 2.0 (the "License"); -// you may not use this file except in compliance with the License. -// You may obtain a copy of the License at -// https://www.apache.org/licenses/LICENSE-2.0 -// -// Unless required by applicable law or agreed to in writing, software -// distributed under the License is distributed on an "AS IS" BASIS, -// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -// See the License for the specific language governing permissions and -// limitations under the License. -// ------------------------------------------------------------------------ - - - -//! \addtogroup fn_eps -//! @{ - - - -template -arma_warn_unused -inline -const eOp -eps(const Base& X, const typename arma_not_cx::result* junk = nullptr) - { - arma_debug_sigprint(); - arma_ignore(junk); - - return eOp(X.get_ref()); - } - - - -template -arma_warn_unused -inline -Mat< typename T1::pod_type > -eps(const Base< std::complex, T1>& X, const typename arma_cx_only::result* junk = nullptr) - { - arma_debug_sigprint(); - arma_ignore(junk); - - typedef typename T1::pod_type T; - typedef typename T1::elem_type eT; - - const unwrap tmp(X.get_ref()); - const Mat& A = tmp.M; - - Mat out(A.n_rows, A.n_cols, arma_nozeros_indicator()); - - T* out_mem = out.memptr(); - const eT* A_mem = A.memptr(); - - const uword n_elem = A.n_elem; - - for(uword i=0; i -arma_warn_unused -arma_inline -typename arma_integral_only::result -eps(const eT& x) - { - arma_ignore(x); - - return eT(0); - } - - - -template -arma_warn_unused -arma_inline -typename arma_real_only::result -eps(const eT& x) - { - return eop_aux::direct_eps(x); - } - - - -template -arma_warn_unused -arma_inline -typename arma_real_only::result -eps(const std::complex& x) - { - return eop_aux::direct_eps(x); - } - - - -//! @} diff --git a/include/armadillo_bits/op_misc_bones.hpp b/include/armadillo_bits/op_misc_bones.hpp index 1a83b676..a62b1f6b 100644 --- a/include/armadillo_bits/op_misc_bones.hpp +++ b/include/armadillo_bits/op_misc_bones.hpp @@ -79,4 +79,18 @@ struct op_replace +struct op_eps + : public traits_op_passthru + { + template inline static typename get_pod_type::result direct_eps(const eT& x); + + template + inline static void apply(Mat& out, const mtOp& in); + + template + inline static void apply_noalias(Mat& out, const Mat& X); + }; + + + //! @} diff --git a/include/armadillo_bits/op_misc_meat.hpp b/include/armadillo_bits/op_misc_meat.hpp index d06daa74..c4e0e719 100644 --- a/include/armadillo_bits/op_misc_meat.hpp +++ b/include/armadillo_bits/op_misc_meat.hpp @@ -109,6 +109,10 @@ op_real::apply( Cube& out, const mtOpCube inline void @@ -207,6 +211,10 @@ op_imag::apply( Cube& out, const mtOpCube inline void @@ -321,6 +329,10 @@ op_abs::apply( Cube& out, const mtOpCube inline void @@ -411,6 +423,10 @@ op_arg::apply( Cube& out, const mtOpCube inline void @@ -445,4 +461,73 @@ op_replace::apply(Cube& out, const mtOpCube& in) +// + + + +template +inline +typename get_pod_type::result +op_eps::direct_eps(const eT& x) + { + typedef typename get_pod_type::result T; + + const T xx = std::abs(x); + + const T yy = std::nextafter(xx, std::numeric_limits::infinity()); + + return (yy - xx); + } + + + +template +inline +void +op_eps::apply(Mat& out, const mtOp& in) + { + arma_debug_sigprint(); + + typedef typename T1::pod_type T; + + const quasi_unwrap U(in.m); + + if(U.is_alias(out)) + { + Mat tmp; + + op_eps::apply_noalias(tmp, U.M); + + out.steal_mem(tmp); + } + else + { + op_eps::apply_noalias(out, U.M); + } + } + + + +template +inline +void +op_eps::apply_noalias(Mat& out, const Mat& X) + { + arma_debug_sigprint(); + + out.set_size(X.n_rows, X.n_cols); + + T* out_mem = out.memptr(); + const eT* X_mem = X.memptr(); + + const uword n_elem = X.n_elem; + + for(uword i=0; i