ntpd: try to avoid using libm. -1.2k if we succeed
uclibc's sqrt(x) is pathetic, 411 bytes? it can be ~100... Signed-off-by: Denys Vlasenko <vda.linux@googlemail.com>
This commit is contained in:
parent
510f56aa6f
commit
d498ff0ac4
@ -299,7 +299,45 @@ static ALWAYS_INLINE double MIND(double a, double b)
|
||||
return a;
|
||||
return b;
|
||||
}
|
||||
#define SQRT(x) (sqrt(x))
|
||||
static NOINLINE double my_SQRT(double X)
|
||||
{
|
||||
union {
|
||||
float f;
|
||||
int32_t i;
|
||||
} v;
|
||||
double invsqrt;
|
||||
double Xhalf = X * 0.5;
|
||||
|
||||
/* Fast and good approximation to 1/sqrt(X), black magic */
|
||||
v.f = X;
|
||||
/*v.i = 0x5f3759df - (v.i >> 1);*/
|
||||
v.i = 0x5f375a86 - (v.i >> 1); /* - this constant is slightly better */
|
||||
invsqrt = v.f; /* better than 0.2% accuracy */
|
||||
|
||||
/* Refining it using Newton's method: x1 = x0 - f(x0)/f'(x0)
|
||||
* f(x) = 1/(x*x) - X (f==0 when x = 1/sqrt(X))
|
||||
* f'(x) = -2/(x*x*x)
|
||||
* f(x)/f'(x) = (X - 1/(x*x)) / (2/(x*x*x)) = X*x*x*x/2 - x/2
|
||||
* x1 = x0 - (X*x0*x0*x0/2 - x0/2) = 1.5*x0 - X*x0*x0*x0/2 = x0*(1.5 - (X/2)*x0*x0)
|
||||
*/
|
||||
invsqrt = invsqrt * (1.5 - Xhalf * invsqrt * invsqrt); /* ~0.05% accuracy */
|
||||
/* invsqrt = invsqrt * (1.5 - Xhalf * invsqrt * invsqrt); 2nd iter: ~0.0001% accuracy */
|
||||
/* With 4 iterations, more than half results will be exact,
|
||||
* at 6th iterations result stabilizes with about 72% results exact.
|
||||
* We are well satisfied with 0.05% accuracy.
|
||||
*/
|
||||
|
||||
return X * invsqrt; /* X * 1/sqrt(X) ~= sqrt(X) */
|
||||
}
|
||||
static ALWAYS_INLINE double SQRT(double X)
|
||||
{
|
||||
/* If this arch doesn't use IEEE 754 floats, fall back to using libm */
|
||||
if (sizeof(float) != 4)
|
||||
return sqrt(X);
|
||||
|
||||
/* This avoids needing libm, saves about 1.2k on x86-32 */
|
||||
return my_SQRT(X);
|
||||
}
|
||||
|
||||
static double
|
||||
gettime1900d(void)
|
||||
|
Loading…
Reference in New Issue
Block a user