From c4e44e97f8562254d9da898f6ed7e6cb4d8a3ce4 Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Sun, 6 Mar 2005 07:11:53 +0000 Subject: Trim off whitespace --- libm/Makefile | 10 ++-- libm/README | 4 +- libm/e_acos.c | 10 ++-- libm/e_acosh.c | 12 ++-- libm/e_asin.c | 20 +++---- libm/e_atan2.c | 20 +++---- libm/e_atanh.c | 8 +-- libm/e_cosh.c | 18 +++--- libm/e_exp.c | 30 +++++----- libm/e_fmod.c | 12 ++-- libm/e_gamma.c | 2 +- libm/e_gamma_r.c | 6 +- libm/e_hypot.c | 22 ++++---- libm/e_j0.c | 26 ++++----- libm/e_j1.c | 30 +++++----- libm/e_jn.c | 42 +++++++------- libm/e_lgamma.c | 2 +- libm/e_lgamma_r.c | 32 +++++------ libm/e_log.c | 38 ++++++------- libm/e_log10.c | 10 ++-- libm/e_pow.c | 42 +++++++------- libm/e_rem_pio2.c | 50 ++++++++--------- libm/e_remainder.c | 10 ++-- libm/e_scalb.c | 4 +- libm/e_sinh.c | 14 ++--- libm/e_sqrt.c | 94 +++++++++++++++---------------- libm/fp_private.h | 4 +- libm/fpmacros.c | 92 +++++++++++++++---------------- libm/k_cos.c | 20 +++---- libm/k_rem_pio2.c | 60 ++++++++++---------- libm/k_sin.c | 18 +++--- libm/k_standard.c | 34 ++++++------ libm/k_tan.c | 20 +++---- libm/math_private.h | 40 +++++++------- libm/powerpc/Makefile | 8 +-- libm/powerpc/s_ceil.c | 16 +++--- libm/powerpc/s_copysign.c | 6 +- libm/powerpc/s_floor.c | 16 +++--- libm/powerpc/s_frexp.c | 8 +-- libm/powerpc/s_ldexp.c | 12 ++-- libm/powerpc/s_logb.c | 26 ++++----- libm/powerpc/s_modf.c | 138 +++++++++++++++++++++++----------------------- libm/powerpc/s_rint.c | 44 +++++++-------- libm/powerpc/w_scalb.c | 32 +++++------ libm/s_asinh.c | 16 +++--- libm/s_atan.c | 16 +++--- libm/s_cbrt.c | 12 ++-- libm/s_ceil.c | 6 +- libm/s_ceilf.c | 4 +- libm/s_copysign.c | 2 +- libm/s_cos.c | 8 +-- libm/s_erf.c | 34 ++++++------ libm/s_expm1.c | 48 ++++++++-------- libm/s_fabs.c | 2 +- libm/s_finite.c | 2 +- libm/s_floor.c | 6 +- libm/s_floorf.c | 4 +- libm/s_frexp.c | 8 +-- libm/s_ilogb.c | 4 +- libm/s_ldexp.c | 2 +- libm/s_lib_version.c | 2 +- libm/s_log1p.c | 40 +++++++------- libm/s_logb.c | 6 +- libm/s_matherr.c | 2 +- libm/s_modf.c | 4 +- libm/s_nextafter.c | 10 ++-- libm/s_rint.c | 6 +- libm/s_scalbn.c | 14 ++--- libm/s_significand.c | 2 +- libm/s_sin.c | 8 +-- libm/s_tan.c | 8 +-- libm/s_tanh.c | 4 +- libm/w_acos.c | 2 +- libm/w_acosh.c | 4 +- libm/w_asin.c | 4 +- libm/w_atan2.c | 4 +- libm/w_atanh.c | 6 +- libm/w_cabs.c | 2 +- libm/w_cosh.c | 6 +- libm/w_drem.c | 2 +- libm/w_exp.c | 6 +- libm/w_fmod.c | 4 +- libm/w_gamma.c | 4 +- libm/w_gamma_r.c | 6 +- libm/w_hypot.c | 2 +- libm/w_j0.c | 2 +- libm/w_j1.c | 6 +- libm/w_jn.c | 6 +- libm/w_lgamma.c | 4 +- libm/w_lgamma_r.c | 6 +- libm/w_log.c | 4 +- libm/w_log10.c | 6 +- libm/w_pow.c | 14 ++--- libm/w_remainder.c | 6 +- libm/w_scalb.c | 8 +-- libm/w_sinh.c | 6 +- libm/w_sqrt.c | 4 +- libm/w_sqrtf.c | 2 +- 98 files changed, 784 insertions(+), 784 deletions(-) diff --git a/libm/Makefile b/libm/Makefile index 7a0edbea5..e8aaddc12 100644 --- a/libm/Makefile +++ b/libm/Makefile @@ -6,14 +6,14 @@ # math library for Apple's MacOS X/Darwin math library, which was # itself swiped from FreeBSD. The original copyright information # is as follows: -# +# # 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. -# +# # It has been ported to work with uClibc and generally behave # by Erik Andersen # @@ -36,7 +36,7 @@ include $(TOPDIR)Rules.mak CFLAGS+=$(SSP_ALL_CFLAGS) -DIRS = +DIRS = ifeq ($(strip $(HAS_FPU)),y) ifeq ($(TARGET_ARCH),$(wildcard $(TARGET_ARCH))) DIRS = $(TARGET_ARCH) @@ -71,7 +71,7 @@ else CSRC = w_acos.c w_asin.c s_atan.c w_atan2.c s_ceil.c s_cos.c \ w_cosh.c w_exp.c s_fabs.c s_floor.c w_fmod.c s_frexp.c \ s_ldexp.c w_log.c w_log10.c s_modf.c w_pow.c s_sin.c \ - w_sinh.c w_sqrt.c s_tan.c s_tanh.c + w_sinh.c w_sqrt.c s_tan.c s_tanh.c CSRC+= s_expm1.c s_scalbn.c s_copysign.c e_acos.c e_asin.c e_atan2.c \ k_cos.c e_cosh.c e_exp.c e_fmod.c e_log.c e_log10.c e_pow.c \ k_sin.c e_sinh.c e_sqrt.c k_tan.c e_rem_pio2.c k_rem_pio2.c \ diff --git a/libm/README b/libm/README index c275d1b9a..2c69a54af 100644 --- a/libm/README +++ b/libm/README @@ -7,10 +7,10 @@ is as follows: Developed at SunPro, a Sun Microsystems, Inc. business. Permission to use, copy, modify, and distribute this - software is freely granted, provided that this notice + software is freely granted, provided that this notice is preserved. -It has been ported to work with uClibc and generally behave +It has been ported to work with uClibc and generally behave by Erik Andersen 22 May, 2001 diff --git a/libm/e_acos.c b/libm/e_acos.c index 78bdae9f8..9f6305e1e 100644 --- a/libm/e_acos.c +++ b/libm/e_acos.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -15,14 +15,14 @@ static char rcsid[] = "$NetBSD: e_acos.c,v 1.9 1995/05/12 04:57:13 jtc Exp $"; #endif /* __ieee754_acos(x) - * Method : + * Method : * acos(x) = pi/2 - asin(x) * acos(-x) = pi/2 + asin(x) * For |x|<=0.5 * acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c) * For x>0.5 * acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2))) - * = 2asin(sqrt((1-x)/2)) + * = 2asin(sqrt((1-x)/2)) * = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z) * = 2f + (2c + 2s*z*R(z)) * where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term @@ -42,9 +42,9 @@ static char rcsid[] = "$NetBSD: e_acos.c,v 1.9 1995/05/12 04:57:13 jtc Exp $"; #include "math_private.h" #ifdef __STDC__ -static const double +static const double #else -static double +static double #endif one= 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ diff --git a/libm/e_acosh.c b/libm/e_acosh.c index 8383519df..3b3c6251e 100644 --- a/libm/e_acosh.c +++ b/libm/e_acosh.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -16,7 +16,7 @@ static char rcsid[] = "$NetBSD: e_acosh.c,v 1.9 1995/05/12 04:57:18 jtc Exp $"; /* __ieee754_acosh(x) * Method : - * Based on + * Based on * acosh(x) = log [ x + sqrt(x*x-1) ] * we have * acosh(x) := log(x)+ln2, if x is large; else @@ -32,9 +32,9 @@ static char rcsid[] = "$NetBSD: e_acosh.c,v 1.9 1995/05/12 04:57:18 jtc Exp $"; #include "math_private.h" #ifdef __STDC__ -static const double +static const double #else -static double +static double #endif one = 1.0, ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ @@ -45,7 +45,7 @@ ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ double __ieee754_acosh(x) double x; #endif -{ +{ double t; int32_t hx; u_int32_t lx; @@ -55,7 +55,7 @@ ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ } else if(hx >=0x41b00000) { /* x > 2**28 */ if(hx >=0x7ff00000) { /* x is inf of NaN */ return x+x; - } else + } else return __ieee754_log(x)+ln2; /* acosh(huge)=log(2x) */ } else if(((hx-0x3ff00000)|lx)==0) { return 0.0; /* acosh(1) = 0 */ diff --git a/libm/e_asin.c b/libm/e_asin.c index 210fdbf3b..d7a785350 100644 --- a/libm/e_asin.c +++ b/libm/e_asin.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -15,12 +15,12 @@ static char rcsid[] = "$NetBSD: e_asin.c,v 1.9 1995/05/12 04:57:22 jtc Exp $"; #endif /* __ieee754_asin(x) - * Method : + * Method : * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ... * we approximate asin(x) on [0,0.5] by * asin(x) = x + x*x^2*R(x^2) * where - * R(x^2) is a rational approximation of (asin(x)-x)/x^3 + * R(x^2) is a rational approximation of (asin(x)-x)/x^3 * and its remez error is bounded by * |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75) * @@ -49,9 +49,9 @@ static char rcsid[] = "$NetBSD: e_asin.c,v 1.9 1995/05/12 04:57:22 jtc Exp $"; #include "math_private.h" #ifdef __STDC__ -static const double +static const double #else -static double +static double #endif one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ huge = 1.000e+300, @@ -86,12 +86,12 @@ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ GET_LOW_WORD(lx,x); if(((ix-0x3ff00000)|lx)==0) /* asin(1)=+-pi/2 with inexact */ - return x*pio2_hi+x*pio2_lo; - return (x-x)/(x-x); /* asin(|x|>1) is NaN */ + return x*pio2_hi+x*pio2_lo; + return (x-x)/(x-x); /* asin(|x|>1) is NaN */ } else if (ix<0x3fe00000) { /* |x|<0.5 */ if(ix<0x3e400000) { /* if |x| < 2**-27 */ if(huge+x>one) return x;/* return x with inexact if x!=0*/ - } else { + } else { t = x*x; p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); @@ -116,6 +116,6 @@ qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */ p = 2.0*s*r-(pio2_lo-2.0*c); q = pio4_hi-2.0*w; t = pio4_hi-(p-q); - } - if(hx>0) return t; else return -t; + } + if(hx>0) return t; else return -t; } diff --git a/libm/e_atan2.c b/libm/e_atan2.c index 2b5edeacf..8bfc79232 100644 --- a/libm/e_atan2.c +++ b/libm/e_atan2.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -17,7 +17,7 @@ static char rcsid[] = "$NetBSD: e_atan2.c,v 1.8 1995/05/10 20:44:51 jtc Exp $"; /* __ieee754_atan2(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): + * 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, * @@ -35,9 +35,9 @@ static char rcsid[] = "$NetBSD: e_atan2.c,v 1.8 1995/05/10 20:44:51 jtc Exp $"; * 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 + * 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. */ @@ -45,9 +45,9 @@ static char rcsid[] = "$NetBSD: e_atan2.c,v 1.8 1995/05/10 20:44:51 jtc Exp $"; #include "math_private.h" #ifdef __STDC__ -static const double +static const double #else -static double +static double #endif tiny = 1.0e-300, zero = 0.0, @@ -62,7 +62,7 @@ pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ double __ieee754_atan2(y,x) double y,x; #endif -{ +{ double z; int32_t k,m,hx,hy,ix,iy; u_int32_t lx,ly; @@ -80,7 +80,7 @@ pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ /* when y = 0 */ if((iy|ly)==0) { switch(m) { - case 0: + 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 */ @@ -88,7 +88,7 @@ pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */ } /* when x = 0 */ if((ix|lx)==0) return (hy<0)? -pi_o_2-tiny: pi_o_2+tiny; - + /* when x is INF */ if(ix==0x7ff00000) { if(iy==0x7ff00000) { diff --git a/libm/e_atanh.c b/libm/e_atanh.c index 559e8f158..d199a05cc 100644 --- a/libm/e_atanh.c +++ b/libm/e_atanh.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -21,7 +21,7 @@ static char rcsid[] = "$NetBSD: e_atanh.c,v 1.8 1995/05/10 20:44:55 jtc Exp $"; * 1 2x x * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) * 2 1 - x 1 - x - * + * * For x<0.5 * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) * @@ -61,14 +61,14 @@ static double zero = 0.0; ix = hx&0x7fffffff; if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */ return (x-x)/(x-x); - if(ix==0x3ff00000) + if(ix==0x3ff00000) return x/zero; if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */ SET_HIGH_WORD(x,ix); if(ix<0x3fe00000) { /* x < 0.5 */ t = x+x; t = 0.5*log1p(t+t*x/(one-x)); - } else + } else t = 0.5*log1p((x+x)/(one-x)); if(hx>=0) return t; else return -t; } diff --git a/libm/e_cosh.c b/libm/e_cosh.c index 5d73fb5c4..88499684a 100644 --- a/libm/e_cosh.c +++ b/libm/e_cosh.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -15,18 +15,18 @@ static char rcsid[] = "$NetBSD: e_cosh.c,v 1.7 1995/05/10 20:44:58 jtc Exp $"; #endif /* __ieee754_cosh(x) - * Method : + * Method : * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2 - * 1. Replace x by |x| (cosh(x) = cosh(-x)). - * 2. - * [ exp(x) - 1 ]^2 + * 1. Replace x by |x| (cosh(x) = cosh(-x)). + * 2. + * [ exp(x) - 1 ]^2 * 0 <= x <= ln2/2 : cosh(x) := 1 + ------------------- * 2*exp(x) * * exp(x) + 1/exp(x) * ln2/2 <= x <= 22 : cosh(x) := ------------------- * 2 - * 22 <= x <= lnovft : cosh(x) := exp(x)/2 + * 22 <= x <= lnovft : cosh(x) := exp(x)/2 * lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2) * ln2ovft < x : cosh(x) := huge*huge (overflow) * @@ -50,7 +50,7 @@ static double one = 1.0, half=0.5, huge = 1.0e300; double __ieee754_cosh(x) double x; #endif -{ +{ double t,w; int32_t ix; u_int32_t lx; @@ -60,7 +60,7 @@ static double one = 1.0, half=0.5, huge = 1.0e300; ix &= 0x7fffffff; /* x is INF or NaN */ - if(ix>=0x7ff00000) return x*x; + if(ix>=0x7ff00000) return x*x; /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ if(ix<0x3fd62e43) { @@ -81,7 +81,7 @@ static double one = 1.0, half=0.5, huge = 1.0e300; /* |x| in [log(maxdouble), overflowthresold] */ GET_LOW_WORD(lx,x); - if (ix<0x408633CE || + if (ix<0x408633CE || ((ix==0x408633ce)&&(lx<=(u_int32_t)0x8fb9f87d))) { w = __ieee754_exp(half*fabs(x)); t = half*w; diff --git a/libm/e_exp.c b/libm/e_exp.c index f92910e85..f4d832bbb 100644 --- a/libm/e_exp.c +++ b/libm/e_exp.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -22,36 +22,36 @@ static char rcsid[] = "$NetBSD: e_exp.c,v 1.8 1995/05/10 20:45:03 jtc Exp $"; * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. * Given x, find r and integer k such that * - * x = k*ln2 + r, |r| <= 0.5*ln2. + * x = k*ln2 + r, |r| <= 0.5*ln2. * - * Here r will be represented as r = hi-lo for better + * Here r will be represented as r = hi-lo for better * accuracy. * * 2. Approximation of exp(r) by a special rational function on * the interval [0,0.34658]: * Write * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... - * We use a special Reme algorithm on [0,0.34658] to generate - * a polynomial of degree 5 to approximate R. The maximum error + * We use a special Reme algorithm on [0,0.34658] to generate + * a polynomial of degree 5 to approximate R. The maximum error * of this polynomial approximation is bounded by 2**-59. In * other words, * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 * (where z=r*r, and the values of P1 to P5 are listed below) * and * | 5 | -59 - * | 2.0+P1*z+...+P5*z - R(z) | <= 2 + * | 2.0+P1*z+...+P5*z - R(z) | <= 2 * | | * The computation of exp(r) thus becomes * 2*r * exp(r) = 1 + ------- * R - r - * r*R1(r) + * r*R1(r) * = 1 + r + ----------- (for better accuracy) * 2 - R1(r) * where * 2 4 10 * R1(r) = r - (P1*r + P2*r + ... + P5*r ). - * + * * 3. Scale back to obtain exp(x): * From step 1, we have * exp(x) = 2^k * exp(r) @@ -66,13 +66,13 @@ static char rcsid[] = "$NetBSD: e_exp.c,v 1.8 1995/05/10 20:45:03 jtc Exp $"; * 1 ulp (unit in the last place). * * Misc. info. - * For IEEE double + * For IEEE double * if x > 7.09782712893383973096e+02 then exp(x) overflow * if x < -7.45133219101941108420e+02 then exp(x) underflow * * Constants: - * The hexadecimal values are the intended ones for the following - * constants. The decimal values may be used, provided that the + * 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. */ @@ -128,7 +128,7 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ if(hx>=0x7ff00000) { u_int32_t lx; GET_LOW_WORD(lx,x); - if(((hx&0xfffff)|lx)!=0) + if(((hx&0xfffff)|lx)!=0) return x+x; /* NaN */ else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */ } @@ -137,7 +137,7 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ } /* argument reduction */ - if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ + if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb; } else { @@ -147,7 +147,7 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ lo = t*ln2LO[0]; } x = hi - lo; - } + } else if(hx < 0x3e300000) { /* when |x|<2**-28 */ if(huge+x>one) return one+x;/* trigger inexact */ } @@ -156,7 +156,7 @@ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ /* x is now in primary range */ t = x*x; c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5)))); - if(k==0) return one-((x*c)/(c-2.0)-x); + if(k==0) return one-((x*c)/(c-2.0)-x); else y = one-((lo-(x*c)/(2.0-c))-hi); if(k >= -1021) { u_int32_t hy; diff --git a/libm/e_fmod.c b/libm/e_fmod.c index 2ce613574..565c3772b 100644 --- a/libm/e_fmod.c +++ b/libm/e_fmod.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -14,7 +14,7 @@ static char rcsid[] = "$NetBSD: e_fmod.c,v 1.8 1995/05/10 20:45:07 jtc Exp $"; #endif -/* +/* * __ieee754_fmod(x,y) * Return x mod y in exact arithmetic * Method: shift and subtract @@ -51,7 +51,7 @@ static double one = 1.0, Zero[] = {0.0, -0.0,}; return (x*y)/(x*y); if(hx<=hy) { if((hx>31]; /* |x|=|y| return x*0*/ } @@ -74,7 +74,7 @@ static double one = 1.0, Zero[] = {0.0, -0.0,}; } else iy = (hy>>20)-1023; /* set up {hx,lx}, {hy,ly} and align y to x */ - if(ix >= -1022) + if(ix >= -1022) hx = 0x00100000|(0x000fffff&hx); else { /* subnormal x, shift x to normal */ n = -1022-ix; @@ -86,7 +86,7 @@ static double one = 1.0, Zero[] = {0.0, -0.0,}; lx = 0; } } - if(iy >= -1022) + if(iy >= -1022) hy = 0x00100000|(0x000fffff&hy); else { /* subnormal y, shift y to normal */ n = -1022-iy; @@ -115,7 +115,7 @@ static double one = 1.0, Zero[] = {0.0, -0.0,}; /* convert back to floating value and restore the sign */ if((hx|lx)==0) /* return sign(x)*0 */ - return Zero[(u_int32_t)sx>>31]; + return Zero[(u_int32_t)sx>>31]; while(hx<0x00100000) { /* normalize x */ hx = hx+hx+(lx>>31); lx = lx+lx; iy -= 1; diff --git a/libm/e_gamma.c b/libm/e_gamma.c index c4ea7a903..a82261370 100644 --- a/libm/e_gamma.c +++ b/libm/e_gamma.c @@ -6,7 +6,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== * diff --git a/libm/e_gamma_r.c b/libm/e_gamma_r.c index 909c4203d..039eba49b 100644 --- a/libm/e_gamma_r.c +++ b/libm/e_gamma_r.c @@ -6,15 +6,15 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== * */ /* __ieee754_gamma_r(x, signgamp) - * Reentrant version of the logarithm of the Gamma function - * with user provide pointer for the sign of Gamma(x). + * Reentrant version of the logarithm of the Gamma function + * with user provide pointer for the sign of Gamma(x). * * Method: See __ieee754_lgamma_r */ diff --git a/libm/e_hypot.c b/libm/e_hypot.c index 24c8ae452..2052e215b 100644 --- a/libm/e_hypot.c +++ b/libm/e_hypot.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -16,12 +16,12 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; /* __ieee754_hypot(x,y) * - * Method : - * If (assume round-to-nearest) z=x*x+y*y - * has error less than sqrt(2)/2 ulp, than + * 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 + * So, compute sqrt(x*x+y*y) with some care as * follows to get the error below 1 ulp: * * Assume x>y>0; @@ -31,10 +31,10 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; * 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, + * 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 + * + * NOTE: scaling may be necessary if some argument is too * large or too tiny * * Special cases: @@ -42,8 +42,8 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; * 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) + * hypot(x,y) returns sqrt(x^2+y^2) with error less + * than 1 ulps (units in the last place) */ #include "math.h" @@ -84,7 +84,7 @@ static char rcsid[] = "$NetBSD: e_hypot.c,v 1.9 1995/05/12 04:57:27 jtc Exp $"; SET_HIGH_WORD(b,hb); } if(hb < 0x20b00000) { /* b < 2**-500 */ - if(hb <= 0x000fffff) { /* subnormal b or 0 */ + if(hb <= 0x000fffff) { /* subnormal b or 0 */ u_int32_t low; GET_LOW_WORD(low,b); if((hb|low)==0) return a; diff --git a/libm/e_j0.c b/libm/e_j0.c index 46fa16656..b347481ec 100644 --- a/libm/e_j0.c +++ b/libm/e_j0.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -33,20 +33,20 @@ static char rcsid[] = "$NetBSD: e_j0.c,v 1.8 1995/05/10 20:45:23 jtc Exp $"; * (To avoid cancellation, use * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) * to compute the worse one.) - * + * * 3 Special cases * j0(nan)= nan * j0(0) = 1 * j0(inf) = 0 - * + * * Method -- y0(x): * 1. For x<2. - * Since + * Since * y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...) * therefore y0(x)-2/pi*j0(x)*ln(x) is an even function. * We use the following function to approximate y0, * y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2 - * where + * where * U(z) = u00 + u01*z + ... + u06*z^6 * V(z) = 1 + v01*z + ... + v04*z^4 * with absolute approximation error bounded by 2**-72. @@ -69,9 +69,9 @@ static double pzero(), qzero(); #endif #ifdef __STDC__ -static const double +static const double #else -static double +static double #endif huge = 1e300, one = 1.0, @@ -94,9 +94,9 @@ static double zero = 0.0; #endif #ifdef __STDC__ - double __ieee754_j0(double x) + double __ieee754_j0(double x) #else - double __ieee754_j0(x) + double __ieee754_j0(x) double x; #endif { @@ -163,9 +163,9 @@ v03 = 2.59150851840457805467e-07, /* 0x3E91642D, 0x7FF202FD */ v04 = 4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */ #ifdef __STDC__ - double __ieee754_y0(double x) + double __ieee754_y0(double x) #else - double __ieee754_y0(x) + double __ieee754_y0(x) double x; #endif { @@ -175,7 +175,7 @@ v04 = 4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */ EXTRACT_WORDS(hx,lx,x); ix = 0x7fffffff&hx; /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */ - if(ix>=0x7ff00000) return one/(x+x*x); + if(ix>=0x7ff00000) return one/(x+x*x); if((ix|lx)==0) return -one/zero; if(hx<0) return zero/zero; if(ix >= 0x40000000) { /* |x| >= 2.0 */ @@ -349,7 +349,7 @@ static double pS2[5] = { s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); return one+ r/s; } - + /* For x >= 8, the asymptotic expansions of qzero is * -1/8 s + 75/1024 s^3 - ..., where s = 1/x. diff --git a/libm/e_j1.c b/libm/e_j1.c index 5eb81ee73..e57049e84 100644 --- a/libm/e_j1.c +++ b/libm/e_j1.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -34,16 +34,16 @@ static char rcsid[] = "$NetBSD: e_j1.c,v 1.8 1995/05/10 20:45:27 jtc Exp $"; * (To avoid cancellation, use * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) * to compute the worse one.) - * + * * 3 Special cases * j1(nan)= nan * j1(0) = 0 * j1(inf) = 0 - * + * * Method -- y1(x): - * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN + * 1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN * 2. For x<2. - * Since + * Since * y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...) * therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function. * We use the following function to approximate y1, @@ -69,9 +69,9 @@ static double pone(), qone(); #endif #ifdef __STDC__ -static const double +static const double #else -static double +static double #endif huge = 1e300, one = 1.0, @@ -95,9 +95,9 @@ static double zero = 0.0; #endif #ifdef __STDC__ - double __ieee754_j1(double x) + double __ieee754_j1(double x) #else - double __ieee754_j1(x) + double __ieee754_j1(x) double x; #endif { @@ -164,9 +164,9 @@ static double V0[5] = { }; #ifdef __STDC__ - double __ieee754_y1(double x) + double __ieee754_y1(double x) #else - double __ieee754_y1(x) + double __ieee754_y1(x) double x; #endif { @@ -176,7 +176,7 @@ static double V0[5] = { EXTRACT_WORDS(hx,lx,x); ix = 0x7fffffff&hx; /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */ - if(ix>=0x7ff00000) return one/(x+x*x); + if(ix>=0x7ff00000) return one/(x+x*x); if((ix|lx)==0) return -one/zero; if(hx<0) return zero/zero; if(ix >= 0x40000000) { /* |x| >= 2.0 */ @@ -206,10 +206,10 @@ static double V0[5] = { z = invsqrtpi*(u*ss+v*cc)/sqrt(x); } return z; - } + } if(ix<=0x3c900000) { /* x < 2**-54 */ return(-tpi/x); - } + } z = x*x; u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4]))); v = one+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4])))); @@ -347,7 +347,7 @@ static double ps2[5] = { s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4])))); return one+ r/s; } - + /* For x >= 8, the asymptotic expansions of qone is * 3/8 s - 105/1024 s^3 - ..., where s = 1/x. diff --git a/libm/e_jn.c b/libm/e_jn.c index 870824cf8..857c4a3f5 100644 --- a/libm/e_jn.c +++ b/libm/e_jn.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -18,7 +18,7 @@ static char rcsid[] = "$NetBSD: e_jn.c,v 1.9 1995/05/10 20:45:34 jtc Exp $"; * __ieee754_jn(n, x), __ieee754_yn(n, x) * floating point Bessel's function of the 1st and 2nd kind * of order n - * + * * Special cases: * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal; * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. @@ -37,7 +37,7 @@ static char rcsid[] = "$NetBSD: e_jn.c,v 1.9 1995/05/10 20:45:34 jtc Exp $"; * yn(n,x) is similar in all respects, except * that forward recursion is used for all * values of n>1. - * + * */ #include "math.h" @@ -76,7 +76,7 @@ static double zero = 0.00000000000000000000e+00; ix = 0x7fffffff&hx; /* if J(n,NaN) is NaN */ if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x; - if(n<0){ + if(n<0){ n = -n; x = -x; hx ^= 0x80000000; @@ -87,13 +87,13 @@ static double zero = 0.00000000000000000000e+00; x = fabs(x); if((ix|lx)==0||ix>=0x7ff00000) /* if x is 0 or inf */ b = zero; - else if((double)n<=x) { + else if((double)n<=x) { /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */ if(ix>=0x52D00000) { /* x > 2**302 */ - /* (x >> n**2) + /* (x >> n**2) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Let s=sin(x), c=cos(x), + * Let s=sin(x), c=cos(x), * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then * * n sin(xn)*sqt2 cos(xn)*sqt2 @@ -110,7 +110,7 @@ static double zero = 0.00000000000000000000e+00; case 3: temp = cos(x)-sin(x); break; } b = invsqrtpi*temp/sqrt(x); - } else { + } else { a = __ieee754_j0(x); b = __ieee754_j1(x); for(i=1;i33) /* underflow */ @@ -136,14 +136,14 @@ static double zero = 0.00000000000000000000e+00; } } else { /* use backward recurrence */ - /* x x^2 x^2 + /* x x^2 x^2 * J(n,x)/J(n-1,x) = ---- ------ ------ ..... * 2n - 2(n+1) - 2(n+2) * - * 1 1 1 + * 1 1 1 * (for large x) = ---- ------ ------ ..... * 2n 2(n+1) 2(n+2) - * -- - ------ - ------ - + * -- - ------ - ------ - * x x x * * Let w = 2n/x and h=2/x, then the above quotient @@ -159,9 +159,9 @@ static double zero = 0.00000000000000000000e+00; * To determine how many terms needed, let * Q(0) = w, Q(1) = w(w+h) - 1, * Q(k) = (w+k*h)*Q(k-1) - Q(k-2), - * When Q(k) > 1e4 good for single - * When Q(k) > 1e9 good for double - * When Q(k) > 1e17 good for quadruple + * When Q(k) > 1e4 good for single + * When Q(k) > 1e9 good for double + * When Q(k) > 1e17 good for quadruple */ /* determine k */ double t,v; @@ -183,7 +183,7 @@ static double zero = 0.00000000000000000000e+00; * single 8.8722839355e+01 * double 7.09782712893383973096e+02 * long double 1.1356523406294143949491931077970765006170e+04 - * then recurrent value may overflow and the result is + * then recurrent value may overflow and the result is * likely underflow to zero */ tmp = n; @@ -219,9 +219,9 @@ static double zero = 0.00000000000000000000e+00; } #ifdef __STDC__ - double __ieee754_yn(int n, double x) + double __ieee754_yn(int n, double x) #else - double __ieee754_yn(n,x) + double __ieee754_yn(n,x) int n; double x; #endif { @@ -244,10 +244,10 @@ static double zero = 0.00000000000000000000e+00; if(n==1) return(sign*__ieee754_y1(x)); if(ix==0x7ff00000) return zero; if(ix>=0x52D00000) { /* x > 2**302 */ - /* (x >> n**2) + /* (x >> n**2) * Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi) * Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi) - * Let s=sin(x), c=cos(x), + * Let s=sin(x), c=cos(x), * xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then * * n sin(xn)*sqt2 cos(xn)*sqt2 @@ -270,7 +270,7 @@ static double zero = 0.00000000000000000000e+00; b = __ieee754_y1(x); /* quit if b is -inf */ GET_HIGH_WORD(high,b); - for(i=1;i= 0x7ff00000) return x+x; k += (hx>>20)-1023; hx &= 0x000fffff; @@ -126,14 +126,14 @@ static double zero = 0.0; if(k==0) return f-R; else {dk=(double)k; return dk*ln2_hi-((R-dk*ln2_lo)-f);} } - s = f/(2.0+f); + s = f/(2.0+f); dk = (double)k; z = s*s; i = hx-0x6147a; w = z*z; j = 0x6b851-hx; - t1= w*(Lg2+w*(Lg4+w*Lg6)); - t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); + t1= w*(Lg2+w*(Lg4+w*Lg6)); + t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7))); i |= j; R = t2+t1; if(i>0) { diff --git a/libm/e_log10.c b/libm/e_log10.c index 5d004ac4e..2716a72c1 100644 --- a/libm/e_log10.c +++ b/libm/e_log10.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -16,26 +16,26 @@ static char rcsid[] = "$NetBSD: e_log10.c,v 1.9 1995/05/10 20:45:51 jtc Exp $"; /* __ieee754_log10(x) * Return the base 10 logarithm of x - * + * * Method : * Let log10_2hi = leading 40 bits of log10(2) and * log10_2lo = log10(2) - log10_2hi, * ivln10 = 1/log(10) rounded. * Then - * n = ilogb(x), + * n = ilogb(x), * if(n<0) n = n+1; * x = scalbn(x,-n); * log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x)) * * Note 1: - * To guarantee log10(10**n)=n, where 10**n is normal, the rounding + * To guarantee log10(10**n)=n, where 10**n is normal, the rounding * mode must set to Round-to-Nearest. * Note 2: * [1/log(10)] rounded to 53 bits has error .198 ulps; * log10 is monotonic at all binary break points. * * Special cases: - * log10(x) is NaN with signal if x < 0; + * log10(x) is NaN with signal if x < 0; * log10(+INF) is +INF with no signal; log10(0) is -INF with signal; * log10(NaN) is that NaN with no signal; * log10(10**N) = N for N=0,1,...,22. diff --git a/libm/e_pow.c b/libm/e_pow.c index 4f6a44f20..b970775cd 100644 --- a/libm/e_pow.c +++ b/libm/e_pow.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -21,7 +21,7 @@ static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $"; * 1. Compute and return log2(x) in two pieces: * log2(x) = w1 + w2, * where w1 has 53-24 = 29 bit trailing zeros. - * 2. Perform y*log2(x) = n+y' by simulating muti-precision + * 2. Perform y*log2(x) = n+y' by simulating muti-precision * arithmetic, where |y'|<=0.5. * 3. Return x**y = 2**n*exp(y'*log2) * @@ -49,13 +49,13 @@ static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $"; * Accuracy: * pow(x,y) returns x**y nearly rounded. In particular * pow(integer,integer) - * always returns the correct integer provided it is + * always returns the correct integer provided it is * representable. * * 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 + * 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. */ @@ -63,9 +63,9 @@ static char rcsid[] = "$NetBSD: e_pow.c,v 1.9 1995/05/12 04:57:32 jtc Exp $"; #include "math_private.h" #ifdef __STDC__ -static const double +static const double #else -static double +static double #endif bp[] = {1.0, 1.5,}, dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ @@ -117,12 +117,12 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ ix = hx&0x7fffffff; iy = hy&0x7fffffff; /* y==zero: x**0 = 1 */ - if((iy|ly)==0) return one; + if((iy|ly)==0) return one; /* +-NaN return x+y */ if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) || - iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) - return x+y; + iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0))) + return x+y; /* determine if y is an odd int when x < 0 * yisint = 0 ... y is not an integer @@ -130,7 +130,7 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ * yisint = 2 ... y is an even int */ yisint = 0; - if(hx<0) { + if(hx<0) { if(iy>=0x43400000) yisint = 2; /* even integer y */ else if(iy>=0x3ff00000) { k = (iy>>20)-0x3ff; /* exponent */ @@ -141,11 +141,11 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ j = iy>>(20-k); if((j<<(20-k))==iy) yisint = 2-(j&1); } - } - } + } + } /* special value of y */ - if(ly==0) { + if(ly==0) { if (iy==0x7ff00000) { /* y is +-inf */ if(((ix-0x3ff00000)|lx)==0) return y - y; /* inf**+-1 is NaN */ @@ -153,14 +153,14 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ return (hy>=0)? y: zero; else /* (|x|<1)**-,+inf = inf,0 */ return (hy<0)?-y: zero; - } + } if(iy==0x3ff00000) { /* y is +-1 */ if(hy<0) return one/x; else return x; } if(hy==0x40000000) return x*x; /* y is 2 */ if(hy==0x3fe00000) { /* y is 0.5 */ if(hx>=0) /* x >= +0 */ - return __ieee754_sqrt(x); + return __ieee754_sqrt(x); } } @@ -173,13 +173,13 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ if(hx<0) { if(((ix-0x3ff00000)|yisint)==0) { z = (z-z)/(z-z); /* (-1)**non-int is NaN */ - } else if(yisint==1) + } else if(yisint==1) z = -z; /* (x<0)**odd = -(|x|**odd) */ } return z; } } - + /* (x<0)**(non-int) is NaN */ if(((((u_int32_t)hx>>31)-1)|yisint)==0) return (x-x)/(x-x); @@ -192,7 +192,7 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ /* over/underflow if x is not close to one */ if(ix<0x3fefffff) return (hy<0)? huge*huge:tiny*tiny; if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny; - /* now |1-x| is tiny <= 2**-20, suffice to compute + /* now |1-x| is tiny <= 2**-20, suffice to compute log(x) by x-x^2/2+x^3/3-x^4/4 */ t = x-1; /* t has 20 trailing zeros */ w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25)); @@ -289,7 +289,7 @@ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ n = ((n&0x000fffff)|0x00100000)>>(20-k); if(j<0) n = -n; p_h -= t; - } + } t = p_l+p_h; SET_LOW_WORD(t,0); u = t*lg2_h; diff --git a/libm/e_rem_pio2.c b/libm/e_rem_pio2.c index e98c699cd..3dd7f7b49 100644 --- a/libm/e_rem_pio2.c +++ b/libm/e_rem_pio2.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -15,8 +15,8 @@ static char rcsid[] = "$NetBSD: e_rem_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $ #endif /* __ieee754_rem_pio2(x,y) - * - * return the remainder of x rem pi/2 in y[0]+y[1] + * + * return the remainder of x rem pi/2 in y[0]+y[1] * use __kernel_rem_pio2() */ @@ -24,24 +24,24 @@ static char rcsid[] = "$NetBSD: e_rem_pio2.c,v 1.8 1995/05/10 20:46:02 jtc Exp $ #include "math_private.h" /* - * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi + * Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi */ #ifdef __STDC__ static const int32_t two_over_pi[] = { #else static int32_t two_over_pi[] = { #endif -0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, -0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, -0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, -0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, -0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, -0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, -0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, -0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, -0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, -0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, -0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, +0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, +0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, +0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, +0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, +0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, +0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, +0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, +0xF17B3D, 0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, +0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, +0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, +0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B, }; #ifdef __STDC__ @@ -68,9 +68,9 @@ static int32_t npio2_hw[] = { */ #ifdef __STDC__ -static const double +static const double #else -static double +static double #endif zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */ half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ @@ -100,7 +100,7 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ if(ix<=0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */ {y[0] = x; y[1] = 0; return 0;} if(ix<0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */ - if(hx>0) { + if(hx>0) { z = x - pio2_1; if(ix!=0x3ff921fb) { /* 33+53 bit pi is good enough */ y[0] = z - pio2_1t; @@ -130,27 +130,27 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ fn = (double)n; r = t-fn*pio2_1; w = fn*pio2_1t; /* 1st round good to 85 bit */ - if(n<32&&ix!=npio2_hw[n-1]) { + if(n<32&&ix!=npio2_hw[n-1]) { y[0] = r-w; /* quick check no cancellation */ } else { u_int32_t high; j = ix>>20; - y[0] = r-w; + y[0] = r-w; GET_HIGH_WORD(high,y[0]); i = j-((high>>20)&0x7ff); if(i>16) { /* 2nd iteration needed, good to 118 */ t = r; - w = fn*pio2_2; + w = fn*pio2_2; r = t-w; - w = fn*pio2_2t-((t-r)-w); + w = fn*pio2_2t-((t-r)-w); y[0] = r-w; GET_HIGH_WORD(high,y[0]); i = j-((high>>20)&0x7ff); if(i>49) { /* 3rd iteration need, 151 bits acc */ t = r; /* will cover all possible cases */ - w = fn*pio2_3; + w = fn*pio2_3; r = t-w; - w = fn*pio2_3t-((t-r)-w); + w = fn*pio2_3t-((t-r)-w); y[0] = r-w; } } @@ -159,7 +159,7 @@ pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ if(hx<0) {y[0] = -y[0]; y[1] = -y[1]; return -n;} else return n; } - /* + /* * all other (large) arguments */ if(ix>=0x7ff00000) { /* x is inf or NaN */ diff --git a/libm/e_remainder.c b/libm/e_remainder.c index 641808118..52e8a28d3 100644 --- a/libm/e_remainder.c +++ b/libm/e_remainder.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -15,11 +15,11 @@ static char rcsid[] = "$NetBSD: e_remainder.c,v 1.8 1995/05/10 20:46:05 jtc Exp #endif /* __ieee754_remainder(x,p) - * Return : - * returns x REM p = x - [x/p]*p as if in infinite - * precise arithmetic, where [x/p] is the (infinite bit) + * 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 : + * Method : * Based on fmod() return x-[x/p]chopped*p exactlp. */ diff --git a/libm/e_scalb.c b/libm/e_scalb.c index 7f66ec773..9a366f330 100644 --- a/libm/e_scalb.c +++ b/libm/e_scalb.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -16,7 +16,7 @@ static char rcsid[] = "$NetBSD: e_scalb.c,v 1.6 1995/05/10 20:46:09 jtc Exp $"; /* * __ieee754_scalb(x, fn) is provide for - * passing various standard test suite. One + * passing various standard test suite. One * should use scalbn() instead. */ diff --git a/libm/e_sinh.c b/libm/e_sinh.c index e1e614c9e..c63ff8ade 100644 --- a/libm/e_sinh.c +++ b/libm/e_sinh.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -15,15 +15,15 @@ static char rcsid[] = "$NetBSD: e_sinh.c,v 1.7 1995/05/10 20:46:13 jtc Exp $"; #endif /* __ieee754_sinh(x) - * Method : + * Method : * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2 - * 1. Replace x by |x| (sinh(-x) = -sinh(x)). - * 2. + * 1. Replace x by |x| (sinh(-x) = -sinh(x)). + * 2. * E + E/(E+1) * 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x) * 2 * - * 22 <= x <= lnovft : sinh(x) := exp(x)/2 + * 22 <= x <= lnovft : sinh(x) := exp(x)/2 * lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2) * ln2ovft < x : sinh(x) := x*shuge (overflow) * @@ -47,7 +47,7 @@ static double one = 1.0, shuge = 1.0e307; double __ieee754_sinh(x) double x; #endif -{ +{ double t,w,h; int32_t ix,jx; u_int32_t lx; @@ -57,7 +57,7 @@ static double one = 1.0, shuge = 1.0e307; ix = jx&0x7fffffff; /* x is INF or NaN */ - if(ix>=0x7ff00000) return x+x; + if(ix>=0x7ff00000) return x+x; h = 0.5; if (jx<0) h = -h; diff --git a/libm/e_sqrt.c b/libm/e_sqrt.c index 15fba001d..fae52c842 100644 --- a/libm/e_sqrt.c +++ b/libm/e_sqrt.c @@ -5,7 +5,7 @@ * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice + * software is freely granted, provided that this notice * is preserved. * ==================================================== */ @@ -19,10 +19,10 @@ static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $"; * ------------------------------------------ * | Use the hardware sqrt if you have one | * ------------------------------------------ - * Method: - * Bit by bit method using integer arithmetic. (Slow, but portable) + * Method: + * Bit by bit method using integer arithmetic. (Slow, but portable) * 1. Normalization - * Scale x to y in [1,4) with even powers of 2: + * Scale x to y in [1,4) with even powers of 2: * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then * sqrt(x) = 2^k * sqrt(y) * 2. Bit by bit computation @@ -31,9 +31,9 @@ static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $"; * i+1 2 * s = 2*q , and y = 2 * ( y - q ). (1) * i i i i - * - * To compute q from q , one checks whether - * i+1 i + * + * To compute q from q , one checks whether + * i+1 i * * -(i+1) 2 * (q + 2 ) <= y. (2) @@ -43,12 +43,12 @@ static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $"; * i+1 i i+1 i * * With some algebric manipulation, it is not difficult to see - * that (2) is equivalent to + * that (2) is equivalent to * -(i+1) * s + 2 <= y (3) * i i * - * The advantage of (3) is that s and y can be computed by + * The advantage of (3) is that s and y can be computed by * i i * the following recurrence formula: * if (3) is false @@ -60,10 +60,10 @@ static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $"; * -i -(i+1) * s = s + 2 , y = y - s - 2 (5) * i+1 i i+1 i i - * - * One may easily use induction to prove (4) and (5). + * + * One may easily use induction to prove (4) and (5). * Note. Since the left hand side of (3) contain only i+2 bits, - * it does not necessary to do a full (53-bit) comparison + * it does not necessary to do a full (53-bit) comparison * in (3). * 3. Final rounding * After generating the 53 bits result, we compute one more bit. @@ -73,7 +73,7 @@ static char rcsid[] = "$NetBSD: e_sqrt.c,v 1.8 1995/05/10 20:46:17 jtc Exp $"; * The rounding mode can be detected by checking whether * huge + tiny is equal to huge, and whether huge - tiny is * equal to huge for some floating point number "huge" and "tiny". - * + * * Special cases: * sqrt(+-0) = +-0 ... exact * sqrt(inf) = inf @@ -101,17 +101,17 @@ static double one = 1.0, tiny=1.0e-300; #endif { double z; - int32_t sign = (int)0x80000000; + int32_t sign = (int)0x80000000; int32_t ix0,s0,q,m,t,i; u_int32_t r,t1,s1,ix1,q1; EXTRACT_WORDS(ix0,ix1,x); /* take care of Inf and NaN */ - if((ix0&0x7ff00000)==0x7ff00000) { + if((ix0&0x7ff00000)==0x7ff00000) { return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf sqrt(-inf)=sNaN */ - } + } /* take care of zero */ if(ix0<=0) { if(((ix0&(~sign))|ix1)==0) return x;/* sqrt(+-0) = +-0 */ @@ -145,12 +145,12 @@ static double one = 1.0, tiny=1.0e-300; r = 0x00200000; /* r = moving bit from right to left */ while(r!=0) { - t = s0+r; - if(t<=ix0) { - s0 = t+r; - ix0 -= t; - q += r; - } + t = s0+r; + if(t<=ix0) { + s0 = t+r; + ix0 -= t; + q += r; + } ix0 += ix0 + ((ix1&sign)>>31); ix1 += ix1; r>>=1; @@ -158,9 +158,9 @@ static double one = 1.0, tiny=1.0e-300; r = sign; while(r!=0) { - t1 = s1+r; + t1 = s1+r; t = s0; - if((tone) { if (q1==(u_int32_t)0xfffffffe) q+=1; - q1+=2; + q1+=2; } else q1 += (q1&1); } @@ -197,18 +197,18 @@ static double one = 1.0, tiny=1.0e-300; /* Other methods (use floating-point arithmetic) ------------- -(This is a copy of a drafted paper by Prof W. Kahan +(This is a copy of a drafted paper by Prof W. Kahan and K.C. Ng, written in May, 1986) - Two algorithms are given here to implement sqrt(x) + Two algorithms are given here to implement sqrt(x) (IEEE double precision arithmetic) in software. Both supply sqrt(x) correctly rounded. The first algorithm (in Section A) uses newton iterations and involves four divisions. The second one uses reciproot iterations to avoid division, but requires more multiplications. Both algorithms need the ability - to chop results of arithmetic operations instead of round them, + to chop results of arithmetic operations instead of round them, and the INEXACT flag to indicate when an arithmetic operation - is executed exactly with no roundoff error, all