math: use double_t for temporaries to avoid stores on i386

When FLT_EVAL_METHOD!=0 (only i386 with x87 fp) the excess
precision of an expression must be removed in an assignment.
(gcc needs -fexcess-precision=standard or -std=c99 for this)

This is done by extra load/store instructions which adds code
bloat when lot of temporaries are used and it makes the result
less precise in many cases.
Using double_t and float_t avoids these issues on i386 and
it makes no difference on other archs.

For now only a few functions are modified where the excess
precision is clearly beneficial (mostly polynomial evaluations
with temporaries).

object size differences on i386, gcc-4.8:
             old   new
__cosdf.o    123    95
__cos.o      199   169
__sindf.o    131    95
__sin.o      225   203
__tandf.o    207   151
__tan.o      605   499
erff.o      1470  1416
erf.o       1703  1649
j0f.o       1779  1745
j0.o        2308  2274
j1f.o       1602  1568
j1.o        2286  2252
tgamma.o    1431  1424
math/*.o   64164 63635

diff --git a/src/math/__cos.c b/src/math/__cos.c
index 8699c1d51933172b1010b7d65b769930414ff9a5..46cefb38134398a4e0e7b40621eb6ea6e514ca38 100644 (file)
--- a/src/math/__cos.c
+++ b/src/math/__cos.c
@@ -60,7 +60,7 @@ C6  = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
 
 double __cos(double x, double y)
 {
-       double hz,z,r,w;
+       double_t hz,z,r,w;
 
        z  = x*x;
        w  = z*z;

diff --git a/src/math/__cosdf.c b/src/math/__cosdf.c
index a65f7f21481643ffb2aaa40ed83b062162a338e2..2124989b3299e53893f2fc9267567e2e7738509b 100644 (file)
--- a/src/math/__cosdf.c
+++ b/src/math/__cosdf.c
@@ -25,7 +25,7 @@ C3  =  0x199342e0ee5069.0p-68; /*  0.0000243904487962774090654 */
 
 float __cosdf(double x)
 {
-       double r, w, z;
+       double_t r, w, z;
 
        /* Try to optimize for parallel evaluation as in __tandf.c. */
        z = x*x;

diff --git a/src/math/__log1p.h b/src/math/__log1p.h
index ec2c77b9318e2113297ff0cabb7df732f27b6efd..57187115c3a4bc1406b8d95d0b0a2635ffa1ac44 100644 (file)
--- a/src/math/__log1p.h
+++ b/src/math/__log1p.h
@@ -81,7 +81,7 @@ Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
  */
 static inline double __log1p(double f)
 {
-       double hfsq,s,z,R,w,t1,t2;
+       double_t hfsq,s,z,R,w,t1,t2;
 
        s = f/(2.0+f);
        z = s*s;

diff --git a/src/math/__log1pf.h b/src/math/__log1pf.h
index 99492c5a1ac9b892628d79b61f38df37eb1a85de..f2fbef291fb443c52c41332880128be63278eac6 100644 (file)
--- a/src/math/__log1pf.h
+++ b/src/math/__log1pf.h
@@ -22,7 +22,7 @@ Lg4 = 0xf89e26.0p-26; /* 0.24279078841 */
 
 static inline float __log1pf(float f)
 {
-       float hfsq,s,z,R,w,t1,t2;
+       float_t hfsq,s,z,R,w,t1,t2;
 
        s = f/(2.0f + f);
        z = s*s;

diff --git a/src/math/__sin.c b/src/math/__sin.c
index 9aead04becfd4e68177cd18cf9d4b50ab2208533..40309496646ea67bc247a46b8cd363a73ca9801e 100644 (file)
--- a/src/math/__sin.c
+++ b/src/math/__sin.c
@@ -51,7 +51,7 @@ S6  =  1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
 
 double __sin(double x, double y, int iy)
 {
-       double z,r,v,w;
+       double_t z,r,v,w;
 
        z = x*x;
        w = z*z;

diff --git a/src/math/__sindf.c b/src/math/__sindf.c
index 83c0d7a5f8c277513fd09d6d64a7ea5b9a74079c..8fec2a3f660c627b9dd575f96601c9530dc9a67e 100644 (file)
--- a/src/math/__sindf.c
+++ b/src/math/__sindf.c
@@ -25,7 +25,7 @@ S4 =  0x16cd878c3b46a7.0p-71; /*  0.0000027183114939898219064 */
 
 float __sindf(double x)
 {
-       double r, s, w, z;
+       double_t r, s, w, z;
 
        /* Try to optimize for parallel evaluation as in __tandf.c. */
        z = x*x;

diff --git a/src/math/__tan.c b/src/math/__tan.c
index 01e3fe48baccec982e262ea65b1fb9c8162fce38..fc739f95fd3a26c91799809633f0d2ec459305cb 100644 (file)
--- a/src/math/__tan.c
+++ b/src/math/__tan.c
@@ -65,7 +65,7 @@ pio4lo =     3.06161699786838301793e-17; /* 3C81A626, 33145C07 */
 
 double __tan(double x, double y, int iy)
 {
-       double z, r, v, w, s, sign;
+       double_t z, r, v, w, s, sign;
        int32_t ix, hx;
 
        GET_HIGH_WORD(hx,x);
@@ -106,7 +106,8 @@ double __tan(double x, double y, int iy)
                 * -1.0 / (x+r) here
                 */
                /* compute -1.0 / (x+r) accurately */
-               double a, t;
+               double_t a;
+               double z, t;
                z = w;
                SET_LOW_WORD(z,0);
                v = r - (z - x);        /* z+v = r+x */

diff --git a/src/math/__tandf.c b/src/math/__tandf.c
index 36a8214ed36a3f2ab0ee1bdd97e3a3e1f83540e1..3e632fdf4abce4558b7850b288d9f45c52833449 100644 (file)
--- a/src/math/__tandf.c
+++ b/src/math/__tandf.c
@@ -27,7 +27,7 @@ static const double T[] = {
 
 float __tandf(double x, int iy)
 {
-       double z,r,w,s,t,u;
+       double_t z,r,w,s,t,u;
 
        z = x*x;
        /*

diff --git a/src/math/acos.c b/src/math/acos.c
index cd5d06a632921a668eaa7e7bf63ab449de423478..ea9c87bf087809f9a6b1fa6a2ed2f065d8184f72 100644 (file)
--- a/src/math/acos.c
+++ b/src/math/acos.c
@@ -51,7 +51,7 @@ qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
 
 static double R(double z)
 {
-       double p, q;
+       double_t p, q;
        p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
        q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
        return p/q;

diff --git a/src/math/acosf.c b/src/math/acosf.c
index 5d7c027089897272d3df3c49fcc4f42b7ea11f72..8ee1a71d0cc2965e08e170be1dc4529d0165f3b4 100644 (file)
--- a/src/math/acosf.c
+++ b/src/math/acosf.c
@@ -25,7 +25,7 @@ qS1 = -7.0662963390e-01;
 
 static float R(float z)
 {
-       float p, q;
+       float_t p, q;
        p = z*(pS0+z*(pS1+z*pS2));
        q = 1.0f+z*qS1;
        return p/q;

diff --git a/src/math/asin.c b/src/math/asin.c
index d61c04b4b86ef0452a76238663991747024335e4..3e8f99ed558fb8df591c662648e0ab38335fd182 100644 (file)
--- a/src/math/asin.c
+++ b/src/math/asin.c
@@ -58,7 +58,7 @@ qS4 =  7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
 
 static double R(double z)
 {
-       double p, q;
+       double_t p, q;
        p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
        q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
        return p/q;

diff --git a/src/math/asinf.c b/src/math/asinf.c
index 462bf0431e7e610d39d4f456b631120f86c2d4af..51fe6c617cc03569c759cd56acb1a0ecc2e55b49 100644 (file)
--- a/src/math/asinf.c
+++ b/src/math/asinf.c
@@ -26,7 +26,7 @@ qS1 = -7.0662963390e-01;
 
 static float R(float z)
 {
-       float p, q;
+       float_t p, q;
        p = z*(pS0+z*(pS1+z*pS2));
        q = 1.0f+z*qS1;
        return p/q;

diff --git a/src/math/atan.c b/src/math/atan.c
index 3c9a59ff3067abdeca4a1c6270f89dc9882e7e9d..5a1d33e6a05c644610933b1bea8ea8e39f3f99f1 100644 (file)
--- a/src/math/atan.c
+++ b/src/math/atan.c
@@ -62,7 +62,7 @@ static const double aT[] = {
 
 double atan(double x)
 {
-       double w,s1,s2,z;
+       double_t w,s1,s2,z;
        uint32_t ix,sign;
        int id;
 

diff --git a/src/math/atanf.c b/src/math/atanf.c
index 4b59509ad4a47aafc60127129c8edf1380d641ff..ac8bfd0665f08547e444939c540dd6eafe461fd0 100644 (file)
--- a/src/math/atanf.c
+++ b/src/math/atanf.c
@@ -40,7 +40,7 @@ static const float aT[] = {
 
 float atanf(float x)
 {
-       float w,s1,s2,z;
+       float_t w,s1,s2,z;
        uint32_t ix,sign;
        int id;
 

diff --git a/src/math/erf.c b/src/math/erf.c
index c0fc41db3c95740edb493f4617834c805219dcf5..2f30a298f997cde0316af8681bb1447f085de54d 100644 (file)
--- a/src/math/erf.c
+++ b/src/math/erf.c
@@ -176,7 +176,7 @@ sb7  = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
 
 static double erfc1(double x)
 {
-       double s,P,Q;
+       double_t s,P,Q;
 
        s = fabs(x) - 1;
        P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
@@ -186,7 +186,8 @@ static double erfc1(double x)
 
 static double erfc2(uint32_t ix, double x)
 {
-       double s,z,R,S;
+       double_t s,R,S;
+       double z;
 
        if (ix < 0x3ff40000)  /* |x| < 1.25 */
                return erfc1(x);

diff --git a/src/math/erff.c b/src/math/erff.c
index e2cfc984306d4b52e84b1709977e74f33c617f17..ed5f3975741595057cbc981d13fd6c43ce11270a 100644 (file)
--- a/src/math/erff.c
+++ b/src/math/erff.c
@@ -86,7 +86,7 @@ sb7  = -2.2440952301e+01; /* 0xc1b38712 */
 
 static float erfc1(float x)
 {
-       float s,P,Q;
+       float_t s,P,Q;
 
        s = fabsf(x) - 1;
        P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
@@ -96,7 +96,8 @@ static float erfc1(float x)
 
 static float erfc2(uint32_t ix, float x)
 {
-       float s,z,R,S;
+       float_t s,R,S;
+       float z;
 
        if (ix < 0x3fa00000)  /* |x| < 1.25 */
                return erfc1(x);

diff --git a/src/math/j0.c b/src/math/j0.c
index b281e136cd9307cccc083262fd3d227eb62326c1..d722d942804524e65c08f0bf1cea8d5c8c738d77 100644 (file)
--- a/src/math/j0.c
+++ b/src/math/j0.c
@@ -263,7 +263,7 @@ static const double pS2[5] = {
 static double pzero(double x)
 {
        const double *p,*q;
-       double z,r,s;
+       double_t z,r,s;
        uint32_t ix;
 
        GET_HIGH_WORD(ix, x);
@@ -359,7 +359,7 @@ static const double qS2[6] = {
 static double qzero(double x)
 {
        const double *p,*q;
-       double s,r,z;
+       double_t s,r,z;
        uint32_t ix;
 
        GET_HIGH_WORD(ix, x);

diff --git a/src/math/j0f.c b/src/math/j0f.c
index 79bab62a60fa474ae7aeb144ac15f0cc59fbd82a..4b0ee3b712e818d95240d74e58b66099ee643fd6 100644 (file)
--- a/src/math/j0f.c
+++ b/src/math/j0f.c
@@ -201,7 +201,7 @@ static const float pS2[5] = {
 static float pzerof(float x)
 {
        const float *p,*q;
-       float z,r,s;
+       float_t z,r,s;
        uint32_t ix;
 
        GET_FLOAT_WORD(ix, x);
@@ -297,7 +297,7 @@ static const float qS2[6] = {
 static float qzerof(float x)
 {
        const float *p,*q;
-       float s,r,z;
+       float_t s,r,z;
        uint32_t ix;
 
        GET_FLOAT_WORD(ix, x);

diff --git a/src/math/j1.c b/src/math/j1.c
index ac7bb1ebefc5261b15201c7bc00ae6cac160d0a0..df724d172ea0064ea919d5eeee4627243f7df1d5 100644 (file)
--- a/src/math/j1.c
+++ b/src/math/j1.c
@@ -250,7 +250,7 @@ static const double ps2[5] = {
 static double pone(double x)
 {
        const double *p,*q;
-       double z,r,s;
+       double_t z,r,s;
        uint32_t ix;
 
        GET_HIGH_WORD(ix, x);
@@ -346,7 +346,7 @@ static const double qs2[6] = {
 static double qone(double x)
 {
        const double *p,*q;
-       double  s,r,z;
+       double_t s,r,z;
        uint32_t ix;
 
        GET_HIGH_WORD(ix, x);

diff --git a/src/math/j1f.c b/src/math/j1f.c
index 5a760f712181b06c6ea41b96c40b62afceccfb4e..6abde34927da2998ac0dfa94a651182612b98582 100644 (file)
--- a/src/math/j1f.c
+++ b/src/math/j1f.c
@@ -198,7 +198,7 @@ static const float ps2[5] = {
 static float ponef(float x)
 {
        const float *p,*q;
-       float z,r,s;
+       float_t z,r,s;
        uint32_t ix;
 
        GET_FLOAT_WORD(ix, x);
@@ -294,7 +294,7 @@ static const float qs2[6] = {
 static float qonef(float x)
 {
        const float *p,*q;
-       float s,r,z;
+       float_t s,r,z;
        uint32_t ix;
 
        GET_FLOAT_WORD(ix, x);

diff --git a/src/math/tgamma.c b/src/math/tgamma.c
index a3f203c174468f89c04061cf53cfbaf2efbcdb69..691e86a4bbf504373780e8a4ca24fa523cbfd676 100644 (file)
--- a/src/math/tgamma.c
+++ b/src/math/tgamma.c
@@ -89,7 +89,7 @@ static const double fact[] = {
 /* S(x) rational function for positive x */
 static double S(double x)
 {
-       double num = 0, den = 0;
+       double_t num = 0, den = 0;
        int i;
 
        /* to avoid overflow handle large x differently */