008-03-03 Steven Munroe * sysdeps/powerpc/powerpc32/power4/fpu/e_log.c: New file. * sysdeps/powerpc/powerpc32/power4/fpu/w_log.c: New file. * sysdeps/powerpc/powerpc32/power6/fpu/w_log.c: New file. * sysdeps/powerpc/powerpc64/power4/fpu/e_log.c: New file. * sysdeps/powerpc/powerpc64/power4/fpu/w_log.c: New file. * sysdeps/powerpc/powerpc64/power6/fpu/w_log.c: New file. diff -urN libc25-cvstip-20070919/sysdeps/powerpc/powerpc32/power4/fpu/e_log.c libc25/sysdeps/powerpc/powerpc32/power4/fpu/e_log.c --- /dev/null 2007-05-04 05:56:50.000000000 -0500 +++ libc25/sysdeps/powerpc/powerpc32/power4/fpu/e_log.c 2008-01-04 12:46:04.000000000 -0600 @@ -0,0 +1,208 @@ +/* + * IBM Accurate Mathematical Library + * written by International Business Machines Corp. + * Copyright (C) 2001 Free Software Foundation + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation; either version 2.1 of the License, or + * (at your option) any later version. + * + * This program 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ +/*********************************************************************/ +/* */ +/* MODULE_NAME:ulog.c */ +/* */ +/* FUNCTION:ulog */ +/* */ +/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h ulog.h */ +/* mpexp.c mplog.c mpa.c */ +/* ulog.tbl */ +/* */ +/* An ultimate log routine. Given an IEEE double machine number x */ +/* it computes the correctly rounded (to nearest) value of log(x). */ +/* Assumption: Machine arithmetic operations are performed in */ +/* round to nearest mode of IEEE 754 standard. */ +/* */ +/*********************************************************************/ + + +#include "endian.h" +#include "dla.h" +#include "mpa.h" +#include "MathLib.h" +#include "math_private.h" + +#undef ABS +#define ABS(x) __builtin_fabs(x) + +void __mplog(mp_no *, mp_no *, int); + +/*********************************************************************/ +/* An ultimate log routine. Given an IEEE double machine number x */ +/* it computes the correctly rounded (to nearest) value of log(x). */ +/*********************************************************************/ +double __ieee754_log(double x) { +#define M 4 + static const int pr[M]={8,10,18,32}; + int i,j,n,ux,dx,p; +#if 0 + int k; +#endif + double dbl_n,u,p0,q,r0,w,nln2a,luai,lubi,lvaj,lvbj, + sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb, + t1,t2,t3,t4,t5,t6,t7,t8,t,ra,rb,ww, + a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c; + number num; + mp_no mpx,mpy,mpy1,mpy2,mperr; + +#include "ulog.tbl" +#include "ulog.h" + + /* Treating special values of x ( x<=0, x=INF, x=NaN etc.). */ + + num.d = x; + + w = x-ONE; + t8 = MHALF*w; + n=0; + + ux = num.i[HIGH_HALF]; dx = num.i[LOW_HALF]; + if (ux < 0x00100000) { + if (((ux & 0x7fffffff) | dx) == 0) return MHALF/ZERO; /* return -INF */ + if (ux < 0) return (x-x)/ZERO; /* return NaN */ + n -= 54; x *= two54.d; /* scale x */ + num.d = x; + } + if (ux >= 0x7ff00000) return x+x; /* INF or NaN */ + + /* Regular values of x */ + + if (ABS(w) > U03) { goto case_03; } + + + /*--- Stage I, the case abs(x-1) < 0.03 */ + EMULV(t8,w,a,aa,t1,t2,t3,t4,t5) + EADD(w,a,b,bb) + + /* Evaluate polynomial II */ + polII = (b0.d+w*(b1.d+w*(b2.d+w*(b3.d+w*(b4.d+ + w*(b5.d+w*(b6.d+w*(b7.d+w*b8.d))))))))*w*w*w; + c = (aa+bb)+polII; + + /* End stage I, case abs(x-1) < 0.03 */ + if ((y=b+(c+b*E2)) == b+(c-b*E2)) return y; + + /*--- Stage II, the case abs(x-1) < 0.03 */ + + a = d11.d+w*(d12.d+w*(d13.d+w*(d14.d+w*(d15.d+w*(d16.d+ + w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d)))))))); + EMULV(w,a,s2,ss2,t1,t2,t3,t4,t5) + ADD2(d10.d,dd10.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d9.d,dd9.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d8.d,dd8.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d7.d,dd7.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d6.d,dd6.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d5.d,dd5.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d4.d,dd4.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d3.d,dd3.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d2.d,dd2.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + MUL2(w,ZERO,s2,ss2,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(w,ZERO, s3,ss3, b, bb,t1,t2) + + /* End stage II, case abs(x-1) < 0.03 */ + if ((y=b+(bb+b*E4)) == b+(bb-b*E4)) return y; + goto stage_n; + + /*--- Stage I, the case abs(x-1) > 0.03 */ + case_03: + + /* Find n,u such that x = u*2**n, 1/sqrt(2) < u < sqrt(2) */ + n += (num.i[HIGH_HALF] >> 20) - 1023; + num.i[HIGH_HALF] = (num.i[HIGH_HALF] & 0x000fffff) | 0x3ff00000; + if (num.d > SQRT_2) { num.d *= HALF; n++; } + u = num.d; dbl_n = (double) n; + + /* Find i such that ui=1+(i-75)/2**8 is closest to u (i= 0,1,2,...,181) */ + num.d += h1.d; + i = (num.i[HIGH_HALF] & 0x000fffff) >> 12; + + /* Find j such that vj=1+(j-180)/2**16 is closest to v=u/ui (j= 0,...,361) */ + num.d = u*Iu[i].d + h2.d; + j = (num.i[HIGH_HALF] & 0x000fffff) >> 4; + + /* Compute w=(u-ui*vj)/(ui*vj) */ + p0=(ONE+(i-75)*DEL_U)*(ONE+(j-180)*DEL_V); + q=u-p0; r0=Iu[i].d*Iv[j].d; w=q*r0; + + /* Evaluate polynomial I */ + polI = w+(a2.d+a3.d*w)*w*w; + + /* Add up everything */ + nln2a = dbl_n*LN2A; + luai = Lu[i][0].d; lubi = Lu[i][1].d; + lvaj = Lv[j][0].d; lvbj = Lv[j][1].d; + EADD(luai,lvaj,sij,ssij) + EADD(nln2a,sij,A ,ttij) + B0 = (((lubi+lvbj)+ssij)+ttij)+dbl_n*LN2B; + B = polI+B0; + + /* End stage I, case abs(x-1) >= 0.03 */ + if ((y=A+(B+E1)) == A+(B-E1)) return y; + + + /*--- Stage II, the case abs(x-1) > 0.03 */ + + /* Improve the accuracy of r0 */ + EMULV(p0,r0,sa,sb,t1,t2,t3,t4,t5) + t=r0*((ONE-sa)-sb); + EADD(r0,t,ra,rb) + + /* Compute w */ + MUL2(q,ZERO,ra,rb,w,ww,t1,t2,t3,t4,t5,t6,t7,t8) + + EADD(A,B0,a0,aa0) + + /* Evaluate polynomial III */ + s1 = (c3.d+(c4.d+c5.d*w)*w)*w; + EADD(c2.d,s1,s2,ss2) + MUL2(s2,ss2,w,ww,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8) + MUL2(s3,ss3,w,ww,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(s2,ss2,w,ww,s3,ss3,t1,t2) + ADD2(s3,ss3,a0,aa0,a1,aa1,t1,t2) + + /* End stage II, case abs(x-1) >= 0.03 */ + if ((y=a1+(aa1+E3)) == a1+(aa1-E3)) return y; + + + /* Final stages. Use multi-precision arithmetic. */ + stage_n: + + for (i=0; i +#include +#include "math_private.h" + + +#ifdef __STDC__ + double __log(double x) /* wrapper log */ +#else + double __log(x) /* wrapper log */ + double x; +#endif +{ +#ifdef _IEEE_LIBM + return __ieee754_log(x); +#else + double z; + + z = __ieee754_log(x); + if(_LIB_VERSION == _IEEE_ || (x != x) || x > 0.0) return z; + + if(x==0.0) + return __kernel_standard(x,x,16); /* log(0) */ + else + return __kernel_standard(x,x,17); /* log(x<0) */ +#endif +} +weak_alias (__log, log) +#ifdef NO_LONG_DOUBLE +strong_alias (__log, __logl) +weak_alias (__log, logl) +#endif +#if LONG_DOUBLE_COMPAT(libm, GLIBC_2_0) +compat_symbol (libm, __log, logl, GLIBC_2_0); +#endif diff -urN libc25-cvstip-20070919/sysdeps/powerpc/powerpc32/power6/fpu/w_log.c libc25/sysdeps/powerpc/powerpc32/power6/fpu/w_log.c --- /dev/null 2007-05-04 05:56:50.000000000 -0500 +++ libc25/sysdeps/powerpc/powerpc32/power6/fpu/w_log.c 2008-02-29 14:28:12.000000000 -0600 @@ -0,0 +1,75 @@ +/* @(#)w_log.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: w_log.c,v 1.6 1995/05/10 20:49:33 jtc Exp $"; +#endif + +/* + * wrapper log(x) + */ + +#include +#include +#include "math_private.h" + +union double_long +{ + double d; + long long ll; +}; + +#ifdef __STDC__ + double __log(double x) /* wrapper log */ +#else + double __log(x) /* wrapper log */ + double x; +#endif +{ +#ifdef _IEEE_LIBM + return __ieee754_log(x); +#else + double z; + volatile union double_long tmp; + long long x_ll; + + /* Copy the double value x to a long to test for NaN without the + overhead of calling isnan or triggering VXSNAN. + Mark the union volatile and store x before calling __ieee754_log() + so we can access the value as a 64-bit long long after + __ieee754_log(). This avoids the load-hit-store hazard when the + store queue is deep. */ + tmp.d = x; + + z = __ieee754_log(x); + + x_ll = tmp.ll; + if(_LIB_VERSION == _IEEE_ + /* equivalent to isnan(x) */ + || ((x_ll & 0x7fffffffffffffffLL) > 0x7ff0000000000000LL) + /* equivalent to (x > 0.0) */ + || x_ll > 0LL) return z; + + if(x==0.0) + return __kernel_standard(x,x,16); /* log(0) */ + else + return __kernel_standard(x,x,17); /* log(x<0) */ +#endif +} +weak_alias (__log, log) +#ifdef NO_LONG_DOUBLE +strong_alias (__log, __logl) +weak_alias (__log, logl) +#endif +#if LONG_DOUBLE_COMPAT(libm, GLIBC_2_0) +compat_symbol (libm, __log, logl, GLIBC_2_0); +#endif diff -urN libc25-cvstip-20070919/sysdeps/powerpc/powerpc64/power4/fpu/e_log.c libc25/sysdeps/powerpc/powerpc64/power4/fpu/e_log.c --- /dev/null 2007-05-04 05:56:50.000000000 -0500 +++ libc25/sysdeps/powerpc/powerpc64/power4/fpu/e_log.c 2008-01-04 12:46:14.000000000 -0600 @@ -0,0 +1,208 @@ +/* + * IBM Accurate Mathematical Library + * written by International Business Machines Corp. + * Copyright (C) 2001 Free Software Foundation + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation; either version 2.1 of the License, or + * (at your option) any later version. + * + * This program 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ +/*********************************************************************/ +/* */ +/* MODULE_NAME:ulog.c */ +/* */ +/* FUNCTION:ulog */ +/* */ +/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h ulog.h */ +/* mpexp.c mplog.c mpa.c */ +/* ulog.tbl */ +/* */ +/* An ultimate log routine. Given an IEEE double machine number x */ +/* it computes the correctly rounded (to nearest) value of log(x). */ +/* Assumption: Machine arithmetic operations are performed in */ +/* round to nearest mode of IEEE 754 standard. */ +/* */ +/*********************************************************************/ + + +#include "endian.h" +#include "dla.h" +#include "mpa.h" +#include "MathLib.h" +#include "math_private.h" + +#undef ABS +#define ABS(x) __builtin_fabs(x) + +void __mplog(mp_no *, mp_no *, int); + +/*********************************************************************/ +/* An ultimate log routine. Given an IEEE double machine number x */ +/* it computes the correctly rounded (to nearest) value of log(x). */ +/*********************************************************************/ +double __ieee754_log(double x) { +#define M 4 + static const int pr[M]={8,10,18,32}; + long i,j,n,ux,p; + + double dbl_n,u,p0,q,r0,w,nln2a,luai,lubi,lvaj,lvbj, + sij,ssij,ttij,A,B,B0,y,y1,y2,polI,polII,sa,sb, + t1,t2,t3,t4,t5,t6,t7,t8,t,ra,rb,ww, + a0,aa0,s1,s2,ss2,s3,ss3,a1,aa1,a,aa,b,bb,c; + + typedef union { long l; double d; } longnum; + longnum num; + mp_no mpx,mpy,mpy1,mpy2,mperr; + +#include "ulog.tbl" +#include "ulog.h" + + /* Treating special values of x ( x<=0, x=INF, x=NaN etc.). */ + + num.d = x; + + n=0; + + ux = num.l; + if (ux < 0x0010000000000000L) { + if ((ux & 0x7fffffffffffffffL) == 0) return MHALF/ZERO; /* return -INF */ + if (ux < 0) return (x-x)/ZERO; /* return NaN */ + n -= 54; x *= two54.d; /* scale x */ + num.d = x; + } + if (ux >= 0x7ff0000000000000L) return x+x; /* INF or NaN */ + + /* Regular values of x */ + w = x-ONE; + if (ABS(w) > U03) { goto case_03; } + + + /*--- Stage I, the case abs(x-1) < 0.03 */ + + t8 = MHALF*w; + EMULV(t8,w,a,aa,t1,t2,t3,t4,t5) + EADD(w,a,b,bb) + + /* Evaluate polynomial II */ + polII = (b0.d+w*(b1.d+w*(b2.d+w*(b3.d+w*(b4.d+ + w*(b5.d+w*(b6.d+w*(b7.d+w*b8.d))))))))*w*w*w; + c = (aa+bb)+polII; + + /* End stage I, case abs(x-1) < 0.03 */ + if ((y=b+(c+b*E2)) == b+(c-b*E2)) return y; + + /*--- Stage II, the case abs(x-1) < 0.03 */ + + a = d11.d+w*(d12.d+w*(d13.d+w*(d14.d+w*(d15.d+w*(d16.d+ + w*(d17.d+w*(d18.d+w*(d19.d+w*d20.d)))))))); + EMULV(w,a,s2,ss2,t1,t2,t3,t4,t5) + ADD2(d10.d,dd10.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d9.d,dd9.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d8.d,dd8.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d7.d,dd7.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d6.d,dd6.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d5.d,dd5.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d4.d,dd4.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d3.d,dd3.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(d2.d,dd2.d,s2,ss2,s3,ss3,t1,t2) + MUL2(w,ZERO,s3,ss3,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + MUL2(w,ZERO,s2,ss2,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(w,ZERO, s3,ss3, b, bb,t1,t2) + + /* End stage II, case abs(x-1) < 0.03 */ + if ((y=b+(bb+b*E4)) == b+(bb-b*E4)) return y; + goto stage_n; + + /*--- Stage I, the case abs(x-1) > 0.03 */ + case_03: + + /* Find n,u such that x = u*2**n, 1/sqrt(2) < u < sqrt(2) */ + n += (num.l >> 52) - 1023; + num.l = (num.l & 0x000fffffffffffffL) | 0x3ff0000000000000L; + if (num.d > SQRT_2) { num.d *= HALF; n++; } + u = num.d; dbl_n = (double) n; + + /* Find i such that ui=1+(i-75)/2**8 is closest to u (i= 0,1,2,...,181) */ + num.d += h1.d; + i = (num.l & 0x000fffff00000000L) >> (12 + 32); + + /* Find j such that vj=1+(j-180)/2**16 is closest to v=u/ui (j= 0,...,361) */ + num.d = u*Iu[i].d + h2.d; + j = (num.l & 0x000fffff00000000L) >> (4 + 32); + + /* Compute w=(u-ui*vj)/(ui*vj) */ + p0=(ONE+(i-75)*DEL_U)*(ONE+(j-180)*DEL_V); + q=u-p0; r0=Iu[i].d*Iv[j].d; w=q*r0; + + /* Evaluate polynomial I */ + polI = w+(a2.d+a3.d*w)*w*w; + + /* Add up everything */ + nln2a = dbl_n*LN2A; + luai = Lu[i][0].d; lubi = Lu[i][1].d; + lvaj = Lv[j][0].d; lvbj = Lv[j][1].d; + EADD(luai,lvaj,sij,ssij) + EADD(nln2a,sij,A ,ttij) + B0 = (((lubi+lvbj)+ssij)+ttij)+dbl_n*LN2B; + B = polI+B0; + + /* End stage I, case abs(x-1) >= 0.03 */ + if ((y=A+(B+E1)) == A+(B-E1)) return y; + + + /*--- Stage II, the case abs(x-1) > 0.03 */ + + /* Improve the accuracy of r0 */ + EMULV(p0,r0,sa,sb,t1,t2,t3,t4,t5) + t=r0*((ONE-sa)-sb); + EADD(r0,t,ra,rb) + + /* Compute w */ + MUL2(q,ZERO,ra,rb,w,ww,t1,t2,t3,t4,t5,t6,t7,t8) + + EADD(A,B0,a0,aa0) + + /* Evaluate polynomial III */ + s1 = (c3.d+(c4.d+c5.d*w)*w)*w; + EADD(c2.d,s1,s2,ss2) + MUL2(s2,ss2,w,ww,s3,ss3,t1,t2,t3,t4,t5,t6,t7,t8) + MUL2(s3,ss3,w,ww,s2,ss2,t1,t2,t3,t4,t5,t6,t7,t8) + ADD2(s2,ss2,w,ww,s3,ss3,t1,t2) + ADD2(s3,ss3,a0,aa0,a1,aa1,t1,t2) + + /* End stage II, case abs(x-1) >= 0.03 */ + if ((y=a1+(aa1+E3)) == a1+(aa1-E3)) return y; + + + /* Final stages. Use multi-precision arithmetic. */ + stage_n: + + for (i=0; i +#include +#include "math_private.h" + +#ifdef __STDC__ + double __log(double x) /* wrapper log */ +#else + double __log(x) /* wrapper log */ + double x; +#endif +{ +#ifdef _IEEE_LIBM + return __ieee754_log(x); +#else + double z; + z = __ieee754_log(x); + if(_LIB_VERSION == _IEEE_ || (x != x) || x > 0.0) return z; + if(x==0.0) + + return __kernel_standard(x,x,16); /* log(0) */ + else + return __kernel_standard(x,x,17); /* log(x<0) */ +#endif +} +weak_alias (__log, log) +#ifdef NO_LONG_DOUBLE +strong_alias (__log, __logl) +weak_alias (__log, logl) +#endif +#if LONG_DOUBLE_COMPAT(libm, GLIBC_2_0) +compat_symbol (libm, __log, logl, GLIBC_2_0); +#endif diff -urN libc25-cvstip-20070919/sysdeps/powerpc/powerpc64/power6/fpu/w_log.c libc25/sysdeps/powerpc/powerpc64/power6/fpu/w_log.c --- /dev/null 2007-05-04 05:56:50.000000000 -0500 +++ libc25/sysdeps/powerpc/powerpc64/power6/fpu/w_log.c 2008-02-29 14:27:28.000000000 -0600 @@ -0,0 +1,76 @@ +/* @(#)w_log.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#if defined(LIBM_SCCS) && !defined(lint) +static char rcsid[] = "$NetBSD: w_log.c,v 1.6 1995/05/10 20:49:33 jtc Exp $"; +#endif + +/* + * wrapper log(x) + */ + +#include +#include +#include "math_private.h" + + union double2long + { + double d; + long l; + }; + +#ifdef __STDC__ + double __log(double x) /* wrapper log */ +#else + double __log(x) /* wrapper log */ + double x; +#endif +{ +#ifdef _IEEE_LIBM + return __ieee754_log(x); +#else + double z; + volatile union double2long dx; + long lx_l; + + + /* Copy the double value x to a long to test for NaN without the + overhead of calling isnan or triggering VXSNAN. + Mark the union volatile and store x before calling __ieee754_log() + so we can access the value as a 64-bit long after + __ieee754_log(). This avoids the load-hit-store hazard when the + store queue is deep. */ + dx.d = x; + + z = __ieee754_log(x); + + lx_l = dx.l; + if(_LIB_VERSION == _IEEE_ + /* equivalent to isnan(x) */ + || ((lx_l & 0x7fffffffffffffff) > 0x7ff0000000000000) + /* equivalent to (x > 0.0) */ + || (lx_l > 0)) return z; + + if(x == 0.0) + return __kernel_standard(x,x,16); /* log(0) */ + else + return __kernel_standard(x,x,17); /* log(x<0) */ +#endif +} +weak_alias (__log, log) +#ifdef NO_LONG_DOUBLE +strong_alias (__log, __logl) +weak_alias (__log, logl) +#endif +#if LONG_DOUBLE_COMPAT(libm, GLIBC_2_0) +compat_symbol (libm, __log, logl, GLIBC_2_0); +#endif