netlist: Completed __float128 support. [Couriersud]

Both compiling the core and the shaders with __float128 now work.
The support was added to be ready to deal with academic edge cases.

Performance drops to 10% of double - thus disabled by default.
This commit is contained in:
couriersud 2019-11-03 15:08:41 +01:00
parent 36c17abc79
commit 34ccb11c53
20 changed files with 283 additions and 59 deletions

View File

@ -193,7 +193,7 @@ public:
// FIXME: make this a parameter // FIXME: make this a parameter
// avoid calls due to noise // avoid calls due to noise
if (std::fabs(cur - m_last) > 1e-6) if (plib::abs(cur - m_last) > 1e-6)
{ {
m_cpu_device->update_icount(exec().time()); m_cpu_device->update_icount(exec().time());
(*m_callback)(cur, m_cpu_device->local_time()); (*m_callback)(cur, m_cpu_device->local_time());
@ -446,7 +446,7 @@ public:
nl_fptype val = m_in() * m_mult() + m_offset(); nl_fptype val = m_in() * m_mult() + m_offset();
sound_update(exec().time()); sound_update(exec().time());
/* ignore spikes */ /* ignore spikes */
if (std::abs(val) < 32767.0) if (plib::abs(val) < 32767.0)
m_cur = val; m_cur = val;
else if (val > 0.0) else if (val > 0.0)
m_cur = 32767.0; m_cur = 32767.0;

View File

@ -37,9 +37,9 @@ namespace analog
m_VT = nlconst::magic(0.0258) * n; m_VT = nlconst::magic(0.0258) * n;
m_VT_inv = plib::reciprocal(m_VT); m_VT_inv = plib::reciprocal(m_VT);
} }
nl_fptype I(const nl_fptype V) const { return m_Is * std::exp(V * m_VT_inv) - m_Is; } nl_fptype I(const nl_fptype V) const { return m_Is * plib::exp(V * m_VT_inv) - m_Is; }
nl_fptype g(const nl_fptype V) const { return m_Is * m_VT_inv * std::exp(V * m_VT_inv); } nl_fptype g(const nl_fptype V) const { return m_Is * m_VT_inv * plib::exp(V * m_VT_inv); }
nl_fptype V(const nl_fptype I) const { return std::log1p(I / m_Is) * m_VT; } // log1p(x)=log(1.0 + x) nl_fptype V(const nl_fptype I) const { return plib::log1p(I / m_Is) * m_VT; } // log1p(x)=log(1.0 + x)
nl_fptype gI(const nl_fptype I) const { return m_VT_inv * (I + m_Is); } nl_fptype gI(const nl_fptype I) const { return m_VT_inv * (I + m_Is); }
private: private:

View File

@ -163,8 +163,8 @@ namespace analog
if (nVd > m_Vcrit) if (nVd > m_Vcrit)
{ {
const nl_fptype d = std::min(+fp_constants<nl_fptype>::DIODE_MAXDIFF(), nVd - m_Vd); const nl_fptype d = std::min(+fp_constants<nl_fptype>::DIODE_MAXDIFF(), nVd - m_Vd);
const nl_fptype a = std::abs(d) * m_VtInv; const nl_fptype a = plib::abs(d) * m_VtInv;
m_Vd = m_Vd + nlconst::magic(d < 0 ? -1.0 : 1.0) * std::log1p(a) * m_Vt; m_Vd = m_Vd + nlconst::magic(d < 0 ? -1.0 : 1.0) * plib::log1p(a) * m_Vt;
} }
else else
m_Vd = std::max(-fp_constants<nl_fptype>::DIODE_MAXDIFF(), nVd); m_Vd = std::max(-fp_constants<nl_fptype>::DIODE_MAXDIFF(), nVd);
@ -177,7 +177,7 @@ namespace analog
} }
else else
{ {
IseVDVt = std::exp(m_logIs + m_Vd * m_VtInv); IseVDVt = plib::exp(m_logIs + m_Vd * m_VtInv);
m_Id = IseVDVt - m_Is; m_Id = IseVDVt - m_Is;
m_G = IseVDVt * m_VtInv + m_gmin; m_G = IseVDVt * m_VtInv + m_gmin;
} }
@ -193,7 +193,7 @@ namespace analog
else /* log stepping should already be done in mosfet */ else /* log stepping should already be done in mosfet */
{ {
m_Vd = nVd; m_Vd = nVd;
IseVDVt = std::exp(std::min(+fp_constants<nl_fptype>::DIODE_MAXVOLT(), m_logIs + m_Vd * m_VtInv)); IseVDVt = plib::exp(std::min(+fp_constants<nl_fptype>::DIODE_MAXVOLT(), m_logIs + m_Vd * m_VtInv));
m_Id = IseVDVt - m_Is; m_Id = IseVDVt - m_Is;
m_G = IseVDVt * m_VtInv + m_gmin; m_G = IseVDVt * m_VtInv + m_gmin;
} }
@ -203,7 +203,7 @@ namespace analog
void set_param(const nl_fptype Is, const nl_fptype n, nl_fptype gmin, nl_fptype temp) void set_param(const nl_fptype Is, const nl_fptype n, nl_fptype gmin, nl_fptype temp)
{ {
m_Is = Is; m_Is = Is;
m_logIs = std::log(Is); m_logIs = plib::log(Is);
m_n = n; m_n = n;
m_gmin = gmin; m_gmin = gmin;
@ -211,7 +211,7 @@ namespace analog
m_Vmin = nlconst::magic(-5.0) * m_Vt; m_Vmin = nlconst::magic(-5.0) * m_Vt;
m_Vcrit = m_Vt * std::log(m_Vt / m_Is / nlconst::sqrt2()); m_Vcrit = m_Vt * plib::log(m_Vt / m_Is / nlconst::sqrt2());
m_VtInv = nlconst::one() / m_Vt; m_VtInv = nlconst::one() / m_Vt;
//printf("%g %g\n", m_Vmin, m_Vcrit); //printf("%g %g\n", m_Vmin, m_Vcrit);
} }

View File

@ -263,7 +263,7 @@ namespace analog
else if (m_model.m_NSUB > nlconst::zero()) else if (m_model.m_NSUB > nlconst::zero())
{ {
nl_assert_always(m_model.m_NSUB * nlconst::magic(1e6) >= constants::NiSi(), "Error calculating phi for model " + m_model.name()); nl_assert_always(m_model.m_NSUB * nlconst::magic(1e6) >= constants::NiSi(), "Error calculating phi for model " + m_model.name());
m_phi = nlconst::two() * Vt * std::log (m_model.m_NSUB * nlconst::magic(1e6) / constants::NiSi()); m_phi = nlconst::two() * Vt * plib::log (m_model.m_NSUB * nlconst::magic(1e6) / constants::NiSi());
} }
else else
m_phi = nlconst::magic(0.6); m_phi = nlconst::magic(0.6);
@ -274,7 +274,7 @@ namespace analog
else else
{ {
if (Cox > nlconst::zero() && m_model.m_NSUB > nlconst::zero()) if (Cox > nlconst::zero() && m_model.m_NSUB > nlconst::zero())
m_gamma = std::sqrt (nlconst::two() m_gamma = plib::sqrt (nlconst::two()
* constants::Q_e() * constants::eps_Si() * constants::eps_0() * constants::Q_e() * constants::eps_Si() * constants::eps_0()
* m_model.m_NSUB * nlconst::magic(1e6)) / Cox; * m_model.m_NSUB * nlconst::magic(1e6)) / Cox;
else else
@ -418,8 +418,8 @@ namespace analog
else else
{ {
// linear // linear
const nl_fptype Sqr1 = static_cast<nl_fptype>(std::pow(Vdsat - Vds, 2)); const nl_fptype Sqr1 = static_cast<nl_fptype>(plib::pow(Vdsat - Vds, 2));
const nl_fptype Sqr2 = static_cast<nl_fptype>(std::pow(nlconst::two() * Vdsat - Vds, 2)); const nl_fptype Sqr2 = static_cast<nl_fptype>(plib::pow(nlconst::two() * Vdsat - Vds, 2));
Cgb = 0; Cgb = 0;
Cgs = m_CoxWL * (nlconst::one() - Sqr1 / Sqr2) * nlconst::magic(2.0 / 3.0); Cgs = m_CoxWL * (nlconst::one() - Sqr1 / Sqr2) * nlconst::magic(2.0 / 3.0);
Cgd = m_CoxWL * (nlconst::one() - Vdsat * Vdsat / Sqr2) * nlconst::magic(2.0 / 3.0); Cgd = m_CoxWL * (nlconst::one() - Vdsat * Vdsat / Sqr2) * nlconst::magic(2.0 / 3.0);
@ -452,9 +452,9 @@ namespace analog
const nl_fptype k = nlconst::magic(3.5); // see "Circuit Simulation", page 185 const nl_fptype k = nlconst::magic(3.5); // see "Circuit Simulation", page 185
nl_fptype d = (Vgs - m_Vgs); nl_fptype d = (Vgs - m_Vgs);
Vgs = m_Vgs + plib::reciprocal(k) * nlconst::magic(d < 0 ? -1.0 : 1.0) * std::log1p(k * std::abs(d)); Vgs = m_Vgs + plib::reciprocal(k) * nlconst::magic(d < 0 ? -1.0 : 1.0) * plib::log1p(k * plib::abs(d));
d = (Vgd - m_Vgd); d = (Vgd - m_Vgd);
Vgd = m_Vgd + plib::reciprocal(k) * nlconst::magic(d < 0 ? -1.0 : 1.0) * std::log1p(k * std::abs(d)); Vgd = m_Vgd + plib::reciprocal(k) * nlconst::magic(d < 0 ? -1.0 : 1.0) * plib::log1p(k * plib::abs(d));
m_Vgs = Vgs; m_Vgs = Vgs;
m_Vgd = Vgd; m_Vgd = Vgd;
@ -475,13 +475,13 @@ namespace analog
// calculate Vth // calculate Vth
const nl_fptype Vbulk = is_forward ? Vbs : Vbd; const nl_fptype Vbulk = is_forward ? Vbs : Vbd;
const nl_fptype phi_m_Vbulk = (m_phi > Vbulk) ? std::sqrt(m_phi - Vbulk) : nlconst::zero(); const nl_fptype phi_m_Vbulk = (m_phi > Vbulk) ? plib::sqrt(m_phi - Vbulk) : nlconst::zero();
const nl_fptype Vth = m_vto * m_polarity + m_gamma * (phi_m_Vbulk - std::sqrt(m_phi)); const nl_fptype Vth = m_vto * m_polarity + m_gamma * (phi_m_Vbulk - plib::sqrt(m_phi));
const nl_fptype Vctrl = (is_forward ? Vgs : Vgd) - Vth; const nl_fptype Vctrl = (is_forward ? Vgs : Vgd) - Vth;
nl_fptype Ids(0), gm(0), gds(0), gmb(0); nl_fptype Ids(0), gm(0), gds(0), gmb(0);
const nl_fptype absVds = std::abs(Vds); const nl_fptype absVds = plib::abs(Vds);
if (Vctrl <= nlconst::zero()) if (Vctrl <= nlconst::zero())
{ {

View File

@ -199,7 +199,7 @@ namespace netlist
{ {
const nl_fptype cVt = nlconst::magic(0.0258 * 1.0); // * m_n; const nl_fptype cVt = nlconst::magic(0.0258 * 1.0); // * m_n;
const nl_fptype cId = m_model.m_DAB; // 3 mA const nl_fptype cId = m_model.m_DAB; // 3 mA
const nl_fptype cVd = cVt * std::log(cId / nlconst::magic(1e-15) + nlconst::one()); const nl_fptype cVd = cVt * plib::log(cId / nlconst::magic(1e-15) + nlconst::one());
m_VH.push(m_VCC() - m_model.m_VLH - cVd); m_VH.push(m_VCC() - m_model.m_VLH - cVd);
m_VL.push(m_GND() + m_model.m_VLL + cVd); m_VL.push(m_GND() + m_model.m_VLL + cVd);

View File

@ -69,13 +69,13 @@ NETLIB_UPDATE_TERMINALS(LVCCS)
const nl_fptype vi = m_IP.net().Q_Analog() - m_IN.net().Q_Analog(); const nl_fptype vi = m_IP.net().Q_Analog() - m_IN.net().Q_Analog();
const auto c1(nlconst::magic(0.2)); const auto c1(nlconst::magic(0.2));
if (std::abs(m_mult / m_cur_limit() * vi) > nlconst::half()) if (plib::abs(m_mult / m_cur_limit() * vi) > nlconst::half())
m_vi = m_vi + c1 * std::tanh((vi - m_vi) / c1); m_vi = m_vi + c1 * plib::tanh((vi - m_vi) / c1);
else else
m_vi = vi; m_vi = vi;
const nl_fptype x = m_mult / m_cur_limit() * m_vi; const nl_fptype x = m_mult / m_cur_limit() * m_vi;
const nl_fptype X = std::tanh(x); const nl_fptype X = plib::tanh(x);
const nl_fptype beta = m_mult * (nlconst::one() - X*X); const nl_fptype beta = m_mult * (nlconst::one() - X*X);
const nl_fptype I = m_cur_limit() * X - beta * m_vi; const nl_fptype I = m_cur_limit() * X - beta * m_vi;

View File

@ -64,7 +64,7 @@ NETLIB_RESET(POT)
{ {
nl_fptype v = m_Dial(); nl_fptype v = m_Dial();
if (m_DialIsLog()) if (m_DialIsLog())
v = (std::exp(v) - nlconst::one()) / (std::exp(nlconst::one()) - nlconst::one()); v = (plib::exp(v) - nlconst::one()) / (plib::exp(nlconst::one()) - nlconst::one());
m_R1.set_R(std::max(m_R() * v, exec().gmin())); m_R1.set_R(std::max(m_R() * v, exec().gmin()));
m_R2.set_R(std::max(m_R() * (nlconst::one() - v), exec().gmin())); m_R2.set_R(std::max(m_R() * (nlconst::one() - v), exec().gmin()));
@ -77,7 +77,7 @@ NETLIB_UPDATE_PARAM(POT)
nl_fptype v = m_Dial(); nl_fptype v = m_Dial();
if (m_DialIsLog()) if (m_DialIsLog())
v = (std::exp(v) - nlconst::one()) / (std::exp(nlconst::one()) - nlconst::one()); v = (plib::exp(v) - nlconst::one()) / (plib::exp(nlconst::one()) - nlconst::one());
if (m_Reverse()) if (m_Reverse())
v = nlconst::one() - v; v = nlconst::one() - v;
m_R1.set_R(std::max(m_R() * v, exec().gmin())); m_R1.set_R(std::max(m_R() * v, exec().gmin()));
@ -94,7 +94,7 @@ NETLIB_RESET(POT2)
nl_fptype v = m_Dial(); nl_fptype v = m_Dial();
if (m_DialIsLog()) if (m_DialIsLog())
v = (std::exp(v) - nlconst::one()) / (std::exp(nlconst::one()) - nlconst::one()); v = (plib::exp(v) - nlconst::one()) / (plib::exp(nlconst::one()) - nlconst::one());
if (m_Reverse()) if (m_Reverse())
v = nlconst::one() - v; v = nlconst::one() - v;
m_R1.set_R(std::max(m_R() * v, exec().gmin())); m_R1.set_R(std::max(m_R() * v, exec().gmin()));
@ -108,7 +108,7 @@ NETLIB_UPDATE_PARAM(POT2)
nl_fptype v = m_Dial(); nl_fptype v = m_Dial();
if (m_DialIsLog()) if (m_DialIsLog())
v = (std::exp(v) - nlconst::one()) / (std::exp(nlconst::one()) - nlconst::one()); v = (plib::exp(v) - nlconst::one()) / (plib::exp(nlconst::one()) - nlconst::one());
if (m_Reverse()) if (m_Reverse())
v = nlconst::one() - v; v = nlconst::one() - v;
m_R1.set_R(std::max(m_R() * v, exec().gmin())); m_R1.set_R(std::max(m_R() * v, exec().gmin()));

View File

@ -52,7 +52,7 @@ endif
CFLAGS = $(LTO) -g -O3 -std=c++14 -I$(CURDIR)/.. -I$(CURDIR)/../.. $(CEXTRAFLAGS) CFLAGS = $(LTO) -g -O3 -std=c++14 -I$(CURDIR)/.. -I$(CURDIR)/../.. $(CEXTRAFLAGS)
LDFLAGS = $(LTO) -g -O3 -std=c++14 $(LDEXTRAFLAGS) LDFLAGS = $(LTO) -g -O3 -std=c++14 $(LDEXTRAFLAGS)
LIBS = -lpthread -ldl LIBS = -lpthread -ldl $(EXTRALIBS)
CC = g++ CC = g++
LD = @g++ LD = @g++
@ -235,7 +235,7 @@ native:
$(MAKE) CEXTRAFLAGS="-march=native -msse4.2 -Wall -Wpedantic -Wsign-compare -Wextra " $(MAKE) CEXTRAFLAGS="-march=native -msse4.2 -Wall -Wpedantic -Wsign-compare -Wextra "
gcc9: gcc9:
$(MAKE) CC=g++-9 LD=g++-9 CEXTRAFLAGS="-march=native -msse4.2 -Wall -Wpedantic -Wsign-compare -Wextra " $(MAKE) CC=g++-9 LD=g++-9 CEXTRAFLAGS="-march=native -msse4.2 -Wall -Wpedantic -Wsign-compare -Wextra" EXTRALIBS="-lquadmath"
clang: clang:
$(MAKE) CC=clang++-10 LD=clang++-10 CEXTRAFLAGS="-march=native -msse4.2 -Weverything -Werror -Wno-padded -Wno-weak-vtables -Wno-unused-template -Wno-missing-variable-declarations -Wno-float-equal -Wconversion -Wno-c++98-compat -Wno-c++98-compat-pedantic -Wno-format-nonliteral -Wno-exit-time-destructors" $(MAKE) CC=clang++-10 LD=clang++-10 CEXTRAFLAGS="-march=native -msse4.2 -Weverything -Werror -Wno-padded -Wno-weak-vtables -Wno-unused-template -Wno-missing-variable-declarations -Wno-float-equal -Wconversion -Wno-c++98-compat -Wno-c++98-compat-pedantic -Wno-format-nonliteral -Wno-exit-time-destructors"

View File

@ -255,7 +255,7 @@ namespace netlist
NETLIB_RESET(74123) NETLIB_RESET(74123)
{ {
m_KP = plib::reciprocal(nlconst::one() + std::exp(m_K())); m_KP = plib::reciprocal(nlconst::one() + plib::exp(m_K()));
m_RP.reset(); m_RP.reset();
m_RN.reset(); m_RN.reset();

View File

@ -107,13 +107,17 @@
static constexpr const auto NETLIST_INTERNAL_RES = 1000000000; static constexpr const auto NETLIST_INTERNAL_RES = 1000000000;
static constexpr const auto NETLIST_CLOCK = NETLIST_INTERNAL_RES; static constexpr const auto NETLIST_CLOCK = NETLIST_INTERNAL_RES;
//============================================================ /*! Floating point types used
// Floating point types used *
// * nl_fptype is the floating point type used throughout the
// Don't change this. Simple analog circuits like pong * netlist core.
// work with float. Kidniki just doesn't work at all *
// due to numeric issues * Don't change this! Simple analog circuits like pong
//============================================================ * work with float. Kidniki just doesn't work at all.
*
* FIXME: More work needed. Review magic numbers.
*
*/
using nl_fptype = double; using nl_fptype = double;
@ -162,8 +166,8 @@ namespace netlist
template <> template <>
struct fp_constants<float> struct fp_constants<float>
{ {
static inline constexpr float DIODE_MAXDIFF() noexcept { return 1e2; } static inline constexpr float DIODE_MAXDIFF() noexcept { return 1e20; }
static inline constexpr float DIODE_MAXVOLT() noexcept { return 20.0; } static inline constexpr float DIODE_MAXVOLT() noexcept { return 90.0; }
static inline constexpr float TIMESTEP_MAXDIFF() noexcept { return 1e30f; } static inline constexpr float TIMESTEP_MAXDIFF() noexcept { return 1e30f; }
static inline constexpr float TIMESTEP_MINDIV() noexcept { return 1e-8f; } static inline constexpr float TIMESTEP_MINDIV() noexcept { return 1e-8f; }

View File

@ -121,8 +121,8 @@ namespace netlist
void nlparse_t::register_param_x(const pstring &param, const nl_fptype value) void nlparse_t::register_param_x(const pstring &param, const nl_fptype value)
{ {
if (std::abs(value - std::floor(value)) > nlconst::magic(1e-30) if (plib::abs(value - plib::floor(value)) > nlconst::magic(1e-30)
|| std::abs(value) > nlconst::magic(1e9)) || plib::abs(value) > nlconst::magic(1e9))
register_param(param, plib::pfmt("{1:.9}").e(value)); register_param(param, plib::pfmt("{1:.9}").e(value));
else else
register_param(param, plib::pfmt("{1}")(static_cast<long>(value))); register_param(param, plib::pfmt("{1}")(static_cast<long>(value)));
@ -1126,7 +1126,7 @@ void setup_t::prepare_to_run()
//FIXME: check for errors ... //FIXME: check for errors ...
bool err(false); bool err(false);
auto v = plib::pstonum_ne<nl_fptype>(p->second, err); auto v = plib::pstonum_ne<nl_fptype>(p->second, err);
if (err || std::abs(v - std::floor(v)) > nlconst::magic(1e-6) ) if (err || plib::abs(v - plib::floor(v)) > nlconst::magic(1e-6) )
log().fatal(MF_HND_VAL_NOT_SUPPORTED(p->second)); log().fatal(MF_HND_VAL_NOT_SUPPORTED(p->second));
// FIXME comparison with zero // FIXME comparison with zero
d.second->set_hint_deactivate(v == nlconst::zero()); d.second->set_hint_deactivate(v == nlconst::zero());

View File

@ -247,6 +247,7 @@ namespace netlist
// FIXME: quick hack // FIXME: quick hack
void register_param_x(const pstring &param, const nl_fptype value); void register_param_x(const pstring &param, const nl_fptype value);
template <typename T> template <typename T>
typename std::enable_if<std::is_floating_point<T>::value || std::is_integral<T>::value>::type typename std::enable_if<std::is_floating_point<T>::value || std::is_integral<T>::value>::type
register_param(const pstring &param, T value) register_param(const pstring &param, T value)
@ -254,6 +255,13 @@ namespace netlist
register_param_x(param, static_cast<nl_fptype>(value)); register_param_x(param, static_cast<nl_fptype>(value));
} }
#if PUSE_FLOAT128
void register_param(const pstring &param, __float128 value)
{
register_param_x(param, static_cast<nl_fptype>(value));
}
#endif
void register_lib_entry(const pstring &name, const pstring &sourcefile); void register_lib_entry(const pstring &name, const pstring &sourcefile);
void register_frontier(const pstring &attach, const nl_fptype r_IN, const nl_fptype r_OUT); void register_frontier(const pstring &attach, const nl_fptype r_IN, const nl_fptype r_OUT);

View File

@ -192,6 +192,11 @@ public:
typename std::enable_if<std::is_floating_point<T>::value, pfmt &>::type typename std::enable_if<std::is_floating_point<T>::value, pfmt &>::type
e(const T &x) {return format_element('e', x); } e(const T &x) {return format_element('e', x); }
#if PUSE_FLOAT128
// FIXME: not what we want
pfmt & e(const __float128 &x) {return format_element('e', static_cast<long double>(x)); }
#endif
template <typename T> template <typename T>
typename std::enable_if<std::is_floating_point<T>::value, pfmt &>::type typename std::enable_if<std::is_floating_point<T>::value, pfmt &>::type
g(const T &x) {return format_element('g', x); } g(const T &x) {return format_element('g', x); }

View File

@ -206,10 +206,10 @@ namespace plib {
OP(MULT, 1, ST2 * ST1) OP(MULT, 1, ST2 * ST1)
OP(SUB, 1, ST2 - ST1) OP(SUB, 1, ST2 - ST1)
OP(DIV, 1, ST2 / ST1) OP(DIV, 1, ST2 / ST1)
OP(POW, 1, std::pow(ST2, ST1)) OP(POW, 1, plib::pow(ST2, ST1))
OP(SIN, 0, std::sin(ST2)) OP(SIN, 0, plib::sin(ST2))
OP(COS, 0, std::cos(ST2)) OP(COS, 0, plib::cos(ST2))
OP(TRUNC, 0, std::trunc(ST2)) OP(TRUNC, 0, plib::trunc(ST2))
case RAND: case RAND:
stack[ptr++] = lfsr_random(); stack[ptr++] = lfsr_random();
break; break;
@ -227,5 +227,8 @@ namespace plib {
template class pfunction<float>; template class pfunction<float>;
template class pfunction<double>; template class pfunction<double>;
template class pfunction<long double>; template class pfunction<long double>;
#if (PUSE_FLOAT128)
template class pfunction<__float128>;
#endif
} // namespace plib } // namespace plib

View File

@ -118,7 +118,9 @@ namespace plib {
extern template class pfunction<float>; extern template class pfunction<float>;
extern template class pfunction<double>; extern template class pfunction<double>;
extern template class pfunction<long double>; extern template class pfunction<long double>;
#if (PUSE_FLOAT128)
extern template class pfunction<__float128>;
#endif
} // namespace plib } // namespace plib
#endif /* PEXCEPTION_H_ */ #endif /* PEXCEPTION_H_ */

View File

@ -100,13 +100,21 @@ namespace plib
constexpr internal_type as_raw() const noexcept { return m_time; } constexpr internal_type as_raw() const noexcept { return m_time; }
template <typename FT> template <typename FT, typename = std::enable_if<std::is_floating_point<FT>::value, FT>>
constexpr typename std::enable_if<std::is_floating_point<FT>::value, FT>::type constexpr FT
as_fp() const noexcept as_fp() const noexcept
{ {
return static_cast<FT>(m_time) * inv_res<FT>(); return static_cast<FT>(m_time) * inv_res<FT>();
} }
#if PUSE_FLOAT128
constexpr __float128
as_fp() const noexcept
{
return static_cast<__float128>(m_time) * inv_res<__float128>();
}
#endif
constexpr double as_double() const noexcept { return as_fp<double>(); } constexpr double as_double() const noexcept { return as_fp<double>(); }
constexpr double as_float() const noexcept { return as_fp<float>(); } constexpr double as_float() const noexcept { return as_fp<float>(); }
constexpr double as_long_double() const noexcept { return as_fp<long double>(); } constexpr double as_long_double() const noexcept { return as_fp<long double>(); }
@ -122,8 +130,12 @@ namespace plib
static constexpr const ptime from_raw(const internal_type raw) noexcept { return ptime(raw); } static constexpr const ptime from_raw(const internal_type raw) noexcept { return ptime(raw); }
template <typename FT> template <typename FT>
static constexpr const typename std::enable_if<std::is_floating_point<FT>::value, ptime>::type static constexpr const typename std::enable_if<std::is_floating_point<FT>::value
from_fp(const FT t) noexcept { return ptime(static_cast<internal_type>(std::floor(t * static_cast<FT>(RES) + static_cast<FT>(0.5))), RES); } #if PUSE_FLOAT128
|| std::is_same<FT, __float128>::value
#endif
, ptime>::type
from_fp(const FT t) noexcept { return ptime(static_cast<internal_type>(plib::floor(t * static_cast<FT>(RES) + static_cast<FT>(0.5))), RES); }
static constexpr const ptime from_double(const double t) noexcept static constexpr const ptime from_double(const double t) noexcept
{ return from_fp<double>(t); } { return from_fp<double>(t); }

View File

@ -314,9 +314,10 @@ namespace plib
/*! hypot function /*! hypot function
* *
* @tparam T type of the argument * @tparam T type of the arguments
* @param v argument * @param v1 first argument
* @return absolute value of argument * @param v1 second argument
* @return sqrt(v1*v1+v2*v2)
*/ */
template <typename T> template <typename T>
static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type
@ -325,6 +326,127 @@ namespace plib
return std::hypot(v1, v2); return std::hypot(v1, v2);
} }
/*! exp function
*
* @tparam T type of the argument
* @param v argument
* @return exp(v)
*/
template <typename T>
static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type
exp(T v) noexcept
{
return std::exp(v);
}
/*! log function
*
* @tparam T type of the argument
* @param v argument
* @return log(v)
*/
template <typename T>
static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type
log(T v) noexcept
{
return std::log(v);
}
/*! tanh function
*
* @tparam T type of the argument
* @param v argument
* @return tanh(v)
*/
template <typename T>
static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type
tanh(T v) noexcept
{
return std::tanh(v);
}
/*! floor function
*
* @tparam T type of the argument
* @param v argument
* @return floor(v)
*/
template <typename T>
static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type
floor(T v) noexcept
{
return std::floor(v);
}
/*! log1p function
*
* @tparam T type of the argument
* @param v argument
* @return log(1 + v)
*/
template <typename T>
static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type
log1p(T v) noexcept
{
return std::log1p(v);
}
/*! sin function
*
* @tparam T type of the argument
* @param v argument
* @return sin(v)
*/
template <typename T>
static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type
sin(T v) noexcept
{
return std::sin(v);
}
/*! cos function
*
* @tparam T type of the argument
* @param v argument
* @return cos(v)
*/
template <typename T>
static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type
cos(T v) noexcept
{
return std::cos(v);
}
/*! trunc function
*
* @tparam T type of the argument
* @param v argument
* @return trunc(v)
*/
template <typename T>
static inline constexpr typename std::enable_if<std::is_floating_point<T>::value, T>::type
trunc(T v) noexcept
{
return std::trunc(v);
}
/*! pow function
*
* @tparam T1 type of the first argument
* @tparam T2 type of the second argument
* @param v argument
* @param p power
* @return v^p
*
* FIXME: limited implementation
*/
template <typename T1, typename T2>
static inline T1
pow(T1 v, T2 p) noexcept
{
return std::pow(v, p);
}
#if (PUSE_FLOAT128) #if (PUSE_FLOAT128)
static inline constexpr __float128 reciprocal(__float128 v) noexcept static inline constexpr __float128 reciprocal(__float128 v) noexcept
{ {
@ -345,6 +467,61 @@ namespace plib
{ {
return hypotq(v1, v2); return hypotq(v1, v2);
} }
static inline __float128 exp(__float128 v) noexcept
{
return expq(v);
}
static inline __float128 log(__float128 v) noexcept
{
return logq(v);
}
static inline __float128 tanh(__float128 v) noexcept
{
return tanhq(v);
}
static inline __float128 floor(__float128 v) noexcept
{
return floorq(v);
}
static inline __float128 log1p(__float128 v) noexcept
{
return log1pq(v);
}
static inline __float128 sin(__float128 v) noexcept
{
return sinq(v);
}
static inline __float128 cos(__float128 v) noexcept
{
return cosq(v);
}
static inline __float128 trunc(__float128 v) noexcept
{
return truncq(v);
}
template <typename T>
static inline __float128 pow(__float128 v, T p) noexcept
{
return powq(v, static_cast<__float128>(p));
}
static inline __float128 pow(__float128 v, int p) noexcept
{
if (p==2)
return v*v;
else
return powq(v, static_cast<__float128>(p));
}
#endif #endif
static_assert(noexcept(constants<double>::one()) == true, "Not evaluated as constexpr"); static_assert(noexcept(constants<double>::one()) == true, "Not evaluated as constexpr");
@ -430,6 +607,19 @@ namespace plib
} }
}; };
#if PUSE_FLOAT128
template<>
struct pstonum_helper<__float128>
{
// FIXME: use strtoflt128 from quadmath.h
template <typename S>
__float128 operator()(std::locale loc, const S &arg, std::size_t *idx)
{
return static_cast<__float128>(pstonum_locale<long double>(loc, arg, idx));
}
};
#endif
template<typename T, typename S> template<typename T, typename S>
T pstonum(const S &arg, const std::locale &loc = std::locale::classic()) T pstonum(const S &arg, const std::locale &loc = std::locale::classic())
{ {

View File

@ -514,7 +514,7 @@ namespace solver
if (colu==row) colu = static_cast<unsigned>(diag); if (colu==row) colu = static_cast<unsigned>(diag);
else if (colu==diag) colu = static_cast<unsigned>(row); else if (colu==diag) colu = static_cast<unsigned>(row);
weight = weight + std::abs(static_cast<nl_fptype>(colu) - static_cast<nl_fptype>(diag)); weight = weight + plib::abs(static_cast<nl_fptype>(colu) - static_cast<nl_fptype>(diag));
touched[colu] = true; touched[colu] = true;
} }
} }

View File

@ -476,8 +476,8 @@ namespace solver
m_h_n_m_1[k] = hn; m_h_n_m_1[k] = hn;
m_DD_n_m_1[k] = DD_n; m_DD_n_m_1[k] = DD_n;
if (std::fabs(DD2) > fp_constants<nl_fptype>::TIMESTEP_MINDIV()) // avoid div-by-zero if (plib::abs(DD2) > fp_constants<nl_fptype>::TIMESTEP_MINDIV()) // avoid div-by-zero
new_net_timestep = std::sqrt(m_params.m_dynamic_lte / std::fabs(nlconst::magic(0.5)*DD2)); new_net_timestep = plib::sqrt(m_params.m_dynamic_lte / plib::abs(nlconst::magic(0.5)*DD2));
else else
new_net_timestep = m_params.m_max_timestep; new_net_timestep = m_params.m_max_timestep;

View File

@ -96,7 +96,7 @@ namespace solver
if (this->m_params.m_use_gabs) if (this->m_params.m_use_gabs)
{ {
for (std::size_t i = 0; i < term_count; i++) for (std::size_t i = 0; i < term_count; i++)
gabs_t = gabs_t + std::abs(go[i]); gabs_t = gabs_t + plib::abs(go[i]);
gabs_t *= nlconst::magic(0.5); // derived by try and error gabs_t *= nlconst::magic(0.5); // derived by try and error
if (gabs_t <= gtot_t) if (gabs_t <= gtot_t)