Message-ID: <3D0A53C7.BF363CE0@yahoo.com>
Date: Fri, 14 Jun 2002 16:36:23 -0400
From: CBFalconer <cbfalconer AT yahoo DOT com>
Organization: Ched Research
X-Mailer: Mozilla 4.75 [en] (Win98; U)
X-Accept-Language: en
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>
	 <3D09D6D8 DOT 926022C0 AT phekda DOT freeserve DOT co DOT uk> <3D0A04A1 DOT 5F275E5C AT cyberoptics DOT com>
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit
Reply-To: djgpp-workers AT delorie DOT com

Eric Rudd wrote:
> Richard Dawe wrote:
> 
> > 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.
> 
> My libc routine is wrong -- see below.
> 
> > 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?
> 
> That formula is mathematically correct, but impractical to use for negative a.
> In practice, for negative a one computes
> 
>    a^b = (sign(a)*|a|)^b = (sign(a))^b * exp(b*ln(|a|))
> 
> The exp() term is as easy to compute when a is negative as when it's positive.
> The first term simply corrects the sign: for a>=0, the result is never negative;
> for a<0, if b is even, the result is positive; if b is odd, the result is
> negative; if b is non-integer, pow() returns NaN, since the result is complex.
> 
> Since I wrote the pow() function in libc, I figured I was in as good a position
> as anyone to figure out the cause of the discrepancy.  I took a look at my code,
> and very quickly saw the reason:  I had used the "fistl" instruction to write b
> out to memory as a long int.  I then tested its evenness by and'ing it with 1.
> This test works as long as the number is within the range of a long int, but
> fails for really big numbers such as Moshier was using, since fistl clips (or
> saturates) out-of-range numbers, instead of writing them out modulo 2^32, as I
> wish it would.  I'll take a look at my code, and see if there's a reasonably
> efficient test of evenness for big integers.  There's a "fistq" instruction,
> which writes out 64-bit integers and would presumably solve the problem, but it
> will be slower.  The reason 64-bit integers should solve the problem is that any
> doubles > 2^64 in absolute value are always even integers, since their
> significands are only 53 bits.

Here is something I wrote some time ago.  I think it does not have
that problem.

/* ======= bigpow.c by C.B. Falconer ======== */
/*   Put in public domain - use as you wish   */

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

typedef int bool;

typedef
enum error {okay = 0,
            badnumber,
            badexponent,
            xtrajunk,
            neg} anerror;

/* =========================================
 * Give help message and exit with failure */
void helpexit(char * progname)
{
    printf("Usage: %s number exponent\n", progname);
    printf("Where number and exponent are numeric values\n");
    exit(EXIT_FAILURE);
} /* helpexit */

/* ========================
 * Translate error status */
void showfault(anerror err, char *baditem)
{
    char * msg = NULL;

    switch (err) {
case badnumber:
        msg = "Invalid number value %s\n";
        break;

case badexponent:
        msg = "Invalid exponent value %s\n";
        break;

case xtrajunk:
        msg = "Extra characters in %s\n";
        break;

case neg:
        msg = "Number can not be negative %s\n";
        break;

default:
        break;
    } /* switch (err) */

    if (msg)
        printf(msg, baditem);
} /* showfault */

/* ====================================================
 * Check value is an integer, i.e. no fractional part */
bool nonintegral(double val)
{
    return ((val - (long)val) != 0.0);
} /* nonintegral */

/* =============================================
 * Check whether (known) integral value is odd */
bool odd(double val)
{
    return ((long)fabs(val) & 1);
} /* odd */

/* ============================
 * Parse number out of string */
anerror parsefails(char *arg, double *val)
{
    char * last;

    *val = strtod(arg, &last);
    if (last == arg) return badnumber;
    else if (*last)  return xtrajunk;
    else             return okay;
} /* parsefails */

/* =========================================
 * Compute power and its base 10 logarithm *
 * no must be >= 0.0                       */
void bigpow(double no,    double expo,            /* input */
            double *mant, int    *characteristic) /* output */
{
    double log;

    if (no == 0.0) {
        *mant = no;
        *characteristic = 0;
    }
    else {
        log = log10(no) * expo;
        *characteristic = (int) log;
        log = log - *characteristic;
        *mant = pow(10.0, log);
    }
} /* bigpow */

/* ========================================
 * compute and show large power of number */
int main(int argc, char ** argv)
{
    double  number, exponent;
    double  result;
    int     exp;
    anerror error;
    bool    resultneg = 0;

    if (argc < 3) helpexit(argv[0]);
    else if ((error = parsefails(argv[1], &number))) {
        showfault(error, argv[1]);
        helpexit(argv[0]);
    }
    else if ((error = parsefails(argv[2], &exponent))) {
        if (error == badnumber) error = badexponent;
        showfault(error, argv[2]);
        helpexit(argv[0]);
    }
    else if (number < 0.0) {
        if (nonintegral(exponent)) {
            showfault(neg, argv[1]);
            helpexit(argv[0]);
        }
        else {
            resultneg = odd(exponent);
            number = - number;
        }
    }
    /* input seems to be viable, number is non-negative */
    bigpow(number, exponent, &result, &exp);
    if (resultneg) result = - result;

    printf("%s ** %s = %fe%d\n",
           argv[1],
           argv[2],
           result,
           exp);
    return 0;
} /* main */
/* ============ End of bigpow.c ============= */


gcc 3.1 apparently generates the following, which I believe has
the same problem, but it is guarded by the fact that the value is
known integral.  The definition of integral in the nonintegral
function precludes it being excessively large.  I think.

000001c4 <_odd>:

/* =============================================
 * Check whether (known) integral value is odd */
bool odd(double val)
{
 1c4:   55                      push   %ebp
 1c5:   89 e5                   mov    %esp,%ebp
 1c7:   83 ec 18                sub    $0x18,%esp
 1ca:   8b 45 08                mov    0x8(%ebp),%eax
 1cd:   8b 55 0c                mov    0xc(%ebp),%edx
 1d0:   89 45 f8                mov    %eax,0xfffffff8(%ebp)
 1d3:   89 55 fc                mov    %edx,0xfffffffc(%ebp)
    return ((long)fabs(val) & 1);
 1d6:   83 ec 08                sub    $0x8,%esp
 1d9:   ff 75 fc                pushl  0xfffffffc(%ebp)
 1dc:   ff 75 f8                pushl  0xfffffff8(%ebp)
 1df:   e8 1c fe ff ff          call   0 <.text>
 1e4:   83 c4 10                add    $0x10,%esp
 1e7:   d9 7d f6                fnstcw 0xfffffff6(%ebp)
 1ea:   66 8b 45 f6             mov    0xfffffff6(%ebp),%ax
 1ee:   b4 0c                   mov    $0xc,%ah
 1f0:   66 89 45 f4             mov    %ax,0xfffffff4(%ebp)
 1f4:   d9 6d f4                fldcw  0xfffffff4(%ebp)
 1f7:   db 5d f0                fistpl 0xfffffff0(%ebp)
 1fa:   d9 6d f6                fldcw  0xfffffff6(%ebp)
 1fd:   8b 45 f0                mov    0xfffffff0(%ebp),%eax
 200:   83 e0 01                and    $0x1,%eax
} /* odd */
 203:   c9                      leave
 204:   c3                      ret
 205:   90                      nop

-- 
Chuck F (cbfalconer AT yahoo DOT com) (cbfalconer AT worldnet DOT att DOT net)
   Available for consulting/temporary embedded and systems.
   <http://cbfalconer.home.att.net>  USE worldnet address!