updated gamma and gammaf

copied y0 y1 yn BeOS behaviour (even if the spec accepts both)
added ilogb ilogbf significand significandf (x86 and generic)


git-svn-id: file:///srv/svn/repos/haiku/haiku/trunk@17248 a95241bf-73f2-0310-859d-f6bbb57e9c96
This commit is contained in:
Jérôme Duval
2006-04-27 19:04:30 +00:00
parent f90e454336
commit b7b5cda158
15 changed files with 319 additions and 14 deletions
@@ -1,5 +1,5 @@
/* Implementation of gamma function according to ISO C.
Copyright (C) 1997, 1999, 2001 Free Software Foundation, Inc.
Copyright (C) 1997, 1999, 2001, 2004 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <[email protected]>, 1997.
@@ -35,9 +35,9 @@ __ieee754_gamma_r (double x, int *signgamp)
if (((hx & 0x7fffffff) | lx) == 0)
{
/* Return value for x == 0 is NaN with invalid exception. */
/* Return value for x == 0 is Inf with divide by zero exception. */
*signgamp = 0;
return x / x;
return 1.0 / x;
}
if (hx < 0 && (u_int32_t) hx < 0xfff00000 && __rint (x) == x)
{
@@ -1,5 +1,5 @@
/* Implementation of gamma function according to ISO C.
Copyright (C) 1997, 1999, 2001 Free Software Foundation, Inc.
Copyright (C) 1997, 1999, 2001, 2004 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <[email protected]>, 1997.
@@ -34,9 +34,9 @@ __ieee754_gammaf_r (float x, int *signgamp)
if ((hx & 0x7fffffff) == 0)
{
/* Return value for x == 0 is NaN with invalid exception. */
/* Return value for x == 0 is Inf with divide by zero exception. */
*signgamp = 0;
return x / x;
return 1.0 / x;
}
if (hx < 0 && (u_int32_t) hx < 0xff800000 && __rintf (x) == x)
{
@@ -185,10 +185,14 @@ V[] = {1.27304834834123699328e-02, /* 0x3F8A1270, 0x91C9C71A */
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0, y0(0) is -inf. */
if(ix>=0x7ff00000) return one/(x+x*x);
if((ix|lx)==0) return -one/zero;
if(hx<0) return zero/zero;
if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception. */
#ifdef __BEOS__
if(hx<0) return -HUGE_VAL+x; // also valid says the spec
#else
if(hx<0) return zero/(zero*x);
#endif
if(ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4
@@ -190,8 +190,12 @@ static double V0[5] = {
ix = 0x7fffffff&hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if(ix>=0x7ff00000) return one/(x+x*x);
if((ix|lx)==0) return -one/zero;
if(hx<0) return zero/zero;
if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception. */;
#ifdef __BEOS__
if(hx<0) return -HUGE_VAL+x;
#else
if(hx<0) return zero/(zero*x);
#endif
if(ix >= 0x40000000) { /* |x| >= 2.0 */
__sincos (x, &s, &c);
ss = -s-c;
@@ -20,7 +20,7 @@ static char rcsid[] = "$NetBSD: e_jn.c,v 1.9 1995/05/10 20:45:34 jtc Exp $";
* of order n
*
* Special cases:
* y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
* y0(0)=y1(0)=yn(n,0) = -inf with overflow signal;
* y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
* Note 2. About jn(n,x), yn(n,x)
* For n=0, j0(x) is called,
@@ -236,8 +236,12 @@ static double zero = 0.00000000000000000000e+00;
ix = 0x7fffffff&hx;
/* if Y(n,NaN) is NaN */
if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
if((ix|lx)==0) return -one/zero;
if(hx<0) return zero/zero;
if((ix|lx)==0) return -HUGE_VAL+x; /* -inf and overflow exception. */;
#ifdef __BEOS__
if(hx<0) return -HUGE_VAL+x;
#else
if(hx<0) return zero/(zero*x);
#endif
sign = 1;
if(n<0){
n = -n;
@@ -0,0 +1,64 @@
/* @(#)s_ilogb.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: s_ilogb.c,v 1.9 1995/05/10 20:47:28 jtc Exp $";
#endif
/* ilogb(double x)
* return the binary exponent of non-zero x
* ilogb(0) = FP_ILOGB0
* ilogb(NaN) = FP_ILOGBNAN (no signal is raised)
* ilogb(+-Inf) = INT_MAX (no signal is raised)
*/
#include <limits.h>
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
int __ilogb(double x)
#else
int __ilogb(x)
double x;
#endif
{
int32_t hx,lx,ix;
GET_HIGH_WORD(hx,x);
hx &= 0x7fffffff;
if(hx<0x00100000) {
GET_LOW_WORD(lx,x);
if((hx|lx)==0)
return FP_ILOGB0; /* ilogb(0) = FP_ILOGB0 */
else /* subnormal x */
if(hx==0) {
for (ix = -1043; lx>0; lx<<=1) ix -=1;
} else {
for (ix = -1022,hx<<=11; hx>0; hx<<=1) ix -=1;
}
return ix;
}
else if (hx<0x7ff00000) return (hx>>20)-1023;
else if (FP_ILOGBNAN != INT_MAX) {
/* ISO C99 requires ilogb(+-Inf) == INT_MAX. */
GET_LOW_WORD(lx,x);
if(((hx^0x7ff00000)|lx) == 0)
return INT_MAX;
}
return FP_ILOGBNAN;
}
weak_alias (__ilogb, ilogb)
#ifdef NO_LONG_DOUBLE
strong_alias (__ilogb, __ilogbl)
weak_alias (__ilogb, ilogbl)
#endif
@@ -0,0 +1,50 @@
/* s_ilogbf.c -- float version of s_ilogb.c.
* Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
*/
/*
* ====================================================
* 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.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: s_ilogbf.c,v 1.4 1995/05/10 20:47:31 jtc Exp $";
#endif
#include <limits.h>
#include "math.h"
#include "math_private.h"
#ifdef __STDC__
int __ilogbf(float x)
#else
int __ilogbf(x)
float x;
#endif
{
int32_t hx,ix;
GET_FLOAT_WORD(hx,x);
hx &= 0x7fffffff;
if(hx<0x00800000) {
if(hx==0)
return FP_ILOGB0; /* ilogb(0) = FP_ILOGB0 */
else /* subnormal x */
for (ix = -126,hx<<=8; hx>0; hx<<=1) ix -=1;
return ix;
}
else if (hx<0x7f800000) return (hx>>23)-127;
else if (FP_ILOGBNAN != INT_MAX) {
/* ISO C99 requires ilogbf(+-Inf) == INT_MAX. */
if (hx==0x7f800000)
return INT_MAX;
}
return FP_ILOGBNAN;
}
weak_alias (__ilogbf, ilogbf)
@@ -0,0 +1,39 @@
/* @(#)s_signif.c 5.1 93/09/24 */
/*
* ====================================================
* 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.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: s_significand.c,v 1.6 1995/05/10 20:48:11 jtc Exp $";
#endif
/*
* significand(x) computes just
* scalb(x, (double) -ilogb(x)),
* for exercising the fraction-part(F) IEEE 754-1985 test vector.
*/
#include <math.h>
#include "math_private.h"
#ifdef __STDC__
double __significand(double x)
#else
double __significand(x)
double x;
#endif
{
return __ieee754_scalb(x,(double) -__ilogb(x));
}
weak_alias (__significand, significand)
#ifdef NO_LONG_DOUBLE
strong_alias (__significand, __significandl)
weak_alias (__significand, significandl)
#endif
@@ -0,0 +1,32 @@
/* s_significandf.c -- float version of s_significand.c.
* Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
*/
/*
* ====================================================
* 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.
* ====================================================
*/
#if defined(LIBM_SCCS) && !defined(lint)
static char rcsid[] = "$NetBSD: s_significandf.c,v 1.3 1995/05/10 20:48:13 jtc Exp $";
#endif
#include <math.h>
#include "math_private.h"
#ifdef __STDC__
float __significandf(float x)
#else
float __significandf(x)
float x;
#endif
{
return __ieee754_scalbf(x,(float) -__ilogbf(x));
}
weak_alias (__significandf, significandf)
@@ -70,6 +70,7 @@ local genericSources =
s_floor.c s_floorf.c # s_floorl.c
s_fpclassify.c s_fpclassifyf.c # s_fpclassifyl.c
s_frexp.c s_frexpf.c # s_frexpl.c
s_ilogb.c s_ilogbf.c
s_isinf.c s_isinff.c # s_isinfl.c
s_ldexp.c s_ldexpf.c # s_ldexpl.c
s_log1p.c s_log1pf.c
@@ -79,6 +80,7 @@ local genericSources =
s_round.c s_roundf.c # s_roundl.c
s_scalbn.c s_scalbnf.c # s_scalbnl.c
s_signbit.c s_signbitf.c # s_signbitl.c
s_significand.c s_significandf.c
s_signgam.c
s_sin.c s_sinf.c # s_sinl.c
s_sincos.c s_sincosf.c
@@ -143,6 +143,7 @@ MergeObject posix_gnu_arch_$(TARGET_ARCH)_s.o :
s_finite.S s_finitef.S s_finitel.S
s_floor.S s_floorf.S s_floorl.S
s_frexp.S s_frexpf.S s_frexpl.S
s_ilogb.S s_ilogbf.S
s_isinfl.c
s_isnanl.c
s_llrint.S s_llrintf.S s_llrintl.S
@@ -151,6 +152,7 @@ MergeObject posix_gnu_arch_$(TARGET_ARCH)_s.o :
s_lrint.S s_lrintf.S s_lrintl.S
s_rint.S s_rintf.S s_rintl.c
s_scalbn.S s_scalbnf.S s_scalbnl.S
s_significand.S s_significandf.S
s_sin.S s_sinf.S s_sinl.S
s_sincos.S s_sincosf.S s_sincosl.S
s_tan.S s_tanf.S s_tanl.S
@@ -0,0 +1,36 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*/
#include <machine/asm.h>
RCSID("$NetBSD: s_ilogb.S,v 1.5 1995/10/12 15:53:09 jtc Exp $")
ENTRY(__ilogb)
fldl 4(%esp)
/* I added the following ugly construct because ilogb(+-Inf) is
required to return INT_MAX in ISO C99.
-- jakub@redhat.com. */
fxam /* Is NaN or +-Inf? */
fstsw %ax
movb $0x45, %dh
andb %ah, %dh
cmpb $0x05, %dh
je 1f /* Is +-Inf, jump. */
fxtract
pushl %eax
fstp %st
fistpl (%esp)
fwait
popl %eax
ret
1: fstp %st
movl $0x7fffffff, %eax
ret
END (__ilogb)
weak_alias (__ilogb, ilogb)
@@ -0,0 +1,36 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*/
#include <machine/asm.h>
RCSID("$NetBSD: s_ilogbf.S,v 1.4 1995/10/22 20:32:43 pk Exp $")
ENTRY(__ilogbf)
flds 4(%esp)
/* I added the following ugly construct because ilogb(+-Inf) is
required to return INT_MAX in ISO C99.
-- jakub@redhat.com. */
fxam /* Is NaN or +-Inf? */
fstsw %ax
movb $0x45, %dh
andb %ah, %dh
cmpb $0x05, %dh
je 1f /* Is +-Inf, jump. */
fxtract
pushl %eax
fstp %st
fistpl (%esp)
fwait
popl %eax
ret
1: fstp %st
movl $0x7fffffff, %eax
ret
END (__ilogbf)
weak_alias (__ilogbf, ilogbf)
@@ -0,0 +1,16 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*/
#include <machine/asm.h>
RCSID("$NetBSD: s_significand.S,v 1.4 1995/05/09 00:21:47 jtc Exp $")
ENTRY(__significand)
fldl 4(%esp)
fxtract
fstp %st(1)
ret
END (__significand)
weak_alias (__significand, significand)
@@ -0,0 +1,16 @@
/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*/
#include <machine/asm.h>
RCSID("$NetBSD: s_significandf.S,v 1.3 1995/05/09 00:24:07 jtc Exp $")
ENTRY(__significandf)
flds 4(%esp)
fxtract
fstp %st(1)
ret
END (__significandf)
weak_alias (__significandf, significandf)