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 #1635 and #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 #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 #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.

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

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.

> 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()?

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

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

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

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

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.

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

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

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

The patch implements replacements for copysign, log1p (taken from Mark's
patch #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

Excellent!  I'll take a look.

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?

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

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.

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.

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.

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

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.

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.

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

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.

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

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

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.

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

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

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.

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

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.

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.