diff options
Diffstat (limited to 'sysdeps/ieee754/ldbl-96')
48 files changed, 3622 insertions, 0 deletions
diff --git a/sysdeps/ieee754/ldbl-96/e_acoshl.c b/sysdeps/ieee754/ldbl-96/e_acoshl.c new file mode 100644 index 0000000000..a60704aa29 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/e_acoshl.c @@ -0,0 +1,72 @@ +/* e_acoshl.c -- long double version of e_acosh.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* __ieee754_acoshl(x) + * Method : + * Based on + * acoshl(x) = logl [ x + sqrtl(x*x-1) ] + * we have + * acoshl(x) := logl(x)+ln2, if x is large; else + * acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else + * acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1. + * + * Special cases: + * acoshl(x) is NaN with signal if x<1. + * acoshl(NaN) is NaN without signal. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double +#else +static long double +#endif +one = 1.0, +ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */ + +#ifdef __STDC__ + long double __ieee754_acoshl(long double x) +#else + long double __ieee754_acoshl(x) + long double x; +#endif +{ + long double t; + u_int32_t se,i0,i1; + GET_LDOUBLE_WORDS(se,i0,i1,x); + if(se<0x3fff || se & 0x8000) { /* x < 1 */ + return (x-x)/(x-x); + } else if(se >=0x401b) { /* x > 2**28 */ + if(se >=0x7fff) { /* x is inf of NaN */ + return x+x; + } else + return __ieee754_logl(x)+ln2; /* acoshl(huge)=logl(2x) */ + } else if(((se-0x3fff)|i0|i1)==0) { + return 0.0; /* acosh(1) = 0 */ + } else if (se > 0x4000) { /* 2**28 > x > 2 */ + t=x*x; + return __ieee754_logl(2.0*x-one/(x+__ieee754_sqrtl(t-one))); + } else { /* 1<x<2 */ + t = x-one; + return __log1pl(t+__sqrtl(2.0*t+t*t)); + } +} diff --git a/sysdeps/ieee754/ldbl-96/e_atan2l.c b/sysdeps/ieee754/ldbl-96/e_atan2l.c new file mode 100644 index 0000000000..aff7a3d976 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/e_atan2l.c @@ -0,0 +1,136 @@ +/* e_atan2l.c -- long double version of e_atan2.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* __ieee754_atan2l(y,x) + * Method : + * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). + * 2. Reduce x to positive by (if x and y are unexceptional): + * ARG (x+iy) = arctan(y/x) ... if x > 0, + * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, + * + * Special cases: + * + * ATAN2((anything), NaN ) is NaN; + * ATAN2(NAN , (anything) ) is NaN; + * ATAN2(+-0, +(anything but NaN)) is +-0 ; + * ATAN2(+-0, -(anything but NaN)) is +-pi ; + * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2; + * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ; + * ATAN2(+-(anything but INF and NaN), -INF) is +-pi; + * ATAN2(+-INF,+INF ) is +-pi/4 ; + * ATAN2(+-INF,-INF ) is +-3pi/4; + * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2; + * + * Constants: + * The hexadecimal values are the intended ones for the following + * constants. The decimal values may be used, provided that the + * compiler will convert from decimal to binary accurately enough + * to produce the hexadecimal values shown. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double +#else +static long double +#endif +tiny = 1.0e-4900L, +zero = 0.0, +pi_o_4 = 7.85398163397448309628202E-01L, /* 0x3FFE, 0xC90FDAA2, 0x2168C235 */ +pi_o_2 = 1.5707963267948966192564E+00L, /* 0x3FFF, 0xC90FDAA2, 0x2168C235 */ +pi = 3.14159265358979323851281E+00L, /* 0x4000, 0xC90FDAA2, 0x2168C235 */ +pi_lo = -5.01655761266833202345176e-20L;/* 0xBFBE, 0xECE675D1, 0xFC8F8CBB */ + +#ifdef __STDC__ + long double __ieee754_atan2l(long double y, long double x) +#else + long double __ieee754_atan2l(y,x) + long double y,x; +#endif +{ + long double z; + int32_t k,m,hx,hy,ix,iy; + u_int32_t sx,sy,lx,ly; + + GET_LDOUBLE_WORDS(sx,hx,lx,x); + ix = sx&0x7fff; + lx |= hx & 0x7fffffff; + GET_LDOUBLE_WORDS(sy,hy,ly,y); + iy = sy&0x7fff; + ly |= hy & 0x7fffffff; + if(((2*ix|((lx|-lx)>>31))>0xfffe)|| + ((2*iy|((ly|-ly)>>31))>0xfffe)) /* x or y is NaN */ + return x+y; + if((sx-0x3fff|lx)==0) return __atanl(y); /* x=1.0 */ + m = ((sy>>15)&1)|((sx>>14)&2); /* 2*sign(x)+sign(y) */ + + /* when y = 0 */ + if((iy|ly)==0) { + switch(m) { + case 0: + case 1: return y; /* atan(+-0,+anything)=+-0 */ + case 2: return pi+tiny;/* atan(+0,-anything) = pi */ + case 3: return -pi-tiny;/* atan(-0,-anything) =-pi */ + } + } + /* when x = 0 */ + if((ix|lx)==0) return (sy>=0x8000)? -pi_o_2-tiny: pi_o_2+tiny; + + /* when x is INF */ + if(ix==0x7fff) { + if(iy==0x7fff) { + switch(m) { + case 0: return pi_o_4+tiny;/* atan(+INF,+INF) */ + case 1: return -pi_o_4-tiny;/* atan(-INF,+INF) */ + case 2: return 3.0*pi_o_4+tiny;/*atan(+INF,-INF)*/ + case 3: return -3.0*pi_o_4-tiny;/*atan(-INF,-INF)*/ + } + } else { + switch(m) { + case 0: return zero ; /* atan(+...,+INF) */ + case 1: return -zero ; /* atan(-...,+INF) */ + case 2: return pi+tiny ; /* atan(+...,-INF) */ + case 3: return -pi-tiny ; /* atan(-...,-INF) */ + } + } + } + /* when y is INF */ + if(iy==0x7fff) return (sy>=0x8000)? -pi_o_2-tiny: pi_o_2+tiny; + + /* compute y/x */ + k = sy-sx; + if(k > 70) z=pi_o_2+0.5*pi_lo; /* |y/x| > 2**70 */ + else if(sx>=0x8000&&k<-70) z=0.0; /* |y|/x < -2**70 */ + else z=__atanl(fabsl(y/x)); /* safe to do y/x */ + switch (m) { + case 0: return z ; /* atan(+,+) */ + case 1: { + u_int32_t sz; + GET_LDOUBLE_EXP(sz,z); + SET_LDOUBLE_EXP(z,sz ^ 0x8000); + } + return z ; /* atan(-,+) */ + case 2: return pi-(z-pi_lo);/* atan(+,-) */ + default: /* case 3 */ + return (z-pi_lo)-pi;/* atan(-,-) */ + } +} diff --git a/sysdeps/ieee754/ldbl-96/e_atanhl.c b/sysdeps/ieee754/ldbl-96/e_atanhl.c new file mode 100644 index 0000000000..fdcd1e9fe8 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/e_atanhl.c @@ -0,0 +1,79 @@ +/* s_atanhl.c -- long double version of s_atan.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* __ieee754_atanhl(x) + * Method : + * 1.Reduced x to positive by atanh(-x) = -atanh(x) + * 2.For x>=0.5 + * 1 2x x + * atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) + * 2 1 - x 1 - x + * + * For x<0.5 + * atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x)) + * + * Special cases: + * atanhl(x) is NaN if |x| > 1 with signal; + * atanhl(NaN) is that NaN with no signal; + * atanhl(+-1) is +-INF with signal. + * + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double one = 1.0, huge = 1e4900L; +#else +static long double one = 1.0, huge = 1e4900L; +#endif + +#ifdef __STDC__ +static const long double zero = 0.0; +#else +static double long zero = 0.0; +#endif + +#ifdef __STDC__ + long double __ieee754_atanhl(long double x) +#else + long double __ieee754_atanhl(x) + long double x; +#endif +{ + long double t; + int32_t ix; + u_int32_t se,i0,i1; + GET_LDOUBLE_WORDS(se,i0,i1,x); + ix = se&0x7fff; + if ((ix+((((i0&0x7fffffff)|i1)|(-((i0&0x7fffffff)|i1)))>>31))>0x3fff) + /* |x|>1 */ + return (x-x)/(x-x); + if(ix==0x3fff) + return x/zero; + if(ix<0x3fe3&&(huge+x)>zero) return x; /* x<2**-28 */ + SET_LDOUBLE_EXP(x,ix); + if(ix<0x3ffe) { /* x < 0.5 */ + t = x+x; + t = 0.5*__log1pl(t+t*x/(one-x)); + } else + t = 0.5*__log1pl((x+x)/(one-x)); + if(se<=0x7fff) return t; else return -t; +} diff --git a/sysdeps/ieee754/ldbl-96/e_coshl.c b/sysdeps/ieee754/ldbl-96/e_coshl.c new file mode 100644 index 0000000000..6af846cb2d --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/e_coshl.c @@ -0,0 +1,94 @@ +/* + * ==================================================== + * 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: e_cosh.c,v 1.7 1995/05/10 20:44:58 jtc Exp $"; +#endif + +/* __ieee754_coshl(x) + * Method : + * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2 + * 1. Replace x by |x| (coshl(x) = coshl(-x)). + * 2. + * [ exp(x) - 1 ]^2 + * 0 <= x <= ln2/2 : coshl(x) := 1 + ------------------- + * 2*exp(x) + * + * exp(x) + 1/exp(x) + * ln2/2 <= x <= 22 : coshl(x) := ------------------- + * 2 + * 22 <= x <= lnovft : coshl(x) := expl(x)/2 + * lnovft <= x <= ln2ovft: coshl(x) := expl(x/2)/2 * expl(x/2) + * ln2ovft < x : coshl(x) := huge*huge (overflow) + * + * Special cases: + * coshl(x) is |x| if x is +INF, -INF, or NaN. + * only coshl(0)=1 is exact for finite x. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double one = 1.0, half=0.5, huge = 1.0e4900L; +#else +static long double one = 1.0, half=0.5, huge = 1.0e4900L; +#endif + +#ifdef __STDC__ + long double __ieee754_coshl(long double x) +#else + long double __ieee754_coshl(x) + long double x; +#endif +{ + long double t,w; + int32_t ex; + u_int32_t mx,lx; + + /* High word of |x|. */ + GET_LDOUBLE_WORDS(ex,mx,lx,x); + ex &= 0x7fff; + + /* x is INF or NaN */ + if(ex==0x7fff) return x*x; + + /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */ + if(ex < 0x3ffd || (ex == 0x3ffd && mx < 0xb17217f7u)) { + t = __expm1l(fabsl(x)); + w = one+t; + if (ex<0x3fbc) return w; /* cosh(tiny) = 1 */ + return one+(t*t)/(w+w); + } + + /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */ + if (ex < 0x4003 || (ex == 0x4003 && mx < 0xb0000000u)) { + t = __ieee754_expl(fabsl(x)); + return half*t+half/t; + } + + /* |x| in [22, ln(maxdouble)] return half*exp(|x|) */ + if (ex < 0x400c || (ex == 0x400c && mx < 0xb1700000u)) + return half*__ieee754_expl(fabsl(x)); + + /* |x| in [log(maxdouble), overflowthresold] */ + if (ex < 0x400d + || (ex == 0x400d && (mx < 0xb170b513u + || (mx == 0xb170b513u && lx < 0xa1dfd60cu)))) + { + w = __ieee754_expl(half*fabsl(x)); + t = half*w; + return t*w; + } + + /* |x| > overflowthresold, cosh(x) overflow */ + return huge*huge; +} diff --git a/sysdeps/ieee754/ldbl-96/e_gammal_r.c b/sysdeps/ieee754/ldbl-96/e_gammal_r.c new file mode 100644 index 0000000000..104992450b --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/e_gammal_r.c @@ -0,0 +1,50 @@ +/* Implementation of gamma function according to ISO C. + Copyright (C) 1997, 1999 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> +#include <math_private.h> + + +long double +__ieee754_gammal_r (long double x, int *signgamp) +{ + /* We don't have a real gamma implementation now. We'll use lgamma + and the exp function. But due to the required boundary + conditions we must check some values separately. */ + u_int32_t es, hx, lx; + + GET_LDOUBLE_WORDS (es, hx, lx, x); + + if (((es & 0x7fff) | hx | lx) == 0) + { + /* Return value for x == 0 is NaN with invalid exception. */ + *signgamp = 0; + return x / x; + } + if ((hx & 0x8000) != 0 && (hx & 0x7fff) != 0x7fff && __rintl (x) == x) + { + /* Return value for integer x < 0 is NaN with invalid exception. */ + *signgamp = 0; + return (x - x) / (x - x); + } + + /* XXX FIXME. */ + return __ieee754_expl (__ieee754_lgammal_r (x, signgamp)); +} diff --git a/sysdeps/ieee754/ldbl-96/e_hypotl.c b/sysdeps/ieee754/ldbl-96/e_hypotl.c new file mode 100644 index 0000000000..1a40c556dc --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/e_hypotl.c @@ -0,0 +1,133 @@ +/* e_hypotl.c -- long double version of e_hypot.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* __ieee754_hypotl(x,y) + * + * Method : + * If (assume round-to-nearest) z=x*x+y*y + * has error less than sqrt(2)/2 ulp, than + * sqrt(z) has error less than 1 ulp (exercise). + * + * So, compute sqrt(x*x+y*y) with some care as + * follows to get the error below 1 ulp: + * + * Assume x>y>0; + * (if possible, set rounding to round-to-nearest) + * 1. if x > 2y use + * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y + * where x1 = x with lower 32 bits cleared, x2 = x-x1; else + * 2. if x <= 2y use + * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y)) + * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1, + * y1= y with lower 32 bits chopped, y2 = y-y1. + * + * NOTE: scaling may be necessary if some argument is too + * large or too tiny + * + * Special cases: + * hypot(x,y) is INF if x or y is +INF or -INF; else + * hypot(x,y) is NAN if x or y is NAN. + * + * Accuracy: + * hypot(x,y) returns sqrt(x^2+y^2) with error less + * than 1 ulps (units in the last place) + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + long double __ieee754_hypotl(long double x, long double y) +#else + long double __ieee754_hypotl(x,y) + long double x, y; +#endif +{ + long double a,b,t1,t2,y1,y2,w; + u_int32_t j,k,ea,eb; + + GET_LDOUBLE_EXP(ea,x); + ea &= 0x7fff; + GET_LDOUBLE_EXP(eb,y); + eb &= 0x7fff; + if(eb > ea) {a=y;b=x;j=ea; ea=eb;eb=j;} else {a=x;b=y;} + SET_LDOUBLE_EXP(a,ea); /* a <- |a| */ + SET_LDOUBLE_EXP(b,eb); /* b <- |b| */ + if((ea-eb)>0x46) {return a+b;} /* x/y > 2**70 */ + k=0; + if(ea > 0x5f3f) { /* a>2**8000 */ + if(ea == 0x7fff) { /* Inf or NaN */ + u_int32_t exp,high,low; + w = a+b; /* for sNaN */ + GET_LDOUBLE_WORDS(exp,high,low,a); + if(((high&0x7fffffff)|low)==0) w = a; + GET_LDOUBLE_WORDS(exp,high,low,b); + if(((eb^0x7fff)|(high&0x7fffffff)|low)==0) w = b; + return w; + } + /* scale a and b by 2**-9600 */ + ea -= 0x2580; eb -= 0x2580; k += 9600; + SET_LDOUBLE_EXP(a,ea); + SET_LDOUBLE_EXP(b,eb); + } + if(eb < 0x20bf) { /* b < 2**-8000 */ + if(eb == 0) { /* subnormal b or 0 */ + u_int32_t exp,high,low; + GET_LDOUBLE_WORDS(exp,high,low,b); + if((high|low)==0) return a; + SET_LDOUBLE_WORDS(t1, 0x7ffd, 0, 0); /* t1=2^16382 */ + b *= t1; + a *= t1; + k -= 16382; + } else { /* scale a and b by 2^9600 */ + ea += 0x2580; /* a *= 2^9600 */ + eb += 0x2580; /* b *= 2^9600 */ + k -= 9600; + SET_LDOUBLE_EXP(a,ea); + SET_LDOUBLE_EXP(b,eb); + } + } + /* medium size a and b */ + w = a-b; + if (w>b) { + u_int32_t high; + GET_LDOUBLE_MSW(high,a); + SET_LDOUBLE_WORDS(t1,ea,high,0); + t2 = a-t1; + w = __ieee754_sqrtl(t1*t1-(b*(-b)-t2*(a+t1))); + } else { + u_int32_t high; + GET_LDOUBLE_MSW(high,b); + a = a+a; + SET_LDOUBLE_WORDS(y1,eb,high,0); + y2 = b - y1; + GET_LDOUBLE_MSW(high,a); + SET_LDOUBLE_WORDS(t1,ea+1,high,0); + t2 = a - t1; + w = __ieee754_sqrtl(t1*y1-(w*(-w)-(t1*y2+t2*b))); + } + if(k!=0) { + u_int32_t exp; + t1 = 1.0; + GET_LDOUBLE_EXP(exp,t1); + SET_LDOUBLE_EXP(t1,exp+k); + return t1*w; + } else return w; +} diff --git a/sysdeps/ieee754/ldbl-96/e_remainderl.c b/sysdeps/ieee754/ldbl-96/e_remainderl.c new file mode 100644 index 0000000000..e721a6e8cd --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/e_remainderl.c @@ -0,0 +1,83 @@ +/* e_remainderl.c -- long double version of e_remainder.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* __ieee754_remainderl(x,p) + * Return : + * returns x REM p = x - [x/p]*p as if in infinite + * precise arithmetic, where [x/p] is the (infinite bit) + * integer nearest x/p (in half way case choose the even one). + * Method : + * Based on fmod() return x-[x/p]chopped*p exactlp. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double zero = 0.0; +#else +static long double zero = 0.0; +#endif + + +#ifdef __STDC__ + long double __ieee754_remainderl(long double x, long double p) +#else + long double __ieee754_remainderl(x,p) + long double x,p; +#endif +{ + u_int32_t sx,sex,sep,x0,x1,p0,p1; + long double p_half; + + GET_LDOUBLE_WORDS(sex,x0,x1,x); + GET_LDOUBLE_WORDS(sep,p0,p1,p); + sx = sex&0x8000; + sep &= 0x7fff; + sex &= 0x7fff; + + /* purge off exception values */ + if((sep|p0|p1)==0) return (x*p)/(x*p); /* p = 0 */ + if((sex==0x7fff)|| /* x not finite */ + ((sep==0x7fff)&& /* p is NaN */ + ((p0|p1)!=0))) + return (x*p)/(x*p); + + + if (sep<0x7ffe) x = __ieee754_fmodl(x,p+p); /* now x < 2p */ + if (((sex-sep)|(x0-p0)|(x1-p1))==0) return zero*x; + x = fabsl(x); + p = fabsl(p); + if (sep<0x0002) { + if(x+x>p) { + x-=p; + if(x+x>=p) x -= p; + } + } else { + p_half = 0.5*p; + if(x>p_half) { + x-=p; + if(x>=p_half) x -= p; + } + } + GET_LDOUBLE_EXP(sex,x); + SET_LDOUBLE_EXP(x,sex^sx); + return x; +} diff --git a/sysdeps/ieee754/ldbl-96/e_sinhl.c b/sysdeps/ieee754/ldbl-96/e_sinhl.c new file mode 100644 index 0000000000..4f9cfe2c38 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/e_sinhl.c @@ -0,0 +1,91 @@ +/* e_asinhl.c -- long double version of e_asinh.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* __ieee754_sinhl(x) + * Method : + * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 + * 1. Replace x by |x| (sinhl(-x) = -sinhl(x)). + * 2. + * E + E/(E+1) + * 0 <= x <= 25 : sinhl(x) := --------------, E=expm1l(x) + * 2 + * + * 25 <= x <= lnovft : sinhl(x) := expl(x)/2 + * lnovft <= x <= ln2ovft: sinhl(x) := expl(x/2)/2 * expl(x/2) + * ln2ovft < x : sinhl(x) := x*shuge (overflow) + * + * Special cases: + * sinhl(x) is |x| if x is +INF, -INF, or NaN. + * only sinhl(0)=0 is exact for finite x. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double one = 1.0, shuge = 1.0e4931L; +#else +static long double one = 1.0, shuge = 1.0e4931L; +#endif + +#ifdef __STDC__ + long double __ieee754_sinhl(long double x) +#else + long double __ieee754_sinhl(x) + long double x; +#endif +{ + long double t,w,h; + u_int32_t jx,ix,i0,i1; + + /* Words of |x|. */ + GET_LDOUBLE_WORDS(jx,i0,i1,x); + ix = jx&0x7fff; + + /* x is INF or NaN */ + if(ix==0x7fff) return x+x; + + h = 0.5; + if (jx & 0x8000) h = -h; + /* |x| in [0,25], return sign(x)*0.5*(E+E/(E+1))) */ + if (ix < 0x4003 || (ix == 0x4003 && i0 <= 0xc8000000)) { /* |x|<25 */ + if (ix<0x3fe3) /* |x|<2**-28 */ + if(shuge+x>one) return x;/* sinh(tiny) = tiny with inexact */ + t = __expm1l(fabsl(x)); + if(ix<0x3fff) return h*(2.0*t-t*t/(t+one)); + return h*(t+t/(t+one)); + } + + /* |x| in [25, log(maxdouble)] return 0.5*exp(|x|) */ + if (ix < 0x400c || (ix == 0x400c && i0 < 0xb17217f7)) + return h*__ieee754_expl(fabsl(x)); + + /* |x| in [log(maxdouble), overflowthreshold] */ + if (ix<0x400c || (ix == 0x400c && (i0 < 0xb174ddc0 + || (i0 == 0xb174ddc0 + && i1 <= 0x31aec0ea)))) { + w = __ieee754_expl(0.5*fabsl(x)); + t = h*w; + return t*w; + } + + /* |x| > overflowthreshold, sinhl(x) overflow */ + return x*shuge; +} diff --git a/sysdeps/ieee754/ldbl-96/ldbl2mpn.c b/sysdeps/ieee754/ldbl-96/ldbl2mpn.c new file mode 100644 index 0000000000..78e0b7097f --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/ldbl2mpn.c @@ -0,0 +1,95 @@ +/* Copyright (C) 1995, 1996, 1997, 1998 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "longlong.h" +#include <ieee754.h> +#include <float.h> +#include <math.h> +#include <stdlib.h> + +/* Convert a `long double' in IEEE854 standard double-precision format to a + multi-precision integer representing the significand scaled up by its + number of bits (64 for long double) and an integral power of two + (MPN frexpl). */ + +mp_size_t +__mpn_extract_long_double (mp_ptr res_ptr, mp_size_t size, + int *expt, int *is_neg, + long double value) +{ + union ieee854_long_double u; + u.d = value; + + *is_neg = u.ieee.negative; + *expt = (int) u.ieee.exponent - IEEE854_LONG_DOUBLE_BIAS; + +#if BITS_PER_MP_LIMB == 32 + res_ptr[0] = u.ieee.mantissa1; /* Low-order 32 bits of fraction. */ + res_ptr[1] = u.ieee.mantissa0; /* High-order 32 bits. */ + #define N 2 +#elif BITS_PER_MP_LIMB == 64 + /* Hopefully the compiler will combine the two bitfield extracts + and this composition into just the original quadword extract. */ + res_ptr[0] = ((unsigned long int) u.ieee.mantissa0 << 32) | u.ieee.mantissa1; + #define N 1 +#else + #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" +#endif + + if (u.ieee.exponent == 0) + { + /* A biased exponent of zero is a special case. + Either it is a zero or it is a denormal number. */ + if (res_ptr[0] == 0 && res_ptr[N - 1] == 0) /* Assumes N<=2. */ + /* It's zero. */ + *expt = 0; + else + { + /* It is a denormal number, meaning it has no implicit leading + one bit, and its exponent is in fact the format minimum. */ + int cnt; + + if (res_ptr[N - 1] != 0) + { + count_leading_zeros (cnt, res_ptr[N - 1]); + if (cnt != 0) + { +#if N == 2 + res_ptr[N - 1] = res_ptr[N - 1] << cnt + | (res_ptr[0] >> (BITS_PER_MP_LIMB - cnt)); + res_ptr[0] <<= cnt; +#else + res_ptr[N - 1] <<= cnt; +#endif + } + *expt = LDBL_MIN_EXP - 1 - cnt; + } + else + { + count_leading_zeros (cnt, res_ptr[0]); + res_ptr[N - 1] = res_ptr[0] << cnt; + res_ptr[0] = 0; + *expt = LDBL_MIN_EXP - 1 - BITS_PER_MP_LIMB - cnt; + } + } + } + + return N; +} diff --git a/sysdeps/ieee754/ldbl-96/math_ldbl.h b/sysdeps/ieee754/ldbl-96/math_ldbl.h new file mode 100644 index 0000000000..dccc4a1240 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/math_ldbl.h @@ -0,0 +1,98 @@ +#ifndef _MATH_PRIVATE_H_ +#error "Never use <math_ldbl.h> directly; include <math_private.h> instead." +#endif + +/* A union which permits us to convert between a long double and + three 32 bit ints. */ + +#if __FLOAT_WORD_ORDER == BIG_ENDIAN + +typedef union +{ + long double value; + struct + { + unsigned int sign_exponent:16; + unsigned int empty:16; + u_int32_t msw; + u_int32_t lsw; + } parts; +} ieee_long_double_shape_type; + +#endif + +#if __FLOAT_WORD_ORDER == LITTLE_ENDIAN + +typedef union +{ + long double value; + struct + { + u_int32_t lsw; + u_int32_t msw; + unsigned int sign_exponent:16; + unsigned int empty:16; + } parts; +} ieee_long_double_shape_type; + +#endif + +/* Get three 32 bit ints from a double. */ + +#define GET_LDOUBLE_WORDS(exp,ix0,ix1,d) \ +do { \ + ieee_long_double_shape_type ew_u; \ + ew_u.value = (d); \ + (exp) = ew_u.parts.sign_exponent; \ + (ix0) = ew_u.parts.msw; \ + (ix1) = ew_u.parts.lsw; \ +} while (0) + +/* Set a double from two 32 bit ints. */ + +#define SET_LDOUBLE_WORDS(d,exp,ix0,ix1) \ +do { \ + ieee_long_double_shape_type iw_u; \ + iw_u.parts.sign_exponent = (exp); \ + iw_u.parts.msw = (ix0); \ + iw_u.parts.lsw = (ix1); \ + (d) = iw_u.value; \ +} while (0) + +/* Get the more significant 32 bits of a long double mantissa. */ + +#define GET_LDOUBLE_MSW(v,d) \ +do { \ + ieee_long_double_shape_type sh_u; \ + sh_u.value = (d); \ + (v) = sh_u.parts.msw; \ +} while (0) + +/* Set the more significant 32 bits of a long double mantissa from an int. */ + +#define SET_LDOUBLE_MSW(d,v) \ +do { \ + ieee_long_double_shape_type sh_u; \ + sh_u.value = (d); \ + sh_u.parts.msw = (v); \ + (d) = sh_u.value; \ +} while (0) + +/* Get int from the exponent of a long double. */ + +#define GET_LDOUBLE_EXP(exp,d) \ +do { \ + ieee_long_double_shape_type ge_u; \ + ge_u.value = (d); \ + (exp) = ge_u.parts.sign_exponent; \ +} while (0) + +/* Set exponent of a long double from an int. */ + +#define SET_LDOUBLE_EXP(d,exp) \ +do { \ + ieee_long_double_shape_type se_u; \ + se_u.value = (d); \ + se_u.parts.sign_exponent = (exp); \ + (d) = se_u.value; \ +} while (0) diff --git a/sysdeps/ieee754/ldbl-96/mpn2ldbl.c b/sysdeps/ieee754/ldbl-96/mpn2ldbl.c new file mode 100644 index 0000000000..1f049ba12e --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/mpn2ldbl.c @@ -0,0 +1,47 @@ +/* Copyright (C) 1995, 1996, 1997, 1998 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include <ieee754.h> +#include <float.h> +#include <math.h> + +/* Convert a multi-precision integer of the needed number of bits (64 for + long double) and an integral power of two to a `long double' in IEEE854 + extended-precision format. */ + +long double +__mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign) +{ + union ieee854_long_double u; + + u.ieee.negative = sign; + u.ieee.exponent = expt + IEEE854_LONG_DOUBLE_BIAS; +#if BITS_PER_MP_LIMB == 32 + u.ieee.mantissa1 = frac_ptr[0]; + u.ieee.mantissa0 = frac_ptr[1]; +#elif BITS_PER_MP_LIMB == 64 + u.ieee.mantissa1 = frac_ptr[0] & ((1L << 32) - 1); + u.ieee.mantissa0 = frac_ptr[0] >> 32; +#else + #error "mp_limb size " BITS_PER_MP_LIMB "not accounted for" +#endif + + return u.d; +} diff --git a/sysdeps/ieee754/ldbl-96/printf_fphex.c b/sysdeps/ieee754/ldbl-96/printf_fphex.c new file mode 100644 index 0000000000..8dfa387df5 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/printf_fphex.c @@ -0,0 +1,61 @@ +#ifndef LONG_DOUBLE_DENORM_BIAS +# define LONG_DOUBLE_DENORM_BIAS (IEEE854_LONG_DOUBLE_BIAS - 1) +#endif + +#define PRINT_FPHEX_LONG_DOUBLE \ +do { \ + /* The "strange" 80 bit format on ix86 and m68k has an explicit \ + leading digit in the 64 bit mantissa. */ \ + unsigned long long int num; \ + \ + assert (sizeof (long double) == 12); \ + \ + num = (((unsigned long long int) fpnum.ldbl.ieee.mantissa0) << 32 \ + | fpnum.ldbl.ieee.mantissa1); \ + \ + zero_mantissa = num == 0; \ + \ + if (sizeof (unsigned long int) > 6) \ + numstr = _itoa_word (num, numbuf + sizeof numbuf, 16, \ + info->spec == 'A'); \ + else \ + numstr = _itoa (num, numbuf + sizeof numbuf, 16, info->spec == 'A'); \ + \ + /* Fill with zeroes. */ \ + while (numstr > numbuf + (sizeof numbuf - 64 / 4)) \ + *--numstr = '0'; \ + \ + /* We use a full nibble for the leading digit. */ \ + leading = *numstr++; \ + \ + /* We have 3 bits from the mantissa in the leading nibble. \ + Therefore we are here using `IEEE854_LONG_DOUBLE_BIAS + 3'. */ \ + exponent = fpnum.ldbl.ieee.exponent; \ + \ + if (exponent == 0) \ + { \ + if (zero_mantissa) \ + expnegative = 0; \ + else \ + { \ + /* This is a denormalized number. */ \ + expnegative = 1; \ + /* This is a hook for the m68k long double format, where the \ + exponent bias is the same for normalized and denormalized \ + numbers. */ \ + exponent = LONG_DOUBLE_DENORM_BIAS + 3; \ + } \ + } \ + else if (exponent >= IEEE854_LONG_DOUBLE_BIAS + 3) \ + { \ + expnegative = 0; \ + exponent -= IEEE854_LONG_DOUBLE_BIAS + 3; \ + } \ + else \ + { \ + expnegative = 1; \ + exponent = -(exponent - (IEEE854_LONG_DOUBLE_BIAS + 3)); \ + } \ +} while (0) + +#include <sysdeps/generic/printf_fphex.c> diff --git a/sysdeps/ieee754/ldbl-96/s_asinhl.c b/sysdeps/ieee754/ldbl-96/s_asinhl.c new file mode 100644 index 0000000000..6eb434c44b --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_asinhl.c @@ -0,0 +1,70 @@ +/* s_asinhl.c -- long double version of s_asinh.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* asinhl(x) + * Method : + * Based on + * asinhl(x) = signl(x) * logl [ |x| + sqrtl(x*x+1) ] + * we have + * asinhl(x) := x if 1+x*x=1, + * := signl(x)*(logl(x)+ln2)) for large |x|, else + * := signl(x)*logl(2|x|+1/(|x|+sqrtl(x*x+1))) if|x|>2, else + * := signl(x)*log1pl(|x| + x^2/(1 + sqrtl(1+x^2))) + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double +#else +static long double +#endif +one = 1.000000000000000000000e+00L, /* 0x3FFF, 0x00000000, 0x00000000 */ +ln2 = 6.931471805599453094287e-01L, /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */ +huge= 1.000000000000000000e+4900L; + +#ifdef __STDC__ + long double __asinhl(long double x) +#else + long double __asinhl(x) + long double x; +#endif +{ + long double t,w; + int32_t hx,ix; + GET_LDOUBLE_EXP(hx,x); + ix = hx&0x7fff; + if(ix==0x7fff) return x+x; /* x is inf or NaN */ + if(ix< 0x3fde) { /* |x|<2**-34 */ + if(huge+x>one) return x; /* return x inexact except 0 */ + } + if(ix>0x4020) { /* |x| > 2**34 */ + w = __ieee754_logl(fabsl(x))+ln2; + } else if (ix>0x4000) { /* 2**34 > |x| > 2.0 */ + t = fabsl(x); + w = __ieee754_logl(2.0*t+one/(__ieee754_sqrtl(x*x+one)+t)); + } else { /* 2.0 > |x| > 2**-28 */ + t = x*x; + w =__log1pl(fabsl(x)+t/(one+__ieee754_sqrtl(one+t))); + } + if(hx&0x8000) return -w; else return w; +} +weak_alias (__asinhl, asinhl) diff --git a/sysdeps/ieee754/ldbl-96/s_cbrtl.c b/sysdeps/ieee754/ldbl-96/s_cbrtl.c new file mode 100644 index 0000000000..1d021b7c3c --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_cbrtl.c @@ -0,0 +1,78 @@ +/* Compute cubic root of double value. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Dirk Alboth <dirka@uni-paderborn.de> and + Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include "math.h" +#include "math_private.h" + + +#define CBRT2 1.2599210498948731648 /* 2^(1/3) */ +#define SQR_CBRT2 1.5874010519681994748 /* 2^(2/3) */ + +/* We don't use long double values here since U need not be computed + with full precision. */ +static const double factor[5] = +{ + 1.0 / SQR_CBRT2, + 1.0 / CBRT2, + 1.0, + CBRT2, + SQR_CBRT2 +}; + + +long double +__cbrtl (long double x) +{ + long double xm, ym, u, t2; + int xe; + + /* Reduce X. XM now is an range 1.0 to 0.5. */ + xm = __frexpl (fabs (x), &xe); + + /* If X is not finite or is null return it (with raising exceptions + if necessary. + Note: *Our* version of `frexp' sets XE to zero if the argument is + Inf or NaN. This is not portable but faster. */ + if (xe == 0 && fpclassify (x) <= FP_ZERO) + return x + x; + + u = (0.338058687610520237 + + (1.67595307700780102 + + (-2.82414939754975962 + + (4.09559907378707839 + + (-4.11151425200350531 + + (2.65298938441952296 + + (-0.988553671195413709 + + 0.161617097923756032 * xm) + * xm) + * xm) + * xm) + * xm) + * xm) + *xm); + + t2 = u * u * u; + + ym = u * (t2 + 2.0 * xm) / (2.0 * t2 + xm) * factor[2 + xe % 3]; + + return __ldexpl (x > 0.0 ? ym : -ym, xe / 3); +} +weak_alias (__cbrtl, cbrtl) diff --git a/sysdeps/ieee754/ldbl-96/s_ceill.c b/sysdeps/ieee754/ldbl-96/s_ceill.c new file mode 100644 index 0000000000..d53f3954ba --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_ceill.c @@ -0,0 +1,89 @@ +/* s_ceill.c -- long double version of s_ceil.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* + * ceill(x) + * Return x rounded toward -inf to integral value + * Method: + * Bit twiddling. + * Exception: + * Inexact flag raised if x not equal to ceil(x). + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double huge = 1.0e4930; +#else +static long double huge = 1.0e4930; +#endif + +#ifdef __STDC__ + long double __ceill(long double x) +#else + long double __ceill(x) + long double x; +#endif +{ + int32_t i1,j0; + u_int32_t i,j,se,i0,sx; + GET_LDOUBLE_WORDS(se,i0,i1,x); + sx = (se>>15)&1; + j0 = (se&0x7fff)-0x3fff; + if(j0<31) { + if(j0<0) { /* raise inexact if x != 0 */ + if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ + if(sx) {se=0x8000;i0=0;i1=0;} + else if((i0|i1)!=0) { se=0x3fff;i0=0;i1=0;} + } + } else { + i = (0x7fffffff)>>j0; + if(((i0&i)|i1)==0) return x; /* x is integral */ + if(huge+x>0.0) { /* raise inexact flag */ + if(sx==0) { + if (j0>0) i0 += (0x80000000)>>j0; + else ++se; + } + i0 &= (~i); i1=0; + } + } + } else if (j0>62) { + if(j0==0x4000) return x+x; /* inf or NaN */ + else return x; /* x is integral */ + } else { + i = ((u_int32_t)(0xffffffff))>>(j0-31); + if((i1&i)==0) return x; /* x is integral */ + if(huge+x>0.0) { /* raise inexact flag */ + if(sx==0) { + if(j0==31) i0+=1; + else { + j = i1 + (1<<(63-j0)); + if(j<i1) i0+=1; /* got a carry */ + i1 = j; + } + } + i1 &= (~i); + } + } + SET_LDOUBLE_WORDS(x,se,i0,i1); + return x; +} +weak_alias (__ceill, ceill) diff --git a/sysdeps/ieee754/ldbl-96/s_copysignl.c b/sysdeps/ieee754/ldbl-96/s_copysignl.c new file mode 100644 index 0000000000..9976575b3b --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_copysignl.c @@ -0,0 +1,43 @@ +/* s_copysignl.c -- long double version of s_copysign.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* + * copysignl(long double x, long double y) + * copysignl(x,y) returns a value with the magnitude of x and + * with the sign bit of y. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + long double __copysignl(long double x, long double y) +#else + long double __copysignl(x,y) + long double x,y; +#endif +{ + u_int32_t es1,es2; + GET_LDOUBLE_EXP(es1,x); + GET_LDOUBLE_EXP(es2,y); + SET_LDOUBLE_EXP(x,(es1&0x7fff)|(es2&0x8000)); + return x; +} +weak_alias (__copysignl, copysignl) diff --git a/sysdeps/ieee754/ldbl-96/s_cosl.c b/sysdeps/ieee754/ldbl-96/s_cosl.c new file mode 100644 index 0000000000..9765f7fd4e --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_cosl.c @@ -0,0 +1,88 @@ +/* s_cosl.c -- long double version of s_cos.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* cosl(x) + * Return cosine function of x. + * + * kernel function: + * __kernel_sinl ... sine function on [-pi/4,pi/4] + * __kernel_cosl ... cosine function on [-pi/4,pi/4] + * __ieee754_rem_pio2l ... argument reduction routine + * + * Method. + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * in [-pi/4 , +pi/4], and let n = k mod 4. + * We have + * + * n sin(x) cos(x) tan(x) + * ---------------------------------------------------------- + * 0 S C T + * 1 C -S -1/T + * 2 -S -C T + * 3 -C S -1/T + * ---------------------------------------------------------- + * + * Special cases: + * Let trig be any of sin, cos, or tan. + * trig(+-INF) is NaN, with signals; + * trig(NaN) is that NaN; + * + * Accuracy: + * TRIG(x) returns trig(x) nearly rounded + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + long double __cosl(long double x) +#else + long double __cosl(x) + long double x; +#endif +{ + long double y[2],z=0.0; + int32_t n, se, i0, i1; + + /* High word of x. */ + GET_LDOUBLE_WORDS(se,i0,i1,x); + + /* |x| ~< pi/4 */ + se &= 0x7fff; + if(se < 0x3ffe || (se == 0x3ffe && i0 <= 0xc90fdaa2)) + return __kernel_cosl(x,z); + + /* cos(Inf or NaN) is NaN */ + else if (se==0x7fff) return x-x; + + /* argument reduction needed */ + else { + n = __ieee754_rem_pio2l(x,y); + switch(n&3) { + case 0: return __kernel_cosl(y[0],y[1]); + case 1: return -__kernel_sinl(y[0],y[1],1); + case 2: return -__kernel_cosl(y[0],y[1]); + default: + return __kernel_sinl(y[0],y[1],1); + } + } +} +weak_alias (__cosl, cosl) diff --git a/sysdeps/ieee754/ldbl-96/s_fabsl.c b/sysdeps/ieee754/ldbl-96/s_fabsl.c new file mode 100644 index 0000000000..f7170503fb --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_fabsl.c @@ -0,0 +1,40 @@ +/* s_fabsl.c -- long double version of s_fabs.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* + * fabsl(x) returns the absolute value of x. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + long double __fabsl(long double x) +#else + long double __fabsl(x) + long double x; +#endif +{ + u_int32_t exp; + GET_LDOUBLE_EXP(exp,x); + SET_LDOUBLE_EXP(x,exp&0x7fff); + return x; +} +weak_alias (__fabsl, fabsl) diff --git a/sysdeps/ieee754/ldbl-96/s_finitel.c b/sysdeps/ieee754/ldbl-96/s_finitel.c new file mode 100644 index 0000000000..6e444e90d0 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_finitel.c @@ -0,0 +1,40 @@ +/* s_finitel.c -- long double version of s_finite.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* + * finitel(x) returns 1 is x is finite, else 0; + * no branching! + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + int __finitel(long double x) +#else + int __finitel(x) + long double x; +#endif +{ + int32_t exp; + GET_LDOUBLE_EXP(exp,x); + return (int)((u_int32_t)((exp&0x7fff)-0x7fff)>>31); +} +weak_alias (__finitel, finitel) diff --git a/sysdeps/ieee754/ldbl-96/s_floorl.c b/sysdeps/ieee754/ldbl-96/s_floorl.c new file mode 100644 index 0000000000..fb0c37e801 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_floorl.c @@ -0,0 +1,90 @@ +/* s_floorl.c -- long double version of s_floor.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* + * floorl(x) + * Return x rounded toward -inf to integral value + * Method: + * Bit twiddling. + * Exception: + * Inexact flag raised if x not equal to floor(x). + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double huge = 1.0e4930; +#else +static long double huge = 1.0e4930; +#endif + +#ifdef __STDC__ + long double __floorl(long double x) +#else + long double __floorl(x) + long double x; +#endif +{ + int32_t i1,j0; + u_int32_t i,j,se,i0,sx; + GET_LDOUBLE_WORDS(se,i0,i1,x); + sx = (se>>15)&1; + j0 = (se&0x7fff)-0x3fff; + if(j0<31) { + if(j0<0) { /* raise inexact if x != 0 */ + if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */ + if(sx==0) {se=0;i0=i1=0;} + else if(((se&0x7fff)|i0|i1)!=0) + { se=0xbfff;i0=i1=0;} + } + } else { + i = (0x7fffffff)>>j0; + if(((i0&i)|i1)==0) return x; /* x is integral */ + if(huge+x>0.0) { /* raise inexact flag */ + if(sx) { + if (j0>0) i0 += (0x80000000)>>j0; + else ++se; + } + i0 &= (~i); i1=0; + } + } + } else if (j0>62) { + if(j0==0x4000) return x+x; /* inf or NaN */ + else return x; /* x is integral */ + } else { + i = ((u_int32_t)(0xffffffff))>>(j0-31); + if((i1&i)==0) return x; /* x is integral */ + if(huge+x>0.0) { /* raise inexact flag */ + if(sx) { + if(j0==31) i0+=1; + else { + j = i1+(1<<(63-j0)); + if(j<i1) i0 +=1 ; /* got a carry */ + i1=j; + } + } + i1 &= (~i); + } + } + SET_LDOUBLE_WORDS(x,se,i0,i1); + return x; +} +weak_alias (__floorl, floorl) diff --git a/sysdeps/ieee754/ldbl-96/s_fpclassifyl.c b/sysdeps/ieee754/ldbl-96/s_fpclassifyl.c new file mode 100644 index 0000000000..4df0b44f75 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_fpclassifyl.c @@ -0,0 +1,44 @@ +/* Return classification value corresponding to argument. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + Fixed by Andreas Schwab <schwab@issan.informatik.uni-dortmund.de>. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +int +__fpclassifyl (long double x) +{ + u_int32_t ex, hx, lx, m; + int retval = FP_NORMAL; + + GET_LDOUBLE_WORDS (ex, hx, lx, x); + m = (hx & 0x7fffffff) | lx; + ex &= 0x7fff; + if ((ex | m) == 0) + retval = FP_ZERO; + else if (ex == 0 && (hx & 0x80000000) == 0) + retval = FP_SUBNORMAL; + else if (ex == 0x7fff) + retval = m != 0 ? FP_NAN : FP_INFINITE; + + return retval; +} diff --git a/sysdeps/ieee754/ldbl-96/s_frexpl.c b/sysdeps/ieee754/ldbl-96/s_frexpl.c new file mode 100644 index 0000000000..7a49bc37c8 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_frexpl.c @@ -0,0 +1,69 @@ +/* s_frexpl.c -- long double version of s_frexp.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* + * for non-zero x + * x = frexpl(arg,&exp); + * return a long double fp quantity x such that 0.5 <= |x| <1.0 + * and the corresponding binary exponent "exp". That is + * arg = x*2^exp. + * If arg is inf, 0.0, or NaN, then frexpl(arg,&exp) returns arg + * with *exp=0. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double +#else +static long double +#endif +#if LDBL_MANT_DIG == 64 +two65 = 3.68934881474191032320e+19L; /* 0x4040, 0x80000000, 0x00000000 */ +#else +# error "Cannot handle this MANT_DIG" +#endif + + +#ifdef __STDC__ + long double __frexpl(long double x, int *eptr) +#else + long double __frexpl(x, eptr) + long double x; int *eptr; +#endif +{ + u_int32_t se, hx, ix, lx; + GET_LDOUBLE_WORDS(se,hx,lx,x); + ix = 0x7fff&se; + *eptr = 0; + if(ix==0x7fff||((ix|hx|lx)==0)) return x; /* 0,inf,nan */ + if (ix==0x0000) { /* subnormal */ + x *= two65; + GET_LDOUBLE_EXP(se,x); + ix = se&0x7fff; + *eptr = -65; + } + *eptr += ix-16382; + se = (se & 0x8000) | 0x3ffe; + SET_LDOUBLE_EXP(x,se); + return x; +} +weak_alias (__frexpl, frexpl) diff --git a/sysdeps/ieee754/ldbl-96/s_ilogbl.c b/sysdeps/ieee754/ldbl-96/s_ilogbl.c new file mode 100644 index 0000000000..d44229dcda --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_ilogbl.c @@ -0,0 +1,56 @@ +/* s_ilogbl.c -- long double version of s_ilogb.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* ilogbl(long double x) + * return the binary exponent of non-zero x + * ilogbl(0) = 0x80000001 + * ilogbl(inf/NaN) = 0x7fffffff (no signal is raised) + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + int __ilogbl(long double x) +#else + int __ilogbl(x) + long double x; +#endif +{ + int32_t es,hx,lx,ix; + + GET_LDOUBLE_EXP(es,x); + es &= 0x7fff; + if(es==0) { + GET_LDOUBLE_WORDS(es,hx,lx,x); + if((hx|lx)==0) + return FP_ILOGB0; /* ilogbl(0) = FP_ILOGB0 */ + else /* subnormal x */ + if(hx==0) { + for (ix = -16415; lx>0; lx<<=1) ix -=1; + } else { + for (ix = -16383; hx>0; hx<<=1) ix -=1; + } + return ix; + } + else if (es<0x7fff) return es-0x3fff; + else return FP_ILOGBNAN; +} +weak_alias (__ilogbl, ilogbl) diff --git a/sysdeps/ieee754/ldbl-96/s_isinfl.c b/sysdeps/ieee754/ldbl-96/s_isinfl.c new file mode 100644 index 0000000000..6f7c07c5af --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_isinfl.c @@ -0,0 +1,29 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Change for long double by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: $"; +#endif + +/* + * isinfl(x) returns 1 if x is inf, -1 if x is -inf, else 0; + * no branching! + */ + +#include "math.h" +#include "math_private.h" + +int +__isinfl (long double x) +{ + int32_t se,hx,lx; + GET_LDOUBLE_WORDS(se,hx,lx,x); + lx |= (hx & 0x7fffffff) | ((se & 0x7fff) ^ 0x7fff); + lx |= -lx; + se &= 0x8000; + return ~(lx >> 31) & (1 - (se >> 14)); +} +weak_alias (__isinfl, isinfl) diff --git a/sysdeps/ieee754/ldbl-96/s_isnanl.c b/sysdeps/ieee754/ldbl-96/s_isnanl.c new file mode 100644 index 0000000000..0a7ff38433 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_isnanl.c @@ -0,0 +1,44 @@ +/* s_isnanl.c -- long double version of s_isnan.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* + * isnanl(x) returns 1 is x is nan, else 0; + * no branching! + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + int __isnanl(long double x) +#else + int __isnanl(x) + long double x; +#endif +{ + int32_t se,hx,lx; + GET_LDOUBLE_WORDS(se,hx,lx,x); + se = (se & 0x7fff) << 1; + lx |= hx & 0x7fffffff; + se |= (u_int32_t)(lx|(-lx))>>31; + se = 0xfffe - se; + return (int)(((u_int32_t)(se))>>31); +} +weak_alias (__isnanl, isnanl) diff --git a/sysdeps/ieee754/ldbl-96/s_llrintl.c b/sysdeps/ieee754/ldbl-96/s_llrintl.c new file mode 100644 index 0000000000..2aeaa1e102 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_llrintl.c @@ -0,0 +1,77 @@ +/* Round argument to nearest integral value according to current rounding + direction. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + +static const long double two63[2] = +{ + 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ + -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ +}; + + +long long int +__llrintl (long double x) +{ + int32_t se,j0; + u_int32_t i0, i1; + long long int result; + volatile long double w; + long double t; + int sx; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + + sx = (se >> 15) & 1; + j0 = (se & 0x7fff) - 0x3fff; + + if (j0 < (int32_t) (8 * sizeof (long long int)) - 1) + { + if (j0 < -1) + return 0; + else if (j0 >= 63) + result = (((long long int) i0 << 32) | i1) << (j0 - 63); + else + { + w = two63[sx] + x; + t = w - two63[sx]; + GET_LDOUBLE_WORDS (se, i0, i1, t); + j0 = (se & 0x7fff) - 0x3fff; + + if (j0 < 31) + result = i0 >> (31 - j0); + else + result = ((long long int) i0 << (j0 - 31)) | (i1 >> (63 - j0)); + } + } + else + { + /* The number is too large. It is left implementation defined + what happens. */ + return (long long int) x; + } + + return sx ? -result : result; +} + +weak_alias (__llrintl, llrintl) diff --git a/sysdeps/ieee754/ldbl-96/s_llroundl.c b/sysdeps/ieee754/ldbl-96/s_llroundl.c new file mode 100644 index 0000000000..4a537c871e --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_llroundl.c @@ -0,0 +1,81 @@ +/* Round long double value to long long int. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +long long int +__llroundl (long double x) +{ + int32_t j0; + u_int32_t se, i1, i0; + long long int result; + int sign; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + j0 = (se & 0x7fff) - 0x3fff; + sign = (se & 0x8000) != 0 ? -1 : 1; + + if (j0 < 31) + { + if (j0 < 0) + return j0 < -1 ? 0 : sign; + else + { + u_int32_t j = i0 + (0x40000000 >> j0); + if (j < i0) + { + j >>= 1; + j |= 0x80000000; + ++j0; + } + + result = j >> (31 - j0); + } + } + else if (j0 < (int32_t) (8 * sizeof (long long int)) - 1) + { + if (j0 >= 63) + result = (((long long int) i0 << 32) | i1) << (j0 - 63); + else + { + u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); + if (j < i1) + ++i0; + + if (j0 == 31) + result = (long long int) i0; + else + result = ((long long int) i0 << (j0 - 31)) | (j >> (63 - j0)); + } + } + else + { + /* The number is too large. It is left implementation defined + what happens. */ + return (long long int) x; + } + + return sign * result; +} + +weak_alias (__llroundl, llroundl) diff --git a/sysdeps/ieee754/ldbl-96/s_logbl.c b/sysdeps/ieee754/ldbl-96/s_logbl.c new file mode 100644 index 0000000000..2cd9d105f8 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_logbl.c @@ -0,0 +1,47 @@ +/* s_logbl.c -- long double version of s_logb.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* + * long double logbl(x) + * IEEE 754 logb. Included to pass IEEE test suite. Not recommend. + * Use ilogb instead. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + long double __logbl(long double x) +#else + long double __logbl(x) + long double x; +#endif +{ + int32_t es,lx,ix; + GET_LDOUBLE_WORDS(es,ix,lx,x); + es &= 0x7fff; /* exponent */ + if((es|ix|lx)==0) return -1.0/fabs(x); + if(es==0x7fff) return x*x; + if(es==0) /* IEEE 754 logb */ + return -16382.0; + else + return (long double) (es-0x3fff); +} +weak_alias (__logbl, logbl) diff --git a/sysdeps/ieee754/ldbl-96/s_lrintl.c b/sysdeps/ieee754/ldbl-96/s_lrintl.c new file mode 100644 index 0000000000..673cf3d9db --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_lrintl.c @@ -0,0 +1,86 @@ +/* Round argument to nearest integral value according to current rounding + direction. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + +static const long double two63[2] = +{ + 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ + -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ +}; + + +long int +__lrintl (long double x) +{ + int32_t se,j0; + u_int32_t i0, i1; + long int result; + volatile long double w; + long double t; + int sx; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + + sx = (se >> 15) & 1; + j0 = (se & 0x7fff) - 0x3fff; + + if (j0 < 31) + { + if (j0 < -1) + return 0; + else + { + w = two63[sx] + x; + t = w - two63[sx]; + GET_LDOUBLE_WORDS (se, i0, i1, t); + j0 = (se & 0x7fff) - 0x3fff; + + result = i0 >> (31 - j0); + } + } + else if (j0 < (int32_t) (8 * sizeof (long int)) - 1) + { + if (j0 >= 63) + result = ((long int) i0 << (j0 - 31)) | (i1 << (j0 - 63)); + else + { + w = two63[sx] + x; + t = w - two63[sx]; + GET_LDOUBLE_WORDS (se, i0, i1, t); + j0 = (se & 0x7fff) - 0x3fff; + + result = ((long int) i0 << (j0 - 31)) | (i1 >> (63 - j0)); + } + } + else + { + /* The number is too large. It is left implementation defined + what happens. */ + return (long int) x; + } + + return sx ? -result : result; +} + +weak_alias (__lrintl, lrintl) diff --git a/sysdeps/ieee754/ldbl-96/s_lroundl.c b/sysdeps/ieee754/ldbl-96/s_lroundl.c new file mode 100644 index 0000000000..3bdac830b4 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_lroundl.c @@ -0,0 +1,78 @@ +/* Round long double value to long int. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +long int +__lroundl (long double x) +{ + int32_t j0; + u_int32_t se, i1, i0; + long int result; + int sign; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + j0 = (se & 0x7fff) - 0x3fff; + sign = (se & 0x8000) != 0 ? -1 : 1; + + if (j0 < 31) + { + if (j0 < 0) + return j0 < -1 ? 0 : sign; + else + { + u_int32_t j = i0 + (0x40000000 >> j0); + if (j < i0) + { + j >>= 1; + j |= 0x80000000; + ++j0; + } + + result = j >> (31 - j0); + } + } + else if (j0 < (int32_t) (8 * sizeof (long int)) - 1) + { + if (j0 >= 63) + result = ((long int) i0 << (j0 - 31)) | (i1 << (j0 - 63)); + else + { + u_int32_t j = i1 + (0x80000000 >> (j0 - 31)); + if (j < i1) + ++i0; + + result = ((long int) i0 << (j0 - 31)) | (j >> (63 - j0)); + } + } + else + { + /* The number is too large. It is left implementation defined + what happens. */ + return (long int) x; + } + + return sign * result; +} + +weak_alias (__lroundl, lroundl) diff --git a/sysdeps/ieee754/ldbl-96/s_modfl.c b/sysdeps/ieee754/ldbl-96/s_modfl.c new file mode 100644 index 0000000000..fb1b3acf30 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_modfl.c @@ -0,0 +1,85 @@ +/* s_modfl.c -- long double version of s_modf.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* + * modfl(long double x, long double *iptr) + * return fraction part of x, and return x's integral part in *iptr. + * Method: + * Bit twiddling. + * + * Exception: + * No exception. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double one = 1.0; +#else +static long double one = 1.0; +#endif + +#ifdef __STDC__ + long double __modfl(long double x, long double *iptr) +#else + long double __modfl(x, iptr) + long double x,*iptr; +#endif +{ + int32_t i0,i1,j0; + u_int32_t i,se; + GET_LDOUBLE_WORDS(se,i0,i1,x); + j0 = (se&0x7fff)-0x3fff; /* exponent of x */ + if(j0<32) { /* integer part in high x */ + if(j0<0) { /* |x|<1 */ + SET_LDOUBLE_WORDS(*iptr,se&0x8000,0,0); /* *iptr = +-0 */ + return x; + } else { + i = (0x7fffffff)>>j0; + if(((i0&i)|i1)==0) { /* x is integral */ + *iptr = x; + SET_LDOUBLE_WORDS(x,se&0x8000,0,0); /* return +-0 */ + return x; + } else { + SET_LDOUBLE_WORDS(*iptr,se,i0&(~i),0); + return x - *iptr; + } + } + } else if (j0>63) { /* no fraction part */ + *iptr = x*one; + /* We must handle NaNs separately. */ + if (j0 == 0x4000 && ((i0 & 0x7fffffff) | i1)) + return x*one; + SET_LDOUBLE_WORDS(x,se&0x8000,0,0); /* return +-0 */ + return x; + } else { /* fraction part in low x */ + i = ((u_int32_t)(0xffffffff))>>(j0-32); + if((i1&i)==0) { /* x is integral */ + *iptr = x; + SET_LDOUBLE_WORDS(x,se&0x8000,0,0); /* return +-0 */ + return x; + } else { + SET_LDOUBLE_WORDS(*iptr,se,i0,i1&(~i)); + return x - *iptr; + } + } +} +weak_alias (__modfl, modfl) diff --git a/sysdeps/ieee754/ldbl-96/s_nearbyintl.c b/sysdeps/ieee754/ldbl-96/s_nearbyintl.c new file mode 100644 index 0000000000..92c3ebf368 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_nearbyintl.c @@ -0,0 +1,95 @@ +/* s_rintl.c -- long double version of s_rint.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ +/* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. */ + +/* + * ==================================================== + * 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. + * ==================================================== + */ + +/* + * rintl(x) + * Return x rounded to integral value according to the prevailing + * rounding mode. + * Method: + * Using floating addition. + * Exception: + * Inexact flag raised if x not equal to rintl(x). + */ + +#include <fenv.h> +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double +#else +static long double +#endif +TWO63[2]={ + 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ + -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ +}; + +#ifdef __STDC__ + long double __nearbyintl(long double x) +#else + long double __nearbyintl(x) + long double x; +#endif +{ + fenv_t env; + int32_t se,j0,sx; + u_int32_t i,i0,i1; + long double w,t; + GET_LDOUBLE_WORDS(se,i0,i1,x); + sx = (se>>15)&1; + j0 = (se&0x7fff)-0x3fff; + if(j0<31) { + if(j0<0) { + if(((se&0x7fff)|i0|i1)==0) return x; + i1 |= i0; + i0 &= 0xe0000000; + i0 |= (i1|-i1)&0x80000000; + SET_LDOUBLE_MSW(x,i0); + feholdexcept (&env); + w = TWO63[sx]+x; + t = w-TWO63[sx]; + fesetenv (&env); + GET_LDOUBLE_EXP(i0,t); + SET_LDOUBLE_EXP(t,(i0&0x7fff)|(sx<<15)); + return t; + } else { + i = (0x7fffffff)>>j0; + if(((i0&i)|i1)==0) return x; /* x is integral */ + i>>=1; + if(((i0&i)|i1)!=0) { + if (j0==30) i1 = 0x40000000; else + i0 = (i0&(~i))|((0x20000000)>>j0); + } + } + } else if (j0>62) { + if(j0==0x4000) return x+x; /* inf or NaN */ + else return x; /* x is integral */ + } else { + i = ((u_int32_t)(0xffffffff))>>(j0-31); + if((i1&i)==0) return x; /* x is integral */ + i>>=1; + if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-31)); + } + SET_LDOUBLE_WORDS(x,se,i0,i1); + feholdexcept (&env); + w = TWO63[sx]+x; + t = w-TWO63[sx]; + fesetenv (&env); + return t; +} +weak_alias (__nearbyintl, nearbyintl) diff --git a/sysdeps/ieee754/ldbl-96/s_nextafterl.c b/sysdeps/ieee754/ldbl-96/s_nextafterl.c new file mode 100644 index 0000000000..aea57e3086 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_nextafterl.c @@ -0,0 +1,101 @@ +/* s_nextafterl.c -- long double version of s_nextafter.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* IEEE functions + * nextafterl(x,y) + * return the next machine floating-point number of x in the + * direction toward y. + * Special cases: + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + long double __nextafterl(long double x, long double y) +#else + long double __nextafterl(x,y) + long double x,y; +#endif +{ + int32_t hx,hy,ix,iy; + u_int32_t lx,ly,esx,esy; + + GET_LDOUBLE_WORDS(esx,hx,lx,x); + GET_LDOUBLE_WORDS(esy,hy,ly,y); + ix = esx&0x7fff; /* |x| */ + iy = esy&0x7fff; /* |y| */ + + if(((ix==0x7fff)&&((hx|lx)|-(hx|lx))!=0) || /* x is nan */ + ((iy==0x7fff)&&((hy|ly)|-(hy|ly))!=0)) /* y is nan */ + return x+y; + if(x==y) return y; /* x=y, return y */ + if((ix|hx|lx)==0) { /* x == 0 */ + SET_LDOUBLE_WORDS(x,esx&0x8000,0,1);/* return +-minsubnormal */ + y = x*x; + if(y==x) return y; else return x; /* raise underflow flag */ + } + if(esx<0x8000) { /* x > 0 */ + if(ix>iy||((ix==iy) && (hx>hy||((hx==hy)&&(lx>ly))))) { + /* x > y, x -= ulp */ + if(lx==0) { + if (hx==0) esx -= 1; + hx -= 1; + } + lx -= 1; + } else { /* x < y, x += ulp */ + lx += 1; + if(lx==0) { + hx += 1; + if (hx==0) + esx += 1; + } + } + } else { /* x < 0 */ + if(esy>=0||(ix>iy||((ix==iy)&&(hx>hy||((hx==hy)&&(lx>ly)))))){ + /* x < y, x -= ulp */ + if(lx==0) { + if (hx==0) esx -= 1; + hx -= 1; + } + lx -= 1; + } else { /* x > y, x += ulp */ + lx += 1; + if(lx==0) { + hx += 1; + if (hx==0) esx += 1; + } + } + } + esy = esx&0x7fff; + if(esy==0x7fff) return x+x; /* overflow */ + if(esy==0) { /* underflow */ + y = x*x; + if(y!=x) { /* raise underflow flag */ + SET_LDOUBLE_WORDS(y,esx,hx,lx); + return y; + } + } + SET_LDOUBLE_WORDS(x,esx,hx,lx); + return x; +} +weak_alias (__nextafterl, nextafterl) +strong_alias (__nextafterl, __nexttowardl) +weak_alias (__nextafterl, nexttowardl) diff --git a/sysdeps/ieee754/ldbl-96/s_nexttoward.c b/sysdeps/ieee754/ldbl-96/s_nexttoward.c new file mode 100644 index 0000000000..debbb86a3c --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_nexttoward.c @@ -0,0 +1,98 @@ +/* s_nexttoward.c + * Conversion from s_nextafter.c by Ulrich Drepper, Cygnus Support, + * drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* IEEE functions + * nextafterx(x,y) + * return the next machine floating-point number of x in the + * direction toward y. + * Special cases: + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + double __nexttoward(double x, long double y) +#else + double __nexttoward(x,y) + double x; + long double y; +#endif +{ + int32_t hx,ix,iy; + u_int32_t lx,hy,ly,esy; + + EXTRACT_WORDS(hx,lx,x); + GET_LDOUBLE_WORDS(esy,hy,ly,y); + ix = hx&0x7fffffff; /* |x| */ + iy = esy&0x7fff; /* |y| */ + + if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */ + ((iy>=0x7fff)&&((hy|ly)|-(hy|ly))!=0)) /* y is nan */ + return x+y; + if((long double) x==y) return x; /* x=y, return x */ + if((ix|lx)==0) { /* x == 0 */ + double x2; + INSERT_WORDS(x,esy&0x8000?0x80000000:0,1);/* return +-minsub */ + x2 = x*x; + if(x2==x) return x2; else return x; /* raise underflow flag */ + } + if(hx>=0) { /* x > 0 */ + if (esy>=0x8000||((ix>>20)&0x7ff)>iy-0x3c00 + || (((ix>>20)&0x7ff)==iy-0x3c00 + && (((hx<<11)|(lx>>21))>(hy&0x7fffffff) + || (((hx<<11)|(lx>>21))==(hy&0x7fffffff) + && (lx<<11)>ly)))) { /* x > y, x -= ulp */ + if(lx==0) hx -= 1; + lx -= 1; + } else { /* x < y, x += ulp */ + lx += 1; + if(lx==0) hx += 1; + } + } else { /* x < 0 */ + if (esy<0x8000||((ix>>20)&0x7ff)>iy-0x3c00 + || (((ix>>20)&0x7ff)==iy-0x3c00 + && (((hx<<11)|(lx>>21))>(hy&0x7fffffff) + || (((hx<<11)|(lx>>21))==(hy&0x7fffffff) + && (lx<<11)>ly)))) {/* x < y, x -= ulp */ + if(lx==0) hx -= 1; + lx -= 1; + } else { /* x > y, x += ulp */ + lx += 1; + if(lx==0) hx += 1; + } + } + hy = hx&0x7ff00000; + if(hy>=0x7ff00000) return x+x; /* overflow */ + if(hy<0x00100000) { /* underflow */ + double x2 = x*x; + if(x2!=x) { /* raise underflow flag */ + INSERT_WORDS(x2,hx,lx); + return x2; + } + } + INSERT_WORDS(x,hx,lx); + return x; +} +weak_alias (__nexttoward, nexttoward) +#ifdef NO_LONG_DOUBLE +strong_alias (__nexttoward, __nexttowardl) +weak_alias (__nexttoward, nexttowardl) +#endif diff --git a/sysdeps/ieee754/ldbl-96/s_nexttowardf.c b/sysdeps/ieee754/ldbl-96/s_nexttowardf.c new file mode 100644 index 0000000000..b7e9f00185 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_nexttowardf.c @@ -0,0 +1,78 @@ +/* s_nexttowardf.c -- float version of s_nextafter.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + float __nexttowardf(float x, long double y) +#else + float __nexttowardf(x,y) + float x; + long double y; +#endif +{ + int32_t hx,ix,iy; + u_int32_t hy,ly,esy; + + GET_FLOAT_WORD(hx,x); + GET_LDOUBLE_WORDS(esy,hy,ly,y); + ix = hx&0x7fffffff; /* |x| */ + iy = esy&0x7fff; /* |y| */ + + if((ix>0x7f800000) || /* x is nan */ + (iy>=0x7fff&&((hy|ly)|-(hy|ly))!=0)) /* y is nan */ + return x+y; + if((long double) x==y) return y; /* x=y, return y */ + if(ix==0) { /* x == 0 */ + float x2; + SET_FLOAT_WORD(x,(esy&0x8000?0x80000000:0)|1);/* return +-minsub*/ + x2 = x*x; + if(x2==x) return x2; else return x; /* raise underflow flag */ + } + if(hx>=0) { /* x > 0 */ + if(esy>=0x8000||((ix>>23)&0xff)>iy-0x3f80 + || (((ix>>23)&0xff)==iy-0x3f80 + && ((ix&0x7fffff)<<8)>(hy&0x7fffffff))) {/* x > y, x -= ulp */ + hx -= 1; + } else { /* x < y, x += ulp */ + hx += 1; + } + } else { /* x < 0 */ + if(esy<0x8000||((ix>>23)&0xff)>iy-0x3f80 + || (((ix>>23)&0xff)==iy-0x3f80 + && ((ix&0x7fffff)<<8)>(hy&0x7fffffff))) {/* x < y, x -= ulp */ + hx -= 1; + } else { /* x > y, x += ulp */ + hx += 1; + } + } + hy = hx&0x7f800000; + if(hy>=0x7f800000) return x+x; /* overflow */ + if(hy<0x00800000) { /* underflow */ + float x2 = x*x; + if(x2!=x) { /* raise underflow flag */ + SET_FLOAT_WORD(x2,hx); + return x2; + } + } + SET_FLOAT_WORD(x,hx); + return x; +} +weak_alias (__nexttowardf, nexttowardf) diff --git a/sysdeps/ieee754/ldbl-96/s_remquol.c b/sysdeps/ieee754/ldbl-96/s_remquol.c new file mode 100644 index 0000000000..88ff298eb6 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_remquol.c @@ -0,0 +1,109 @@ +/* Compute remainder and a congruent to the quotient. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +static const long double zero = 0.0; + + +long double +__remquol (long double x, long double p, int *quo) +{ + int32_t ex,ep,hx,hp; + u_int32_t sx,lx,lp; + int cquo,qs; + + GET_LDOUBLE_WORDS (ex, hx, lx, x); + GET_LDOUBLE_WORDS (ep, hp, lp, p); + sx = ex & 0x8000; + qs = (sx ^ (ep & 0x8000)) >> 15; + ep &= 0x7fff; + ex &= 0x7fff; + + /* Purge off exception values. */ + if ((ep | hp | lp) == 0) + return (x * p) / (x * p); /* p = 0 */ + if ((ex == 0x7fff) /* x not finite */ + || ((ep == 0x7fff) /* p is NaN */ + && ((hp | lp) != 0))) + return (x * p) / (x * p); + + if (ep <= 0x7ffb) + x = __ieee754_fmodl (x, 8 * p); /* now x < 8p */ + + if (((ex - ep) | (hx - hp) | (lx - lp)) == 0) + { + *quo = qs ? -1 : 1; + return zero * x; + } + + x = fabsl (x); + p = fabsl (p); + cquo = 0; + + if (x >= 4 * p) + { + x -= 4 * p; + cquo += 4; + } + if (x >= 2 * p) + { + x -= 2 * p; + cquo += 2; + } + + if (ep < 0x0002) + { + if (x + x > p) + { + x -= p; + ++cquo; + if (x + x >= p) + { + x -= p; + ++cquo; + } + } + } + else + { + long double p_half = 0.5 * p; + if (x > p_half) + { + x -= p; + ++cquo; + if (x >= p_half) + { + x -= p; + ++cquo; + } + } + } + + *quo = qs ? -cquo : cquo; + + if (sx) + x = -x; + return x; +} +weak_alias (__remquol, remquol) diff --git a/sysdeps/ieee754/ldbl-96/s_rintl.c b/sysdeps/ieee754/ldbl-96/s_rintl.c new file mode 100644 index 0000000000..9d4777dcc4 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_rintl.c @@ -0,0 +1,91 @@ +/* s_rintl.c -- long double version of s_rint.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* + * rintl(x) + * Return x rounded to integral value according to the prevailing + * rounding mode. + * Method: + * Using floating addition. + * Exception: + * Inexact flag raised if x not equal to rintl(x). + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double +#else +static long double +#endif +TWO63[2]={ + 9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */ + -9.223372036854775808000000e+18 /* 0xC03E, 0x00000000, 0x00000000 */ +}; + +#ifdef __STDC__ + long double __rintl(long double x) +#else + long double __rintl(x) + long double x; +#endif +{ + int32_t se,j0,sx; + u_int32_t i,i0,i1; + long double w,t; + GET_LDOUBLE_WORDS(se,i0,i1,x); + sx = (se>>15)&1; + j0 = (se&0x7fff)-0x3fff; + if(j0<31) { + if(j0<0) { + if(((se&0x7fff)|i0|i1)==0) return x; + i1 |= i0; + i0 &= 0xe0000000; + i0 |= (i1|-i1)&0x80000000; + SET_LDOUBLE_MSW(x,i0); + w = TWO63[sx]+x; + t = w-TWO63[sx]; + GET_LDOUBLE_EXP(i0,t); + SET_LDOUBLE_EXP(t,(i0&0x7fff)|(sx<<15)); + return t; + } else { + i = (0x7fffffff)>>j0; + if(((i0&i)|i1)==0) return x; /* x is integral */ + i>>=1; + if(((i0&i)|i1)!=0) { + if(j0==30) i1 = 0x40000000; else + i0 = (i0&(~i))|((0x20000000)>>j0); + } + } + } else if (j0>62) { + if(j0==0x4000) return x+x; /* inf or NaN */ + else return x; /* x is integral */ + } else { + i = ((u_int32_t)(0xffffffff))>>(j0-31); + if((i1&i)==0) return x; /* x is integral */ + i>>=1; + if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-31)); + } + SET_LDOUBLE_WORDS(x,se,i0,i1); + w = TWO63[sx]+x; + return w-TWO63[sx]; +} +weak_alias (__rintl, rintl) diff --git a/sysdeps/ieee754/ldbl-96/s_roundl.c b/sysdeps/ieee754/ldbl-96/s_roundl.c new file mode 100644 index 0000000000..d7482b9b7c --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_roundl.c @@ -0,0 +1,106 @@ +/* Round long double to integer away from zero. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +static const long double huge = 1.0e4930; + + +long double +__roundl (long double x) +{ + int32_t j0; + u_int32_t se, i1, i0; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + j0 = (se & 0x7fff) - 0x3fff; + if (j0 < 31) + { + if (j0 < 0) + { + if (huge + x > 0.0) + { + se &= 0x8000; + i0 = i1 = 0; + if (j0 == -1) + { + se |= 0x3fff; + i0 = 0x80000000; + } + } + } + else + { + u_int32_t i = 0x7fffffff >> j0; + if (((i0 & i) | i1) == 0) + /* X is integral. */ + return x; + if (huge + x > 0.0) + { + /* Raise inexact if x != 0. */ + u_int32_t j = i0 + (0x40000000 >> j0); + if (j < i0) + se += 1; + i0 = (j & ~i) | 0x80000000; + i1 = 0; + } + } + } + else if (j0 > 62) + { + if (j0 == 0x4000) + /* Inf or NaN. */ + return x + x; + else + return x; + } + else + { + u_int32_t i = 0xffffffff >> (j0 - 31); + if ((i1 & i) == 0) + /* X is integral. */ + return x; + + if (huge + x > 0.0) + { + /* Raise inexact if x != 0. */ + u_int32_t j = i1 + (1 << (62 - j0)); + if (j < i1) + { + u_int32_t k = i0 + 1; + if (k < i0) + { + se += 1; + k |= 0x80000000; + } + i0 = k; + } + i1 = j; + } + i1 &= ~i; + } + + SET_LDOUBLE_WORDS (x, se, i0, i1); + return x; +} +weak_alias (__roundl, roundl) diff --git a/sysdeps/ieee754/ldbl-96/s_scalblnl.c b/sysdeps/ieee754/ldbl-96/s_scalblnl.c new file mode 100644 index 0000000000..8e556fabe1 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_scalblnl.c @@ -0,0 +1,71 @@ +/* s_scalbnl.c -- long double version of s_scalbn.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* + * scalbnl (long double x, int n) + * scalbnl(x,n) returns x* 2**n computed by exponent + * manipulation rather than by actually performing an + * exponentiation or a multiplication. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double +#else +static long double +#endif +two63 = 4.50359962737049600000e+15, +twom63 = 1.08420217248550443400e-19, +huge = 1.0e+4900L, +tiny = 1.0e-4900L; + +#ifdef __STDC__ + long double __scalblnl (long double x, long int n) +#else + long double __scalblnl (x,n) + long double x; long int n; +#endif +{ + int32_t k,es,hx,lx; + GET_LDOUBLE_WORDS(es,hx,lx,x); + k = es&0x7fff; /* extract exponent */ + if (k==0) { /* 0 or subnormal x */ + if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ + x *= two63; + GET_LDOUBLE_EXP(es,x); + k = (hx&0x7fff) - 63; + } + if (k==0x7fff) return x+x; /* NaN or Inf */ + k = k+n; + if (n> 50000 || k > 0x7ffe) + return huge*__copysignl(huge,x); /* overflow */ + if (n< -50000) + return tiny*__copysignl(tiny,x); + if (k > 0) /* normal result */ + {SET_LDOUBLE_EXP(x,(es&0x8000)|k); return x;} + if (k <= -63) + return tiny*__copysignl(tiny,x); /*underflow*/ + k += 63; /* subnormal result */ + SET_LDOUBLE_EXP(x,(es&0x8000)|k); + return x*twom63; +} +weak_alias (__scalblnl, scalblnl) diff --git a/sysdeps/ieee754/ldbl-96/s_scalbnl.c b/sysdeps/ieee754/ldbl-96/s_scalbnl.c new file mode 100644 index 0000000000..34c52e773a --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_scalbnl.c @@ -0,0 +1,71 @@ +/* s_scalbnl.c -- long double version of s_scalbn.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* + * scalbnl (long double x, int n) + * scalbnl(x,n) returns x* 2**n computed by exponent + * manipulation rather than by actually performing an + * exponentiation or a multiplication. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double +#else +static long double +#endif +two63 = 4.50359962737049600000e+15, +twom63 = 1.08420217248550443400e-19, +huge = 1.0e+4900L, +tiny = 1.0e-4900L; + +#ifdef __STDC__ + long double __scalbnl (long double x, int n) +#else + long double __scalbnl (x,n) + long double x; int n; +#endif +{ + int32_t k,es,hx,lx; + GET_LDOUBLE_WORDS(es,hx,lx,x); + k = es&0x7fff; /* extract exponent */ + if (k==0) { /* 0 or subnormal x */ + if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ + x *= two63; + GET_LDOUBLE_EXP(es,x); + k = (hx&0x7fff) - 63; + } + if (k==0x7fff) return x+x; /* NaN or Inf */ + k = k+n; + if (n> 50000 || k > 0x7ffe) + return huge*__copysignl(huge,x); /* overflow */ + if (n< -50000) + return tiny*__copysignl(tiny,x); + if (k > 0) /* normal result */ + {SET_LDOUBLE_EXP(x,(es&0x8000)|k); return x;} + if (k <= -63) + return tiny*__copysignl(tiny,x); /*underflow*/ + k += 63; /* subnormal result */ + SET_LDOUBLE_EXP(x,(es&0x8000)|k); + return x*twom63; +} +weak_alias (__scalbnl, scalbnl) diff --git a/sysdeps/ieee754/ldbl-96/s_signbitl.c b/sysdeps/ieee754/ldbl-96/s_signbitl.c new file mode 100644 index 0000000000..b12fdefff4 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_signbitl.c @@ -0,0 +1,32 @@ +/* Return nonzero value if number is negative. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + +int +__signbitl (long double x) +{ + int32_t e; + + GET_LDOUBLE_EXP (e, x); + return e & 0x8000; +} diff --git a/sysdeps/ieee754/ldbl-96/s_sincosl.c b/sysdeps/ieee754/ldbl-96/s_sincosl.c new file mode 100644 index 0000000000..78c78d5619 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_sincosl.c @@ -0,0 +1,74 @@ +/* Compute sine and cosine of argument. + Copyright (C) 1997 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +void +__sincosl (long double x, long double *sinx, long double *cosx) +{ + int32_t se, i0, i1; + + /* High word of x. */ + GET_LDOUBLE_WORDS (se, i0, i1, x); + + /* |x| ~< pi/4 */ + se &= 0x7fff; + if (se < 0x3ffe || (se == 0x3ffe && i0 <= 0xc90fdaa2)) + { + *sinx = __kernel_sinl (x, 0.0, 0); + *cosx = __kernel_cosl (x, 0.0); + } + else if (se == 0x7fff) + { + /* sin(Inf or NaN) is NaN */ + *sinx = *cosx = x - x; + } + else + { + /* Argument reduction needed. */ + long double y[2]; + int n; + + n = __ieee754_rem_pio2l (x, y); + switch (n & 3) + { + case 0: + *sinx = __kernel_sinl (y[0], y[1], 1); + *cosx = __kernel_cosl (y[0], y[1]); + break; + case 1: + *sinx = __kernel_cosl (y[0], y[1]); + *cosx = -__kernel_sinl (y[0], y[1], 1); + break; + case 2: + *sinx = -__kernel_sinl (y[0], y[1], 1); + *cosx = -__kernel_cosl (y[0], y[1]); + break; + default: + *sinx = -__kernel_cosl (y[0], y[1]); + *cosx = __kernel_sinl (y[0], y[1], 1); + break; + } + } +} +weak_alias (__sincosl, sincosl) diff --git a/sysdeps/ieee754/ldbl-96/s_sinl.c b/sysdeps/ieee754/ldbl-96/s_sinl.c new file mode 100644 index 0000000000..4fd48805b4 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_sinl.c @@ -0,0 +1,88 @@ +/* s_sinl.c -- long double version of s_sin.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* sinl(x) + * Return sine function of x. + * + * kernel function: + * __kernel_sinl ... sine function on [-pi/4,pi/4] + * __kernel_cosl ... cose function on [-pi/4,pi/4] + * __ieee754_rem_pio2l ... argument reduction routine + * + * Method. + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * in [-pi/4 , +pi/4], and let n = k mod 4. + * We have + * + * n sin(x) cos(x) tan(x) + * ---------------------------------------------------------- + * 0 S C T + * 1 C -S -1/T + * 2 -S -C T + * 3 -C S -1/T + * ---------------------------------------------------------- + * + * Special cases: + * Let trig be any of sin, cos, or tan. + * trig(+-INF) is NaN, with signals; + * trig(NaN) is that NaN; + * + * Accuracy: + * TRIG(x) returns trig(x) nearly rounded + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + long double __sinl(long double x) +#else + long double __sinl(x) + long double x; +#endif +{ + long double y[2],z=0.0; + int32_t n, se, i0, i1; + + /* High word of x. */ + GET_LDOUBLE_WORDS(se,i0,i1,x); + + /* |x| ~< pi/4 */ + se &= 0x7fff; + if(se < 0x3ffe || (se == 0x3ffe && i0 <= 0xc90fdaa2)) + return __kernel_sinl(x,z,0); + + /* sin(Inf or NaN) is NaN */ + else if (se==0x7fff) return x-x; + + /* argument reduction needed */ + else { + n = __ieee754_rem_pio2l(x,y); + switch(n&3) { + case 0: return __kernel_sinl(y[0],y[1],1); + case 1: return __kernel_cosl(y[0],y[1]); + case 2: return -__kernel_sinl(y[0],y[1],1); + default: + return -__kernel_cosl(y[0],y[1]); + } + } +} +weak_alias (__sinl, sinl) diff --git a/sysdeps/ieee754/ldbl-96/s_tanhl.c b/sysdeps/ieee754/ldbl-96/s_tanhl.c new file mode 100644 index 0000000000..1e3dc3b613 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_tanhl.c @@ -0,0 +1,95 @@ +/* s_tanhl.c -- long double version of s_tanh.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* tanhl(x) + * Return the Hyperbolic Tangent of x + * + * Method : + * x -x + * e - e + * 0. tanhl(x) is defined to be ----------- + * x -x + * e + e + * 1. reduce x to non-negative by tanhl(-x) = -tanhl(x). + * 2. 0 <= x <= 2**-55 : tanhl(x) := x*(one+x) + * -t + * 2**-55 < x <= 1 : tanhl(x) := -----; t = expm1l(-2x) + * t + 2 + * 2 + * 1 <= x <= 23.0 : tanhl(x) := 1- ----- ; t=expm1l(2x) + * t + 2 + * 23.0 < x <= INF : tanhl(x) := 1. + * + * Special cases: + * tanhl(NaN) is NaN; + * only tanhl(0)=0 is exact for finite argument. + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double one=1.0, two=2.0, tiny = 1.0e-4900L; +#else +static long double one=1.0, two=2.0, tiny = 1.0e-4900L; +#endif + +#ifdef __STDC__ + long double __tanhl(long double x) +#else + long double __tanhl(x) + long double x; +#endif +{ + long double t,z; + int32_t se; + u_int32_t j0,j1,ix; + + /* High word of |x|. */ + GET_LDOUBLE_WORDS(se,j0,j1,x); + ix = se&0x7fff; + + /* x is INF or NaN */ + if(ix==0x7fff) { + /* for NaN it's not important which branch: tanhl(NaN) = NaN */ + if (se&0x8000) return one/x-one; /* tanhl(-inf)= -1; */ + else return one/x+one; /* tanhl(+inf)=+1 */ + } + + /* |x| < 23 */ + if (ix < 0x4003 || (ix == 0x4003 && j0 < 0xb8000000u)) {/* |x|<23 */ + if ((ix|j0|j1) == 0) + return x; /* x == +- 0 */ + if (ix<0x3fc8) /* |x|<2**-55 */ + return x*(one+x); /* tanh(small) = small */ + if (ix>=0x3fff) { /* |x|>=1 */ + t = __expm1l(two*fabsl(x)); + z = one - two/(t+two); + } else { + t = __expm1l(-two*fabsl(x)); + z= -t/(t+two); + } + /* |x| > 23, return +-1 */ + } else { + z = one - tiny; /* raised inexact flag */ + } + return (se>0x7fff)? -z: z; +} +weak_alias (__tanhl, tanhl) diff --git a/sysdeps/ieee754/ldbl-96/s_tanl.c b/sysdeps/ieee754/ldbl-96/s_tanl.c new file mode 100644 index 0000000000..97a0b27f32 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_tanl.c @@ -0,0 +1,81 @@ +/* s_tanl.c -- long double version of s_tan.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* tanl(x) + * Return tangent function of x. + * + * kernel function: + * __kernel_tanl ... tangent function on [-pi/4,pi/4] + * __ieee754_rem_pio2l ... argument reduction routine + * + * Method. + * Let S,C and T denote the sin, cos and tan respectively on + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * in [-pi/4 , +pi/4], and let n = k mod 4. + * We have + * + * n sin(x) cos(x) tan(x) + * ---------------------------------------------------------- + * 0 S C T + * 1 C -S -1/T + * 2 -S -C T + * 3 -C S -1/T + * ---------------------------------------------------------- + * + * Special cases: + * Let trig be any of sin, cos, or tan. + * trig(+-INF) is NaN, with signals; + * trig(NaN) is that NaN; + * + * Accuracy: + * TRIG(x) returns trig(x) nearly rounded + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ + long double __tanl(long double x) +#else + long double __tanl(x) + long double x; +#endif +{ + long double y[2],z=0.0; + int32_t n, se; + + /* High word of x. */ + GET_LDOUBLE_EXP(se,x); + + /* |x| ~< pi/4 */ + se &= 0x7fff; + if(se <= 0x3ffe) return __kernel_tanl(x,z,1); + + /* tan(Inf or NaN) is NaN */ + else if (se==0x7fff) return x-x; /* NaN */ + + /* argument reduction needed */ + else { + n = __ieee754_rem_pio2l(x,y); + return __kernel_tanl(y[0],y[1],1-((n&1)<<1)); /* 1 -- n even + -1 -- n odd */ + } +} +weak_alias (__tanl, tanl) diff --git a/sysdeps/ieee754/ldbl-96/s_truncl.c b/sysdeps/ieee754/ldbl-96/s_truncl.c new file mode 100644 index 0000000000..59c3b9c173 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/s_truncl.c @@ -0,0 +1,57 @@ +/* Truncate argument to nearest integral value not larger than the argument. + Copyright (C) 1997, 1999 Free Software Foundation, Inc. + This file is part of the GNU C Library. + Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +#include "math_private.h" + + +long double +__truncl (long double x) +{ + int32_t i0, j0; + u_int32_t se, i1; + int sx; + + GET_LDOUBLE_WORDS (se, i0, i1, x); + sx = se & 0x8000; + j0 = (se & 0x7fff) - 0x3fff; + if (j0 < 31) + { + if (j0 < 0) + /* The magnitude of the number is < 1 so the result is +-0. */ + SET_LDOUBLE_WORDS (x, sx, 0, 0); + else + SET_LDOUBLE_WORDS (x, se, i0 & ~(0x7fffffff >> j0), 0); + } + else if (j0 > 63) + { + if (j0 == 0x4000) + /* x is inf or NaN. */ + return x + x; + } + else + { + SET_LDOUBLE_WORDS (x, se, i0, i1 & ~(0xffffffffu >> (j0 - 31))); + } + + return x; +} +weak_alias (__truncl, truncl) diff --git a/sysdeps/ieee754/ldbl-96/strtold.c b/sysdeps/ieee754/ldbl-96/strtold.c new file mode 100644 index 0000000000..cabb787d4e --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/strtold.c @@ -0,0 +1,42 @@ +/* Copyright (C) 1999 Free Software Foundation, Inc. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Library General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Library General Public License for more details. + + You should have received a copy of the GNU Library General Public + License along with the GNU C Library; see the file COPYING.LIB. If not, + write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + Boston, MA 02111-1307, USA. */ + +#include <math.h> + +/* The actual implementation for all floating point sizes is in strtod.c. + These macros tell it to produce the `long double' version, `strtold'. */ + +# define FLOAT long double +# define FLT LDBL +# ifdef USE_IN_EXTENDED_LOCALE_MODEL +# define STRTOF __strtold_l +# else +# define STRTOF strtold +# endif +# define MPN2FLOAT __mpn_construct_long_double +# define FLOAT_HUGE_VAL HUGE_VALL +# define SET_MANTISSA(flt, mant) \ + do { union ieee854_long_double u; \ + u.d = (flt); \ + if ((mant & 0x7fffffffffffffffULL) == 0) \ + mant = 0x4000000000000000ULL; \ + u.ieee.mantissa0 = (((mant) >> 32) & 0x7fffffff) | 0x80000000; \ + u.ieee.mantissa1 = (mant) & 0xffffffff; \ + (flt) = u.d; \ + } while (0) + +# include "strtod.c" diff --git a/sysdeps/ieee754/ldbl-96/w_expl.c b/sysdeps/ieee754/ldbl-96/w_expl.c new file mode 100644 index 0000000000..b8152cea65 --- /dev/null +++ b/sysdeps/ieee754/ldbl-96/w_expl.c @@ -0,0 +1,60 @@ +/* w_expl.c -- long double version of w_exp.c. + * Conversion to long double by Ulrich Drepper, + * Cygnus Support, drepper@cygnus.com. + */ + +/* + * ==================================================== + * 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: $"; +#endif + +/* + * wrapper expl(x) + */ + +#include "math.h" +#include "math_private.h" + +#ifdef __STDC__ +static const long double +#else +static long double +#endif +o_threshold= 1.135652340629414394949193107797076489134e4, + /* 0x400C, 0xB17217F7, 0xD1CF79AC */ +u_threshold= -1.140019167866942050398521670162263001513e4; + /* 0x400C, 0xB220C447, 0x69C201E8 */ + +#ifdef __STDC__ + long double __expl(long double x) /* wrapper exp */ +#else + long double __expl(x) /* wrapper exp */ + long double x; +#endif +{ +#ifdef _IEEE_LIBM + return __ieee754_expl(x); +#else + long double z; + z = __ieee754_expl(x); + if(_LIB_VERSION == _IEEE_) return z; + if(__finitel(x)) { + if(x>o_threshold) + return __kernel_standard(x,x,206); /* exp overflow */ + else if(x<u_threshold) + return __kernel_standard(x,x,207); /* exp underflow */ + } + return z; +#endif +} +weak_alias (__expl, expl) |