From 7d7c2c1e432d3e8aed926bd65c8fac994c2696f0 Mon Sep 17 00:00:00 2001 From: cshapiro <cshapiro> Date: Fri, 11 Apr 2008 07:49:46 +0000 Subject: [PATCH] Retrieve the constituent words of a double through a union instead of casts and pointer arithmetic. This code now compiles correctly on GCC version 4 without disabling the default assumption of strict aliasing rules. Also, remove lots of dead code associated with parts of fdlibm that have not been imported into the runtime. --- lisp/e_rem_pio2.c | 35 +++------ lisp/fdlibm.h | 189 ++-------------------------------------------- lisp/k_rem_pio2.c | 19 +---- 3 files changed, 21 insertions(+), 222 deletions(-) diff --git a/lisp/e_rem_pio2.c b/lisp/e_rem_pio2.c index 7242bb232..d15af174b 100644 --- a/lisp/e_rem_pio2.c +++ b/lisp/e_rem_pio2.c @@ -23,11 +23,7 @@ /* * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi */ -#ifdef __STDC__ static const int two_over_pi[] = { -#else -static int two_over_pi[] = { -#endif 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, @@ -41,11 +37,7 @@ static int two_over_pi[] = { 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, }; -#ifdef __STDC__ static const int npio2_hw[] = { -#else -static int npio2_hw[] = { -#endif 0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, @@ -64,11 +56,7 @@ static int npio2_hw[] = { * pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3) */ -#ifdef __STDC__ static const double -#else -static double -#endif zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ @@ -80,18 +68,15 @@ pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */ pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ -#ifdef __STDC__ - int __ieee754_rem_pio2(double x, double *y) -#else - int __ieee754_rem_pio2(x,y) - double x,y[]; -#endif +int __ieee754_rem_pio2(double x, double *y) { double z,w,t,r,fn; double tx[3]; int e0,i,j,nx,n,ix,hx; + union { int i[2]; double d; } ux,uy0,uz; - hx = __HI(x); /* high word of x */ + ux.d = x; + hx = ux.i[HIWORD]; /* high word of x */ ix = hx&0x7fffffff; if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ {y[0] = x; y[1] = 0; return 0;} @@ -131,14 +116,16 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ } else { j = ix>>20; y[0] = r-w; - i = j-(((__HI(y[0]))>>20)&0x7ff); + uy0.d = y[0]; + i = j-((uy0.i[HIWORD]>>20)&0x7ff); if(i>16) { /* 2nd iteration needed, good to 118 */ t = r; w = fn*pio2_2; r = t-w; w = fn*pio2_2t-((t-r)-w); y[0] = r-w; - i = j-(((__HI(y[0]))>>20)&0x7ff); + uy0.d = y[0]; + i = j-((uy0.i[HIWORD]>>20)&0x7ff); if(i>49) { /* 3rd iteration need, 151 bits acc */ t = r; /* will cover all possible cases */ w = fn*pio2_3; @@ -159,9 +146,11 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ y[0]=y[1]=x-x; return 0; } /* set z = scalbn(|x|,ilogb(x)-23) */ - __LO(z) = __LO(x); + ux.d = x; + uz.i[LOWORD] = ux.i[LOWORD]; e0 = (ix>>20)-1046; /* e0 = ilogb(z)-23; */ - __HI(z) = ix - (e0<<20); + uz.i[HIWORD] = ix - (e0<<20); + z = uz.d; for(i=0;i<2;i++) { tx[i] = (double)((int)(z)); z = (z-tx[i])*two24; diff --git a/lisp/fdlibm.h b/lisp/fdlibm.h index 02ca7fb2d..3a9fcdb93 100644 --- a/lisp/fdlibm.h +++ b/lisp/fdlibm.h @@ -20,197 +20,24 @@ #endif #ifdef __LITTLE_ENDIAN -#define __HI(x) *(1+(int*)&x) -#define __LO(x) *(int*)&x -#define __HIp(x) *(1+(int*)x) -#define __LOp(x) *(int*)x +enum { HIWORD = 1, LOWORD = 0 }; #else -#define __HI(x) *(int*)&x -#define __LO(x) *(1+(int*)&x) -#define __HIp(x) *(int*)x -#define __LOp(x) *(1+(int*)x) -#endif - -#ifdef __STDC__ -#define __P(p) p -#else -#define __P(p) () +enum { HIWORD = 0, LOWORD = 1 }; #endif /* * ANSI/POSIX */ +extern double fabs(double); +extern double floor(double); -extern int signgam; - -#define MAXFLOAT ((float)3.40282346638528860e+38) - -enum fdversion {fdlibm_ieee = -1, fdlibm_svid, fdlibm_xopen, fdlibm_posix}; - -#define _LIB_VERSION_TYPE enum fdversion -#define _LIB_VERSION _fdlib_version - -/* if global variable _LIB_VERSION is not desirable, one may - * change the following to be a constant by: - * #define _LIB_VERSION_TYPE const enum version - * In that case, after one initializes the value _LIB_VERSION (see - * s_lib_version.c) during compile time, it cannot be modified - * in the middle of a program - */ -extern _LIB_VERSION_TYPE _LIB_VERSION; - -#define _IEEE_ fdlibm_ieee -#define _SVID_ fdlibm_svid -#define _XOPEN_ fdlibm_xopen -#define _POSIX_ fdlibm_posix - -struct exception { - int type; - char *name; - double arg1; - double arg2; - double retval; -}; - -#define HUGE MAXFLOAT - -/* - * set X_TLOSS = pi*2**52, which is possibly defined in <values.h> - * (one may replace the following line by "#include <values.h>") - */ - -#define X_TLOSS 1.41484755040568800000e+16 - -#define DOMAIN 1 -#define SING 2 -#define OVERFLOW 3 -#define UNDERFLOW 4 -#define TLOSS 5 -#define PLOSS 6 - -/* - * ANSI/POSIX - */ -extern double acos __P((double)); -extern double asin __P((double)); -extern double atan __P((double)); -extern double atan2 __P((double, double)); -extern double cos __P((double)); -extern double sin __P((double)); -extern double tan __P((double)); - -extern double cosh __P((double)); -extern double sinh __P((double)); -extern double tanh __P((double)); - -extern double exp __P((double)); -extern double frexp __P((double, int *)); -extern double ldexp __P((double, int)); -extern double log __P((double)); -extern double log10 __P((double)); -extern double modf __P((double, double *)); - -extern double pow __P((double, double)); -extern double sqrt __P((double)); - -extern double ceil __P((double)); -extern double fabs __P((double)); -extern double floor __P((double)); -extern double fmod __P((double, double)); - -extern double erf __P((double)); -extern double erfc __P((double)); -extern double gamma __P((double)); -extern double hypot __P((double, double)); -extern int isnan __P((double)); -extern int finite __P((double)); -extern double j0 __P((double)); -extern double j1 __P((double)); -extern double jn __P((int, double)); -extern double lgamma __P((double)); -extern double y0 __P((double)); -extern double y1 __P((double)); -extern double yn __P((int, double)); - -extern double acosh __P((double)); -extern double asinh __P((double)); -extern double atanh __P((double)); -extern double cbrt __P((double)); -extern double logb __P((double)); -extern double nextafter __P((double, double)); -extern double remainder __P((double, double)); -#ifdef _SCALB_INT -extern double scalb __P((double, int)); -#else -extern double scalb __P((double, double)); -#endif - -extern int matherr __P((struct exception *)); - -/* - * IEEE Test Vector - */ -extern double significand __P((double)); +/* ieee style elementary functions */ +extern int __ieee754_rem_pio2(double,double*); /* * Functions callable from C, intended to support IEEE arithmetic. */ -extern double copysign __P((double, double)); -extern int ilogb __P((double)); -extern double rint __P((double)); -extern double scalbn __P((double, int)); - -/* - * BSD math library entry points - */ -extern double expm1 __P((double)); -extern double log1p __P((double)); - -/* - * Reentrant version of gamma & lgamma; passes signgam back by reference - * as the second argument; user must allocate space for signgam. - */ -#ifdef _REENTRANT -extern double gamma_r __P((double, int *)); -extern double lgamma_r __P((double, int *)); -#endif /* _REENTRANT */ - -/* ieee style elementary functions */ -extern double __ieee754_sqrt __P((double)); -extern double __ieee754_acos __P((double)); -extern double __ieee754_acosh __P((double)); -extern double __ieee754_log __P((double)); -extern double __ieee754_atanh __P((double)); -extern double __ieee754_asin __P((double)); -extern double __ieee754_atan2 __P((double,double)); -extern double __ieee754_exp __P((double)); -extern double __ieee754_cosh __P((double)); -extern double __ieee754_fmod __P((double,double)); -extern double __ieee754_pow __P((double,double)); -extern double __ieee754_lgamma_r __P((double,int *)); -extern double __ieee754_gamma_r __P((double,int *)); -extern double __ieee754_lgamma __P((double)); -extern double __ieee754_gamma __P((double)); -extern double __ieee754_log10 __P((double)); -extern double __ieee754_sinh __P((double)); -extern double __ieee754_hypot __P((double,double)); -extern double __ieee754_j0 __P((double)); -extern double __ieee754_j1 __P((double)); -extern double __ieee754_y0 __P((double)); -extern double __ieee754_y1 __P((double)); -extern double __ieee754_jn __P((int,double)); -extern double __ieee754_yn __P((int,double)); -extern double __ieee754_remainder __P((double,double)); -extern int __ieee754_rem_pio2 __P((double,double*)); -#ifdef _SCALB_INT -extern double __ieee754_scalb __P((double,int)); -#else -extern double __ieee754_scalb __P((double,double)); -#endif +extern double scalbn(double, int); /* fdlibm kernel function */ -extern double __kernel_standard __P((double,double,int)); -extern double __kernel_sin __P((double,double,int)); -extern double __kernel_cos __P((double,double)); -extern double __kernel_tan __P((double,double,int)); -extern int __kernel_rem_pio2 __P((double*,double*,int,int,int,const int*)); +extern int __kernel_rem_pio2(double*,double*,int,int,int,const int*); diff --git a/lisp/k_rem_pio2.c b/lisp/k_rem_pio2.c index ec473ac0d..1d1a82524 100644 --- a/lisp/k_rem_pio2.c +++ b/lisp/k_rem_pio2.c @@ -129,17 +129,9 @@ #include "fdlibm.h" -#ifdef __STDC__ static const int init_jk[] = {2,3,4,6}; /* initial value for jk */ -#else -static int init_jk[] = {2,3,4,6}; -#endif -#ifdef __STDC__ static const double PIo2[] = { -#else -static double PIo2[] = { -#endif 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ @@ -150,22 +142,13 @@ static double PIo2[] = { 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */ }; -#ifdef __STDC__ static const double -#else -static double -#endif zero = 0.0, one = 1.0, two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ -#ifdef __STDC__ - int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int *ipio2) -#else - int __kernel_rem_pio2(x,y,e0,nx,prec,ipio2) - double x[], y[]; int e0,nx,prec; int ipio2[]; -#endif +int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int *ipio2) { int jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih; double z,fw,f[20],fq[20],q[20]; -- GitLab