From bb2376a6ed6f6c2c587d45cb2dbec4f47afe10fe Mon Sep 17 00:00:00 2001 From: Olaf Barthel Date: Sun, 9 Oct 2005 10:38:56 +0000 Subject: [PATCH] - Implemented lrintf(), lrint(), lroundf(), lround(), nearbyintf(), nearbyint(), remquof(), remquo(), roundf(), round(), tgammaf(), tgamma(), truncf(), trunc(). Sort of implemented fmaf() and fma(), which really ought to be done in "SIMD" fashion. This completes the "real" floating point math library (ignoring for a moment that the floating point environment code is still not implemented). git-svn-id: file:///Users/olsen/Code/migration-svn-zu-git/logical-line-staging/clib2/trunk@15035 87f5fb63-7c3d-0410-a384-fd976d0f7a62 --- library/TODO | 20 ++-------- library/changes | 9 +++++ library/include/math.h | 18 ++++++++- library/math_fma.c | 7 ++-- library/math_fmaf.c | 7 ++-- library/math_headers.h | 4 +- library/math_lgamma.c | 21 ++++++++++- library/math_lgammaf.c | 21 ++++++++++- library/math_lrint.c | 79 +++++++++++++++++++++++++++++++++++++-- library/math_lrintf.c | 64 +++++++++++++++++++++++++++++-- library/math_lround.c | 59 +++++++++++++++++++++++++++-- library/math_lroundf.c | 40 ++++++++++++++++++-- library/math_nearbyint.c | 6 +-- library/math_nearbyintf.c | 6 +-- library/math_remquo.c | 28 ++++++++++++-- library/math_remquof.c | 28 ++++++++++++-- library/math_round.c | 71 +++++++++++++++++++++++++++++++++-- library/math_roundf.c | 56 +++++++++++++++++++++++++-- library/math_tgamma.c | 21 +++++++++-- library/math_tgammaf.c | 21 +++++++++-- library/math_trunc.c | 57 ++++++++++++++++++++++++++-- library/math_truncf.c | 44 ++++++++++++++++++++-- 22 files changed, 616 insertions(+), 71 deletions(-) diff --git a/library/TODO b/library/TODO index 9650cc0..3cd823a 100755 --- a/library/TODO +++ b/library/TODO @@ -1,19 +1,5 @@ C99 math functions: - (functions generally missing, including their "float" counterparts) - fma - fmaf - lrint - lrintf - lround - lroundf - nearbyint - nearbyintf - remquo - remquof - round - roundf - tgamma - tgammaf - trunc - truncf + fma() and fmaf() should be implemented as a true "fused" multiply + and add function rather than the sequential operation implied in the + current implementation. diff --git a/library/changes b/library/changes index 00574fc..5bba762 100644 --- a/library/changes +++ b/library/changes @@ -1,3 +1,12 @@ +- Implemented lrintf(), lrint(), lroundf(), lround(), nearbyintf(), + nearbyint(), remquof(), remquo(), roundf(), round(), tgammaf(), + tgamma(), truncf(), trunc(). Sort of implemented fmaf() and fma(), + which really ought to be done in "SIMD" fashion. + + This completes the "real" floating point math library (ignoring + for a moment that the floating point environment code is still + not implemented). + - accept() now calls the bsdsocket.library accept() function first and then hooks up the socket with the clib2 data structures. This makes it possible to have several Processes calling the accept() diff --git a/library/include/math.h b/library/include/math.h index 6dff2b5..6cb9a50 100644 --- a/library/include/math.h +++ b/library/include/math.h @@ -1,5 +1,5 @@ /* - * $Id: math.h,v 1.16 2005-06-26 09:06:12 obarthel Exp $ + * $Id: math.h,v 1.17 2005-10-09 10:38:56 obarthel Exp $ * * :ts=4 * @@ -250,17 +250,25 @@ extern float erfcf(float x); extern float erff(float x); extern float expm1f(float x); extern float fdimf(float x,float y); +extern float fmaf(float x,float y,float z); extern float fmaxf(float x,float y); extern float fminf(float x,float y); extern float hypotf(float x, float y); extern float lgammaf(float x); extern float log1pf(float x); extern float logbf(float x); +extern long int lrintf(float x); +extern long int lroundf(float x); extern float nanf(const char *tagp); +extern float nearbyintf(float x); extern float nextafterf(float x,float y); extern float remainderf(float x, float p); +extern float remquof(float x,float y,int * quo); extern float rintf(float x); +extern float roundf(float x); extern float scalbnf (float x, int n); +extern float tgammaf(float x); +extern float truncf(float x); extern int ilogbf(float x); /****************************************************************************/ @@ -273,17 +281,25 @@ extern double erf(double x); extern double erfc(double x); extern double expm1(double x); extern double fdim(double x,double y); +extern double fma(double x,double y,double z); extern double fmax(double x,double y); extern double fmin(double x,double y); extern double hypot(double x,double y); extern double lgamma(double x); extern double log1p(double x); extern double logb(double x); +extern long int lrint(double x); +extern long int lround(double x); extern double nan(const char *tagp); +extern double nearbyint(double x); extern double nextafter(double x,double y); extern double remainder(double x, double p); +extern double remquo(double x,double y,int * quo); extern double rint(double x); +extern double round(double x); extern double scalbn (double x, int n); +extern double tgamma(double x); +extern double trunc(double x); extern int ilogb(double x); /****************************************************************************/ diff --git a/library/math_fma.c b/library/math_fma.c index 14a4ab0..86af270 100644 --- a/library/math_fma.c +++ b/library/math_fma.c @@ -1,5 +1,5 @@ /* - * $Id: math_fma.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_fma.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -44,8 +44,9 @@ double fma(double x,double y,double z) { - /* ZZZ unimplemented */ - return(0); + /* ZZZ this should be a *fused* multiply & add, and + not a sequential operation as declared below! */ + return((x * y) + z); } /****************************************************************************/ diff --git a/library/math_fmaf.c b/library/math_fmaf.c index ce9660a..97d976c 100644 --- a/library/math_fmaf.c +++ b/library/math_fmaf.c @@ -1,5 +1,5 @@ /* - * $Id: math_fmaf.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_fmaf.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -44,8 +44,9 @@ float fmaf(float x,float y,float z) { - /* ZZZ unimplemented */ - return(0); + /* ZZZ this should be a *fused* multiply & add, and + not a sequential operation as declared below! */ + return((x * y) + z); } /****************************************************************************/ diff --git a/library/math_headers.h b/library/math_headers.h index a81ca7f..e935061 100644 --- a/library/math_headers.h +++ b/library/math_headers.h @@ -1,5 +1,5 @@ /* - * $Id: math_headers.h,v 1.11 2005-05-30 08:10:37 obarthel Exp $ + * $Id: math_headers.h,v 1.12 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -220,6 +220,8 @@ extern float __kernel_cosf(float x, float y); extern float __kernel_sinf(float x, float y, int iy); extern LONG __rem_pio2f(float x, float *y); extern float __kernel_tanf(float x, float y, int iy); +extern double __lgamma(double x,int * gamma_sign_ptr); +extern float __lgammaf(float x,int * gamma_sign_ptr); /****************************************************************************/ diff --git a/library/math_lgamma.c b/library/math_lgamma.c index 4bcd182..431be03 100644 --- a/library/math_lgamma.c +++ b/library/math_lgamma.c @@ -1,5 +1,5 @@ /* - * $Id: math_lgamma.c,v 1.2 2005-05-30 08:47:26 obarthel Exp $ + * $Id: math_lgamma.c,v 1.3 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -157,7 +157,7 @@ sin_pi(double x) } double -lgamma(double x) +__lgamma(double x,int * gamma_sign_ptr) { double t,y,z,nadj=0.0,p,p1,p2,p3,q,r,w; LONG i,hx,lx,ix; @@ -165,11 +165,13 @@ lgamma(double x) EXTRACT_WORDS(hx,lx,x); /* purge off +-inf, NaN, +-0, and negative arguments */ + (*gamma_sign_ptr) = 1; ix = hx&0x7fffffff; if(ix>=0x7ff00000) return x*x; if((ix|lx)==0) return one/zero; if(ix<0x3b900000) { /* |x|<2**-70, return -log(|x|) */ if(hx<0) { + (*gamma_sign_ptr) = -1; return -log(-x); } else return -log(x); } @@ -179,6 +181,8 @@ lgamma(double x) t = sin_pi(x); if(t==zero) return one/zero; /* -integer */ nadj = log(pi/fabs(t*x)); + if(t=0x7f800000) return x*x; if(ix==0) return one/zero; if(ix<0x1c800000) { /* |x|<2**-70, return -log(|x|) */ if(hx<0) { + (*gamma_sign_ptr) = -1; return -logf(-x); } else return -logf(x); } @@ -192,6 +194,8 @@ lgammaf(float x) t = sin_pif(x); if(t==zero) return one/zero; /* -integer */ nadj = logf(pi/fabsf(t*x)); + if(t>31)&1; + + /* Extract exponent field. */ + j0 = ((i0 & 0x7ff00000) >> 20) - 1023; + + if(j0 < 20) + { + if(j0 < -1) + return 0; + else + { + w = TWO52[sx] + x; + t = w - TWO52[sx]; + GET_HIGH_WORD(i0, t); + /* Detect the all-zeros representation of plus and + minus zero, which fails the calculation below. */ + if ((i0 & ~(1 << 31)) == 0) + return 0; + j0 = ((i0 & 0x7ff00000) >> 20) - 1023; + i0 &= 0x000fffff; + i0 |= 0x00100000; + result = i0 >> (20 - j0); + } + } + else if (j0 < (int)(8 * sizeof (long int)) - 1) + { + if (j0 >= 52) + result = ((long int) ((i0 & 0x000fffff) | 0x0010000) << (j0 - 20)) | + (i1 << (j0 - 52)); + else + { + w = TWO52[sx] + x; + t = w - TWO52[sx]; + EXTRACT_WORDS (i0, i1, t); + j0 = ((i0 & 0x7ff00000) >> 20) - 1023; + i0 &= 0x000fffff; + i0 |= 0x00100000; + result = ((long int) i0 << (j0 - 20)) | (i1 >> (52 - j0)); + } + } + else + { + return (long int) x; + } + + return sx ? -result : result; } /****************************************************************************/ diff --git a/library/math_lrintf.c b/library/math_lrintf.c index c241bc8..d6c7cb3 100755 --- a/library/math_lrintf.c +++ b/library/math_lrintf.c @@ -1,5 +1,5 @@ /* - * $Id: math_lrintf.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_lrintf.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -29,6 +29,15 @@ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. + * + * + * PowerPC math library based in part on work by Sun Microsystems + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. */ #ifndef _MATH_HEADERS_H @@ -41,11 +50,60 @@ /****************************************************************************/ +/* Adding a float, x, to 2^23 will cause the result to be rounded based on + the fractional part of x, according to the implementation's current rounding + mode. 2^23 is the smallest float that can be represented using all 23 significant + digits. */ +static const float TWO23[2]={ + 8.3886080000e+06, /* 0x4b000000 */ + -8.3886080000e+06, /* 0xcb000000 */ +}; + +/****************************************************************************/ + long int lrintf(float x) { - /* ZZZ unimplemented */ - return(0); + LONG j0,sx; + ULONG i0; + float t; + volatile float w; + long int result; + + GET_FLOAT_WORD(i0,x); + + /* Extract sign bit. */ + sx = (i0 >> 31); + + /* Extract exponent field. */ + j0 = ((i0 & 0x7f800000) >> 23) - 127; + + if (j0 < (int)(sizeof (long int) * 8) - 1) + { + if (j0 < -1) + return 0; + else if (j0 >= 23) + result = (long int) ((i0 & 0x7fffff) | 0x800000) << (j0 - 23); + else + { + w = TWO23[sx] + x; + t = w - TWO23[sx]; + GET_FLOAT_WORD (i0, t); + /* Detect the all-zeros representation of plus and + minus zero, which fails the calculation below. */ + if ((i0 & ~(1 << 31)) == 0) + return 0; + j0 = ((i0 >> 23) & 0xff) - 0x7f; + i0 &= 0x7fffff; + i0 |= 0x800000; + result = i0 >> (23 - j0); + } + } + else + { + return (long int) x; + } + return sx ? -result : result; } /****************************************************************************/ diff --git a/library/math_lround.c b/library/math_lround.c index 7a5c8e3..08994a7 100644 --- a/library/math_lround.c +++ b/library/math_lround.c @@ -1,5 +1,5 @@ /* - * $Id: math_lround.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_lround.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -29,6 +29,15 @@ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. + * + * + * PowerPC math library based in part on work by Sun Microsystems + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. */ #ifndef _MATH_HEADERS_H @@ -44,8 +53,52 @@ long int lround(double x) { - /* ZZZ unimplemented */ - return(0); + LONG sign, exponent_less_1023; + /* Most significant word, least significant word. */ + ULONG msw, lsw; + long int result; + + EXTRACT_WORDS(msw, lsw, x); + + /* Extract sign. */ + sign = ((msw & 0x80000000) ? -1 : 1); + /* Extract exponent field. */ + exponent_less_1023 = ((msw & 0x7ff00000) >> 20) - 1023; + msw &= 0x000fffff; + msw |= 0x00100000; + + if (exponent_less_1023 < 20) + { + if (exponent_less_1023 < 0) + { + if (exponent_less_1023 < -1) + return 0; + else + return sign; + } + else + { + msw += 0x80000 >> exponent_less_1023; + result = msw >> (20 - exponent_less_1023); + } + } + else if (exponent_less_1023 < (8 * sizeof (long int)) - 1) + { + if (exponent_less_1023 >= 52) + result = ((long int) msw << (exponent_less_1023 - 20)) | (lsw << (exponent_less_1023 - 52)); + else + { + unsigned int tmp = lsw + (0x80000000 >> (exponent_less_1023 - 20)); + if (tmp < lsw) + ++msw; + result = ((long int) msw << (exponent_less_1023 - 20)) | (tmp >> (52 - exponent_less_1023)); + } + } + else + /* Result is too large to be represented by a long int. */ + return (long int)x; + + return sign * result; } /****************************************************************************/ diff --git a/library/math_lroundf.c b/library/math_lroundf.c index 1f4ad51..e7a212a 100755 --- a/library/math_lroundf.c +++ b/library/math_lroundf.c @@ -1,5 +1,5 @@ /* - * $Id: math_lroundf.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_lroundf.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -29,6 +29,15 @@ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. + * + * + * PowerPC math library based in part on work by Sun Microsystems + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. */ #ifndef _MATH_HEADERS_H @@ -44,8 +53,33 @@ long int lroundf(float x) { - /* ZZZ unimplemented */ - return(0); + LONG exponent_less_127; + ULONG w; + long int result; + LONG sign; + + GET_FLOAT_WORD (w, x); + exponent_less_127 = ((w & 0x7f800000) >> 23) - 127; + sign = (w & 0x80000000) != 0 ? -1 : 1; + w &= 0x7fffff; + w |= 0x800000; + + if (exponent_less_127 < (8 * sizeof (long int)) - 1) + { + if (exponent_less_127 < 0) + return exponent_less_127 < -1 ? 0 : sign; + else if (exponent_less_127 >= 23) + result = (long int) w << (exponent_less_127 - 23); + else + { + w += 0x400000 >> exponent_less_127; + result = w >> (23 - exponent_less_127); + } + } + else + return (long int) x; + + return sign * result; } /****************************************************************************/ diff --git a/library/math_nearbyint.c b/library/math_nearbyint.c index 17e2f02..91f1895 100644 --- a/library/math_nearbyint.c +++ b/library/math_nearbyint.c @@ -1,5 +1,5 @@ /* - * $Id: math_nearbyint.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_nearbyint.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -44,8 +44,8 @@ double nearbyint(double x) { - /* ZZZ unimplemented */ - return(0); + /* ZZZ is this such a good idea? */ + return(rint(x)); } /****************************************************************************/ diff --git a/library/math_nearbyintf.c b/library/math_nearbyintf.c index 9e56b85..041361c 100644 --- a/library/math_nearbyintf.c +++ b/library/math_nearbyintf.c @@ -1,5 +1,5 @@ /* - * $Id: math_nearbyintf.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_nearbyintf.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -44,8 +44,8 @@ float nearbyintf(float x) { - /* ZZZ unimplemented */ - return(0); + /* ZZZ is this such a good idea? */ + return(rintf(x)); } /****************************************************************************/ diff --git a/library/math_remquo.c b/library/math_remquo.c index ea7fb2a..54d1549 100644 --- a/library/math_remquo.c +++ b/library/math_remquo.c @@ -1,5 +1,5 @@ /* - * $Id: math_remquo.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_remquo.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -29,6 +29,12 @@ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. + * + * + * Copyright (C) 2002 by Red Hat, Incorporated. All rights reserved. + * + * Permission to use, copy, modify, and distribute this software + * is freely granted, provided that this notice is preserved. */ #ifndef _MATH_HEADERS_H @@ -44,8 +50,24 @@ double remquo(double x,double y,int * quo) { - /* ZZZ unimplemented */ - return(0); + int signx, signy, signres; + int mswx; + int mswy; + double x_over_y; + + GET_HIGH_WORD(mswx, x); + GET_HIGH_WORD(mswy, y); + + signx = (mswx & 0x80000000) >> 31; + signy = (mswy & 0x80000000) >> 31; + + signres = (signx ^ signy) ? -1 : 1; + + x_over_y = fabs(x / y); + + *quo = signres * (lrint(x_over_y) & 0x7f); + + return remainder(x,y); } /****************************************************************************/ diff --git a/library/math_remquof.c b/library/math_remquof.c index e212d4c..e4ee67f 100644 --- a/library/math_remquof.c +++ b/library/math_remquof.c @@ -1,5 +1,5 @@ /* - * $Id: math_remquof.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_remquof.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -29,6 +29,12 @@ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. + * + * + * Copyright (C) 2002 by Red Hat, Incorporated. All rights reserved. + * + * Permission to use, copy, modify, and distribute this software + * is freely granted, provided that this notice is preserved. */ #ifndef _MATH_HEADERS_H @@ -44,8 +50,24 @@ float remquof(float x,float y,int * quo) { - /* ZZZ unimplemented */ - return(0); + int signx, signy, signres; + int wx; + int wy; + float x_over_y; + + GET_FLOAT_WORD(wx, x); + GET_FLOAT_WORD(wy, y); + + signx = (wx & 0x80000000) >> 31; + signy = (wy & 0x80000000) >> 31; + + signres = (signx ^ signy) ? -1 : 1; + + x_over_y = fabsf(x / y); + + *quo = signres * (lrintf(x_over_y) & 0x7f); + + return remainderf(x,y); } /****************************************************************************/ diff --git a/library/math_round.c b/library/math_round.c index a6fdc6b..f93afe2 100644 --- a/library/math_round.c +++ b/library/math_round.c @@ -1,5 +1,5 @@ /* - * $Id: math_round.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_round.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -29,6 +29,15 @@ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. + * + * + * PowerPC math library based in part on work by Sun Microsystems + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. */ #ifndef _MATH_HEADERS_H @@ -44,8 +53,64 @@ double round(double x) { - /* ZZZ unimplemented */ - return(0); + /* Most significant word, least significant word. */ + LONG msw, exponent_less_1023; + ULONG lsw; + + EXTRACT_WORDS(msw, lsw, x); + + /* Extract exponent field. */ + exponent_less_1023 = ((msw & 0x7ff00000) >> 20) - 1023; + + if (exponent_less_1023 < 20) + { + if (exponent_less_1023 < 0) + { + msw &= 0x80000000; + if (exponent_less_1023 == -1) + /* Result is +1.0 or -1.0. */ + msw |= (1023 << 20); + lsw = 0; + } + else + { + ULONG exponent_mask = 0x000fffff >> exponent_less_1023; + if ((msw & exponent_mask) == 0 && lsw == 0) + /* x in an integral value. */ + return x; + + msw += 0x00080000 >> exponent_less_1023; + msw &= ~exponent_mask; + lsw = 0; + } + } + else if (exponent_less_1023 > 51) + { + if (exponent_less_1023 == 1024) + /* x is NaN or infinite. */ + return x + x; + else + return x; + } + else + { + ULONG exponent_mask = 0xffffffff >> (exponent_less_1023 - 20); + ULONG tmp; + + if ((lsw & exponent_mask) == 0) + /* x is an integral value. */ + return x; + + tmp = lsw + (1 << (51 - exponent_less_1023)); + if (tmp < lsw) + msw += 1; + lsw = tmp; + + lsw &= ~exponent_mask; + } + INSERT_WORDS(x, msw, lsw); + + return x; } /****************************************************************************/ diff --git a/library/math_roundf.c b/library/math_roundf.c index c179470..e796313 100644 --- a/library/math_roundf.c +++ b/library/math_roundf.c @@ -1,5 +1,5 @@ /* - * $Id: math_roundf.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_roundf.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -29,6 +29,15 @@ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. + * + * + * PowerPC math library based in part on work by Sun Microsystems + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. */ #ifndef _MATH_HEADERS_H @@ -44,8 +53,49 @@ float roundf(float x) { - /* ZZZ unimplemented */ - return(0); + int signbit; + ULONG w; + /* Most significant word, least significant word. */ + int exponent_less_127; + + GET_FLOAT_WORD(w, x); + + /* Extract sign bit. */ + signbit = w & 0x80000000; + + /* Extract exponent field. */ + exponent_less_127 = (int)((w & 0x7f800000) >> 23) - 127; + + if (exponent_less_127 < 23) + { + if (exponent_less_127 < 0) + { + w &= 0x80000000; + if (exponent_less_127 == -1) + /* Result is +1.0 or -1.0. */ + w |= (127 << 23); + } + else + { + unsigned int exponent_mask = 0x007fffff >> exponent_less_127; + if ((w & exponent_mask) == 0) + /* x has an integral value. */ + return x; + + w += 0x00400000 >> exponent_less_127; + w &= ~exponent_mask; + } + } + else + { + if (exponent_less_127 == 128) + /* x is NaN or infinite. */ + return x + x; + else + return x; + } + SET_FLOAT_WORD(x, w); + return x; } /****************************************************************************/ diff --git a/library/math_tgamma.c b/library/math_tgamma.c index e53f0e5..4b48b80 100644 --- a/library/math_tgamma.c +++ b/library/math_tgamma.c @@ -1,5 +1,5 @@ /* - * $Id: math_tgamma.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_tgamma.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -29,6 +29,15 @@ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. + * + * + * PowerPC math library based in part on work by Sun Microsystems + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. */ #ifndef _MATH_HEADERS_H @@ -44,8 +53,14 @@ double tgamma(double x) { - /* ZZZ unimplemented */ - return(0); + int gamma_sign; + double y; + + y = __lgamma(x,&gamma_sign); + if (gamma_sign < 0) + y = -y; + + return y; } /****************************************************************************/ diff --git a/library/math_tgammaf.c b/library/math_tgammaf.c index a277e21..1c1d236 100644 --- a/library/math_tgammaf.c +++ b/library/math_tgammaf.c @@ -1,5 +1,5 @@ /* - * $Id: math_tgammaf.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_tgammaf.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -29,6 +29,15 @@ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. + * + * + * PowerPC math library based in part on work by Sun Microsystems + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. */ #ifndef _MATH_HEADERS_H @@ -44,8 +53,14 @@ float tgammaf(float x) { - /* ZZZ unimplemented */ - return(0); + int gamma_sign; + float y; + + y = __lgammaf(x,&gamma_sign); + if (gamma_sign < 0) + y = -y; + + return y; } /****************************************************************************/ diff --git a/library/math_trunc.c b/library/math_trunc.c index 9d409b9..ac7a820 100644 --- a/library/math_trunc.c +++ b/library/math_trunc.c @@ -1,5 +1,5 @@ /* - * $Id: math_trunc.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_trunc.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -29,6 +29,15 @@ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. + * + * + * PowerPC math library based in part on work by Sun Microsystems + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. */ #ifndef _MATH_HEADERS_H @@ -44,8 +53,50 @@ double trunc(double x) { - /* ZZZ unimplemented */ - return(0); + int signbit; + /* Most significant word, least significant word. */ + int msw; + unsigned int lsw; + int exponent_less_1023; + + EXTRACT_WORDS(msw, lsw, x); + + /* Extract sign bit. */ + signbit = msw & 0x80000000; + + /* Extract exponent field. */ + exponent_less_1023 = ((msw & 0x7ff00000) >> 20) - 1023; + + if (exponent_less_1023 < 20) + { + /* All significant digits are in msw. */ + if (exponent_less_1023 < 0) + { + /* -1 < x < 1, so result is +0 or -0. */ + INSERT_WORDS(x, signbit, 0); + } + else + { + /* All relevant fraction bits are in msw, so lsw of the result is 0. */ + INSERT_WORDS(x, signbit | (msw & ~(0x000fffff >> exponent_less_1023)), 0); + } + } + else if (exponent_less_1023 > 51) + { + if (exponent_less_1023 == 1024) + { + /* x is infinite, or not a number, so trigger an exception. */ + return x + x; + } + /* All bits in the fraction fields of the msw and lsw are needed in the result. */ + } + else + { + /* All fraction bits in msw are relevant. Truncate irrelevant + bits from lsw. */ + INSERT_WORDS(x, msw, lsw & ~(0xffffffffu >> (exponent_less_1023 - 20))); + } + return x; } /****************************************************************************/ diff --git a/library/math_truncf.c b/library/math_truncf.c index dc7a56a..5de750b 100644 --- a/library/math_truncf.c +++ b/library/math_truncf.c @@ -1,5 +1,5 @@ /* - * $Id: math_truncf.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $ + * $Id: math_truncf.c,v 1.2 2005-10-09 10:38:55 obarthel Exp $ * * :ts=4 * @@ -29,6 +29,15 @@ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * POSSIBILITY OF SUCH DAMAGE. + * + * + * PowerPC math library based in part on work by Sun Microsystems + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. */ #ifndef _MATH_HEADERS_H @@ -44,8 +53,37 @@ float truncf(float x) { - /* ZZZ unimplemented */ - return(0); + LONG signbit, w, exponent_less_127; + + GET_FLOAT_WORD(w,x); + + /* Extract sign bit. */ + signbit = w & 0x80000000; + + /* Extract exponent field. */ + exponent_less_127 = ((w & 0x7f800000) >> 23) - 127; + + if (exponent_less_127 < 23) + { + if (exponent_less_127 < 0) + { + /* -1 < x < 1, so result is +0 or -0. */ + SET_FLOAT_WORD(x, signbit); + } + else + { + SET_FLOAT_WORD(x, signbit | (w & ~(0x007fffff >> exponent_less_127))); + } + } + else + { + if (exponent_less_127 == 255) + /* x is NaN or infinite. */ + return x + x; + + /* All bits in the fraction field are relevant. */ + } + return x; } /****************************************************************************/