Вот код exp из студии (кол-во операций считай сам)
| Код |
/* coefficients */
static const double p[] = { /* courtesy Dr. Tim Prince */
1.0,
420.30235984910635,
15132.70094680474802};
static const double q[] = { /* courtesy Dr. Tim Prince */
30.01511290683317,
3362.72154416553028,
30265.40189360949691};
static const double c1 = 22713.0 / 32768.0;
static const double c2 = 1.4286068203094172321214581765680755e-6;
static const double hugexp = HUGE_EXP;
static const double invln2 = 1.4426950408889634073599246810018921;
_CRTIMP2_PURE short __CLRCALL_PURE_OR_CDECL _Exp(double *px, double y, short eoff)
{ /* compute y*e^(*px), (*px) finite, |y| not huge */
if (*px < -hugexp || y == 0.0)
{ /* certain underflow */
*px = 0.0;
return (0);
}
else if (hugexp < *px)
{ /* certain overflow */
*px = _Inf._Double;
return (_INFCODE);
}
else
{ /* xexp won't overflow */
double g = *px * invln2;
short xexp = (short)(g + (g < 0.0 ? - 0.5 : + 0.5));
g = xexp;
g = (*px - g * c1) - g * c2;
if (-_Eps._Double < g && g < _Eps._Double)
*px = y;
else
{ /* g*g worth computing */
const double z = g * g;
const double w = (q[0] * z + q[1]) * z + q[2];
g *= (z + p[1]) * z + p[2];
*px = (w + g) / (w - g) * 2.0 * y;
--xexp;
}
return (_Dscale(px, (long)xexp + eoff));
}
}
|
а вот ф-ция _Dscale
| Код |
_CRTIMP2_PURE short __CLRCALL_PURE_OR_CDECL _Dscale(double *px, long lexp)
{ /* scale *px by 2^xexp with checking */
unsigned short *ps = (unsigned short *)(char *)px;
short xchar = (short)((ps[_D0] & _DMASK) >> _DOFF);
if (xchar == _DMAX)
return ((short)((ps[_D0] & _DFRAC) != 0 || ps[_D1] != 0
|| ps[_D2] != 0 || ps[_D3] != 0 ? _NANCODE : _INFCODE));
else if (xchar == 0 && 0 < (xchar = _Dnorm(ps)))
return (0);
if (0 < lexp && _DMAX - xchar <= lexp)
{ /* overflow, return +/-INF */
*px = ps[_D0] & _DSIGN ? -_Inf._Double : _Inf._Double;
return (_INFCODE);
}
else if (-xchar < lexp)
{ /* finite result, repack */
ps[_D0] = (unsigned short)(ps[_D0] & ~_DMASK
| (lexp + xchar) << _DOFF);
return (_FINITE);
}
else
{ /* denormalized, scale */
unsigned short sign = (unsigned short)(ps[_D0] & _DSIGN);
ps[_D0] = (unsigned short)(1 << _DOFF | ps[_D0] & _DFRAC);
lexp += xchar - 1;
if (lexp < -(48 + 1 + _DOFF) || 0 <= lexp)
{ /* certain underflow, return +/-0 */
ps[_D0] = sign;
ps[_D1] = 0;
ps[_D2] = 0;
ps[_D3] = 0;
return (0);
}
else
{ /* nonzero, align fraction */
short xexp = (short)lexp;
unsigned short psx = 0;
for (; xexp <= -16; xexp += 16)
{ /* scale by words */
psx = ps[_D3] | (psx != 0 ? 1 : 0);
ps[_D3] = ps[_D2];
ps[_D2] = ps[_D1];
ps[_D1] = ps[_D0];
ps[_D0] = 0;
}
if ((xexp = (short)-xexp) != 0)
{ /* scale by bits */
psx = (ps[_D3] << (16 - xexp)) | (psx != 0 ? 1 : 0);
ps[_D3] = (unsigned short)(ps[_D3] >> xexp
| ps[_D2] << (16 - xexp));
ps[_D2] = (unsigned short)(ps[_D2] >> xexp
| ps[_D1] << (16 - xexp));
ps[_D1] = (unsigned short)(ps[_D1] >> xexp
| ps[_D0] << (16 - xexp));
ps[_D0] >>= xexp;
}
ps[_D0] |= sign;
if ((0x8000 < psx || 0x8000 == psx && (ps[_D3] & 0x0001) != 0)
&& (++ps[_D3] & 0xffff) == 0
&& (++ps[_D2] & 0xffff) == 0
&& (++ps[_D1] & 0xffff) == 0)
++ps[_D0]; /* round up */
else if (ps[_D0] == sign && ps[_D1] == 0
&& ps[_D2] == 0 && ps[_D3] == 0)
return (0);
return (_FINITE);
}
}
}
|
|
)