headmyshoulder / odeint-v2

odeint - solving ordinary differential equations in c++ v2
http://headmyshoulder.github.com/odeint-v2/
Other
337 stars 102 forks source link

stepper algorithm makes an assumption that value_types are of type IEEE floating point, failing compile with alternative number systems #250

Open Ravenwater opened 4 years ago

Ravenwater commented 4 years ago

When trying to study the Lorenz system with very high precision arithmetic, the ODE stepper code fails to compile because it makes an assumption that the value_type is double:

time_type increase_step(time_type dt, value_type error, const int stepper_order) const
    {
        // returns the increased time step
        BOOST_USING_STD_MIN();
        BOOST_USING_STD_MAX();
        using std::pow;                     <--------- assumes your value_type is an IEEE floating point

        // adjust the size if dt is smaller than max_dt (providede max_dt is not zero)
        if(error < 0.5)
        {
            // error should be > 0
            error = max BOOST_PREVENT_MACRO_SUBSTITUTION (
                    static_cast<value_type>( pow( static_cast<value_type>(5.0) , -static_cast<value_type>(stepper_order) ) ) ,
                    error);
            // time_type dt_old = dt;   unused variable warning 
            //error too small - increase dt and keep the evolution and limit scaling factor to 5.0
            dt *= static_cast<value_type>(9)/static_cast<value_type>(10) *
                  pow(error, static_cast<value_type>(-1) / stepper_order);
            if(m_max_dt != static_cast<time_type >(0))
                // limit to maximal stepsize
                dt = detail::min_abs(dt, m_max_dt);
        }
        return dt;
    }

dt is double, and pow() is forced to use the standard library instead of the high precision arithmetic pow() function, so compilation fails:

1>lorenz2.cpp
1>Unknown compiler version - please run the configure tests and report the results
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(148): error C2679: binary '*=': no operator found which takes a right-hand operand of type 'sw::unum::posit<256,5>' (or there is no acceptable conversion)
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(149): note: could be 'built-in C++ operator*=(double, unsigned short)'
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(149): note: or       'built-in C++ operator*=(double, unsigned int)'
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(149): note: or       'built-in C++ operator*=(double, unsigned long)'
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(149): note: or       'built-in C++ operator*=(double, unsigned __int64)'
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(149): note: or       'built-in C++ operator*=(double, short)'
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(149): note: or       'built-in C++ operator*=(double, int)'
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(149): note: or       'built-in C++ operator*=(double, long)'
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(149): note: or       'built-in C++ operator*=(double, __int64)'
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(149): note: or       'built-in C++ operator*=(double, float)'
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(149): note: or       'built-in C++ operator*=(double, double)'
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(149): note: or       'built-in C++ operator*=(double, long double)'
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(149): note: while trying to match the argument list '(double, sw::unum::posit<256,5>)'
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(133): note: while compiling class template member function 'double boost::numeric::odeint::default_step_adjuster<sw::unum::posit<256,5>,double>::increase_step(double,sw::unum::posit<256,5>,const int) const'
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(778): note: see reference to function template instantiation 'double boost::numeric::odeint::default_step_adjuster<sw::unum::posit<256,5>,double>::increase_step(double,sw::unum::posit<256,5>,const int) const' being compiled
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/stepper/controlled_runge_kutta.hpp(905): note: see reference to class template instantiation 'boost::numeric::odeint::default_step_adjuster<sw::unum::posit<256,5>,double>' being compiled
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/numeric/odeint/integrate/integrate.hpp(48): note: see reference to class template instantiation 'boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_dopri5<State,sw::unum::posit<256,5>,State,Time,boost::numeric::odeint::algebra_dispatcher<State>::algebra_type,boost::numeric::odeint::operations_dispatcher_sfinae<StateType,void>::operations_type,boost::numeric::odeint::initially_resizer>,boost::numeric::odeint::default_error_checker<sw::unum::posit<256,5>,boost::numeric::odeint::array_algebra,boost::numeric::odeint::default_operations>,boost::numeric::odeint::default_step_adjuster<sw::unum::posit<256,5>,Time>,boost::numeric::odeint::initially_resizer,boost::numeric::odeint::explicit_error_stepper_fsal_base<boost::numeric::odeint::runge_kutta_dopri5<State,sw::unum::posit<256,5>,State,Time,boost::numeric::odeint::algebra_dispatcher<StateType>::algebra_type,boost::numeric::odeint::operations_dispatcher_sfinae<StateType,void>::operations_type,boost::numeric::odeint::initially_resizer>,5,5,4,State,Value,Deriv,Time,Algebra,Operations,Resizer>::stepper_category>' being compiled
1>        with
1>        [
1>            State=state_type,
1>            Time=double,
1>            StateType=state_type,
1>            Value=sw::unum::posit<256,5>,
1>            Deriv=state_type,
1>            Algebra=boost::numeric::odeint::algebra_dispatcher<state_type>::algebra_type,
1>            Operations=boost::numeric::odeint::operations_dispatcher_sfinae<state_type,void>::operations_type,
1>            Resizer=boost::numeric::odeint::initially_resizer
1>        ]
1>C:\Users\tomtz\Documents\dev\clones\odeint-v2\examples\posit\lorenz2.cpp(51): note: see reference to function template instantiation 'unsigned __int64 boost::numeric::odeint::integrate<void(__cdecl *)(const state_type &,state_type &,double),state_type,double,void(__cdecl *)(const state_type &,double)>(System,State &,Time,Time,Time,Observer)' being compiled
1>        with
1>        [
1>            System=void (__cdecl *)(const state_type &,state_type &,double),
1>            State=state_type,
1>            Time=double,
1>            Observer=void (__cdecl *)(const state_type &,double)
1>        ]
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/math/constants/constants.hpp(265): note: see reference to class template instantiation 'boost::math::constants::detail::constant_half<T>' being compiled
1>        with
1>        [
1>            T=double
1>        ]
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/math/constants/constants.hpp(265): note: see reference to function template instantiation 'double boost::math::constants::half<T,boost::math::policies::policy<boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy>>(void) noexcept' being compiled
1>        with
1>        [
1>            T=double
1>        ]
1>C:\Users\tomtz\Documents\dev\boost_1_65_0\boost/math/special_functions/detail/bernoulli_details.hpp(138): note: see reference to function template instantiation 'double boost::math::constants::half<double>(void) noexcept' being compiled

The program compiled is the Lorenz example of the wiki, with a value_type of posit<256,5>