The following patch fixes some of the special case
math conditions to be a bit more sane. There are also
test/example programs for both float and double. This
was done on cygwin with an older version of newlib,
but I have verified that the patch is compatible with
yesterdays CVS. Here is a description and my rational
for each change. Most of these changes make the
library more compatible with gcc constant
optimizations and other math libraries.
asin() and acos() should produce NaN when given
something outside the -1 to 1 range. The documentation
said that was what was supposed to happen, but it was
returning 0.0 instead. Other libraries also produce
NaN.
log() and log10() of negative numbers should return
NaN. The existing code was returning minus infinity,
but NaN is a better choice when you consider that NaN
often represents Not A real Number. The documentation
in w_log.c was also changed to describe the new
behavior. Other libraries produce NaN.
The pow() function had a number of defects though some
of them may be subjective.
pow(NaN, 0.0) now produces 1.0 (the wrapper catches
this case and was setting it to NaN). The underlying
code/gcc constant optimizations/other libraries
produce 1.0.
pow(1.0, NaN) now produces 1.0. instead of NaN. Other
libraries and gcc constant optimizations produce 1.0.
I split special case 3 (e_pow.c) into a and b parts to
clarify this behavior.
pow(1.0, +-infinity) now produce 1.0 instead of NaN.
Other libraries and gcc constant optimizations produce
1.0. I split special case 9 (e_pow.c) into a and b
parts to clarify this behavior.
pow(-1.0, +-infinity) are unchanged (produce NaN)
since infinity is neither odd nor even. Some libraries
produce 1.0 which is correct from a magnitude point of
view, but the indeterminate sign is not addressed. To
me it is wrong to return 1.0 without something to
indicate that there is uncertainty in the exact value.
In the spirit of the previous paragraph pow(-x<-1.0,
+infinity) now also produces NaN. Some libraries
produce infinity which is the correct magnitude. The
problem is the indeterminate sign.
I hope this is useful,
Cary
#!/bin/sh
# This is a shell archive (produced by GNU sharutils
4.6.3).
# To extract the files from this archive, save it to
some FILE, remove
# everything before the `#!/bin/sh' line above, then
type `sh FILE'.
#
lock_dir=_sh01984
#
# Existing files will *not* be overwritten, unless
`-c' is specified.
#
# This shar contains:
# length mode name
# ------ ----------
------------------------------------------
# 6984 -rw-r--r-- newlib-math.patch
# 2361 -rw-r--r-- chk_double.c
# 2406 -rw-r--r-- chk_float.c
#
MD5SUM=${MD5SUM-md5sum}
f=`${MD5SUM} --version | egrep '^md5sum
.*(core|text)utils'`
test -n "${f}" && md5check=true || md5check=false
${md5check} || \
echo 'Note: not verifying md5sums. Consider
installing GNU coreutils.'
save_IFS="${IFS}"
IFS="${IFS}:"
gettext_dir=FAILED
locale_dir=FAILED
first_param="$1"
for dir in $PATH
do
if test "$gettext_dir" = FAILED && test -f
$dir/gettext \
&& ($dir/gettext --version >/dev/null 2>&1)
then
case `$dir/gettext --version 2>&1 | sed 1q` in
*GNU*) gettext_dir=$dir ;;
esac
fi
if test "$locale_dir" = FAILED && test -f $dir/shar
\
&& ($dir/shar --print-text-domain-dir >/dev/null
2>&1)
then
locale_dir=`$dir/shar --print-text-domain-dir`
fi
done
IFS="$save_IFS"
if test "$locale_dir" = FAILED || test "$gettext_dir"
= FAILED
then
echo=echo
else
TEXTDOMAINDIR=$locale_dir
export TEXTDOMAINDIR
TEXTDOMAIN=sharutils
export TEXTDOMAIN
echo="$gettext_dir/gettext -s"
fi
if (echo "testing\c"; echo 1,2,3) | grep c >/dev/null
then if (echo -n test; echo 1,2,3) | grep n >/dev/null
then shar_n= shar_c='
'
else shar_n=-n shar_c= ; fi
else shar_n= shar_c='\c' ; fi
f=shar-touch.$$
st1=200112312359.59
st2=123123592001.59
st2tr=123123592001.5 # old SysV 14-char limit
st3=1231235901
if touch -am -t ${st1} ${f} >/dev/null 2>&1 && \
test ! -f ${st1} && test -f ${f}; then
shar_touch='touch -am -t $1$2$3$4$5$6.$7 "$8"'
elif touch -am ${st2} ${f} >/dev/null 2>&1 && \
test ! -f ${st2} && test ! -f ${st2tr} && test -f
${f}; then
shar_touch='touch -am $3$4$5$6$1$2.$7 "$8"'
elif touch -am ${st3} ${f} >/dev/null 2>&1 && \
test ! -f ${st3} && test -f ${f}; then
shar_touch='touch -am $3$4$5$6$2 "$8"'
else
shar_touch=:
echo
${echo} 'WARNING: not restoring timestamps.
Consider getting and'
${echo} 'installing GNU `touch'\'', distributed in
GNU coreutils...'
echo
fi
rm -f ${st1} ${st2} ${st2tr} ${st3} ${f}
#
if test ! -d ${lock_dir}
then : ; else ${echo} 'lock directory '${lock_dir}'
exists'
exit 1
fi
if mkdir ${lock_dir}
then ${echo} 'x - created lock directory
`'${lock_dir}\''.'
else ${echo} 'x - failed to create lock directory
`'${lock_dir}\''.'
exit 1
fi
# ============= newlib-math.patch ==============
if test -f 'newlib-math.patch' && test "$first_param"
!= -c; then
${echo} 'x -SKIPPING newlib-math.patch (file already
exists)'
else
${echo} 'x - extracting newlib-math.patch (text)'
sed 's/^X//' << 'SHAR_EOF' > 'newlib-math.patch' &&
--- newlib/libm/math/e_pow.c.orig 2003-10-20
11:46:37.001000000 -0700
+++ newlib/libm/math/e_pow.c 2007-04-25
11:23:11.737375000 -0700
@@ -25,13 +25,16 @@
X * Special cases:
X * 1. (anything) ** 0 is 1
X * 2. (anything) ** 1 is itself
- * 3. (anything) ** NAN is NAN
+ * 3a. (anything) ** NAN is NAN except
+ * 3b. +1 ** NAN is 1
X * 4. NAN ** (anything except 0) is NAN
- * 5. +-(|x| > 1) ** +INF is +INF
+ * 5a. +(|x| > 1) ** +INF is +INF
+ * 5b. -(|x| > 1) ** +INF is NAN
X * 6. +-(|x| > 1) ** -INF is +0
X * 7. +-(|x| < 1) ** +INF is +0
X * 8. +-(|x| < 1) ** -INF is +INF
- * 9. +-1 ** +-INF is NAN
+ * 9a. +1 ** +-INF is 1
+ * 9b. -1 ** +-INF is NAN
X * 10. +0 ** (+anything except 0, NAN)
is +0
X * 11. -0 ** (+anything except 0, NAN, odd integer)
is +0
X * 12. +0 ** (-anything except 0, NAN)
is +INF
@@ -117,10 +120,11 @@
X /* y==zero: x**0 = 1 */
X if((iy|ly)==0) return one;
X
- /* +-NaN return x+y */
+ /* x|y==NaN return NaN unless x==1 then return 1
*/
X if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0))
||
X iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
- return x+y;
+ if(((ix-0x3ff00000)|lx)==0) return one;
+ else return nan("");
X
X /* determine if y is an odd int when x < 0
X * yisint = 0 ... y is not an integer
@@ -146,9 +150,10 @@
X if(ly==0) {
X if (iy==0x7ff00000) { /* y is +-inf */
X if(((ix-0x3ff00000)|lx)==0)
- return y - y; /* inf**+-1 is NaN */
- else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf
= inf,0 */
- return (hy>=0)? y: zero;
+ if(hx>=0) return one; /* +1**+-inf = 1 */
+ else return nan(""); /* -1**+-inf = NaN */
+ else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf
= +x=inf;-x=NaN,0 */
+ return (hy>=0)? (hx>=0)? y: nan(""): zero;
X else /* (|x|<1)**-,+inf = inf,0 */
X return (hy<0)?-y: zero;
X }
--- newlib/libm/math/ef_pow.c.orig 2002-06-13
16:03:01.001000000 -0700
+++ newlib/libm/math/ef_pow.c 2007-04-25
11:23:47.987375000 -0700
@@ -75,10 +75,11 @@
X /* y==zero: x**0 = 1 */
X if(FLT_UWORD_IS_ZERO(iy)) return one;
X
- /* +-NaN return x+y */
+ /* x|y==NaN return NaN unless x==1 then return 1
*/
X if(FLT_UWORD_IS_NAN(ix) ||
X FLT_UWORD_IS_NAN(iy))
- return x+y;
+ if(ix==0x3f800000) return one;
+ else return nanf("");
X
X /* determine if y is an odd int when x < 0
X * yisint = 0 ... y is not an integer
@@ -98,9 +99,10 @@
X /* special value of y */
X if (FLT_UWORD_IS_INFINITE(iy)) { /* y is +-inf */
X if (ix==0x3f800000)
- return y - y; /* inf**+-1 is NaN */
- else if (ix > 0x3f800000)/* (|x|>1)**+-inf =
inf,0 */
- return (hy>=0)? y: zero;
+ if(hx>=0) return one; /* +1**+-inf = 1 */
+ else return nanf(""); /* -1**+-inf = NaN */
+ else if (ix > 0x3f800000)/* (|x|>1)**+-inf =
+x=inf;-x=NaN,0 */
+ return (hy>=0)? (hx>=0)? y: nanf(""): zero;
X else /* (|x|<1)**-,+inf = inf,0 */
X return (hy<0)?-y: zero;
X }
--- newlib/libm/math/w_acos.c.orig 2003-10-20
11:46:37.001000000 -0700
+++ newlib/libm/math/w_acos.c 2007-04-24
18:22:15.659250000 -0700
@@ -101,7 +101,7 @@
X exc.name = "acos";
X exc.err = 0;
X exc.arg1 = exc.arg2 = x;
- exc.retval = 0.0;
+ exc.retval = nan("");
X if (_LIB_VERSION == _POSIX_)
X errno = EDOM;
X else if (!matherr(&exc)) {
--- newlib/libm/math/w_asin.c.orig 2007-04-24
18:08:55.753000000 -0700
+++ newlib/libm/math/w_asin.c 2007-04-24
18:09:58.034250000 -0700
@@ -104,7 +104,7 @@
X exc.name = "asin";
X exc.err = 0;
X exc.arg1 = exc.arg2 = x;
- exc.retval = 0.0;
+ exc.retval = nan("");
X if(_LIB_VERSION == _POSIX_)
X errno = EDOM;
X else if (!matherr(&exc)) {
--- newlib/libm/math/w_log.c.orig 2007-04-24
18:28:28.753000000 -0700
+++ newlib/libm/math/w_log.c 2007-04-24
18:33:35.893625000 -0700
@@ -44,7 +44,7 @@
X RETURNS
X Normally, returns the calculated value. When <[x]>
is zero, the
X returned value is <<-HUGE_VAL>> and <<errno>> is set
to <<ERANGE>>.
-When <[x]> is negative, the returned value is
<<-HUGE_VAL>> and
+When <[x]> is negative, the returned value is NaN
(not a number) and
X <<errno>> is set to <<EDOM>>. You can control the
error behavior via
X <<matherr>>.
X
@@ -105,6 +105,7 @@
X else if (!matherr(&exc)) {
X errno = EDOM;
X }
+ exc.retval = nan("");
X }
X if (exc.err != 0)
X errno = exc.err;
--- newlib/libm/math/w_log10.c.orig 2007-04-24
18:28:51.190500000 -0700
+++ newlib/libm/math/w_log10.c 2007-04-24
18:29:28.799875000 -0700
@@ -103,6 +103,7 @@
X else if (!matherr(&exc)) {
X errno = EDOM;
X }
+ exc.retval = nan("");
X }
X if (exc.err != 0)
X errno = exc.err;
--- newlib/libm/math/w_pow.c.orig 2003-10-20
11:46:37.001000000 -0700
+++ newlib/libm/math/w_pow.c 2007-04-24
19:13:19.503000000 -0700
@@ -93,7 +93,7 @@
X exc.err = 0;
X exc.arg1 = x;
X exc.arg2 = y;
- exc.retval = x;
+ exc.retval = 1.0;
X if (_LIB_VERSION == _IEEE_ ||
X _LIB_VERSION == _POSIX_) exc.retval = 1.0;
X else if (!matherr(&exc)) {
--- newlib/libm/math/wf_acos.c.orig 2000-02-17
11:39:51.001000000 -0800
+++ newlib/libm/math/wf_acos.c 2007-04-24
18:22:38.331125000 -0700
@@ -40,7 +40,7 @@
X exc.name = "acosf";
X exc.err = 0;
X exc.arg1 = exc.arg2 = (double)x;
- exc.retval = 0.0;
+ exc.retval = nan("");
X if (_LIB_VERSION == _POSIX_)
X errno = EDOM;
X else if (!matherr(&exc)) {
--- newlib/libm/math/wf_asin.c.orig 2000-02-17
11:39:51.001000000 -0800
+++ newlib/libm/math/wf_asin.c 2007-04-24
18:11:11.315500000 -0700
@@ -42,7 +42,7 @@
X exc.name = "asinf";
X exc.err = 0;
X exc.arg1 = exc.arg2 = (double)x;
- exc.retval = 0.0;
+ exc.retval = nan("");
X if(_LIB_VERSION == _POSIX_)
X errno = EDOM;
X else if (!matherr(&exc)) {
--- newlib/libm/math/wf_log.c.orig 2007-04-24
18:28:37.456125000 -0700
+++ newlib/libm/math/wf_log.c 2007-04-24
18:29:37.284250000 -0700
@@ -63,6 +63,7 @@
X else if (!matherr(&exc)) {
X errno = EDOM;
X }
+ exc.retval = nan("");
X }
X if (exc.err != 0)
X errno = exc.err;
--- newlib/libm/math/wf_log10.c.orig 2007-04-24
18:28:59.315500000 -0700
+++ newlib/libm/math/wf_log10.c 2007-04-24
18:29:46.596750000 -0700
@@ -64,6 +64,7 @@
X else if (!matherr(&exc)) {
X errno = EDOM;
X }
+ exc.retval = nan("");
X }
X if (exc.err != 0)
X errno = exc.err;
--- newlib/libm/math/wf_pow.c.orig 2000-02-17
11:39:51.001000000 -0800
+++ newlib/libm/math/wf_pow.c 2007-04-24
19:14:46.346750000 -0700
@@ -43,7 +43,7 @@
X exc.err = 0;
X exc.arg1 = (double)x;
X exc.arg2 = (double)y;
- exc.retval = x;
+ exc.retval = 1.0;
X if (_LIB_VERSION == _IEEE_ ||
X _LIB_VERSION == _POSIX_) exc.retval = 1.0;
X else if (!matherr(&exc)) {
SHAR_EOF
(set 20 07 04 25 12 01 51 'newlib-math.patch'; eval
"$shar_touch") &&
chmod 0644 'newlib-math.patch'
if test $? -ne 0
then ${echo} 'restore of newlib-math.patch failed'
fi
if ${md5check}
then (
${MD5SUM} -c >/dev/null 2>&1 || ${echo}
'newlib-math.patch: MD5 check failed'
) << SHAR_EOF
c518ebfd6bbf7c42f2b033b24ebc8d0b newlib-math.patch
SHAR_EOF
else
test `LC_ALL=C wc -c < 'newlib-math.patch'` -ne 6984
&& \
${echo} 'restoration warning: size of
newlib-math.patch is not 6984'
fi
fi
# ============= chk_double.c ==============
if test -f 'chk_double.c' && test "$first_param" !=
-c; then
${echo} 'x -SKIPPING chk_double.c (file already
exists)'
else
${echo} 'x - extracting chk_double.c (text)'
sed 's/^X//' << 'SHAR_EOF' > 'chk_double.c' &&
#include <math.h>
#include <stdio.h>
X
int main()
{
X double var, nan, inf, minf, one, zero;
X inf = 1.0/0.0;
X minf = -1.0 * inf;
X nan = sqrt(-1.0);
X one = 1.0;
X zero = 0.0;
X
X printf("Using nan = %f, inf = %f and -inf =
%f.\n\n", nan, inf, minf);
X
X // Check the 1 ** inf case.
X var = pow(1.1, inf);
X printf("1.1 ** inf is %f.\n", var);
X
X var = pow(1.0, inf);
X printf("1.0 ** inf is %f", var);
X var = pow(one, inf);
X printf(", %f.\n", var);
X
X var = pow(0.9, inf);
X printf("0.9 ** inf is %f.\n\n", var);
X
X
X // Check the 1 ** -inf case.
X var = pow(1.1, minf);
X printf("1.1 ** -inf is %f.\n", var);
X
X var = pow(1.0, minf);
X printf("1.0 ** -inf is %f", var);
X var = pow(one, minf);
X printf(", %f.\n", var);
X
X var = pow(0.9, minf);
X printf("0.9 ** -inf is %f.\n\n", var);
X
X
X // Check the -1 ** inf case.
X var = pow(-1.1, inf);
X printf("-1.1 ** inf is %f.\n", var);
X
X var = pow(-1.0, inf);
X printf("-1.0 ** inf is %f", var);
X var = pow(-one, inf);
X printf(", %f.\n", var);
X
X var = pow(-0.9, inf);
X printf("-0.9 ** inf is %f.\n\n", var);
X
X
X // Check the -1 ** -inf case.
X var = pow(-1.1, minf);
X printf("-1.1 ** -inf is %f.\n", var);
X
X var = pow(-1.0, minf);
X printf("-1.0 ** -inf is %f", var);
X var = pow(-one, minf);
X printf(", %f.\n", var);
X
X var = pow(-0.9, minf);
X printf("-0.9 ** -inf is %f.\n\n", var);
X
X
X // Check the nan cases.
X var = pow(1.0, nan);
X printf("1.0 ** nan is %f", var);
X var = pow(one, nan);
X printf(", %f.\n", var);
X
X var = pow(nan, 0.0);
X printf("nan ** 0.0 is %f", var);
X var = pow(nan, zero);
X printf(", %f.\n\n", var);
X
X // Check log and log10.
X var = log(one);
X printf("log of 1.0 is %f\n", var);
X var = log(zero);
X printf("log of 0.0 is %f\n", var);
X var = log(-one);
X printf("log of -1.0 is %f\n\n", var);
X
X var = log10(one);
X printf("log10 of 1.0 is %f\n", var);
X var = log10(zero);
X printf("log10 of 0.0 is %f\n", var);
X var = log10(-one);
X printf("log10 of -1.0 is %f\n\n", var);
X
X // Check asin and acos.
X var = asin(zero);
X printf("asin of 0.0 is %f\n", var);
X var = asin(one);
X printf("asin of 1.0 is %f\n", var);
X var = asin(2*one);
X printf("asin of 2.0 is %f\n\n", var);
X
X var = acos(zero);
X printf("acos of 0.0 is %f\n", var);
X var = acos(one);
X printf("acos of 1.0 is %f\n", var);
X var = acos(2*one);
X printf("acos of 2.0 is %f\n", var);
X
X return 0;
}
SHAR_EOF
(set 20 07 04 25 11 49 07 'chk_double.c'; eval
"$shar_touch") &&
chmod 0644 'chk_double.c'
if test $? -ne 0
then ${echo} 'restore of chk_double.c failed'
fi
if ${md5check}
then (
${MD5SUM} -c >/dev/null 2>&1 || ${echo}
'chk_double.c: MD5 check failed'
) << SHAR_EOF
031055eb310f4ce425ef60700e01e2ff chk_double.c
SHAR_EOF
else
test `LC_ALL=C wc -c < 'chk_double.c'` -ne 2361 && \
${echo} 'restoration warning: size of chk_double.c
is not 2361'
fi
fi
# ============= chk_float.c ==============
if test -f 'chk_float.c' && test "$first_param" != -c;
then
${echo} 'x -SKIPPING chk_float.c (file already
exists)'
else
${echo} 'x - extracting chk_float.c (text)'
sed 's/^X//' << 'SHAR_EOF' > 'chk_float.c' &&
#include <math.h>
#include <stdio.h>
X
int main()
{
X double var, nan, inf, minf, one, zero;
X inf = 1.0/0.0;
X minf = -1.0 * inf;
X nan = sqrt(-1.0);
X one = 1.0;
X zero = 0.0;
X
X printf("Using nan = %f, inf = %f and -inf =
%f.\n\n", nan, inf, minf);
X
X // Check the 1 ** inf case.
X var = powf(1.1, inf);
X printf("1.1 ** inf is %f.\n", var);
X
X var = powf(1.0, inf);
X printf("1.0 ** inf is %f", var);
X var = powf(one, inf);
X printf(", %f.\n", var);
X
X var = powf(0.9, inf);
X printf("0.9 ** inf is %f.\n\n", var);
X
X
X // Check the 1 ** -inf case.
X var = powf(1.1, minf);
X printf("1.1 ** -inf is %f.\n", var);
X
X var = powf(1.0, minf);
X printf("1.0 ** -inf is %f", var);
X var = powf(one, minf);
X printf(", %f.\n", var);
X
X var = powf(0.9, minf);
X printf("0.9 ** -inf is %f.\n\n", var);
X
X
X // Check the -1 ** inf case.
X var = powf(-1.1, inf);
X printf("-1.1 ** inf is %f.\n", var);
X
X var = powf(-1.0, inf);
X printf("-1.0 ** inf is %f", var);
X var = powf(-one, inf);
X printf(", %f.\n", var);
X
X var = powf(-0.9, inf);
X printf("-0.9 ** inf is %f.\n\n", var);
X
X
X // Check the -1 ** -inf case.
X var = powf(-1.1, minf);
X printf("-1.1 ** -inf is %f.\n", var);
X
X var = powf(-1.0, minf);
X printf("-1.0 ** -inf is %f", var);
X var = powf(-one, minf);
X printf(", %f.\n", var);
X
X var = powf(-0.9, minf);
X printf("-0.9 ** -inf is %f.\n\n", var);
X
X
X // Check the nan cases.
X var = powf(1.0, nan);
X printf("1.0 ** nan is %f", var);
X var = powf(one, nan);
X printf(", %f.\n", var);
X
X var = powf(nan, 0.0);
X printf("nan ** 0.0 is %f", var);
X var = powf(nan, zero);
X printf(", %f.\n\n", var);
X
X // Check logf and log10f.
X var = logf(one);
X printf("log of 1.0 is %f\n", var);
X var = logf(zero);
X printf("log of 0.0 is %f\n", var);
X var = logf(-one);
X printf("log of -1.0 is %f\n\n", var);
X
X var = log10f(one);
X printf("log10f of 1.0 is %f\n", var);
X var = log10f(zero);
X printf("log10f of 0.0 is %f\n", var);
X var = log10f(-one);
X printf("log10f of -1.0 is %f\n\n", var);
X
X // Check asinf and acosf.
X var = asinf(zero);
X printf("asinf of 0.0 is %f\n", var);
X var = asinf(one);
X printf("asinf of 1.0 is %f\n", var);
X var = asinf(2*one);
X printf("asinf of 2.0 is %f\n\n", var);
X
X var = acosf(zero);
X printf("acosf of 0.0 is %f\n", var);
X var = acosf(one);
X printf("acosf of 1.0 is %f\n", var);
X var = acosf(2*one);
X printf("acosf of 2.0 is %f\n", var);
X
X return 0;
}
SHAR_EOF
(set 20 07 04 25 11 49 51 'chk_float.c'; eval
"$shar_touch") &&
chmod 0644 'chk_float.c'
if test $? -ne 0
then ${echo} 'restore of chk_float.c failed'
fi
if ${md5check}
then (
${MD5SUM} -c >/dev/null 2>&1 || ${echo}
'chk_float.c: MD5 check failed'
) << SHAR_EOF
eb034af9ab404965a895ae166f42ba2a chk_float.c
SHAR_EOF
else
test `LC_ALL=C wc -c < 'chk_float.c'` -ne 2406 && \
${echo} 'restoration warning: size of chk_float.c
is not 2406'
fi
fi
if rm -fr ${lock_dir}
then ${echo} 'x - removed lock directory
`'${lock_dir}\''.'
else ${echo} 'x - failed to remove lock directory
`'${lock_dir}\''.'
exit 1
fi
exit 0
__________________________________________________
Do You Yahoo!?
Tired of spam? Yahoo! Mail has the best spam protection around
http://mail.yahoo.com