1
0
mirror of https://github.com/adtools/clib2.git synced 2025-12-08 14:59:05 +00:00

- 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
This commit is contained in:
Olaf Barthel
2005-10-09 10:38:56 +00:00
parent 159e55f1e6
commit bb2376a6ed
22 changed files with 616 additions and 71 deletions

View File

@ -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.

View File

@ -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()

View File

@ -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);
/****************************************************************************/

View File

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

View File

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

View File

@ -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);
/****************************************************************************/

View File

@ -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<zero)
(*gamma_sign_ptr) = -1;
x = -x;
}
@ -250,4 +254,17 @@ lgamma(double x)
/****************************************************************************/
double
lgamma(double x)
{
double result;
int gamma_sign;
result = __lgamma(x,&gamma_sign);
return(result);
}
/****************************************************************************/
#endif /* FLOATING_POINT_SUPPORT */

View File

@ -1,5 +1,5 @@
/*
* $Id: math_lgammaf.c,v 1.2 2005-05-30 08:47:26 obarthel Exp $
* $Id: math_lgammaf.c,v 1.3 2005-10-09 10:38:55 obarthel Exp $
*
* :ts=4
*
@ -170,7 +170,7 @@ sin_pif(float x)
float
lgammaf(float x)
__lgammaf(float x,int * gamma_sign_ptr)
{
float t,y,z,nadj=0.0,p,p1,p2,p3,q,r,w;
LONG i,hx,ix;
@ -178,11 +178,13 @@ lgammaf(float x)
GET_FLOAT_WORD(hx,x);
/* purge off +-inf, NaN, +-0, and negative arguments */
(*gamma_sign_ptr) = 1;
ix = hx&0x7fffffff;
if(ix>=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<zero)
(*gamma_sign_ptr) = -1;
x = -x;
}
@ -263,4 +267,17 @@ lgammaf(float x)
/****************************************************************************/
float
lgammaf(float x)
{
int gamma_sign;
float result;
result = __lgammaf(x,&gamma_sign);
return(result);
}
/****************************************************************************/
#endif /* FLOATING_POINT_SUPPORT */

View File

@ -1,5 +1,5 @@
/*
* $Id: math_lrint.c,v 1.1 2005-05-29 11:19:01 obarthel Exp $
* $Id: math_lrint.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,75 @@
/****************************************************************************/
/* Adding a double, x, to 2^52 will cause the result to be rounded based on
the fractional part of x, according to the implementation's current rounding
mode. 2^52 is the smallest double that can be represented using all 52 significant
digits. */
static const double TWO52[2]={
4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
-4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
};
/****************************************************************************/
long int
lrint(double x)
{
/* ZZZ unimplemented */
return(0);
LONG i0,j0,sx;
ULONG i1;
double t;
volatile double w;
long int result;
EXTRACT_WORDS(i0,i1,x);
/* Extract sign bit. */
sx = (i0>>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;
}
/****************************************************************************/

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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

View File

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