Sender: rich AT phekda DOT freeserve DOT co DOT uk
Message-ID: <3D09D6D8.926022C0@phekda.freeserve.co.uk>
Date: Fri, 14 Jun 2002 12:43:20 +0100
From: Richard Dawe <rich AT phekda DOT freeserve DOT co DOT uk>
X-Mailer: Mozilla 4.77 [en] (X11; U; Linux 2.2.19 i586)
X-Accept-Language: de,fr
MIME-Version: 1.0
To: djgpp-workers AT delorie DOT com
Subject: Re: [steve AT moshier DOT net: bug in fdlibm e_pow.c]
References: <200206131949 DOT g5DJndk20542 AT envy DOT delorie DOT com>
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit
Reply-To: djgpp-workers AT delorie DOT com

Hello.

DJ Delorie wrote:
> Do we have this problem also?

Yes, it looks like it. I saved the sources in pow-bug.c, adding includes for
stdlib.h, stdio.h and math.h. Then I built versions with and without -lm
against CVS using this Makefile:

---Start Makefile---
CFLAGS += -Wall -g -I/develop/djgpp.rw/include
LDFLAGS += -g -L/develop/djgpp.rw/lib

default:        pow-bug pow-bug-m

%.o:    %.c
        $(CC) $(CFLAGS) -o $@ -c $<

pow-bug:        pow-bug.o
        $(CC) $(LDFLAGS) -o $@ $<.

pow-bug-m:      pow-bug.o
        $(CC) $(LDFLAGS) -o $@ $< -lm
---End Makefile---

Here's what I get:

bash-2.04$ ./pow-bug
pow(1.0000000000000002e+00 4.5035996273704970e+15) = 2.7182818284590455e+00
pow(-1.0000000000000002e+00 4.5035996273704970e+15) = 2.7182818284590455e+00
pow(9.9999999999999978e-01 4.5035996273704970e+15) = 3.6787944117144222e-01
pow(-9.9999999999999978e-01 4.5035996273704970e+15) = 3.6787944117144222e-01
bash-2.04$ emacs Makefile
bash-2.04$ make
gcc -g -L/develop/djgpp.rw/lib -o pow-bug-m pow-bug.o -lm
bash-2.04$ ./pow-bug-m
pow(1.0000000000000002e+00 4.5035996273704970e+15) = 2.7182818284590455e+00
pow(-1.0000000000000002e+00 4.5035996273704970e+15) = 0.0000000000000000e+00
pow(9.9999999999999978e-01 4.5035996273704970e+15) = 3.6787944117144222e-01
pow(-9.9999999999999978e-01 4.5035996273704970e+15) = 0.0000000000000000e+00

So libm is affected.

After applying the patch I get:

bash-2.04$ ./pow-bug
pow(1.0000000000000002e+00 4.5035996273704970e+15) = 2.7182818284590455e+00
pow(-1.0000000000000002e+00 4.5035996273704970e+15) = 2.7182818284590455e+00
pow(9.9999999999999978e-01 4.5035996273704970e+15) = 3.6787944117144222e-01
pow(-9.9999999999999978e-01 4.5035996273704970e+15) = 3.6787944117144222e-01
bash-2.04$ ./pow-bug-m
pow(1.0000000000000002e+00 4.5035996273704970e+15) = 2.7182818284590455e+00
pow(-1.0000000000000002e+00 4.5035996273704970e+15) = -2.7182818284590455e+00
pow(9.9999999999999978e-01 4.5035996273704970e+15) = 3.6787944117144222e-01
pow(-9.9999999999999978e-01 4.5035996273704970e+15) = -3.6787944117144222e-01

Note that the results are not consistent - the sign differs in libc vs. libm.

Now, if I can still remember my maths correctly 8) :

    a ** b == exp(b * ln(a))

So if a is negative then ln(a) gives a complex number. How does a function
that deals with real numbers cope with the imaginary part of that?

Looking at the draft of the C99 standard, it has the following in section
"7.12.7.4  The pow functions" it says:

"[#2]  The  pow functions compute x raised to the power y.  A
       domain error occurs if x is negative and y is finite and not
       an  integer  value.   A  domain  error  occurs if the result
       cannot be represented when x is zero and y is less  than  or
       equal to zero.  A range error may occur."

So I don't think this test is entirely correct.

For comparison, here's what I get on RedHat Linux 6.2 with glibc-2.1.3-23:

iolanthe:~/src/tmp =] ./pow-bug 
pow(1.0000000000000002e+00 4.5035996273704970e+15) = 2.7182818278260532e+00
pow(-1.0000000000000002e+00 4.5035996273704970e+15) = -2.7182818278260532e+00
pow(9.9999999999999978e-01 4.5035996273704970e+15) = 3.6787944112862875e-01
pow(-9.9999999999999978e-01 4.5035996273704970e+15) = -3.6787944112862875e-01

I had to link in libm, to get it to link on RH6.2.

So what do we do now? I think we should commit the patch, to be consistent
with other libraries, even if the operation is outside the domain. (If there
is agreement, I'll submit a patch for libm with a changelog entry.) Also, it
looks like there may be a bug in libc. Maybe someone else could look at that?
(I'm not at all familiar with FP.)

FWIW this fix to pow() doesn't make any difference to the results of running
the Cygnus test suite.

Bye, Rich =]

-- 
Richard Dawe [ http://www.phekda.freeserve.co.uk/richdawe/ ]