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; } /****************************************************************************/