Mailing-List: contact cygwin-help AT cygwin DOT com; run by ezmlm List-Subscribe: List-Archive: List-Post: List-Help: , Sender: cygwin-owner AT cygwin DOT com Mail-Followup-To: cygwin AT cygwin DOT com Delivered-To: mailing list cygwin AT cygwin DOT com Message-ID: <4319C016.2958BC77@dessent.net> Date: Sat, 03 Sep 2005 08:24:06 -0700 From: Brian Dessent MIME-Version: 1.0 To: cygwin AT cygwin DOT com, newlib AT sources DOT redhat DOT com Subject: [patch] lrint/lrintf oddity References: <20050903185145 DOT 5bb6c57e DOT cygwin-erikd AT mega-nerd DOT com> Content-Type: multipart/mixed; boundary="------------7CF6A2CE1526A4D41737FB5D" X-IsSubscribed: yes Reply-To: cygwin AT cygwin DOT com --------------7CF6A2CE1526A4D41737FB5D Content-Type: text/plain; charset=us-ascii Content-Transfer-Encoding: 7bit Erik de Castro Lopo wrote: > $ gcc -Wall -O3 lrintf_test.c -o lrintf_test > $ ./lrintf_test.exe > lrintf_test ( 3333.000000) : ok > lrint_test ( 3333.000000) : ok > lrintf_test (33333.000000) : > > Line 70: float : Incorrect sample (#252 : -8399916.000000 => 889181140). Yes, this does in fact seem to be a bug in lrintf. A reduced testcase is: #include #include int main () { float f = -8399916.0; printf ("f = %f, k = %ld\n", f, lrintf (f)); } The value -8399916.0 in binary is 1 10010110 00000000010110000101100 The exponent (j0) is 23 (150 - 127 = 23). In lrintf, the problem is the part of the algorithm that acts on the case of 23 <= j0 < 31. In this case, the mantissa is simply extracted and shifted left by (j0 - 23). This is handled by line 66 of sf_lrint.c: result = (long int) i0 << (j0 - 23); This variable 'result' is then returned, optionally negated if the sign bit was set. However, this calculation makes two errors: First, it must mask off only the mantissa from i0, removing the sign bit and exponent. Second, it neglects to add the implicit MSB "1" which is always present to the left of the decimal but not actually included in the mantissa. The attached patch fixes these two problems and fixes the poster's test case. It appears that lrint() also might suffer from the same problem: else if (j0 < (int)(8 * sizeof (long int)) - 1) { if (j0 >= 52) result = ((long int) i0 << (j0 - 20)) | (i1 << (j0 - 52)); However, this section of code will never be reached on 32 bit hardware since the first 'if' will guarantee that j0 < 31. However this looks like it's still broken for 64 bit machines, and someone should probably audit this. Brian 2005-09-03 Brian Dessent * sf_lrint.c (lrintf): Mask 'i0' correctly when extracting mantissa. --------------7CF6A2CE1526A4D41737FB5D Content-Type: text/plain; charset=us-ascii; name="newlib-lrintf.patch" Content-Transfer-Encoding: 7bit Content-Disposition: inline; filename="newlib-lrintf.patch" Index: sf_lrint.c =================================================================== RCS file: /cvs/src/src/newlib/libm/common/sf_lrint.c,v retrieving revision 1.2 diff -u -p -r1.2 sf_lrint.c --- sf_lrint.c 28 Jun 2005 17:03:18 -0000 1.2 +++ sf_lrint.c 3 Sep 2005 15:13:06 -0000 @@ -63,7 +63,7 @@ TWO23[2]={ if (j0 < -1) return 0; else if (j0 >= 23) - result = (long int) i0 << (j0 - 23); + result = (long int) (i0 & 0x7fffff | 0x800000) << (j0 - 23); else { w = TWO23[sx] + x; --------------7CF6A2CE1526A4D41737FB5D Content-Type: text/plain; charset=us-ascii -- Unsubscribe info: http://cygwin.com/ml/#unsubscribe-simple Problem reports: http://cygwin.com/problems.html Documentation: http://cygwin.com/docs.html FAQ: http://cygwin.com/faq/ --------------7CF6A2CE1526A4D41737FB5D--