From owner-freebsd-current Fri Mar 21 16:29:31 2003 Delivered-To: freebsd-current@freebsd.org Received: from mx1.FreeBSD.org (mx1.freebsd.org [216.136.204.125]) by hub.freebsd.org (Postfix) with ESMTP id D9FAC37B401; Fri, 21 Mar 2003 16:29:29 -0800 (PST) Received: from mailman.zeta.org.au (mailman.zeta.org.au [203.26.10.16]) by mx1.FreeBSD.org (Postfix) with ESMTP id 0323143FBD; Fri, 21 Mar 2003 16:29:26 -0800 (PST) (envelope-from bde@zeta.org.au) Received: from katana.zip.com.au (katana.zip.com.au [61.8.7.246]) by mailman.zeta.org.au (8.9.3/8.8.7) with ESMTP id LAA13411; Sat, 22 Mar 2003 11:29:19 +1100 Date: Sat, 22 Mar 2003 11:29:18 +1100 (EST) From: Bruce Evans X-X-Sender: bde@gamplex.bde.org To: Till Riedel Cc: freebsd-current@FreeBSD.ORG, David Schultz Subject: Re: libm problem In-Reply-To: <20030321235237.GA8097@f111.hadiko.de> Message-ID: <20030322111233.F4471@gamplex.bde.org> References: <20030318173051.GA2322@f111.hadiko.de> <20030319131317.GA670@HAL9000.homeunix.com> <20030321235237.GA8097@f111.hadiko.de> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII Sender: owner-freebsd-current@FreeBSD.ORG Precedence: bulk List-ID: List-Archive: (Web Archive) List-Help: (List Instructions) List-Subscribe: List-Unsubscribe: X-Loop: FreeBSD.ORG 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