Enhancements for mathmodule

[tiran]

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: ```python assignee = None closed_at = created_at = labels = ['extension-modules', 'type-feature'] title = 'Enhancements for mathmodule' updated_at = user = 'https://github.com/tiran' ``` bugs.python.org fields: ```python activity = actor = 'christian.heimes' assignee = 'none' closed = True closed_date = closer = 'christian.heimes' components = ['Extension Modules'] creation = 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'] ```

[tiran]

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.

[mdickinson]

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.

[8d8a8db7-faf5-4c09-a2a3-2697dbaf0735]

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

[tiran]

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

[tiran]

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

[tiran]

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.

[mdickinson]

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

[gvanrossum]

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.

[tiran]

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

[tiran]

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).

[gvanrossum]

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.

[tiran]

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

[gvanrossum]

> 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.

[tiran]

Guido van Rossum wrote: > 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'm fine with a isneg() function but I wouldn't "fix" it for NaN. It has probably some kind of obscure meaning. The best explanation I was able to find, is http://www.cisl.ucar.edu/docs/trap.error/errortypes.html

"Note: Since NaN is "not a number," there really isn't a "+" or "-" sign associated with it."

Christian

[mdickinson]

Why not implement copysign? It's a standard, familiar function with well- documented meaning all over the web, and it can easily be used to create whatever sign test is necessary, letting the user decide what results (s)he wants for +/-0, +/-nan, etc.

[gvanrossum]

Why not implement copysign? It's a standard, familiar function with well- documented meaning all over the web, and it can easily be used to create whatever sign test is necessary, letting the user decide what results (s)he wants for +/-0, +/-nan, etc.

Good idea. Since you seem to like providing patches, can you create one for math.copysign()?

[tiran]

Guido van Rossum wrote:

Good idea. Since you seem to like providing patches, can you create one for math.copysign()?

Don't forget it's copysign() on Unix but _copysign() on Windows.

#if defined(MS_WINDOWS) || defined(HAVE_COPYSIGN)
#ifdef MS_WINDOWS
FUNC2(copysign, _copysign,
#else
FUNC2(copysign, copysign,
#endif
    "doc");
#endif

should work on all systems.

Christian

[gvanrossum]

Guido van Rossum wrote: > Good idea. Since you seem to like providing patches, can you create > one for math.copysign()?

Don't forget it's copysign() on Unix but _copysign() on Windows.

if defined(MS_WINDOWS) || defined(HAVE_COPYSIGN)

ifdef MS_WINDOWS

FUNC2(copysign, _copysign,

else

FUNC2(copysign, copysign,

endif

"doc");

endif

should work on all systems.

Well, the Python API in the math module should always be called copysign(). :-)

And what to do if neither is present? Are there any systems without it?

[tiran]

Guido van Rossum wrote:

Well, the Python API in the math module should always be called copysign(). :-)

And what to do if neither is present? Are there any systems without it?

takes care of it. It's a macro to define a function which accepts two floats and returns a float: FUNC2(funcname, func, docstring).

On Windows _copysign is always defined and on other systems we can use HAVE_COPYSIGN. I added it to configure.in a while ago.

Christian

[gvanrossum]

OK, just check it in then, but do add docs!

On Jan 3, 2008 2:03 PM, Christian Heimes \report@bugs.python.org\ wrote:

Christian Heimes added the comment:

Guido van Rossum wrote: > Well, the Python API in the math module should always be called copysign(). :-) > > And what to do if neither is present? Are there any systems without it?

takes care of it. It's a macro to define a function which accepts two floats and returns a float: FUNC2(funcname, func, docstring).

On Windows _copysign is always defined and on other systems we can use HAVE_COPYSIGN. I added it to configure.in a while ago.

Christian

---

Tracker \report@bugs.python.org\ \http://bugs.python.org/issue1640\

---

[tim-one]

The functionality of what was called (and I agree confusingly so) "sign()" here is supplied by "signbit()" in C99. That standard isn't free, but the relevant part is incorporated in the free Open Group standards:

http://www.opengroup.org/onlinepubs/000095399/functions/signbit.html

The 754 standard says NaNs have sign bits, but assigns no meaning to them. In particular, the value of the sign bit of a NaN resulting from a string->double routine applied to the string "nan" isn't defined by 754, or (AFAICT) by C99 either.

[3af7f971-76e6-491d-9967-91f2c806ac0d]

Is there still interest in implementing the inverse hyperbolic trig functions for real numbers? I would be willing to explore this if there is.

[mdickinson]

George: I'm certainly still interested in having asinh, acosh and atanh in math---I'm not sure about anyone else, but I consider these three functions to be basic ingredients in any math library. They even appear in most calculus texts. And the complex versions are already in cmath.

The right thing to do would be to wrap the libm functions when they exist, but provide fallbacks for the platforms where they don't. I think these functions are present on OS X and in the GNU libm, but they might be missing on Windows.

There's a fallback version of asinh already provided in the patch to issue bpo-1381, which I believe to be reasonably sound (but since I wrote it, a second opinion would be good).

For acosh and atanh, you'd have to find the right way to deal with domain errors (e.g. acosh(0.5), atanh(2.0)) and singularities (atanh(1.0), atanh(-1.0)). I'd suggest trying to follow the details in Appendix F (section F.9.2) of the C99 standard where possible.

[tiran]

I'm +1 in adding fallbacks for important functions like copysign, asinh, acosh and atanh. expm1 and log1p may be worth adding, too. Windows doesn't have any of the functions except of _copysign().

But why write our own version if we can reuse existing implementations? I found a set of very well written and documented implementations in the uclibc sources of libm. The sources are under a BSDish license:

• \====================================================
• Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. *
• Developed at SunPro, a Sun Microsystems, Inc. business.
• Permission to use, copy, modify, and distribute this
• software is freely granted, provided that this notice
• is preserved.
• \====================================================

In bpo-1381 I suggested two new files Python/pymath.c and Include/pymath.h. We could stick all the replacement implementations in the files and maybe move some of the code from Include/pyport.h and Python/hypot.c to the new files.

[tiran]

The patch implements replacements for copysign, log1p (taken from Mark's patch bpo-1381) asinh, acosh, atanh (taken from uclibc's libm). I modified the a*h functions for our needs. They set an errno and use our macros. The patch also adds the three hyberbolic arc functions to the math module. Math related code is now in Include/pymath.h and Python/pymath.c.

copysign is now defined on Windows, too.

I also found a bug in configure.in (my fault). The C99 function is called finite(), not isfinite().

Christian

[mdickinson]

Excellent! I'll take a look.

[mdickinson]

One question: is there a policy on what should happen for singularities and domain errors? If not, I think it would be useful to have one. Following the policy may not be achievable on all platforms, but it should be easy to do on major platforms, and at least we'll know what we're aiming for in general. Maybe it already exists, and I missed it :)

For domain errors (e.g. sqrt(-1), log(-1), acosh(0)), the obvious two options are:

• raise an exception (e.g. OverflowError), or
• return a NaN

For singularities (e.g. log(0), atanh(1)), the options are basically the same:

• raise an exception (preferably something different from OverflowError), or
• return the IEEE-754 recommended value (usually +/-Inf)

I suspect there are use-cases for both types of behaviour here.

Of course, the *right* thing to do, in some sense, would be to have a thread-local floating-point environment that makes it possible for the user to choose whether he/she wants an exception or a special value, much like the way Decimal behaves at the moment. But that would be a big change, almost certainly requiring a PEP and a lot of work.

A few months ago I definitely would have said that an exception should be raised in both cases, as already happens (mostly); but since then NaNs and Infinities have acquired greater legitimacy within Python.

Tim, if you're listening: any words of wisdom?

Should I ask this on python-dev?

[mdickinson]

whoops. OverflowError should be something else in the previous post; of course, OverflowError is inappropriate for domain errors (but entirely appropriate for something like exp(1000)).

Currently, on Linux I get:

• overflow (exp(1000)) -> OverflowError
• domain error (sqrt(-1)) -> ValueError
• singularity (log(0)) -> OverflowError

On OS X I get:

• overflow -> OverflowError
• domain error -> NaN
• singularity -> OverflowError

[tiran]

Mark Dickinson wrote:

Currently, on Linux I get:

• overflow (exp(1000)) -> OverflowError
• domain error (sqrt(-1)) -> ValueError
• singularity (log(0)) -> OverflowError

Windows raises the same exceptions as Linux.

[mdickinson]

Okay: for now, I guess we just follow the pattern that already exists on Linux and Windows.

I think the OS X sqrt(-1) behaviour is a bug.

[3af7f971-76e6-491d-9967-91f2c806ac0d]

Just a quick addition here regarding the singularities to these functions. The atanh(x) is only defined for |x| \< 1, so atanh(1) or atanh(-1) isn't singular there so much as simply isn't defined. So, even though the function approaches infinite as x -> 1, it wouldn't really be correct to return a value at |x| = 1. I think raising an exception at those points would be more correct.

[mdickinson]

No: IEEE-754r and the C99 standard both say clearly that atanh(1) should be infinity and atanh(-1) should be -infinity, and furthermore that the 'divide-by-zero' exception should be raised rather than the 'invalid' exception. It's a singularity, just like log(0). (This makes even more sense viewed from the perspective of complex arithmetic, where atanh is defined at all points in the complex plane except -1 and 1, where it has log-type singularities.)

The general idea is that it's meaningful to set atanh(1) = infinity because that's what the limit of atanh(x) is as x approaches 1 from below; similarly for atanh(-1) and log(0).

[3af7f971-76e6-491d-9967-91f2c806ac0d]

I misunderstood the rationale for the function returning infinite at those points - I didn't realize that C99 was the governing force behind the implementation of these functions, rather than mathematical rigor. Thanks for pointing it out. In that case, I agree with you that, in order to conform, these functions would need to return the values required by those documents.

[mdickinson]

George: I think my last post was a bit rude. I apologize if it came across that way.

Mathematical rigor and IEEE-754 recommendations aren't necessarily in conflict here, though. For example, the natural log function from (0, infinity) to (-infinity, infinity) extends naturally and uniquely to a continuous function on the closed subset [0, infinity] of the extended real line---i.e., the real line together with the two extra points -infinity and infinity. With the appropriate topology, the extended real line is a perfectly well-defined and well-behaved mathematical object, though of course it's no longer a field. Since IEEE-754 floats include infinities, it's reasonable, and sometimes useful, to regard the set of IEEE-floats as a computational model of the extended real line instead of the reals.

At any rate, I agree with you that log(0) and atanh(1) should raise Python exceptions, at least for now. But these calculations are qualitatively different from log(-1) and atanh(2), and it wouldn't be at all unreasonable if they raised a different exception--- e.g. ZeroDivisionError instead of ValueError.

[tim-one]

Mark, these are the "spirit of 754" rules:

1. If the mathematical result is a real number, but of magnitude too large to approximate by a machine float, overflow is signaled and the result is an infinity (with the appropriate sign).

2. If the mathematical result is a real number, but of magnitude too small to approximate by a machine float, underflow is signaled and the result is a zero (with the appropriate sign).

3. At a singularity (a value x such that the limit of f(y) as y approaches x exists and is an infinity), "divide by zero" is signaled and the result is an infinity (with the appropriate sign). This is complicated a little by that the left-side and right-side limits may not be the same; e.g., 1/x approaches +inf or -inf as x approaches 0 from the positive or negative directions. In that specific case, the sign of the zero determines the result of 1/0.

4. At a point where a function has no defined result in the extended reals (i.e., the reals plus an infinity or two), invalid operation is signaled and a NaN is returned.

And these are what Python has historically /tried/ to do (but not always successfully, as platform libm behavior varies a lot):

For #1, raise OverflowError.

For #2, return a zero (with the appropriate sign if that happens by accident ;-)).

For #3 and #4, raise ValueError. It may have made sense to raise Python's ZeroDivisionError in #3, but historically that's only been raised for division by zero and mod by zero.

You're right that you can't make everyone happy, so may as well stick with historical consistency, and wait for a PEP to introduce a more flexible mechanism. As W. Kahan wrote long ago, the reason they're called "exceptions" is that no matter which fixed policy you adopt, someone will take strident exception to it ;-)

[mdickinson]

Thanks, Tim!

Dare I suggest extending these rules to encompass things like sqrt(NaN), log(inf), etc., as follows:

• return a special value in Python where IEEE-754r/C99 specifies a special value, but doesn't raise any of the three divide-by-zero, invalid, or overflow exceptions, and
• raise OverflowError or ValueError as appropriate where IEEE-754r specifies raising one of these exceptions.

So e.g. cos(infinity) should give a ValueError, while log(infinity) and exp(infinity) should not raise any Python exception, but should return an infinity instead. And most single variable operations should return an input NaN unaltered, without raising an exception.

[tiran]

Mark Dickinson wrote:

So e.g. cos(infinity) should give a ValueError, while log(infinity) and exp(infinity) should not raise any Python exception, but should return an infinity instead. And most single variable operations should return an input NaN unaltered, without raising an exception.

The matter should be discussed in a proper PEP and targeted for Python 3.0. Python 3.0 is the right place for subtle changes which may break code. For Python 2.6 we must not change the exception or outcome of a function and new functions should be as consistent with existing ones as possible.

I still don't like the idea of math.atanh(1) == inf. Why? See for yourself:

18.714973875118524
>>> math.atanh(.99999999999999999)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: math domain error

(Linux)

Christian

[tiran]

The mail interface screwed up the message:

>>> math.atanh(.9999999999999999)
18.714973875118524
>>> math.atanh(.99999999999999999)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: math domain error

[tiran]

I've added your complex patch and its tests to my patch. The tests are showing some flaws in the atanh (or log1p) function on Windows:

AssertionError: atanh0022:atanh(-1.000000) returned -18.36840028483855, expected -18.714973875118524

On Linux the tests are passing just fine.

[tiran]

It's most probably the fault of log1p():

AssertionError: atanh0022:atanh(-0.99999999999999989) returned -18.36840028483855, expected -18.714973875118524

[mdickinson]

Christian: I'm definitely not proposing atanh(1) = inf: it should raise ValueError. I'm proposing that we follow Tim's rules for now; this means no change for finite inputs.

The new thing here is that since you've made inf and nan more accessible and consistent across platforms, I think we should make sure that the math functions do the right thing for an *input* of +/-inf or nan. I'm almost sure that the current behavior of e.g. exp(float("inf")) is more-or-less accidental rather than designed.

I think I'm missing the point of your math.atanh(.999...) example. .99999999999999999 *is* already exactly equal to 1.0, so you're just proving that math.atanh(1.0) currently gives a ValueError. (Which, again, I think is the right thing to do.)

>>> x = .99999999999999999
>>> x == 1.0
True

The atanh0022 result is definitely a bug: it looks like either asinh or log1p is buggy. I'll try to figure it out.

[mdickinson]

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).

[mdickinson]

Okay: here's an attempted guess at what's happening on Windows:

Near the end of c_atanh, there's a line that reads:

r.real = log1p(4.*z.real/((1-z.real)*(1-z.real) + ay*ay))/4.;

Somehow, when z.real is 0.99999999999999989 and ay is 0, the argument to log1p is ending up one-quarter of the size that it should be. I suspect that the 1-z.real calculation is producing, for reasons that are beyond me, the float 2-52 instead of the correct value of 2- 53.

Christian: if you have any time to play with this, could you try replacing this line with something like:

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

and see if either the problem goes away or if you can confirm that temp is coming out to be 2.2204460492503131e-16 rather than 1.1102230246251565e-16.

[mdickinson]

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.;

[tiran]

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

[mdickinson]

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.

[tiran]

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

[tiran]

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

[mdickinson]

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.

[mdickinson]

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.