Date: Sat, 22 Mar 2003 11:29:18 +1100 (EST) From: Bruce Evans <bde@zeta.org.au> To: Till Riedel <till@f111.hadiko.de> Cc: freebsd-current@FreeBSD.ORG, David Schultz <das@FreeBSD.ORG> Subject: Re: libm problem Message-ID: <20030322111233.F4471@gamplex.bde.org> In-Reply-To: <20030321235237.GA8097@f111.hadiko.de> References: <20030318173051.GA2322@f111.hadiko.de> <20030319131317.GA670@HAL9000.homeunix.com> <20030321235237.GA8097@f111.hadiko.de>
next in thread | previous in thread | raw e-mail | index | archive | help
On Sat, 22 Mar 2003, Till Riedel wrote: > > > res=pow((float)base,(float)dim); > this is actually not a smart thing. it was cut and paste from libvorbis. > pow is a function for doubles. if you i use powf everything works fine. This should work OK. The casts to float should have no effect except to lose precision is the operands don't fit in about 24 bits. OTOH, powf() is known to be broken. I haven't committed a fix like the following for too long: %%% Index: e_powf.c =================================================================== RCS file: /home/ncvs/src/lib/msun/src/e_powf.c,v retrieving revision 1.9 diff -u -2 -r1.9 e_powf.c --- e_powf.c 17 Jun 2002 15:28:59 -0000 1.9 +++ e_powf.c 17 Jun 2002 15:41:06 -0000 @@ -45,5 +45,9 @@ lg2 = 6.9314718246e-01, /* 0x3f317218 */ lg2_h = 6.93145752e-01, /* 0x3f317200 */ +#if 0 +lg2_l = 1.42860677e-06, /* 0x35bfbe8e */ +#else lg2_l = 1.42860654e-06, /* 0x35bfbe8c */ +#endif ovt = 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */ cp = 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */ @@ -160,7 +164,7 @@ s_h = s; GET_FLOAT_WORD(is,s_h); - SET_FLOAT_WORD(s_h,is&0xfffff000); + SET_FLOAT_WORD(s_h,is&0xfffc0000); /* t_h=ax+bp[k] High */ - SET_FLOAT_WORD(t_h,((ix>>1)|0x20000000)+0x0040000+(k<<21)); + SET_FLOAT_WORD(t_h,((ix>>1)|0x20000000)+0x00400000+(k<<21)); t_l = ax - (t_h-bp[k]); s_l = v*((u-s_h*t_h)-s_h*t_l); @@ -229,5 +233,5 @@ t = p_l+p_h; GET_FLOAT_WORD(is,t); - SET_FLOAT_WORD(t,is&0xfffff000); + SET_FLOAT_WORD(t,is&0xfffffff8); u = t*lg2_h; v = (p_l-(t-p_h))*lg2+t*lg2_l; %%% The change to lg2_l just fixes a tiny error. We gain precision by representing log(2) as approximately (lg2_h + lg2_l) in infinite precision. The bits in lg2_l are relatively uncritical since most of them are to give more than float precision. But they sometimes matter. The Cygnus extensions to support float precision have a number of off by 1 or 2 bit errors converting the low values. I fixed a couple of these, but IIRC there was only one that caused an error that was reported by ucbtest (not this one). The changes in the magic numbers are too fix larger errors. I don't rememebr anything else about them except that the errors were reported by ucbtest and the changes made ucbtest happy. Bruce To Unsubscribe: send mail to majordomo@FreeBSD.org with "unsubscribe freebsd-current" in the body of the message
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20030322111233.F4471>