Bug in printf ?

Bug in printf ?

[haro]

[-- Attachment #1: Type: Text/Plain, Size: 1037 bytes --]

Hi all,

I'm seeing small bug in printf implementation.
May be a bug in newlib ??

    % uname -a
    CYGWIN_NT-5.1 wbbrown 1.5.17(0.129/4/2) 2005-05-25 19:38 i686 unknown unknown Cygwin
    % /bin/printf "%.2f\n" 0.105
    0.10
    % /bin/printf "%.2f\n" 0.115
    0.12
    % /bin/printf "%.2f\n" 0.125
    0.12
    % /bin/printf "%.2f\n" 0.135
    0.14
    % /bin/printf "%.2f\n" 0.145
    0.14
    % /bin/printf "%.2f\n" 0.155
    0.15
    %
    % cat aaa.c
    #include <stdio.h>
    main()
    {
        printf("%0.2f\n", 0.105);
        printf("%0.2f\n", 0.115);
        printf("%0.2f\n", 0.125);
        printf("%0.2f\n", 0.135);
    }
    % make aaa
    gcc     aaa.c   -o aaa
    % ./aaa
    0.10
    0.12
    0.12
    0.14
    %

Thanks in advance,
  Haro
=-----------------------------------------------------------------------
           _ _    Munehiro (haro) Matsuda
 -|- /_\  |_|_|   KGT Inc.
 /|\ |_|  |_|_|   2-8-8 Shinjuku, Shinjuku-ku Tokyo 160-0022, Japan
                  Tel: +81-3-3225-0767  Fax: +81-3-3225-0740

[-- Attachment #2: cygcheck.out --]
[-- Type: Text/Plain, Size: 18522 bytes --]


Cygwin Configuration Diagnostics
Current System Time: Thu Jun 30 19:40:44 2005

Windows XP Professional Ver 5.1 Build 2600 Service Pack 1

Path:	D:\cygwin\usr\local\bin
	E:\home\haro\bin
	D:\cygwin\bin
	D:\cygwin\usr\local\bin
	D:\cygwin\bin
	D:\cygwin\bin
	D:\cygwin\usr\X11R6\bin
	D:\cygwin\usr\local\sbin
	D:\cygwin\usr\sbin
	D:\cygwin\sbin
	C:\WINDOWS\system32
	C:\WINDOWS
	C:\WINDOWS\System32\Wbem
	D:\cygwin\bin
	.\
	D:\cygwin\usr\local\Meadow\1.15\bin
	D:\cygwin\usr\X11R6\bin

Output from D:\cygwin\bin\id.exe (nontsec)
UID: 1003(haro)     GID: 513(なし)
0(root)             513(なし)           544(Administrators) 545(Users)

Output from D:\cygwin\bin\id.exe (ntsec)
UID: 1003(haro)     GID: 513(なし)
0(root)             513(なし)           544(Administrators) 545(Users)

SysDir: C:\WINDOWS\System32
WinDir: C:\WINDOWS

PWD = `/tmp'
CYGWIN = `binmode tty ntsec'
HOME = `/home/haro'
USER = `haro'
MAKE_MODE = `unix'

HOMEPATH = `\Documents and Settings\haro'
APPDATA = `C:\Documents and Settings\haro\Application Data'
PROCESSOR_IDENTIFIER = `x86 Family 6 Model 8 Stepping 6, GenuineIntel'
TERM = `cygwin'
WINDIR = `C:\WINDOWS'
USERDOMAIN = `WBBROWN'
ALLUSERSPROFILE = `C:\Documents and Settings\All Users'
OS = `Windows_NT'
COMMONPROGRAMFILES = `C:\Program Files\Common Files'
TEMP = `/tmp'
PROCESSOR_LEVEL = `6'
SYSTEMDRIVE = `C:'
CLIENTNAME = `Console'
USERPROFILE = `C:\Documents and Settings\haro'
TZ = `JST-9'
LOGONSERVER = `\\WBBROWN'
PROCESSOR_ARCHITECTURE = `x86'
SHLVL = `2'
PATHEXT = `.COM;.EXE;.BAT;.CMD;.VBS;.VBE;.JS;.JSE;.WSF;.WSH'
HOMEDRIVE = `C:'
COMSPEC = `C:\WINDOWS\system32\cmd.exe'
SYSTEMROOT = `C:\WINDOWS'
TMP = `/tmp'
PROCESSOR_REVISION = `0806'
PROGRAMFILES = `C:\Program Files'
NUMBER_OF_PROCESSORS = `1'
SESSIONNAME = `Console'
COMPUTERNAME = `WBBROWN'
LOGNAME = `haro'
OLDPWD = `/home/haro'
MANPATH = `/usr/local/man:/usr/man:/usr/share/man:/usr/autotool/devel/man::/usr/ssl/man'
INFOPATH = `/usr/local/info:/usr/info:/usr/share/info:/usr/autotool/devel/info:/usr/autotool/stable/info:'
SHELL = `/bin/zsh'
CVS_RSH = `/bin/ssh'
PS1 = `%{

[-- Attachment #3: Type: text/plain, Size: 218 bytes --]

--
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/

RE: Bug in printf ?

[Dave Korn]

----Original Message----
>From: haro@kgt.co.jp
>Sent: 30 June 2005 11:43

> Hi all,
> 
> I'm seeing small bug in printf implementation.
> May be a bug in newlib ??

  Confirmed.  I'm on it!

    cheers,
      DaveK
-- 
Can't think of a witty .sigline today....


--
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/

SV: Bug in printf ?

[Peter J. Acklam]

[-- Attachment #1: Type: text/plain, Size: 125 bytes --]

haro@kgt.co.jp wrote:

> I'm seeing small bug in printf implementation.

What bug? I didn't see anything unexpected.

Peter


[-- Attachment #2: Type: text/plain, Size: 218 bytes --]

--
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/

Re: SV: Bug in printf ?

[haro]

From: "Peter J. Acklam" <pjacklam at online no>
Date: Thu, 30 Jun 2005 13:06:51 +0200 (CEST)
::
::> I'm seeing small bug in printf implementation.
::
::What bug? I didn't see anything unexpected.
::
::Peter

Hi Peter,

How come "0.125" gets printed as "0.12", and not "1.3"?

Haro



--
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/

RE: SV: Bug in printf ?

[Dave Korn]

----Original Message----
>From: haro@kgt.co.jp
>Sent: 30 June 2005 12:21

> From: "Peter J. Acklam" <pjacklam at online no>
> Date: Thu, 30 Jun 2005 13:06:51 +0200 (CEST)
>>> 
>>>> I'm seeing small bug in printf implementation.
>>> 
>>> What bug? I didn't see anything unexpected.
>>> 
>>> Peter
> 
> Hi Peter,
> 
> How come "0.125" gets printed as "0.12", and not "1.3"?
> 
> Haro


  Absolutely, there's a rounding error of some sort.  Compare the difference
when compiling the testcase with -mno-cygwin (i.e. using mingw maths lib):

dk@mace /test/signed/lr2> cat aaa.c
#include <stdio.h>
int main(int argc, const char **argv)
{
    printf("%0.2f\n", 0.105);
    printf("%0.2f\n", 0.115);
    printf("%0.2f\n", 0.125);
    printf("%0.2f\n", 0.135);
    return 0;
}

dk@mace /test/signed/lr2> gcc -O0 -g aaa.c -o aaa-c
dk@mace /test/signed/lr2> ./aaa-c.exe
0.10
0.12
0.12
0.14
dk@mace /test/signed/lr2> gcc -O0 -g -mno-cygwin aaa.c -o aaa-m
dk@mace /test/signed/lr2> ./aaa-m.exe
0.11
0.12
0.13
0.14
dk@mace /test/signed/lr2>


    cheers,
      DaveK
-- 
Can't think of a witty .sigline today....


--
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/

Re: SV: Bug in printf ?

[haro]

OOps,

::::What bug? I didn't see anything unexpected.
::::
::::Peter
::
::Hi Peter,
::
::How come "0.125" gets printed as "0.12", and not "1.3"?
                                                    ^^^ "0.13", off cource ;-)

::Haro
::



--
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/

Re: Bug in printf ?

[Peter J. Acklam]

[-- Attachment #1: Type: text/plain, Size: 1381 bytes --]

haro@kgt.co.jp wrote:

> ::How come "0.125" gets printed as "0.12", and not "1.3"?
>                                                     ^^^ "0.13", off cource ;-)

Dealing with integers illustrates the matter more clearly. When
the decimal value is exactly 0.5, then printf should round to the
nearest *even* integer, as far as I know, so you should get

   0   -> 0
   0.5 -> 0
   1   -> 1
   1.5 -> 2
   2   -> 2
   2.5 -> 2
   3   -> 3
   3.5 -> 4
   4   -> 4
   4.5 -> 4
   5   -> 5

Now I realize that Cygwin's printf doesn't get it right, because

   seq 0 .5 5 | while read x; do printf '%-3s -> %.0f\n' $x $x; done

gives

   0   -> 0
   0.5 -> 1
   1   -> 1
   1.5 -> 1
   2   -> 2
   2.5 -> 2
   3   -> 3
   3.5 -> 3
   4   -> 4
   4.5 -> 4
   5   -> 5

Both Solaris' /bin/printf and Perl's printf() give the right
output.

In your examples there is also an additional issue, which is how
numerical values are represented.  Here is Perl:

   seq .105 .01 .155 | while read x; do printf '%-5s -> %.2f\n' $x $x; done

   0.105 -> 0.10
   0.115 -> 0.12
   0.125 -> 0.12
   0.135 -> 0.14
   0.145 -> 0.14
   0.155 -> 0.15

They all seem correct expect the last one.  The reason why 0.155
becomes 0.15 and not 0.16 is that 0.115 can not be represented
exactly.  It is represented as a number which is slightly
*smaller* than 0.155, so it becomes 0.15 after rounding.

Peter


[-- Attachment #2: Type: text/plain, Size: 218 bytes --]

--
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/

SV: SV: Bug in printf ?

[Peter J. Acklam]

[-- Attachment #1: Type: text/plain, Size: 558 bytes --]

Dave Korn [dave.korn@artimi.com] wrote:

> haro@kgt.co.jp wrote:
>
> > How come "0.125" gets printed as "0.12", and not "1.3"?
> 
> Absolutely, there's a rounding error of some sort.

For what it's worth: My Sunblade 100 running Solaris 9 has
Solaris' /bin/printf and GNU's printf as /usr/local/bin/printf
and both give "0.12", not "0.13".

I am quite sure it is the printf returning "0.13" that is buggy.
Note that the value "0.125" can be represented exactly with IEEE
double precision arithmetic, so inexact representation is not the
matter here.

Peter


[-- Attachment #2: Type: text/plain, Size: 218 bytes --]

--
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/

Re: Bug in printf ?

[Peter Mueller]

Hi!


Peter J. Acklam wrote:

>haro@kgt.co.jp wrote:
>
>  
>
>>::How come "0.125" gets printed as "0.12", and not "1.3"?
>>                                                    ^^^ "0.13", off cource ;-)
>>    
>>
>
>Dealing with integers illustrates the matter more clearly. When
>the decimal value is exactly 0.5, then printf should round to the
>nearest *even* integer, as far as I know, so you should get
>  
>
Never heard of such a rule!
As I learned it in school, 0.5 should allways be rounded to 1 .


Bye, Peter
WOTLmade

--
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/

RE: Bug in printf ?

[Dave Korn]

----Original Message----
>From: Peter Mueller
>Sent: 30 June 2005 13:18

> Hi!
> 
> 
> Peter J. Acklam wrote:
> 
>> haro@NOTYOUTOO! wrote:

   ^^^^^^^^^^^^^  Please avoid quoting peoples email addresses in your
replies, it causes them to suffer more spamming.

>>>>> How come "0.125" gets printed as "0.12", and not "1.3"?
>>>                                                    ^^^ "0.13", off
>>> cource ;-) 
>>> 
>>> 
>> 
>> Dealing with integers illustrates the matter more clearly. When
>> the decimal value is exactly 0.5, then printf should round to the
>> nearest *even* integer, as far as I know, so you should get
>> 
>> 
> Never heard of such a rule!
> As I learned it in school, 0.5 should allways be rounded to 1 .


  I guess they didn't teach you IEEE-754-compliant floating point maths in
school.  You may find it informative to google some/any/all/none of the
following terms:

dk@mace /artimi/firmware> grep FE_ /usr/include/mingw/fenv.h  | grep 0x0.00
#define FE_TONEAREST    0x0000
#define FE_DOWNWARD     0x0400
#define FE_UPWARD       0x0800
#define FE_TOWARDZERO   0x0c00
dk@mace /artimi/firmware>


    cheers,
      DaveK
-- 
Can't think of a witty .sigline today....


--
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/

RE: SV: Bug in printf ?

[Dave Korn]

----Original Message----
>From: Peter J. Acklam
>Sent: 30 June 2005 13:13

> Dave Korn [dave.korn@NONONOATHOUSANDTIMESNO] wrote:
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  http://cygwin.com/acronyms#PCYMTNQREAIYR please!  TIA!

> 
>> haro@OOPSNOTAGAIN wrote:
>> 
>>> How come "0.125" gets printed as "0.12", and not "1.3"?
>> 
>> Absolutely, there's a rounding error of some sort.
> 
> For what it's worth: My Sunblade 100 running Solaris 9 has
> Solaris' /bin/printf and GNU's printf as /usr/local/bin/printf
> and both give "0.12", not "0.13".
> 
> I am quite sure it is the printf returning "0.13" that is buggy.

  Have you considered that your sunblade might be operating in a different
rounding mode, by default?

> Note that the value "0.125" can be represented exactly with IEEE
> double precision arithmetic, so inexact representation is not the
> matter here.

  I never suggested it was.  Nonetheless, I would imagine that printf may
well work under-the-hood by shifting the desired decimal places above the
point, rounding to integer, and printing out that number.

    cheers,
      DaveK
-- 
Can't think of a witty .sigline today....


--
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/

RE: SV: Bug in printf ?

[Dave Korn]

----Original Message----
>From: Peter J. Acklam
>Sent: 30 June 2005 13:13

> Dave Korn [dave.korn@NONONOATHOUSANDTIMESNO] wrote:
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  http://cygwin.com/acronyms#PCYMTNQREAIYR please!  TIA!

> 
>> haro@OOPSNOTAGAIN wrote:
>> 
>>> How come "0.125" gets printed as "0.12", and not "1.3"?
>> 
>> Absolutely, there's a rounding error of some sort.
> 
> For what it's worth: My Sunblade 100 running Solaris 9 has
> Solaris' /bin/printf and GNU's printf as /usr/local/bin/printf
> and both give "0.12", not "0.13".
> 
> I am quite sure it is the printf returning "0.13" that is buggy.

  Have you considered that your sunblade might be operating in a different
rounding mode, by default?

> Note that the value "0.125" can be represented exactly with IEEE
> double precision arithmetic, so inexact representation is not the
> matter here.

  I never suggested it was.  Nonetheless, I would imagine that printf may
well work under-the-hood by shifting the desired decimal places above the
point, rounding to integer, and printing out that number.

    cheers,
      DaveK
-- 
Can't think of a witty .sigline today....


--
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/

Re: Bug in printf ?

[haro]

::Dealing with integers illustrates the matter more clearly. When
::the decimal value is exactly 0.5, then printf should round to the
::nearest *even* integer, as far as I know, so you should get

Hi Peter,

Thank you for your explanation. I didn't know, until now, that
rounding should be done to the "nearest *even* integer".
Need to learn every day. ;-)

::Now I realize that Cygwin's printf doesn't get it right, because

Wow, what a coincidence. ;-)

Anyway, thanks for Peter and DaveK for looking into the issue.

Thanks
  Haro
=-----------------------------------------------------------------------
           _ _    Munehiro (haro) Matsuda
 -|- /_\  |_|_|   KGT Inc.
 /|\ |_|  |_|_|   2-8-8 Shinjuku, Shinjuku-ku Tokyo 160-0022, Japan
                  Tel: +81-3-3225-0767  Fax: +81-3-3225-0740



--
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/

SV: SV: Bug in printf ?

[Peter J. Acklam]

[-- Attachment #1: Type: text/plain, Size: 887 bytes --]

Dave Korn wrote:

> Have you considered that your sunblade might be operating
> in a different rounding mode, by default?

I didn't know there were different rounding modes.
I thought everyone used so-called "unbiased rounding",
so I'm sorry for adding confusion.

> I would imagine that printf may well work under-the-hood
> by shifting the desired decimal places above the point,
> rounding to integer, and printing out that number.

Could be.  I haven't dug into those matters, but your
suggestion certainly makes sense.

Anyway, the following line

   seq -4.5 1 4.5 | while read x; do printf '%4s -> %3.0f\n' $x $x; done

which gives

   -4.5 ->  -4
   -3.5 ->  -3
   -2.5 ->  -2
   -1.5 ->  -1
   -0.5 ->  -1
    0.5 ->   1
    1.5 ->   1
    2.5 ->   2
    3.5 ->   3
    4.5 ->   4

suggests that printf rounds towards zero, but something
funny happens for -0.5 and 0.5.

Peter


[-- Attachment #2: Type: text/plain, Size: 218 bytes --]

--
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/

Re: SV: Bug in printf ?

[Brian Dessent]

Dave Korn wrote:

>   Absolutely, there's a rounding error of some sort.  Compare the difference
> when compiling the testcase with -mno-cygwin (i.e. using mingw maths lib):

Isn't this just a case of the Cygwin math library choosing "round to
even" and the MSVCRT/mingw library choosing "0.5 always rounds up"?  The
Goldberg paper covers this, and in fact it claims the IEEE standard
requires round to even for all results of addition, subtraction,
multiplication, and division.  Sooo, doesn't look like a bug at all to
me.

Brian

--
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/

Re: Bug in printf ?

[Yitzchak Scott-Thoennes]

On Thu, Jun 30, 2005 at 09:44:57PM +0900, haro@kgt.co.jp wrote:
> ::Dealing with integers illustrates the matter more clearly. When
> ::the decimal value is exactly 0.5, then printf should round to the
> ::nearest *even* integer, as far as I know, so you should get
> 
> Hi Peter,
> 
> Thank you for your explanation. I didn't know, until now, that
> rounding should be done to the "nearest *even* integer".
> Need to learn every day. ;-)

This is done to prevent systemic bias.

However, note that
http://www.opengroup.org/onlinepubs/009695399/functions/printf.html says:
"The low-order digit shall be rounded in an implementation-defined manner."
so while people have come to expect unbiased rounding, it isn't actually
mandatory AIUI under vanilla SUSv3.

Also see "What Every Computer Scientist Should Know About
Floating-Point Arithmetic" at (among other places)
http://docs.sun.com/source/806-3568/ncg_goldberg.html.

--
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/

RE: SV: Bug in printf ?

[Dave Korn]

[-- Attachment #1: Type: text/plain, Size: 843 bytes --]

----Original Message----
>From: Brian Dessent
>Sent: 30 June 2005 16:58

> Dave Korn wrote:
> 
>>   Absolutely, there's a rounding error of some sort.  Compare the
>> difference when compiling the testcase with -mno-cygwin (i.e. using
>> mingw maths lib): 
> 
> Isn't this just a case of the Cygwin math library choosing "round to
> even" and the MSVCRT/mingw library choosing "0.5 always rounds up"?  

  No, I really don't think that's the case, in fact I reckon they both
default to the same roundingmode, and the discrepancy remains the same in
any rounding mode.  Have a go with the attached version of the testcase,
which allows you to choose any rounding mode you like at runtime: just pass
a hex constant from the set (0, 0x400, 0x800, 0xc00) on the command line.

    cheers,
      DaveK
-- 
Can't think of a witty .sigline today....

[-- Attachment #2: aa2.c --]
[-- Type: text/plain, Size: 1610 bytes --]


#include <stdio.h>
#include <math.h>

#ifdef __MINGW32__
#include <fenv.h>
#endif

#ifndef __MINGW32__

#define FE_TONEAREST	0x0000
#define FE_DOWNWARD	0x0400
#define FE_UPWARD	0x0800
#define FE_TOWARDZERO	0x0c00

 /* 7.6.3.2
    The fesetround function establishes the rounding direction
    represented by its argument round. If the argument is not equal
    to the value of a rounding direction macro, the rounding direction
    is not changed.  */

int fesetround (int mode)
{
  unsigned short _cw;
  if ((mode & ~(FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO))
      != 0)
    return -1;
  __asm__ volatile ("fnstcw %0;": "=m" (_cw));
  _cw &= ~(FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO);
  _cw |= mode;
  __asm__ volatile ("fldcw %0;" : : "m" (_cw));
  return 0;
}

/* 7.6.3.1
   The fegetround function returns the value of the rounding direction
   macro representing the current rounding direction.  */

int
fegetround (void)
{
  unsigned short _cw;
  __asm__ ("fnstcw %0;" : "=m" (_cw));
  return _cw
          & (FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO);
}
#endif

int main(int argc, const char **argv)
{
    if (argc >= 2)
    {
    int nrm, n;
        n = sscanf (argv[1], "%i", &nrm);
        if (n == 1) switch (nrm)
        {
            case 0:
            case 0x400:
            case 0x800:
            case 0xc00:
                fesetround (nrm);
        }
    }

    printf ("Rmode $%08x\n", fegetround ());

    printf("%0.2f\n", 0.105);
    printf("%0.2f\n", 0.115);
    printf("%0.2f\n", 0.125);
    printf("%0.2f\n", 0.135);

    return 0;
}



[-- Attachment #3: Type: text/plain, Size: 218 bytes --]

--
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/

RE: SV: Bug in printf ?

[Dave Korn]

----Original Message----
>From: Peter J. Acklam
>Sent: 30 June 2005 13:13

> Dave Korn [dave.korn@artimiOHHHHHNOOOOOOOOO!] wrote:
                      ^^^^^^^^^^^^^^^^^^^^^^^^
> 
>> haro@kgtPLEASETAKEMORECARE! wrote:
       ^^^^^^^^^^^^^^^^^^^^^^

  Peter, please do http://cygwin.com/acronyms#PCYMTNQREAIYR before you get
me loaded down with even more spam than I already suffer, thank you!

>>> How come "0.125" gets printed as "0.12", and not "1.3"?
>> 
>> Absolutely, there's a rounding error of some sort.
> 
> For what it's worth: My Sunblade 100 running Solaris 9 has
> Solaris' /bin/printf and GNU's printf as /usr/local/bin/printf
> and both give "0.12", not "0.13".
> 
> I am quite sure it is the printf returning "0.13" that is buggy.
> Note that the value "0.125" can be represented exactly with IEEE
> double precision arithmetic, so inexact representation is not the
> matter here.


  Progress report:

  Stepping through the code at the weekend, I followed the 0.105 case as far
down as ldtoa_r, where I observed that although there was code that would
round 0.105 up to 0.11 if asked for only two sig.figs, the loop that
generates successive digits was coming up with the sequence
'0.10499999999.....', and because the rounding only looks at one digit
beyond the requested s.f., it 'correctly' rounds 0.104 down to 0.10.

  I say 'correctly' in quotes, because it's the correct way to round 0.104,
but it's not in general correct to round a number by truncating it and then
rounding the truncated version.  That's at least part of the problem.  I
haven't checked yet whether 0.105 should be an exact representation in fp,
but it could also be the case that some accuracy has been lost by the time
it arrives there.

[  I notice that mingw has a maths library in a subdir called mingwex, and
it looks like this is a version of the same code (SLM's cephes lib, but the
ldtoa implementation appears to be something that's been added subsequently
AFAICT) and would show the same problem, but when you actually compile the
testcase with -mno-cygwin the generated exe calls the msvcrt version of
printf (which DTRT, and which is the version you'd expect it to call), so
I'm not sure why it's there or what it's relevant to.  ]

    cheers,
      DaveK
-- 
Can't think of a witty .sigline today....


--
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/

Re: SV: Bug in printf ?

[Lev Bishop]

[-- Attachment #1: Type: text/plain, Size: 2430 bytes --]

I don't think there is any bug here. This is what I've seen from a
little digging:

1) cygwin strtod rounds to even, with about DECIMAL_DIG (==21) digits
precision, as recommended by 7.20.1.3 of WG14/N843. (It acts strange
when the rounding mode is not round to nearest, but since newlib
doesn't provide fesetround(), that's your own problem if you insist on
changing rounding mode by hand, as I understand it).

2) cygwin gcc rounds with *lots* more precision than DECIMAL_DIG,
about 50 digits. It rounds towards zero. These behaviours are a little
different than the recommendation in 6.4.4.2 Ibid., which says
constants should be converted as if by strtod() at runtime, but its
not a big difference.

3) cygwin printf correctly rounds to even, even out to about 41 digits.

4) I have no idea what mingw is doing, but it's different to the
above. Gcc constructs the same double precision constants as on cygwin
but strtod() is different and seems to have less precision, and
printf() seems to work with about 16 digits precision. At a glance I'd
say it rounds to zero, with exactly the precision of

On 05/07/05, Dave Korn wrote:
>  Stepping through the code at the weekend, I followed the 0.105 case as far
> down as ldtoa_r, where I observed that although there was code that would
> round 0.105 up to 0.11 if asked for only two sig.figs, the loop that
> generates successive digits was coming up with the sequence
> '0.10499999999.....', and because the rounding only looks at one digit
> beyond the requested s.f., it 'correctly' rounds 0.104 down to 0.10.

The closest double to 0.105 is actually
0.10499999999999999611421941381195210851728916168212890625, and this
is what printf() is (correctly) handed, and what it (correctly) rounds
to 0.10.
 
>  I say 'correctly' in quotes, because it's the correct way to round 0.104,
> but it's not in general correct to round a number by truncating it and then
> rounding the truncated version.

But I don't think that "truncating it and then rounding the truncated
version" is what printf() does. If we give it 0.125 and 0.375, which
are both representible exactly, then it rounds to 0.12 and 0.38, ie
one goes up and one goes down, to round to even. You can see it's not
simply truncating if you give it 0.1250000000000001, which gets
rounded to 0.13 not 0.12 as it would for truncation.

To summarize: cygwin does the right thing

[-- Attachment #2: test.c --]
[-- Type: text/plain, Size: 3032 bytes --]

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>

#ifdef __MINGW32__
#include <fenv.h>
#endif


// show the hex and decimal interpretations of the double
union U{double d;unsigned char c[8];};
void display(union U x)
{
  for(int b=7;b>=0;b--)printf("%02x",x.c[b]);
  printf(" %.21f",x.d);
}


int main(int argc, const char**argv)
{

#ifdef __MINGW32__
  if (argc >= 2)
    {
    int nrm, n;
        n = sscanf (argv[1], "%i", &nrm);
        if (n == 1) switch (nrm)
        {
            case 0:
            case 0x400:
            case 0x800:
            case 0xc00:
                fesetround (nrm);
        }
    }

  printf ("Rmode $%08x\n", fegetround ());
#endif

  printf("DECIMAL_DIG=%d\n",DECIMAL_DIG);

#define PASTE(x) {.s=#x,.d1=strtod(#x,0),.d2=x}
  struct S {const char*s;union U d1,d2;}q[]=
    {
      // test the rounding with the last digits...
      // first right on the boundary, to test round-to-even
      PASTE(1.0000000000000000000000000000000000000000000000000000000), // 1
      PASTE(1.0000000000000001110223024625156540423631668090820312500), // 1 + 2^-53
      PASTE(1.0000000000000002220446049250313080847263336181640625000), // 1 + 2^-52
      PASTE(1.0000000000000003330669073875469621270895004272460937500), // 1 + 2^-52 +2^-53
      {0},
      // now just above the boundary
      PASTE(1.00000000000000000001),
      PASTE(1.00000000000000011103),
      PASTE(1.00000000000000022205),
      PASTE(1.00000000000000033307),
      {0},
      // now just below the boundary
      PASTE(0.99999999999999999999),
      PASTE(1.00000000000000011101),
      PASTE(1.00000000000000022204),
      PASTE(1.00000000000000033301),
      {0},
      PASTE(0.105),
    };
 
  printf("                 1.23456789012345678901234567890\n"); // for counting sig figs
  for(int a=0;a<sizeof(q)/sizeof(struct S);a++)
    {
      if(!q[a].s){printf("\n");continue;}
      display(q[a].d1);
      printf(" ");
      display(q[a].d2);
      printf("\n");
    }
  
  printf("\n");

  union U w[]={0.105,0.115,0.125,0.135,0.145,0.155,0.375};
  for(int a=0;a<sizeof(w)/sizeof(union U);a++)
    {
      display(w[a]);
      printf("\n");
    }

  printf("%.40f\n",0.105);

  // here's some numbers with an exact representation in double
  printf("%.20f %.2f\n",.125,.125); 
  printf("%.20f %.2f\n",.375,.375); 
  printf("%.20f %.13f\n",0.95367431640625000000,0.95367431640625000000);
  printf("%.20f %.17f\n",3.5096855163574218750,3.5096855163574218750); // 1 + 657899/2^18
  printf("%.20f %.15f\n",-1.5035324096679687500,-1.5035324096679687500);

  // how many digits precision does gcc use when making double constants?
  display((union U){.d=1.000000000000000333066907387546962127089500427246});
}

// Useful Mathematica stuff...
/*
  Table[{x/1000.,
            BaseForm[FromDigits[RealDigits[#2, 2] // First // Rest, 2], 16],
            N[#2*2^(-53 + #[[2]]),
              41]} &[#,
        Ceiling[FromDigits[#[[1]], 2]/
            2]] &[RealDigits[x/1000, 2, 54]], {x, 105, 155, 10}]

*/

[-- Attachment #3: Type: text/plain, Size: 218 bytes --]

--
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/

Re: SV: Bug in printf ?

[Lev Bishop]

On 06/07/05, Lev Bishop wrote:
> 4) I have no idea what mingw is doing, but it's different to the
> above. Gcc constructs the same double precision constants as on cygwin
> but strtod() is different and seems to have less precision, and
> printf() seems to work with about 16 digits precision. At a glance I'd
> say it rounds to zero, with exactly the precision of
... the precision of a raw double.

Ie, it looks like the cygwin build works slightly nicer than mingw
with plenty of extra precision for conversions to and from decimal,
and round-to-even for runtime conversions, which is a good thing.

--
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/

RE: SV: Bug in printf ?

[Dave Korn]

----Original Message----
>From: Lev Bishop
>Sent: 06 July 2005 09:15

> On 06/07/05, Lev Bishop wrote:
>> 4) I have no idea what mingw is doing, but it's different to the
>> above. Gcc constructs the same double precision constants as on cygwin
>> but strtod() is different and seems to have less precision, and
>> printf() seems to work with about 16 digits precision. At a glance I'd
>> say it rounds to zero, with exactly the precision of
> ... the precision of a raw double.
> 
> Ie, it looks like the cygwin build works slightly nicer than mingw
> with plenty of extra precision for conversions to and from decimal,
> and round-to-even for runtime conversions, which is a good thing.


  Thanks a great deal for looking at this; I'm not any kind of an FP expert,
I just have a grasp of the basics.  I've tried the testcase on an RH system
and yes, glibc and newlib agree here.  What mingw does is simple: it calls
M$'s implementation in the msvcrt library, so blame the beast of redmond for
the technically-incorrect-yet-reasonable-seeming behaviour.


    cheers,
      DaveK
-- 
Can't think of a witty .sigline today....


--
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/