Sprinter/mame
2
0
Fork 0
mirror of https://github.com/holub/mame synced 2025-06-01 02:21:48 +03:00
Code Issues Actions 3 Packages Projects Releases Wiki Activity

Introduced a nl_math static class to have all math related

Browse Source 
functions in one place. This is a portability exercise, e.g.
log1p is not supported on all platforms. (nw)
This commit is contained in:
couriersud 2015-05-10 18:34:03 +02:00
parent e923df8137
commit 423def0ccc
11 changed files with 106 additions and 83 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/emu/netlist/analog/nld_bjt.c
View File

@ -25,9 +25,9 @@ public:
m_VT = 0.0258 * n;
m_VT_inv = 1.0 / m_VT;
}
nl_double I(const nl_double V) const { return m_Is * exp(V * m_VT_inv) - m_Is; }
nl_double g(const nl_double V) const { return m_Is * m_VT_inv * exp(V * m_VT_inv); }
nl_double V(const nl_double I) const { return log(1.0 + I / m_Is) * m_VT; }
nl_double I(const nl_double V) const { return m_Is * nl_math::exp(V * m_VT_inv) - m_Is; }
nl_double g(const nl_double V) const { return m_Is * m_VT_inv * nl_math::exp(V * m_VT_inv); }
nl_double V(const nl_double I) const { return nl_math::log1p(I / m_Is) * m_VT; } // log1p(x)=log(1.0 + x)
nl_double gI(const nl_double I) const { return m_VT_inv * (I + m_Is); }
private:

10
src/emu/netlist/analog/nld_ms_direct.h
View File

@ -106,8 +106,8 @@ ATTR_HOT nl_double netlist_matrix_solver_direct_t<m_N, _storage_N>::compute_next
n->m_h_n_m_1 = hn;
n->m_DD_n_m_1 = DD_n;
if (fabs(DD2) > 1e-50) // avoid div-by-zero
new_net_timestep = sqrt(m_params.m_lte / fabs(0.5*DD2));
if (nl_math::abs(DD2) > 1e-50) // avoid div-by-zero
new_net_timestep = nl_math::sqrt(m_params.m_lte / nl_math::abs(0.5*DD2));
else
new_net_timestep = m_params.m_max_timestep;
@ -290,7 +290,7 @@ ATTR_HOT void netlist_matrix_solver_direct_t<m_N, _storage_N>::LE_solve()
int maxrow = i;
for (int j = i + 1; j < kN; j++)
{
if (fabs(m_A[j][i]) > fabs(m_A[maxrow][i]))
if (nl_math::abs(m_A[j][i]) > nl_math::abs(m_A[maxrow][i]))
maxrow = j;
}
@ -359,8 +359,8 @@ ATTR_HOT nl_double netlist_matrix_solver_direct_t<m_N, _storage_N>::delta(
nl_double cerr2 = 0;
for (int i = 0; i < this->N(); i++)
{
const nl_double e = fabs(V[i] - this->m_nets[i]->m_cur_Analog);
const nl_double e2 = fabs(m_RHS[i] - this->m_last_RHS[i]);
const nl_double e = nl_math::abs(V[i] - this->m_nets[i]->m_cur_Analog);
const nl_double e2 = nl_math::abs(m_RHS[i] - this->m_last_RHS[i]);
cerr = (e > cerr ? e : cerr);
cerr2 = (e2 > cerr2 ? e2 : cerr2);
}

2
src/emu/netlist/analog/nld_ms_direct1.h
View File

@ -44,7 +44,7 @@ ATTR_HOT inline int netlist_matrix_solver_direct1_t::vsolve_non_dynamic()
nl_double new_val = m_RHS[0] / m_A[0][0];
nl_double e = (new_val - net->m_cur_Analog);
nl_double cerr = fabs(e);
nl_double cerr = nl_math::abs(e);
net->m_cur_Analog = new_val;

30
src/emu/netlist/analog/nld_ms_jakobi.h
View File

@ -104,7 +104,7 @@ ATTR_HOT nl_double netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve()
this->m_Vdelta[k] = nv;
}
if (sqo > 1e-90)
m_lp_fact = std::min(sqrt(sq/sqo), 2.0);
m_lp_fact = std::min(nl_math::sqrt(sq/sqo), 2.0);
else
m_lp_fact = 0.0;
}
@ -143,7 +143,7 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
for (int i = 0; i < iN; i++)
{
frob += this->m_A[k][i] * this->m_A[k][i];
s = s + fabs(this->m_A[k][i]);
s = s + nl_math::abs(this->m_A[k][i]);
}
if (s<rmin)
@ -152,9 +152,9 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
rmax = s;
}
#if 0
nl_double frobA = sqrt(frob /(iN));
nl_double frobA = nl_math::sqrt(frob /(iN));
if (1 &&frobA < 1.0)
//ws = 2.0 / (1.0 + sqrt(1.0-frobA));
//ws = 2.0 / (1.0 + nl_math::sqrt(1.0-frobA));
ws = 2.0 / (2.0 - frobA);
else
ws = 1.0;
@ -169,7 +169,7 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
// suitable parameter.
nl_double rm = (rmax + rmin) * 0.5;
if (rm < 1.0)
ws = 2.0 / (1.0 + sqrt(1.0-rm));
ws = 2.0 / (1.0 + nl_math::sqrt(1.0-rm));
else
ws = 1.0;
if (ws > 1.02 && rmax > 1.001)
@ -200,13 +200,13 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
{
//if (i < k) frobL += this->m_A[k][i] * this->m_A[k][i] / this->m_A[k][k] /this-> m_A[k][k];
//if (i > k) frobU += this->m_A[k][i] * this->m_A[k][i] / this->m_A[k][k] / this->m_A[k][k];
//norm_t += fabs(this->m_A[k][i]);
//norm_t += nl_math::abs(this->m_A[k][i]);
}
//if (norm_t > norm) norm = norm_t;
const nl_double new_val = (1.0-ws) * new_v[k] + ws * (this->m_RHS[k] - Idrive + this->m_A[k][k] * new_v[k]) / this->m_A[k][k];
const nl_double e = fabs(new_val - new_v[k]);
const nl_double e = nl_math::abs(new_val - new_v[k]);
cerr = (e > cerr ? e : cerr);
new_v[k] = new_val;
}
@ -216,10 +216,10 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
resched = true;
}
resched_cnt++;
//ATTR_UNUSED nl_double frobUL = sqrt((frobU + frobL) / (double) (iN) / (double) (iN));
//ATTR_UNUSED nl_double frobUL = nl_math::sqrt((frobU + frobL) / (double) (iN) / (double) (iN));
} while (resched && (resched_cnt < this->m_params.m_gs_loops));
//printf("Frobenius %f %f %f %f %f\n", sqrt(frobU), sqrt(frobL), frobUL, frobA, norm);
//printf("Omega Estimate1 %f %f\n", 2.0 / (1.0 + sqrt(1-frobUL)), 2.0 / (1.0 + sqrt(1-frobA)) ); // printf("Frobenius %f\n", sqrt(frob / (double) (iN * iN) ));
//printf("Frobenius %f %f %f %f %f\n", nl_math::sqrt(frobU), nl_math::sqrt(frobL), frobUL, frobA, norm);
//printf("Omega Estimate1 %f %f\n", 2.0 / (1.0 + nl_math::sqrt(1-frobUL)), 2.0 / (1.0 + nl_math::sqrt(1-frobA)) ); // printf("Frobenius %f\n", sqrt(frob / (double) (iN * iN) ));
//printf("Omega Estimate2 %f %f\n", 2.0 / (2.0 - frobUL), 2.0 / (2.0 - frobA) ); // printf("Frobenius %f\n", sqrt(frob / (double) (iN * iN) ));
@ -250,7 +250,7 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
*
* and estimate using
*
* omega = 2.0 / (1.0 + sqrt(1-rho))
* omega = 2.0 / (1.0 + nl_math::sqrt(1-rho))
*/
const nl_double ws = this->m_params.m_sor; //1.045; //2.0 / (1.0 + /*sin*/(3.14159 * 5.5 / (double) (m_nets.count()+1)));
@ -284,7 +284,7 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
}
if (USE_GABS)
for (int i = 0; i < term_count; i++)
gabs_t = gabs_t + fabs(go[i]);
gabs_t = gabs_t + nl_math::abs(go[i]);
for (int i = this->m_terms[k]->m_railstart; i < term_count; i++)
RHS_t = RHS_t + go[i] * *other_cur_analog[i];
@ -292,7 +292,7 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
RHS[k] = RHS_t;
//if (fabs(gabs_t - fabs(gtot_t)) > 1e-20)
//if (nl_math::abs(gabs_t - nl_math::abs(gtot_t)) > 1e-20)
// printf("%d %e abs: %f tot: %f\n",k, gabs_t / gtot_t -1.0, gabs_t, gtot_t);
gabs_t *= 0.95; // avoid rounding issues
@ -329,9 +329,9 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
Idrive = Idrive + go[i] * old_V[net_other[i]];
//const nl_double new_val = new_V[k] * one_m_w[k] + (Idrive + RHS[k]) * w[k];
//resched = resched || (std::abs(new_val - new_V[k]) > accuracy);
//resched = resched || (nl_math::abs(new_val - new_V[k]) > accuracy);
const nl_double new_val = old_V[k] * one_m_w[k] + (Idrive + RHS[k]) * w[k];
resched = resched || (std::abs(new_val - old_V[k]) > accuracy);
resched = resched || (nl_math::abs(new_val - old_V[k]) > accuracy);
new_V[k] = new_val;
}
for (int k = 0; k < iN; k++)

6
src/emu/netlist/analog/nld_ms_sor.h
View File

@ -94,7 +94,7 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
*
* and estimate using
*
* omega = 2.0 / (1.0 + sqrt(1-rho))
* omega = 2.0 / (1.0 + nl_math::sqrt(1-rho))
*/
const nl_double ws = this->m_params.m_sor; //1.045; //2.0 / (1.0 + /*sin*/(3.14159 * 5.5 / (double) (m_nets.count()+1)));
@ -128,7 +128,7 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
if (USE_GABS)
for (int i = 0; i < term_count; i++)
gabs_t = gabs_t + fabs(go[i]);
gabs_t = gabs_t + nl_math::abs(go[i]);
for (int i = this->m_terms[k]->m_railstart; i < term_count; i++)
RHS_t = RHS_t + go[i] * *other_cur_analog[i];
@ -174,7 +174,7 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_t<m_N, _storage_N>::vsolve_non_dyn
const nl_double new_val = new_V[k] * one_m_w[k] + (Idrive + RHS[k]) * w[k];
resched = resched || (std::abs(new_val - new_V[k]) > accuracy);
resched = resched || (nl_math::abs(new_val - new_V[k]) > accuracy);
new_V[k] = new_val;
}

28
src/emu/netlist/analog/nld_ms_sor_mat.h
View File

@ -114,7 +114,7 @@ ATTR_HOT nl_double netlist_matrix_solver_SOR_mat_t<m_N, _storage_N>::vsolve()
this->m_Vdelta[k] = nv;
}
if (sqo > 1e-90)
m_lp_fact = std::min(sqrt(sq/sqo), 2.0);
m_lp_fact = std::min(nl_math::sqrt(sq/sqo), (nl_double) 2.0);
else
m_lp_fact = 0.0;
}
@ -154,16 +154,16 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_mat_t<m_N, _storage_N>::vsolve_non
for (int k = 0; k < iN; k++)
{
#if 0
nl_double akk = fabs(this->m_A[k][k]);
nl_double akk = nl_math::abs(this->m_A[k][k]);
if ( akk > lambdaN)
lambdaN = akk;
if (akk < lambda1)
lambda1 = akk;
#else
nl_double akk = fabs(this->m_A[k][k]);
nl_double akk = nl_math::abs(this->m_A[k][k]);
nl_double s = 0.0;
for (int i=0; i<iN; i++)
s = s + fabs(this->m_A[k][i]);
s = s + nl_math::abs(this->m_A[k][i]);
akk = s / akk - 1.0;
//if ( akk > lambdaN)
// lambdaN = akk;
@ -194,7 +194,7 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_mat_t<m_N, _storage_N>::vsolve_non
for (int i = 0; i < iN; i++)
{
frob += this->m_A[k][i] * this->m_A[k][i];
s = s + fabs(this->m_A[k][i]);
s = s + nl_math::abs(this->m_A[k][i]);
}
if (s<rmin)
@ -203,9 +203,9 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_mat_t<m_N, _storage_N>::vsolve_non
rmax = s;
}
#if 0
nl_double frobA = sqrt(frob /(iN));
nl_double frobA = nl_math::sqrt(frob /(iN));
if (1 &&frobA < 1.0)
//ws = 2.0 / (1.0 + sqrt(1.0-frobA));
//ws = 2.0 / (1.0 + nl_math::sqrt(1.0-frobA));
ws = 2.0 / (2.0 - frobA);
else
ws = 1.0;
@ -220,7 +220,7 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_mat_t<m_N, _storage_N>::vsolve_non
// suitable parameter.
nl_double rm = (rmax + rmin) * 0.5;
if (rm < 1.0)
ws = 2.0 / (1.0 + sqrt(1.0-rm));
ws = 2.0 / (1.0 + nl_math::sqrt(1.0-rm));
else
ws = 1.0;
if (ws > 1.02 && rmax > 1.001)
@ -252,19 +252,19 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_mat_t<m_N, _storage_N>::vsolve_non
{
if (i < k) frobL += this->m_A[k][i] * this->m_A[k][i] / this->m_A[k][k] /this-> m_A[k][k];
if (i > k) frobU += this->m_A[k][i] * this->m_A[k][i] / this->m_A[k][k] / this->m_A[k][k];
norm_t += fabs(this->m_A[k][i]);
norm_t += nl_math::abs(this->m_A[k][i]);
}
#endif
//if (norm_t > norm) norm = norm_t;
#if 0
const nl_double new_val = (1.0-ws) * new_v[k] + ws * (this->m_RHS[k] - Idrive + this->m_A[k][k] * new_v[k]) / this->m_A[k][k];
const nl_double e = fabs(new_val - new_v[k]);
const nl_double e = nl_math::abs(new_val - new_v[k]);
cerr = (e > cerr ? e : cerr);
new_v[k] = new_val;
#else
const nl_double delta = m_omega * (this->m_RHS[k] - Idrive) / this->m_A[k][k];
const nl_double adelta = std::abs(delta);
const nl_double adelta = nl_math::abs(delta);
cerr = (adelta > cerr ? adelta : cerr);
new_v[k] += delta;
#endif
@ -275,10 +275,10 @@ ATTR_HOT inline int netlist_matrix_solver_SOR_mat_t<m_N, _storage_N>::vsolve_non
resched = true;
}
resched_cnt++;
//ATTR_UNUSED nl_double frobUL = sqrt((frobU + frobL) / (double) (iN) / (double) (iN));
//ATTR_UNUSED nl_double frobUL = nl_math::sqrt((frobU + frobL) / (double) (iN) / (double) (iN));
} while (resched && (resched_cnt < this->m_params.m_gs_loops));
//printf("Frobenius %f %f %f %f %f\n", sqrt(frobU), sqrt(frobL), frobUL, frobA, norm);
//printf("Omega Estimate1 %f %f\n", 2.0 / (1.0 + sqrt(1-frobUL)), 2.0 / (1.0 + sqrt(1-frobA)) ); // printf("Frobenius %f\n", sqrt(frob / (double) (iN * iN) ));
//printf("Frobenius %f %f %f %f %f\n", nl_math::sqrt(frobU), nl_math::sqrt(frobL), frobUL, frobA, norm);
//printf("Omega Estimate1 %f %f\n", 2.0 / (1.0 + nl_math::sqrt(1-frobUL)), 2.0 / (1.0 + nl_math::sqrt(1-frobA)) ); // printf("Frobenius %f\n", sqrt(frob / (double) (iN * iN) ));
//printf("Omega Estimate2 %f %f\n", 2.0 / (2.0 - frobUL), 2.0 / (2.0 - frobA) ); // printf("Frobenius %f\n", sqrt(frob / (double) (iN * iN) ));

7
src/emu/netlist/analog/nld_twoterm.c
View File

@ -20,13 +20,14 @@ ATTR_COLD netlist_generic_diode::netlist_generic_diode()
ATTR_COLD void netlist_generic_diode::set_param(const nl_double Is, const nl_double n, nl_double gmin)
{
static const int csqrt2 = nl_math::sqrt(2.0);
m_Is = Is;
m_n = n;
m_gmin = gmin;
m_Vt = 0.0258 * m_n;
m_Vcrit = m_Vt * log(m_Vt / m_Is / sqrt(2.0));
m_Vcrit = m_Vt * nl_math::log(m_Vt / m_Is / csqrt2);
m_VtInv = 1.0 / m_Vt;
}
@ -162,7 +163,7 @@ NETLIB_UPDATE_PARAM(POT)
{
nl_double v = m_Dial.Value();
if (m_DialIsLog.Value())
v = (exp(v) - 1.0) / (exp(1.0) - 1.0);
v = (nl_math::exp(v) - 1.0) / (nl_math::exp(1.0) - 1.0);
// FIXME: Only attached nets should be brought up to current time
//netlist().solver()->update_to_current_time(); // bring up current time
@ -208,7 +209,7 @@ NETLIB_UPDATE_PARAM(POT2)
nl_double v = m_Dial.Value();
if (m_DialIsLog.Value())
v = (exp(v) - 1.0) / (exp(1.0) - 1.0);
v = (nl_math::exp(v) - 1.0) / (nl_math::exp(1.0) - 1.0);
if (m_Reverse.Value())
v = 1.0 - v;

38
src/emu/netlist/analog/nld_twoterm.h
View File

@ -226,20 +226,17 @@ public:
{
m_Vd = nVd;
const nl_double eVDVt = exp(m_Vd * m_VtInv);
const nl_double eVDVt = nl_math::exp(m_Vd * m_VtInv);
m_Id = m_Is * (eVDVt - 1.0);
m_G = m_Is * m_VtInv * eVDVt + m_gmin;
}
else
{
#if defined(_MSC_VER) && _MSC_VER < 1800
m_Vd = m_Vd + log((nVd - m_Vd) * m_VtInv + 1.0) * m_Vt;
#else
nl_double a = (nVd - m_Vd) * m_VtInv;
if (a<1e-12 - 1.0) a = 1e-12 - 1.0;
m_Vd = m_Vd + log1p(a) * m_Vt;
#endif
const nl_double eVDVt = exp(m_Vd * m_VtInv);
const nl_double eVDVt = nl_math::exp(m_Vd * m_VtInv);
m_Id = m_Is * (eVDVt - 1.0);
m_G = m_Is * m_VtInv * eVDVt + m_gmin;
@ -277,35 +274,6 @@ private:
// nld_D
// ----------------------------------------------------------------------------------------
// this one has an accuracy of better than 5%. That's enough for our purpose
// add c3 and it'll be better than 1%
#if 0
inline nl_double fastexp_h(const nl_double x)
{
static const nl_double ln2r = 1.442695040888963387;
static const nl_double ln2 = 0.693147180559945286;
//static const nl_double c3 = 0.166666666666666667;
const nl_double y = x * ln2r;
const unsigned int t = y;
const nl_double z = (x - ln2 * (double) t);
const nl_double zz = z * z;
//const nl_double zzz = zz * z;
return (double)(1 << t)*(1.0 + z + 0.5 * zz); // + c3*zzz;
}
inline nl_double fastexp(const nl_double x)
{
if (x<0)
return 1.0 / fastexp_h(-x);
else
return fastexp_h(x);
}
#endif
class NETLIB_NAME(D) : public NETLIB_NAME(twoterm)
{
public:

3
src/emu/netlist/nl_base.h
View File

@ -1305,7 +1305,6 @@ ATTR_HOT inline void netlist_logic_input_t::activate_lh()
ATTR_HOT inline void netlist_net_t::push_to_queue(const netlist_time &delay)
{
//if (UNEXPECTED(m_num_cons == 0 || is_queued()))
if (!is_queued() && (num_cons() > 0))
{
m_time = netlist().time() + delay;
@ -1319,7 +1318,6 @@ ATTR_HOT inline void netlist_net_t::push_to_queue(const netlist_time &delay)
ATTR_HOT inline void netlist_net_t::reschedule_in_queue(const netlist_time &delay)
{
//if (UNEXPECTED(m_num_cons == 0 || is_queued()))
if (is_queued())
netlist().remove_from_queue(*this);
@ -1334,6 +1332,7 @@ ATTR_HOT inline void netlist_net_t::reschedule_in_queue(const netlist_time &dela
ATTR_HOT inline const netlist_sig_t netlist_logic_input_t::Q() const
{
nl_assert(family() == LOGIC);
return net().as_logic().Q();
}

4
src/emu/netlist/nl_dice_compat.h
View File

@ -161,7 +161,7 @@ inline int CAPACITOR_tc_hl(const double c, const double r)
* Vt = (VH-VL)*exp(-t/RC)
* ln(Vt/(VH-VL))*RC = -t
*/
static const double TIME_CONSTANT = -log(2.0 / (3.7-0.3));
static const double TIME_CONSTANT = -nl_math::log(2.0 / (3.7-0.3));
int ret = (int) (TIME_CONSTANT * (130.0 + r) * c * 1e9);
return ret;
}
@ -172,7 +172,7 @@ inline int CAPACITOR_tc_lh(const double c, const double r)
* Vt = (VH-VL)*(1-exp(-t/RC))
* -t=ln(1-Vt/(VH-VL))*RC
*/
static const double TIME_CONSTANT = -log(1.0 - 0.8 / (3.7-0.3));
static const double TIME_CONSTANT = -nl_math::log(1.0 - 0.8 / (3.7-0.3));
int ret = (int) (TIME_CONSTANT * (1.0 + r) * c * 1e9);
return ret;
}

55
src/emu/netlist/nl_util.h
View File

@ -10,6 +10,7 @@
#include "pstring.h"
#include "plists.h"
#include <cmath>
class nl_util
{
@ -40,4 +41,58 @@ public:
}
};
class nl_math
{
// this is purely static
private:
nl_math() {};
public:
ATTR_HOT inline static float exp(const float x) { return std::exp(x); }
ATTR_HOT inline static double abs(const double x) { return std::abs(x); }
ATTR_HOT inline static float abs(const float x) { return std::abs(x); }
ATTR_HOT inline static double log(const double x) { return std::log(x); }
ATTR_HOT inline static float log(const float x) { return std::log(x); }
#if defined(_MSC_VER) && _MSC_VER < 1800
ATTR_HOT inline static double log1p(const double x) { return nl_math::log(1.0 + x); }
ATTR_HOT inline static float log1p(const float x) { return nl_math::log(1.0 + x); }
#else
ATTR_HOT inline static double log1p(const double x) { return log1p(x); }
ATTR_HOT inline static float log1p(const float x) { return log1pf(x); }
#endif
ATTR_HOT inline static double sqrt(const double x) { return std::sqrt(x); }
ATTR_HOT inline static float sqrt(const float x) { return std::sqrt(x); }
// this one has an accuracy of better than 5%. That's enough for our purpose
// add c3 and it'll be better than 1%
#if 0
inline static double fastexp_h(const double x)
{
static const double ln2r = 1.442695040888963387;
static const double ln2 = 0.693147180559945286;
static const double c3 = 0.166666666666666667;
const double y = x * ln2r;
const unsigned int t = y;
const double z = (x - ln2 * (double) t);
const double zz = z * z;
const double zzz = zz * z;
return (double)(1 << t)*(1.0 + z + 0.5 * zz + c3 * zzz);
}
ATTR_HOT inline static double exp(const double x)
{
if (x<0)
return 1.0 / fastexp_h(-x);
else
return fastexp_h(x);
}
#else
ATTR_HOT inline static double exp(const double x) { return std::exp(x); }
#endif
};
#endif /* NL_UTIL_H_ */