Mail Archives: djgpp-workers/2002/06/14/19:56:53
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!
- Raw text -