Enhancements for mathmodule

@tim-one

BPO	1640
Nosy	@tim-one, @mdickinson, @tiran
Dependencies	bpo-1635: Float patch for inf and nan on Windows (and other platforms)
Files	trunk_float_math_combined.patch trunk_float_math_combined2.patch trunk_pymath_hyberbolic_complex2.patch invhyp.c

^{Note: these values reflect the state of the issue at the time it was migrated and might not reflect the current state.}

Show more details

GitHub fields:

assignee = None
closed_at = <Date 2008-05-01.21:18:02.604>
created_at = <Date 2007-12-17.15:48:44.018>
labels = ['extension-modules', 'type-feature']
title = 'Enhancements for mathmodule'
updated_at = <Date 2008-05-01.21:18:02.603>
user = 'https://github.com/tiran'

bugs.python.org fields:

activity = <Date 2008-05-01.21:18:02.603>
actor = 'christian.heimes'
assignee = 'none'
closed = True
closed_date = <Date 2008-05-01.21:18:02.604>
closer = 'christian.heimes'
components = ['Extension Modules']
creation = <Date 2007-12-17.15:48:44.018>
creator = 'christian.heimes'
dependencies = ['1635']
files = ['8978', '8985', '9251', '9252']
hgrepos = []
issue_num = 1640
keywords = ['patch']
message_count = 56.0
messages = ['58694', '58698', '58699', '58702', '58706', '58735', '58749', '58758', '58770', '58786', '59115', '59137', '59138', '59152', '59153', '59154', '59155', '59156', '59165', '59166', '59211', '60257', '60259', '60264', '61290', '61306', '61308', '61309', '61310', '61311', '61366', '61367', '61370', '61371', '61372', '61376', '61385', '61386', '61394', '61398', '61399', '61401', '61402', '61403', '61414', '61419', '61423', '61431', '61441', '61442', '61443', '61444', '61474', '61479', '61498', '66047']
nosy_count = 5.0
nosy_names = ['tim.peters', 'mark.dickinson', 'Rhamphoryncus', 'christian.heimes', 'gmcastil']
pr_nums = []
priority = 'normal'
resolution = 'accepted'
stage = None
status = 'closed'
superseder = None
type = 'enhancement'
url = 'https://bugs.python.org/issue1640'
versions = ['Python 2.6', 'Python 3.0']

The patch adds several small enhancements to the math module and pyport.h.

Py_MATH_PI and Py_MATH_E in long double precision
Py_IS_NAN and Py_IS_INFINITY use isnan() and isinf() functions were
available (checked by configure)
isnan and isinf for the math module
Bessel (1st and 2nd kind), erf, erfc, lgamma function for math module.
They are defined in almost (or all?) math libraries.

Together with http://bugs.python.org/issue1635 it implements most of PEP-754.

Cool! If there's a move to add functions to the math module, there are
some others that are part of C99 (but not C89), would be good to have, and
that I'd consider more fundamental than the Bessel, error, gamma
functions; for example, the inverse hyperbolic trig functions (acosh,
asinh, atanh), log1p, expm1, copysign.

Minor typo. Should be IEEE:
"Return the sign of an int, long or float. On platforms with full IEE
754\n\"

Mark Dickinson wrote:

Mark Dickinson added the comment:

Cool! If there's a move to add functions to the math module, there are
some others that are part of C99 (but not C89), would be good to have, and
that I'd consider more fundamental than the Bessel, error, gamma
functions; for example, the inverse hyperbolic trig functions (acosh,
asinh, atanh), log1p, expm1, copysign.

I've added the inverse hyperbolic, log1p and expm1. copysign is too low
level but I've added sign(x) -> -1/0/+1. It uses copysign() where
available so you can emulate copysign(x, y) with sign(x) * y.

Do you want some more functions?
http://www.dinkumware.com/manuals/?manual=compleat&page=math.html

Feel free to provide a patch :)

Christian

Combined patch for bpo-1635 and bpo-1640.

UPDATES:

Fixed Py_NAN macro
Fixed typo
Added inverse hyperbolic trigonometric functions
Added expm1 and log1p

Missing:

doc updates
some real unit tests for the a*h and other new functions

The patch fixes some compatibility problems on Windows and implements
several methods like acosh in C when the system's math lib doesn't
provide the functions.

Unfortunately, implementing asinh, acosh, atanh correctly is not as simple as just using the formulas
(this is one of the reasons that it's useful to be able to get at the library definitions). For
example, the formula

asinh(x) = log(x + sqrt(1.+(x*x)))

has significant error (i.e. much more than just a few ulps) when x is small (so that the argument of
the log is close to 1), and when x is large and negative (so that the addition involves adding two
numbers of almost equal magnitude but opposite sign). It also overflows unnecessarily at around
1e154 (on an IEEE754 system), as a result of the intermediate calculation x*x overflowing, and it
fails to give the correct signs on zero arguments. (asinh(0.) = 0., asinh(-0.) = -0.)

So either some serious work is required here, or code should be borrowed in an appropriately legal
fashion from some other library, or those few people who don't have asinh, acosh, etc. already in
libm will have to live without.

Mark

i suggest abandoning any attempts at implementing math ourselves. I
also suggest not combining this with bpo-1635 but reviewing and
(eventually) applying the latter first.

Guido van Rossum wrote:

Guido van Rossum added the comment:

i suggest abandoning any attempts at implementing math ourselves. I
also suggest not combining this with bpo-1635 but reviewing and
(eventually) applying the latter first.

How do you like a compromise? I rip out the new functions for the math
module except isnan(), isinf() and sign(). The new patch would contain
the inf and nan roundtrip plus the helper methods.

Mark is right. My naive functions are a numerical nightmare for small
and large numbers. I copied them from my math handbook 'cause I couldn't
find solid implementations in my numerical handbook. It contains all
sort of algorithms for approximations, ODE, PDE, FFT but no sample of
the inverse hyperbolic functions.

We could reuse the functions from the BSD libm. They are used in
numerous C runtime libraries (BSD, Mac OS X, uclibc).
http://packages.e.kth.se/common/src/os/NetBSD/1.3.2/lib/libm/src/

Christian

The trunk_math_sign_inf_nan patch contains just three new method
isnan(), isinf() and sign(). It also fixes a minor numerical issue with
one function that did small / (small / large) instead of small * (large
/ small).

One nit: you added a blank line to the end of test_math.py.
This will cause the checkin to fail. :-)

One bigger issue: the sign() function doesn't seem to work properly for
nans. E.g. on Linux I get:

>>> inf = 1e1000
>>> nan = inf/inf
>>> mnan = -nan
>>> math.sign(nan)
-1
>>> math.sign(mnan)
1

IOW a positive nan is considered negative and vice versa. (This is
probably due to the way nans defy testing, always returning false.)

I'm also curious why math.sign(0.0) returns 1 -- this is going to cause
a lot of confusion.

This is also missing doc patches.

Guido van Rossum wrote:

Guido van Rossum added the comment:

One nit: you added a blank line to the end of test_math.py.
This will cause the checkin to fail. :-)

*grr* stupid white space check hook

One bigger issue: the sign() function doesn't seem to work properly for
nans. E.g. on Linux I get:

>>> inf = 1e1000
>>> nan = inf/inf
>>> mnan = -nan
>>> math.sign(nan)
-1
>>> math.sign(mnan)
1

IOW a positive nan is considered negative and vice versa. (This is
probably due to the way nans defy testing, always returning false.)

If I recall the definition correctly NaNs don't have a sign. The content
of the sign bit is not defined for NaNs. I could fix the sign but it's
just eye candy and a waste of CPU cycles. IMO it would be more
appropriate to return 0 for NaNs instead of +1 or -1.

I'm also curious why math.sign(0.0) returns 1 -- this is going to
cause a lot of confusion.

math.sign(0.0) == 1 and math.sign(-0.0) == -1 is the main purpose of the
sign() function. Most ordinary users are still going to use x > 0.0 or x
< 0.0 instead of math.sign(). math.sign() is for the power users who
need to determinate whether an operation hits the 0 from the left or
right side of the number line.

Because 0.0 == -0.0 it's not possible the distinguish a positive from a
negative zero with comparison operations. if x == 0.0 and
str(x).startswith("-") was the only existing way to detect a negative zero.

Christian

> One nit: you added a blank line to the end of test_math.py.
> This will cause the checkin to fail. :-)

*grr* stupid white space check hook

No, you edited a line that didn't need editing. :-)

> One bigger issue: the sign() function doesn't seem to work properly for
> nans. E.g. on Linux I get:
>
>>>> inf = 1e1000
>>>> nan = inf/inf
>>>> mnan = -nan
>>>> math.sign(nan)
> -1
>>>> math.sign(mnan)
> 1
>
> IOW a positive nan is considered negative and vice versa. (This is
> probably due to the way nans defy testing, always returning false.)

If I recall the definition correctly NaNs don't have a sign. The content
of the sign bit is not defined for NaNs. I could fix the sign but it's
just eye candy and a waste of CPU cycles. IMO it would be more
appropriate to return 0 for NaNs instead of +1 or -1.

Perhaps you recall wrong, as negating the nan returns one with the
opposite sign? This seems to indicate that there are positive and
negative nans.

> I'm also curious why math.sign(0.0) returns 1 -- this is going to
> cause a lot of confusion.

math.sign(0.0) == 1 and math.sign(-0.0) == -1 is the main purpose of the
sign() function. Most ordinary users are still going to use x > 0.0 or x
< 0.0 instead of math.sign(). math.sign() is for the power users who
need to determinate whether an operation hits the 0 from the left or
right side of the number line.

Because 0.0 == -0.0 it's not possible the distinguish a positive from a
negative zero with comparison operations. if x == 0.0 and
str(x).startswith("-") was the only existing way to detect a negative zero.

Hm, OK, but then passing a zero of some other type (e.g. int) should
also return +1 as the sign. I also think the function's name should be
changed, because I (and I assume many others) have grown up with a
sign() function that essentially returns cmp(x, 0.0).

Perhaps it would be better to have a function math.isneg() that
returns True for -0.0 and anything smaller and False for +0.0 and
anything larger. It could also return the proper sign of a nan.

I suggest that you check in the isinf() and isnan() functions and
their tests once you have docs for them, and hold off on the sign()
function until we've got agreement on it.

Sorry: those lines should have been:

double temp = 1-z.real;
printf("temp is %.19e\n", temp);
r.real = log1p(4.*z.real/(temp*temp + ay*ay))/4.;

Mark Dickinson wrote:

Mark Dickinson added the comment:

Christian: would it be possible for you to tell me what the argument of the
log1p call is on Windows in that last branch of c_atanh, for the testcase
atanh0022.

On OS X I get:

Input to log1p is 3.2451855365842669075e+32

It's hard to imagine that there's anything wrong with log1p here, since all
it does for a large input like this is compute log(1+x).

You got me wrong (and I didn't explain it properly). All complex
functions pass the test. math.atanh() fails. I think my implementation
of Python/pymath.c:atanh() doesn't return the right value for arguments
almost 1.0.

Christian

Sorry: I should have read more carefully. So math.atanh works on Linux
but is producing some strange results on Windows.

It's still rather puzzling though. I still suspect that it's the
argument to log1p that's coming out wrong rather than the result.

Mark Dickinson wrote:

Sorry: I should have read more carefully. So math.atanh works on Linux
but is producing some strange results on Windows.

It's still rather puzzling though. I still suspect that it's the
argument to log1p that's coming out wrong rather than the result.

It uses t = 0.5 * log1p((x + x) / (1.0 - x)) for t > 0.5. I presume the
culprit is in "2 / x" where x is almost 0. Do you have an idea how we
can increase the accuracy for values nearly 1.?

Christian

I disabled the tests for atanh0022 and atanh0023.

Other changes:
Added math.log1p + tests
Added SUN license to Doc/licenses.rst
Added docs to Doc/library/math.rst

The problem with atanh is that it should be using absx throughout, rather than x.

So in "if (absx < 0.5)" branch and the following branch, replace all occurrences of x with
absx, and it should work.

One other comment: asinh shouldn't directly set errno for a NaN. It should do the same as
acosh and atanh: return x+x.

This makes asinh(float("nan")) return a nan, which makes it consistent with acosh and
atanh, consistent with the way that Linux and OS X behave, and consistent with the other
single-argument functions in math.

I think asinh should also return x+x for x an infinity. This again should make it
consistent with the way that the libm asinh works on OS X and Linux.

Also, for the C-level routines, atanh(1.0) and atanh(-1.0) should definitely return
infinity and -infinity (and probably set errno as well.)

Note that this is not an argument about what Python should do: Python will still raise a
ValueError for atanh(1.0) and atanh(-1.0). But the atanh is supposed to be a drop-in
replacement for the libm atanh, on those platforms where it's missing. And the C99
standard is clear about return values, even though it's not useful when it comes to
deciding whether to set errno or not.

One more comment: the Py_IS_NAN test in acosh is likely never reached, since the
comparison x >= two_pow_p28 will probably evaluate to 0 when x is a NaN.

Also: -1.0 shouldn't be returned at this level to indicate an error; these are pure C functions
we're writing---the Python wrappers, and Python error conventions, apply further up the chain
somewhere.

Just so I'm doing something a little more constructive than yelling criticisms from the sideline,
I've attached a file invhyp.c that's how I think the C functions should look. Some notes:

it doesn't touch errno, but lets the platform decide how to handle errors (i.e. produce a
special value/set errno/signal a floating-point exception/some combination of these). This will
make the asinh, acosh, atanh functions behave in the same way that the regular libm functions
behave on any platform. So e.g. if a particular platform is used to setting errno for domain
errors like sqrt(-1), it'll do so for atanh/asinh/acosh. And another platform that signals a
floating-point exception for sqrt(-1) will do the same for atanh(3).
I've left in the "huge+x > 1.0" test in asinh; it's kinda pointless, but also fairly harmless.
I think its only point is to make sure that the inexact flag gets set where appropriate, and
Python doesn't much care about the inexact flag.
I'm reasonably sure that the test "x == 1." in acosh is safe.

Mark Dickinson wrote:

it doesn't touch errno, but lets the platform decide how to handle errors (i.e. produce a
special value/set errno/signal a floating-point exception/some combination of these). This will
make the asinh, acosh, atanh functions behave in the same way that the regular libm functions
behave on any platform. So e.g. if a particular platform is used to setting errno for domain
errors like sqrt(-1), it'll do so for atanh/asinh/acosh. And another platform that signals a
floating-point exception for sqrt(-1) will do the same for atanh(3).

I tried your patch. It fixes the problem with atanh0022 and 0023 test
but it breaks in other places. math.acosh(0) returns NaN and does NOT
raise an exception. math.atanh(1) raises OverflowError instead of
ValueError.

The libm of uclibc *does* set errno for signaling NaN. The relevant code
is in w_acos.c and k_standard.c. I don't know how emit a signal for a
NaN but setting errno sounds reasonable for me and it gives the desired
output.

w_acosh.c:
double acosh(double x)          /* wrapper acosh */
{
#ifdef _IEEE_LIBM
        return __ieee754_acosh(x);
#else
        double z;
        z = __ieee754_acosh(x);
        if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
        if(x<1.0) {
                return __kernel_standard(x,x,29); /* acosh(x<1) */
        } else
            return z;
#endif
}

kernel_standard:
            case 129:
                /* acosh(x<1) */
                exc.type = DOMAIN;
                exc.name = type < 100 ? "acosh" : "acoshf";
                exc.retval = zero/zero;
                if (_LIB_VERSION == _POSIX_)
                  errno = EDOM;
                else if (!matherr(&exc)) {
                  if (_LIB_VERSION == _SVID_) {
                    (void) WRITE2("acosh: DOMAIN error\n", 20);
                  }
                  errno = EDOM;
                }
                break;

Hmmm.

For atanh(1): raising OverflowError is actually consistent with what currently happens.
The only singularity that's already present in math is log(0), and it seems that that
raises OverflowError on OS X, Linux and Windows... I wonder whether this is what Tim
meant to say?

For acosh(0): you're right, and I'm wrong---this should definitely return a NaN and set
errno. I guess that dividing 0 by 0 doesn't set errno on Windows. Okay: let's set it
directly there.

I do still think that asinh(nan), atanh(nan) and acosh(nan) should return nan and not
raise any exceptions, just for the sake of consistency with Linux/OS X and with the other
libm functions.

I guess I don't really care about asinh(+/-inf), etc: an infinite return value will be
caught by the stuff in math_1 anyway.

Mark Dickinson wrote:

For atanh(1): raising OverflowError is actually consistent with what currently happens.
The only singularity that's already present in math is log(0), and it seems that that
raises OverflowError on OS X, Linux and Windows... I wonder whether this is what Tim
meant to say?

But on Linux atanh(1) sets errno = EDOM and raises a ValueError. Either
the outcome is not consistent or the libm developers have a different
philosophy than Time and you. *g*

For acosh(0): you're right, and I'm wrong---this should definitely return a NaN and set
errno. I guess that dividing 0 by 0 doesn't set errno on Windows. Okay: let's set it
directly there.

I'm not an expert on IEEE754 and NaN but I guess "Signaling NaN" also
means it sets errno. Personally I'm not interested in reimplementing or
fixing Windows' mathematical library (although it's tempting *g*) but to
get the results right for our purpose. You are right, we should try to
mimic the C99 standard whenever it is possible. But it's more important
to get the functions of our math module right.

I do still think that asinh(nan), atanh(nan) and acosh(nan) should return nan and not
raise any exceptions, just for the sake of consistency with Linux/OS X and with the other
libm functions.

Yes, it makes perfectly sense and glibc's libm already does the right
thing. However Windows screws up again so I've to fix the function call
manually. Damn you Windows!

I inserted the four lines in mathmodule.c:math_1() and something similar
in math_2() to fix Windows.
#ifndef __GNUC__ /* Windows et al */
	if (Py_IS_NAN(x))
		return PyFloat_FromDouble(x+x);
#endif

I guess I don't really care about asinh(+/-inf), etc: an infinite return value will be
caught by the stuff in math_1 anyway.

Are you fine with OverflowError?

I'm going to create a branch for your work. It's getting tedious sending
large patches back and forth.

Christian

The patch was part of the merge of Mark's and my trunk-math branch. It
was merged into the trunk and 3.0 a while ago.

Navigation Menu

Repository navigation

Enhancements for mathmodule #45981

30 remaining items

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Enhancements for mathmodule #45981

Description

Activity

tiran commented on Dec 17, 2007

mdickinson commented on Dec 17, 2007

Rhamphoryncus commented on Dec 17, 2007

tiran commented on Dec 17, 2007

tiran commented on Dec 17, 2007

tiran commented on Dec 18, 2007

mdickinson commented on Dec 18, 2007

gvanrossum commented on Dec 18, 2007

tiran commented on Dec 18, 2007

tiran commented on Dec 18, 2007

gvanrossum commented on Jan 3, 2008

tiran commented on Jan 3, 2008

gvanrossum commented on Jan 3, 2008

30 remaining items

mdickinson commented on Jan 21, 2008

tiran commented on Jan 21, 2008

mdickinson commented on Jan 21, 2008

tiran commented on Jan 21, 2008

tiran commented on Jan 21, 2008

mdickinson commented on Jan 21, 2008

mdickinson commented on Jan 21, 2008

mdickinson commented on Jan 21, 2008

mdickinson commented on Jan 21, 2008

tiran commented on Jan 22, 2008

mdickinson commented on Jan 22, 2008

tiran commented on Jan 22, 2008

tiran commented on May 1, 2008

Metadata

Metadata

Assignees

Labels

Projects

Milestone

Relationships

Development

Issue actions