Mirror
/
cpython
mirror of https://github.com/python/cpython
4
1
Fork 0
Code Issues Releases Wiki Activity

I finally got the time to update and merge Mark's and my trunk-math branch. The patch is collaborated work of Mark Dickinson and me. It was mostly done a few months ago. The patch fixes a lot of loose ends and edge cases related to operations with NaN, INF, very small values and complex math. 

The patch also adds acosh, asinh, atanh, log1p and copysign to all platforms. Finally it fixes differences between platforms like different results or exceptions for edge cases. Have fun :)
This commit is contained in:
Christian Heimes 2008-04-18 23:13:07 +00:00
parent 858a77099e
commit 6f34109384
25 changed files with 5101 additions and 1158 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

120
Doc/library/cmath.rst
View File

@ -14,8 +14,81 @@ method: these methods are used to convert the object to a complex or
floating-point number, respectively, and the function is then applied to the
result of the conversion.
The functions are:
.. note::
On platforms with hardware and system-level support for signed
zeros, functions involving branch cuts are continuous on *both*
sides of the branch cut: the sign of the zero distinguishes one
side of the branch cut from the other. On platforms that do not
support signed zeros the continuity is as specified below.
Complex coordinates
-------------------
Complex numbers can be expressed by two important coordinate systems.
Python's :class:`complex` type uses rectangular coordinates where a number
on the complex plain is defined by two floats, the real part and the imaginary
part.
Definition::
z = x + 1j * y
x := real(z)
y := imag(z)
In engineering the polar coordinate system is popular for complex numbers. In
polar coordinates a complex number is defined by the radius *r* and the phase
angle *φ*. The radius *r* is the absolute value of the complex, which can be
viewed as distance from (0, 0). The radius *r* is always 0 or a positive float.
The phase angle *φ* is the counter clockwise angle from the positive x axis,
e.g. *1* has the angle *0*, *1j* has the angle *π/2* and *-1* the angle *-π*.
.. note::
While :func:`phase` and func:`polar` return *+π* for a negative real they
may return *-π* for a complex with a very small negative imaginary
part, e.g. *-1-1E-300j*.
Definition::
z = r * exp(1j * φ)
z = r * cis(φ)
r := abs(z) := sqrt(real(z)**2 + imag(z)**2)
phi := phase(z) := atan2(imag(z), real(z))
cis(φ) := cos(φ) + 1j * sin(φ)
.. function:: phase(x)
Return phase, also known as the argument, of a complex.
.. versionadded:: 2.6
.. function:: polar(x)
Convert a :class:`complex` from rectangular coordinates to polar
coordinates. The function returns a tuple with the two elements
*r* and *phi*. *r* is the distance from 0 and *phi* the phase
angle.
.. versionadded:: 2.6
.. function:: rect(r, phi)
Convert from polar coordinates to rectangular coordinates and return
a :class:`complex`.
.. versionadded:: 2.6
cmath functions
---------------
.. function:: acos(x)
@ -37,30 +110,35 @@ The functions are:
.. function:: asinh(x)
Return the hyperbolic arc sine of *x*. There are two branch cuts, extending
left from ``±1j`` to ``±∞j``, both continuous from above. These branch cuts
should be considered a bug to be corrected in a future release. The correct
branch cuts should extend along the imaginary axis, one from ``1j`` up to
``∞j`` and continuous from the right, and one from ``-1j`` down to ``-∞j``
and continuous from the left.
Return the hyperbolic arc sine of *x*. There are two branch cuts:
One extends from ``1j`` along the imaginary axis to ``∞j``,
continuous from the right. The other extends from ``-1j`` along
the imaginary axis to ``-∞j``, continuous from the left.
.. versionchanged:: 2.6
branch cuts moved to match those recommended by the C99 standard
.. function:: atan(x)
Return the arc tangent of *x*. There are two branch cuts: One extends from
``1j`` along the imaginary axis to ``∞j``, continuous from the left. The
``1j`` along the imaginary axis to ``∞j``, continuous from the right. The
other extends from ``-1j`` along the imaginary axis to ``-∞j``, continuous
from the left. (This should probably be changed so the upper cut becomes
continuous from the other side.)
from the left.
.. versionchanged:: 2.6
direction of continuity of upper cut reversed
.. function:: atanh(x)
Return the hyperbolic arc tangent of *x*. There are two branch cuts: One
extends from ``1`` along the real axis to ``∞``, continuous from above. The
extends from ``1`` along the real axis to ``∞``, continuous from below. The
other extends from ``-1`` along the real axis to ``-∞``, continuous from
above. (This should probably be changed so the right cut becomes continuous
from the other side.)
above.
.. versionchanged:: 2.6
direction of continuity of right cut reversed
.. function:: cos(x)
@ -78,6 +156,21 @@ The functions are:
Return the exponential value ``e**x``.
.. function:: isinf(x)
Return *True* if the real or the imaginary part of x is positive
or negative infinity.
.. versionadded:: 2.6
.. function:: isnan(x)
Return *True* if the real or imaginary part of x is not a number (NaN).
.. versionadded:: 2.6
.. function:: log(x[, base])
Returns the logarithm of *x* to the given *base*. If the *base* is not
@ -154,3 +247,4 @@ cuts for numerical purposes, a good reference should be the following:
nothing's sign bit. In Iserles, A., and Powell, M. (eds.), The state of the art
in numerical analysis. Clarendon Press (1987) pp165-211.

47
Doc/library/math.rst
View File

@ -139,6 +139,14 @@ Power and logarithmic functions:
*base* argument added.
.. function:: log1p(x)
Return the natural logarithm of *1+x* (base *e*). The
result is calculated in a way which is accurate for *x* near zero.
.. versionadded:: 2.6
.. function:: log10(x)
Return the base-10 logarithm of *x*.
@ -146,7 +154,11 @@ Power and logarithmic functions:
.. function:: pow(x, y)
Return ``x**y``.
Return ``x**y``. ``1.0**y`` returns *1.0*, even for ``1.0**nan``. ``0**y``
returns *0.* for all positive *y*, *0* and *NAN*.
.. versionchanged:: 2.6
The outcome of ``1**nan`` and ``0**nan`` was undefined.
.. function:: sqrt(x)
@ -197,6 +209,13 @@ Trigonometric functions:
Return the sine of *x* radians.
.. function:: asinh(x)
Return the inverse hyperbolic sine of *x*, in radians.
.. versionadded:: 2.6
.. function:: tan(x)
Return the tangent of *x* radians.
@ -221,6 +240,13 @@ Hyperbolic functions:
Return the hyperbolic cosine of *x*.
.. function:: acosh(x)
Return the inverse hyperbolic cosine of *x*, in radians.
.. versionadded:: 2.6
.. function:: sinh(x)
Return the hyperbolic sine of *x*.
@ -230,6 +256,14 @@ Hyperbolic functions:
Return the hyperbolic tangent of *x*.
.. function:: atanh(x)
Return the inverse hyperbolic tangent of *x*, in radians.
.. versionadded:: 2.6
The module also defines two mathematical constants:
@ -242,6 +276,7 @@ The module also defines two mathematical constants:
The mathematical constant *e*.
.. note::
The :mod:`math` module consists mostly of thin wrappers around the platform C
@ -255,9 +290,17 @@ The module also defines two mathematical constants:
:exc:`OverflowError` isn't defined, and in cases where ``math.log(0)`` raises
:exc:`OverflowError`, ``math.log(0L)`` may raise :exc:`ValueError` instead.
All functions return a quite *NaN* if at least one of the args is *NaN*.
Signaling *NaN*s raise an exception. The exception type still depends on the
platform and libm implementation. It's usually :exc:`ValueError` for *EDOM*
and :exc:`OverflowError` for errno *ERANGE*.
..versionchanged:: 2.6
In earlier versions of Python the outcome of an operation with NaN as
input depended on platform and libm implementation.
.. seealso::
Module :mod:`cmath`
Complex number versions of many of these functions.

1
Include/Python.h
View File

@ -73,6 +73,7 @@
#if defined(PYMALLOC_DEBUG) && !defined(WITH_PYMALLOC)
#error "PYMALLOC_DEBUG requires WITH_PYMALLOC"
#endif
#include "pymath.h"
#include "pymem.h"
#include "object.h"

2
Include/complexobject.h
View File

@ -19,6 +19,7 @@ typedef struct {
#define c_prod _Py_c_prod
#define c_quot _Py_c_quot
#define c_pow _Py_c_pow
#define c_abs _Py_c_abs
PyAPI_FUNC(Py_complex) c_sum(Py_complex, Py_complex);
PyAPI_FUNC(Py_complex) c_diff(Py_complex, Py_complex);
@ -26,6 +27,7 @@ PyAPI_FUNC(Py_complex) c_neg(Py_complex);
PyAPI_FUNC(Py_complex) c_prod(Py_complex, Py_complex);
PyAPI_FUNC(Py_complex) c_quot(Py_complex, Py_complex);
PyAPI_FUNC(Py_complex) c_pow(Py_complex, Py_complex);
PyAPI_FUNC(double) c_abs(Py_complex);
/* Complex object interface */

11
Include/floatobject.h
View File

@ -21,6 +21,17 @@ PyAPI_DATA(PyTypeObject) PyFloat_Type;
#define PyFloat_Check(op) PyObject_TypeCheck(op, &PyFloat_Type)
#define PyFloat_CheckExact(op) (Py_TYPE(op) == &PyFloat_Type)
#ifdef Py_NAN
#define Py_RETURN_NAN PyFloat_FromDouble(Py_NAN)
#endif
#define Py_RETURN_INF(sign) do \
if (copysign(1., sign) == 1.) { \
return PyFloat_FromDouble(Py_HUGE_VAL); \
} else { \
return PyFloat_FromDouble(-Py_HUGE_VAL); \
} while(0)
PyAPI_FUNC(double) PyFloat_GetMax(void);
PyAPI_FUNC(double) PyFloat_GetMin(void);
PyAPI_FUNC(PyObject *) PyFloat_GetInfo(void);

182
Include/pymath.h Normal file
View File

@ -0,0 +1,182 @@
#ifndef Py_PYMATH_H
#define Py_PYMATH_H
#include "pyconfig.h" /* include for defines */
#ifdef HAVE_STDINT_H
#include <stdint.h>
#endif
/**************************************************************************
Symbols and macros to supply platform-independent interfaces to mathematical
functions and constants
**************************************************************************/
/* Python provides implementations for copysign, acosh, asinh, atanh,
* log1p and hypot in Python/pymath.c just in case your math library doesn't
* provide the functions.
*
*Note: PC/pyconfig.h defines copysign as _copysign
*/
#ifndef HAVE_COPYSIGN
extern double copysign(doube, double);
#endif
#ifndef HAVE_ACOSH
extern double acosh(double);
#endif
#ifndef HAVE_ASINH
extern double asinh(double);
#endif
#ifndef HAVE_ATANH
extern double atanh(double);
#endif
#ifndef HAVE_LOG1P
extern double log1p(double);
#endif
#ifndef HAVE_HYPOT
extern double hypot(double, double);
#endif
/* extra declarations */
#ifndef _MSC_VER
#ifndef __STDC__
extern double fmod (double, double);
extern double frexp (double, int *);
extern double ldexp (double, int);
extern double modf (double, double *);
extern double pow(double, double);
#endif /* __STDC__ */
#endif /* _MSC_VER */
#ifdef _OSF_SOURCE
/* OSF1 5.1 doesn't make these available with XOPEN_SOURCE_EXTENDED defined */
extern int finite(double);
extern double copysign(double, double);
#endif
/* High precision defintion of pi and e (Euler)
* The values are taken from libc6's math.h.
*/
#ifndef Py_MATH_PIl
#define Py_MATH_PIl 3.1415926535897932384626433832795029L
#endif
#ifndef Py_MATH_PI
#define Py_MATH_PI 3.14159265358979323846
#endif
#ifndef Py_MATH_El
#define Py_MATH_El 2.7182818284590452353602874713526625L
#endif
#ifndef Py_MATH_E
#define Py_MATH_E 2.7182818284590452354
#endif
/* Py_IS_NAN(X)
* Return 1 if float or double arg is a NaN, else 0.
* Caution:
* X is evaluated more than once.
* This may not work on all platforms. Each platform has *some*
* way to spell this, though -- override in pyconfig.h if you have
* a platform where it doesn't work.
* Note: PC/pyconfig.h defines Py_IS_NAN as _isnan
*/
#ifndef Py_IS_NAN
#ifdef HAVE_ISNAN
#define Py_IS_NAN(X) isnan(X)
#else
#define Py_IS_NAN(X) ((X) != (X))
#endif
#endif
/* Py_IS_INFINITY(X)
* Return 1 if float or double arg is an infinity, else 0.
* Caution:
* X is evaluated more than once.
* This implementation may set the underflow flag if |X| is very small;
* it really can't be implemented correctly (& easily) before C99.
* Override in pyconfig.h if you have a better spelling on your platform.
* Note: PC/pyconfig.h defines Py_IS_INFINITY as _isinf
*/
#ifndef Py_IS_INFINITY
#ifdef HAVE_ISINF
#define Py_IS_INFINITY(X) isinf(X)
#else
#define Py_IS_INFINITY(X) ((X) && (X)*0.5 == (X))
#endif
#endif
/* Py_IS_FINITE(X)
* Return 1 if float or double arg is neither infinite nor NAN, else 0.
* Some compilers (e.g. VisualStudio) have intrisics for this, so a special
* macro for this particular test is useful
* Note: PC/pyconfig.h defines Py_IS_FINITE as _finite
*/
#ifndef Py_IS_FINITE
#ifdef HAVE_FINITE
#define Py_IS_FINITE(X) finite(X)
#else
#define Py_IS_FINITE(X) (!Py_IS_INFINITY(X) && !Py_IS_NAN(X))
#endif
#endif
/* HUGE_VAL is supposed to expand to a positive double infinity. Python
* uses Py_HUGE_VAL instead because some platforms are broken in this
* respect. We used to embed code in pyport.h to try to worm around that,
* but different platforms are broken in conflicting ways. If you're on
* a platform where HUGE_VAL is defined incorrectly, fiddle your Python
* config to #define Py_HUGE_VAL to something that works on your platform.
*/
#ifndef Py_HUGE_VAL
#define Py_HUGE_VAL HUGE_VAL
#endif
/* Py_NAN
* A value that evaluates to a NaN. On IEEE 754 platforms INF*0 or
* INF/INF works. Define Py_NO_NAN in pyconfig.h if your platform
* doesn't support NaNs.
*/
#if !defined(Py_NAN) && !defined(Py_NO_NAN)
#define Py_NAN (Py_HUGE_VAL * 0.)
#endif
/* Py_OVERFLOWED(X)
* Return 1 iff a libm function overflowed. Set errno to 0 before calling
* a libm function, and invoke this macro after, passing the function
* result.
* Caution:
* This isn't reliable. C99 no longer requires libm to set errno under
* any exceptional condition, but does require +- HUGE_VAL return
* values on overflow. A 754 box *probably* maps HUGE_VAL to a
* double infinity, and we're cool if that's so, unless the input
* was an infinity and an infinity is the expected result. A C89
* system sets errno to ERANGE, so we check for that too. We're
* out of luck if a C99 754 box doesn't map HUGE_VAL to +Inf, or
* if the returned result is a NaN, or if a C89 box returns HUGE_VAL
* in non-overflow cases.
* X is evaluated more than once.
* Some platforms have better way to spell this, so expect some #ifdef'ery.
*
* OpenBSD uses 'isinf()' because a compiler bug on that platform causes
* the longer macro version to be mis-compiled. This isn't optimal, and
* should be removed once a newer compiler is available on that platform.
* The system that had the failure was running OpenBSD 3.2 on Intel, with
* gcc 2.95.3.
*
* According to Tim's checkin, the FreeBSD systems use isinf() to work
* around a FPE bug on that platform.
*/
#if defined(__FreeBSD__) || defined(__OpenBSD__)
#define Py_OVERFLOWED(X) isinf(X)
#else
#define Py_OVERFLOWED(X) ((X) != 0.0 && (errno == ERANGE || \
(X) == Py_HUGE_VAL || \
(X) == -Py_HUGE_VAL))
#endif
#endif /* Py_PYMATH_H */

126
Include/pyport.h
View File

@ -353,123 +353,6 @@ extern "C" {
#define Py_SAFE_DOWNCAST(VALUE, WIDE, NARROW) (NARROW)(VALUE)
#endif
/* High precision defintion of pi and e (Euler)
* The values are taken from libc6's math.h.
*/
#ifndef Py_MATH_PIl
#define Py_MATH_PIl 3.1415926535897932384626433832795029L
#endif
#ifndef Py_MATH_PI
#define Py_MATH_PI 3.14159265358979323846
#endif
#ifndef Py_MATH_El
#define Py_MATH_El 2.7182818284590452353602874713526625L
#endif
#ifndef Py_MATH_E
#define Py_MATH_E 2.7182818284590452354
#endif
/* Py_IS_NAN(X)
* Return 1 if float or double arg is a NaN, else 0.
* Caution:
* X is evaluated more than once.
* This may not work on all platforms. Each platform has *some*
* way to spell this, though -- override in pyconfig.h if you have
* a platform where it doesn't work.
*/
#ifndef Py_IS_NAN
#ifdef HAVE_ISNAN
#define Py_IS_NAN(X) isnan(X)
#else
#define Py_IS_NAN(X) ((X) != (X))
#endif
#endif
/* Py_IS_INFINITY(X)
* Return 1 if float or double arg is an infinity, else 0.
* Caution:
* X is evaluated more than once.
* This implementation may set the underflow flag if |X| is very small;
* it really can't be implemented correctly (& easily) before C99.
* Override in pyconfig.h if you have a better spelling on your platform.
*/
#ifndef Py_IS_INFINITY
#ifdef HAVE_ISINF
#define Py_IS_INFINITY(X) isinf(X)
#else
#define Py_IS_INFINITY(X) ((X) && (X)*0.5 == (X))
#endif
#endif
/* Py_IS_FINITE(X)
* Return 1 if float or double arg is neither infinite nor NAN, else 0.
* Some compilers (e.g. VisualStudio) have intrisics for this, so a special
* macro for this particular test is useful
*/
#ifndef Py_IS_FINITE
#ifdef HAVE_FINITE
#define Py_IS_FINITE(X) finite(X)
#else
#define Py_IS_FINITE(X) (!Py_IS_INFINITY(X) && !Py_IS_NAN(X))
#endif
#endif
/* HUGE_VAL is supposed to expand to a positive double infinity. Python
* uses Py_HUGE_VAL instead because some platforms are broken in this
* respect. We used to embed code in pyport.h to try to worm around that,
* but different platforms are broken in conflicting ways. If you're on
* a platform where HUGE_VAL is defined incorrectly, fiddle your Python
* config to #define Py_HUGE_VAL to something that works on your platform.
*/
#ifndef Py_HUGE_VAL
#define Py_HUGE_VAL HUGE_VAL
#endif
/* Py_NAN
* A value that evaluates to a NaN. On IEEE 754 platforms INF*0 or
* INF/INF works. Define Py_NO_NAN in pyconfig.h if your platform
* doesn't support NaNs.
*/
#if !defined(Py_NAN) && !defined(Py_NO_NAN)
#define Py_NAN (Py_HUGE_VAL * 0.)
#endif
/* Py_OVERFLOWED(X)
* Return 1 iff a libm function overflowed. Set errno to 0 before calling
* a libm function, and invoke this macro after, passing the function
* result.
* Caution:
* This isn't reliable. C99 no longer requires libm to set errno under
* any exceptional condition, but does require +- HUGE_VAL return
* values on overflow. A 754 box *probably* maps HUGE_VAL to a
* double infinity, and we're cool if that's so, unless the input
* was an infinity and an infinity is the expected result. A C89
* system sets errno to ERANGE, so we check for that too. We're
* out of luck if a C99 754 box doesn't map HUGE_VAL to +Inf, or
* if the returned result is a NaN, or if a C89 box returns HUGE_VAL
* in non-overflow cases.
* X is evaluated more than once.
* Some platforms have better way to spell this, so expect some #ifdef'ery.
*
* OpenBSD uses 'isinf()' because a compiler bug on that platform causes
* the longer macro version to be mis-compiled. This isn't optimal, and
* should be removed once a newer compiler is available on that platform.
* The system that had the failure was running OpenBSD 3.2 on Intel, with
* gcc 2.95.3.
*
* According to Tim's checkin, the FreeBSD systems use isinf() to work
* around a FPE bug on that platform.
*/
#if defined(__FreeBSD__) || defined(__OpenBSD__)
#define Py_OVERFLOWED(X) isinf(X)
#else
#define Py_OVERFLOWED(X) ((X) != 0.0 && (errno == ERANGE || \
(X) == Py_HUGE_VAL || \
(X) == -Py_HUGE_VAL))
#endif
/* Py_SET_ERRNO_ON_MATH_ERROR(x)
* If a libm function did not set errno, but it looks like the result
* overflowed or not-a-number, set errno to ERANGE or EDOM. Set errno
@ -601,15 +484,6 @@ extern int fdatasync(int);
#endif /* 0 */
/************************
* WRAPPER FOR <math.h> *
************************/
#ifndef HAVE_HYPOT
extern double hypot(double, double);
#endif
/* On 4.4BSD-descendants, ctype functions serves the whole range of
* wchar_t character set rather than single byte code points only.
* This characteristic can break some operations of string object

2355
Lib/test/cmath_testcases.txt Normal file
View File

File diff suppressed because it is too large Load Diff

183
Lib/test/ieee754.txt Normal file
View File

@ -0,0 +1,183 @@
======================================
Python IEEE 754 floating point support
======================================
>>> from sys import float_info as FI
>>> from math import *
>>> PI = pi
>>> E = e
You must never compare two floats with == because you are not going to get
what you expect. We treat two floats as equal if the difference between them
is small than epsilon.
>>> EPS = 1E-15
>>> def equal(x, y):
... """Almost equal helper for floats"""
... return abs(x - y) < EPS
NaNs and INFs
=============
In Python 2.6 and newer NaNs (not a number) and infinity can be constructed
from the strings 'inf' and 'nan'.
>>> INF = float('inf')
>>> NINF = float('-inf')
>>> NAN = float('nan')
>>> INF
inf
>>> NINF
-inf
>>> NAN
nan
The math module's ``isnan`` and ``isinf`` functions can be used to detect INF
and NAN:
>>> isinf(INF), isinf(NINF), isnan(NAN)
(True, True, True)
>>> INF == -NINF
True
Infinity
--------
Ambiguous operations like ``0 * inf`` or ``inf - inf`` result in NaN.
>>> INF * 0
nan
>>> INF - INF
nan
>>> INF / INF
nan
However unambigous operations with inf return inf:
>>> INF * INF
inf
>>> 1.5 * INF
inf
>>> 0.5 * INF
inf
>>> INF / 1000
inf
Not a Number
------------
NaNs are never equal to another number, even itself
>>> NAN == NAN
False
>>> NAN < 0
False
>>> NAN >= 0
False
All operations involving a NaN return a NaN except for the power of *0* and *1*.
>>> 1 + NAN
nan
>>> 1 * NAN
nan
>>> 0 * NAN
nan
>>> 1 ** NAN
1.0
>>> 0 ** NAN
0.0
>>> (1.0 + FI.epsilon) * NAN
nan
Misc Functions
==============
The power of 1 raised to x is always 1.0, even for special values like 0,
infinity and NaN.
>>> pow(1, 0)
1.0
>>> pow(1, INF)
1.0
>>> pow(1, -INF)
1.0
>>> pow(1, NAN)
1.0
The power of 0 raised to x is defined as 0, if x is positive. Negative
values are a domain error or zero division error and NaN result in a
silent NaN.
>>> pow(0, 0)
1.0
>>> pow(0, INF)
0.0
>>> pow(0, -INF)
Traceback (most recent call last):
...
ValueError: math domain error
>>> 0 ** -1
Traceback (most recent call last):
...
ZeroDivisionError: 0.0 cannot be raised to a negative power
>>> pow(0, NAN)
nan
Trigonometric Functions
=======================
>>> sin(INF)
Traceback (most recent call last):
...
ValueError: math domain error
>>> sin(NINF)
Traceback (most recent call last):
...
ValueError: math domain error
>>> sin(NAN)
nan
>>> cos(INF)
Traceback (most recent call last):
...
ValueError: math domain error
>>> cos(NINF)
Traceback (most recent call last):
...
ValueError: math domain error
>>> cos(NAN)
nan
>>> tan(INF)
Traceback (most recent call last):
...
ValueError: math domain error
>>> tan(NINF)
Traceback (most recent call last):
...
ValueError: math domain error
>>> tan(NAN)
nan
Neither pi nor tan are exact, but you can assume that tan(pi/2) is a large value
and tan(pi) is a very small value:
>>> tan(PI/2) > 1E10
True
>>> -tan(-PI/2) > 1E10
True
>>> tan(PI) < 1E-15
True
>>> asin(NAN), acos(NAN), atan(NAN)
(nan, nan, nan)
>>> asin(INF), asin(NINF)
Traceback (most recent call last):
...
ValueError: math domain error
>>> acos(INF), acos(NINF)
Traceback (most recent call last):
...
ValueError: math domain error
>>> equal(atan(INF), PI/2), equal(atan(NINF), -PI/2)
(True, True)
Hyberbolic Functions
====================

336
Lib/test/test_cmath.py
View File

@ -1,6 +1,81 @@
from test.test_support import run_unittest
from test.test_math import parse_testfile, test_file
import unittest
import os, sys
import cmath, math
from cmath import phase, polar, rect, pi
INF = float('inf')
NAN = float('nan')
complex_zeros = [complex(x, y) for x in [0.0, -0.0] for y in [0.0, -0.0]]
complex_infinities = [complex(x, y) for x, y in [
(INF, 0.0), # 1st quadrant
(INF, 2.3),
(INF, INF),
(2.3, INF),
(0.0, INF),
(-0.0, INF), # 2nd quadrant
(-2.3, INF),
(-INF, INF),
(-INF, 2.3),
(-INF, 0.0),
(-INF, -0.0), # 3rd quadrant
(-INF, -2.3),
(-INF, -INF),
(-2.3, -INF),
(-0.0, -INF),
(0.0, -INF), # 4th quadrant
(2.3, -INF),
(INF, -INF),
(INF, -2.3),
(INF, -0.0)
]]
complex_nans = [complex(x, y) for x, y in [
(NAN, -INF),
(NAN, -2.3),
(NAN, -0.0),
(NAN, 0.0),
(NAN, 2.3),
(NAN, INF),
(-INF, NAN),
(-2.3, NAN),
(-0.0, NAN),
(0.0, NAN),
(2.3, NAN),
(INF, NAN)
]]
def almostEqualF(a, b, rel_err=2e-15, abs_err = 5e-323):
"""Determine whether floating-point values a and b are equal to within
a (small) rounding error. The default values for rel_err and
abs_err are chosen to be suitable for platforms where a float is
represented by an IEEE 754 double. They allow an error of between
9 and 19 ulps."""
# special values testing
if math.isnan(a):
return math.isnan(b)
if math.isinf(a):
return a == b
# if both a and b are zero, check whether they have the same sign
# (in theory there are examples where it would be legitimate for a
# and b to have opposite signs; in practice these hardly ever
# occur).
if not a and not b:
return math.copysign(1., a) == math.copysign(1., b)
# if a-b overflows, or b is infinite, return False. Again, in
# theory there are examples where a is within a few ulps of the
# max representable float, and then b could legitimately be
# infinite. In practice these examples are rare.
try:
absolute_error = abs(b-a)
except OverflowError:
return False
else:
return absolute_error <= max(abs_err, rel_err * abs(a))
class CMathTests(unittest.TestCase):
# list of all functions in cmath
@ -12,24 +87,52 @@ class CMathTests(unittest.TestCase):
test_functions.append(lambda x : cmath.log(x, 1729. + 0j))
test_functions.append(lambda x : cmath.log(14.-27j, x))
def cAssertAlmostEqual(self, a, b, rel_eps = 1e-10, abs_eps = 1e-100):
"""Check that two complex numbers are almost equal."""
# the two complex numbers are considered almost equal if
# either the relative error is <= rel_eps or the absolute error
# is tiny, <= abs_eps.
if a == b == 0:
return
absolute_error = abs(a-b)
relative_error = absolute_error/max(abs(a), abs(b))
if relative_error > rel_eps and absolute_error > abs_eps:
self.fail("%s and %s are not almost equal" % (a, b))
def setUp(self):
self.test_values = open(test_file)
def tearDown(self):
self.test_values.close()
def rAssertAlmostEqual(self, a, b, rel_err = 2e-15, abs_err = 5e-323):
"""Check that two floating-point numbers are almost equal."""
# special values testing
if math.isnan(a):
if math.isnan(b):
return
self.fail("%s should be nan" % repr(b))
if math.isinf(a):
if a == b:
return
self.fail("finite result where infinity excpected: "
"expected %s, got %s" % (repr(a), repr(b)))
if not a and not b:
if math.atan2(a, -1.) != math.atan2(b, -1.):
self.fail("zero has wrong sign: expected %s, got %s" %
(repr(a), repr(b)))
# test passes if either the absolute error or the relative
# error is sufficiently small. The defaults amount to an
# error of between 9 ulps and 19 ulps on an IEEE-754 compliant
# machine.
try:
absolute_error = abs(b-a)
except OverflowError:
pass
else:
if absolute_error <= max(abs_err, rel_err * abs(a)):
return
self.fail("%s and %s are not sufficiently close" % (repr(a), repr(b)))
def test_constants(self):
e_expected = 2.71828182845904523536
pi_expected = 3.14159265358979323846
self.assertAlmostEqual(cmath.pi, pi_expected, 9,
self.rAssertAlmostEqual(cmath.pi, pi_expected, 9,
"cmath.pi is %s; should be %s" % (cmath.pi, pi_expected))
self.assertAlmostEqual(cmath.e, e_expected, 9,
self.rAssertAlmostEqual(cmath.e, e_expected, 9,
"cmath.e is %s; should be %s" % (cmath.e, e_expected))
def test_user_object(self):
@ -109,13 +212,13 @@ class CMathTests(unittest.TestCase):
for f in self.test_functions:
# usual usage
self.cAssertAlmostEqual(f(MyComplex(cx_arg)), f(cx_arg))
self.cAssertAlmostEqual(f(MyComplexOS(cx_arg)), f(cx_arg))
self.assertEqual(f(MyComplex(cx_arg)), f(cx_arg))
self.assertEqual(f(MyComplexOS(cx_arg)), f(cx_arg))
# other combinations of __float__ and __complex__
self.cAssertAlmostEqual(f(FloatAndComplex()), f(cx_arg))
self.cAssertAlmostEqual(f(FloatAndComplexOS()), f(cx_arg))
self.cAssertAlmostEqual(f(JustFloat()), f(flt_arg))
self.cAssertAlmostEqual(f(JustFloatOS()), f(flt_arg))
self.assertEqual(f(FloatAndComplex()), f(cx_arg))
self.assertEqual(f(FloatAndComplexOS()), f(cx_arg))
self.assertEqual(f(JustFloat()), f(flt_arg))
self.assertEqual(f(JustFloatOS()), f(flt_arg))
# TypeError should be raised for classes not providing
# either __complex__ or __float__, even if they provide
# __int__, __long__ or __index__. An old-style class
@ -138,7 +241,7 @@ class CMathTests(unittest.TestCase):
# functions, by virtue of providing a __float__ method
for f in self.test_functions:
for arg in [2, 2L, 2.]:
self.cAssertAlmostEqual(f(arg), f(arg.__float__()))
self.assertEqual(f(arg), f(arg.__float__()))
# but strings should give a TypeError
for f in self.test_functions:
@ -182,12 +285,201 @@ class CMathTests(unittest.TestCase):
float_fn = getattr(math, fn)
complex_fn = getattr(cmath, fn)
for v in values:
self.cAssertAlmostEqual(float_fn(v), complex_fn(v))
z = complex_fn(v)
self.rAssertAlmostEqual(float_fn(v), z.real)
self.assertEqual(0., z.imag)
# test two-argument version of log with various bases
for base in [0.5, 2., 10.]:
for v in positive:
self.cAssertAlmostEqual(cmath.log(v, base), math.log(v, base))
z = cmath.log(v, base)
self.rAssertAlmostEqual(math.log(v, base), z.real)
self.assertEqual(0., z.imag)
def test_specific_values(self):
if not float.__getformat__("double").startswith("IEEE"):
return
def rect_complex(z):
"""Wrapped version of rect that accepts a complex number instead of
two float arguments."""
return cmath.rect(z.real, z.imag)
def polar_complex(z):
"""Wrapped version of polar that returns a complex number instead of
two floats."""
return complex(*polar(z))
for id, fn, ar, ai, er, ei, flags in parse_testfile(test_file):
arg = complex(ar, ai)
expected = complex(er, ei)
if fn == 'rect':
function = rect_complex
elif fn == 'polar':
function = polar_complex
else:
function = getattr(cmath, fn)
if 'divide-by-zero' in flags or 'invalid' in flags:
try:
actual = function(arg)
except ValueError:
continue
else:
test_str = "%s: %s(complex(%r, %r))" % (id, fn, ar, ai)
self.fail('ValueError not raised in test %s' % test_str)
if 'overflow' in flags:
try:
actual = function(arg)
except OverflowError:
continue
else:
test_str = "%s: %s(complex(%r, %r))" % (id, fn, ar, ai)
self.fail('OverflowError not raised in test %s' % test_str)
actual = function(arg)
if 'ignore-real-sign' in flags:
actual = complex(abs(actual.real), actual.imag)
expected = complex(abs(expected.real), expected.imag)
if 'ignore-imag-sign' in flags:
actual = complex(actual.real, abs(actual.imag))
expected = complex(expected.real, abs(expected.imag))
# for the real part of the log function, we allow an
# absolute error of up to 2e-15.
if fn in ('log', 'log10'):
real_abs_err = 2e-15
else:
real_abs_err = 5e-323
if not (almostEqualF(expected.real, actual.real,
abs_err = real_abs_err) and
almostEqualF(expected.imag, actual.imag)):
error_message = (
"%s: %s(complex(%r, %r))\n" % (id, fn, ar, ai) +
"Expected: complex(%r, %r)\n" %
(expected.real, expected.imag) +
"Received: complex(%r, %r)\n" %
(actual.real, actual.imag) +
"Received value insufficiently close to expected value.")
self.fail(error_message)
def assertCISEqual(self, a, b):
eps = 1E-7
if abs(a[0] - b[0]) > eps or abs(a[1] - b[1]) > eps:
self.fail((a ,b))
def test_polar(self):
self.assertCISEqual(polar(0), (0., 0.))
self.assertCISEqual(polar(1.), (1., 0.))
self.assertCISEqual(polar(-1.), (1., pi))
self.assertCISEqual(polar(1j), (1., pi/2))
self.assertCISEqual(polar(-1j), (1., -pi/2))
def test_phase(self):
self.assertAlmostEqual(phase(0), 0.)
self.assertAlmostEqual(phase(1.), 0.)
self.assertAlmostEqual(phase(-1.), pi)
self.assertAlmostEqual(phase(-1.+1E-300j), pi)
self.assertAlmostEqual(phase(-1.-1E-300j), -pi)
self.assertAlmostEqual(phase(1j), pi/2)
self.assertAlmostEqual(phase(-1j), -pi/2)
# zeros
self.assertEqual(phase(complex(0.0, 0.0)), 0.0)
self.assertEqual(phase(complex(0.0, -0.0)), -0.0)
self.assertEqual(phase(complex(-0.0, 0.0)), pi)
self.assertEqual(phase(complex(-0.0, -0.0)), -pi)
# infinities
self.assertAlmostEqual(phase(complex(-INF, -0.0)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -2.3)), -pi)
self.assertAlmostEqual(phase(complex(-INF, -INF)), -0.75*pi)
self.assertAlmostEqual(phase(complex(-2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(-0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(0.0, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(2.3, -INF)), -pi/2)
self.assertAlmostEqual(phase(complex(INF, -INF)), -pi/4)
self.assertEqual(phase(complex(INF, -2.3)), -0.0)
self.assertEqual(phase(complex(INF, -0.0)), -0.0)
self.assertEqual(phase(complex(INF, 0.0)), 0.0)
self.assertEqual(phase(complex(INF, 2.3)), 0.0)
self.assertAlmostEqual(phase(complex(INF, INF)), pi/4)
self.assertAlmostEqual(phase(complex(2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-0.0, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-2.3, INF)), pi/2)
self.assertAlmostEqual(phase(complex(-INF, INF)), 0.75*pi)
self.assertAlmostEqual(phase(complex(-INF, 2.3)), pi)
self.assertAlmostEqual(phase(complex(-INF, 0.0)), pi)
# real or imaginary part NaN
for z in complex_nans:
self.assert_(math.isnan(phase(z)))
def test_abs(self):
# zeros
for z in complex_zeros:
self.assertEqual(abs(z), 0.0)
# infinities
for z in complex_infinities:
self.assertEqual(abs(z), INF)
# real or imaginary part NaN
self.assertEqual(abs(complex(NAN, -INF)), INF)
self.assert_(math.isnan(abs(complex(NAN, -2.3))))
self.assert_(math.isnan(abs(complex(NAN, -0.0))))
self.assert_(math.isnan(abs(complex(NAN, 0.0))))
self.assert_(math.isnan(abs(complex(NAN, 2.3))))
self.assertEqual(abs(complex(NAN, INF)), INF)
self.assertEqual(abs(complex(-INF, NAN)), INF)
self.assert_(math.isnan(abs(complex(-2.3, NAN))))
self.assert_(math.isnan(abs(complex(-0.0, NAN))))
self.assert_(math.isnan(abs(complex(0.0, NAN))))
self.assert_(math.isnan(abs(complex(2.3, NAN))))
self.assertEqual(abs(complex(INF, NAN)), INF)
self.assert_(math.isnan(abs(complex(NAN, NAN))))
# result overflows
if float.__getformat__("double").startswith("IEEE"):
self.assertRaises(OverflowError, abs, complex(1.4e308, 1.4e308))
def assertCEqual(self, a, b):
eps = 1E-7
if abs(a.real - b[0]) > eps or abs(a.imag - b[1]) > eps:
self.fail((a ,b))
def test_rect(self):
self.assertCEqual(rect(0, 0), (0, 0))
self.assertCEqual(rect(1, 0), (1., 0))
self.assertCEqual(rect(1, -pi), (-1., 0))
self.assertCEqual(rect(1, pi/2), (0, 1.))
self.assertCEqual(rect(1, -pi/2), (0, -1.))
def test_isnan(self):
self.failIf(cmath.isnan(1))
self.failIf(cmath.isnan(1j))
self.failIf(cmath.isnan(INF))
self.assert_(cmath.isnan(NAN))
self.assert_(cmath.isnan(complex(NAN, 0)))
self.assert_(cmath.isnan(complex(0, NAN)))
self.assert_(cmath.isnan(complex(NAN, NAN)))
self.assert_(cmath.isnan(complex(NAN, INF)))
self.assert_(cmath.isnan(complex(INF, NAN)))
def test_isinf(self):
self.failIf(cmath.isinf(1))
self.failIf(cmath.isinf(1j))
self.failIf(cmath.isinf(NAN))
self.assert_(cmath.isinf(INF))
self.assert_(cmath.isinf(complex(INF, 0)))
self.assert_(cmath.isinf(complex(0, INF)))
self.assert_(cmath.isinf(complex(INF, INF)))
self.assert_(cmath.isinf(complex(NAN, INF)))
self.assert_(cmath.isinf(complex(INF, NAN)))
def test_main():
run_unittest(CMathTests)

21
Lib/test/test_float.py
View File

@ -2,12 +2,12 @@
import unittest, struct
import os
from test import test_support
import math
from math import isinf, isnan
import operator
def isinf(x):
return x * 0.5 == x
def isnan(x):
return x != x
INF = float("inf")
NAN = float("nan")
class FormatFunctionsTestCase(unittest.TestCase):
@ -206,6 +206,17 @@ class InfNanTest(unittest.TestCase):
self.assertEqual(str(1e300 * 1e300 * 0), "nan")
self.assertEqual(str(-1e300 * 1e300 * 0), "nan")
def notest_float_nan(self):
self.assert_(NAN.is_nan())
self.failIf(INF.is_nan())
self.failIf((0.).is_nan())
def notest_float_inf(self):
self.assert_(INF.is_inf())
self.failIf(NAN.is_inf())
self.failIf((0.).is_inf())
def test_main():
test_support.run_unittest(
FormatFunctionsTestCase,

241
Lib/test/test_math.py
View File

@ -4,9 +4,45 @@
from test.test_support import run_unittest, verbose
import unittest
import math
import os
import sys
seps='1e-05'
eps = eval(seps)
eps = 1E-05
NAN = float('nan')
INF = float('inf')
NINF = float('-inf')
# locate file with test values
if __name__ == '__main__':
file = sys.argv[0]
else:
file = __file__
test_dir = os.path.dirname(file) or os.curdir
test_file = os.path.join(test_dir, 'cmath_testcases.txt')
def parse_testfile(fname):
"""Parse a file with test values
Empty lines or lines starting with -- are ignored
yields id, fn, arg_real, arg_imag, exp_real, exp_imag
"""
with open(fname) as fp:
for line in fp:
# skip comment lines and blank lines
if line.startswith('--') or not line.strip():
continue
lhs, rhs = line.split('->')
id, fn, arg_real, arg_imag = lhs.split()
rhs_pieces = rhs.split()
exp_real, exp_imag = rhs_pieces[0], rhs_pieces[1]
flags = rhs_pieces[2:]
yield (id, fn,
float(arg_real), float(arg_imag),
float(exp_real), float(exp_imag),
flags
)
class MathTests(unittest.TestCase):
@ -28,18 +64,57 @@ class MathTests(unittest.TestCase):
self.ftest('acos(-1)', math.acos(-1), math.pi)
self.ftest('acos(0)', math.acos(0), math.pi/2)
self.ftest('acos(1)', math.acos(1), 0)
self.assertRaises(ValueError, math.acos, INF)
self.assertRaises(ValueError, math.acos, NINF)
self.assert_(math.isnan(math.acos(NAN)))
def testAcosh(self):
self.assertRaises(TypeError, math.acosh)
self.ftest('acosh(1)', math.acosh(1), 0)
self.ftest('acosh(2)', math.acosh(2), 1.3169578969248168)
self.assertRaises(ValueError, math.acosh, 0)
self.assertRaises(ValueError, math.acosh, -1)
self.assertEquals(math.acosh(INF), INF)
self.assertRaises(ValueError, math.acosh, NINF)
self.assert_(math.isnan(math.acosh(NAN)))
def testAsin(self):
self.assertRaises(TypeError, math.asin)
self.ftest('asin(-1)', math.asin(-1), -math.pi/2)
self.ftest('asin(0)', math.asin(0), 0)
self.ftest('asin(1)', math.asin(1), math.pi/2)
self.assertRaises(ValueError, math.asin, INF)
self.assertRaises(ValueError, math.asin, NINF)
self.assert_(math.isnan(math.asin(NAN)))
def testAsinh(self):
self.assertRaises(TypeError, math.asinh)
self.ftest('asinh(0)', math.asinh(0), 0)
self.ftest('asinh(1)', math.asinh(1), 0.88137358701954305)
self.ftest('asinh(-1)', math.asinh(-1), -0.88137358701954305)
self.assertEquals(math.asinh(INF), INF)
self.assertEquals(math.asinh(NINF), NINF)
self.assert_(math.isnan(math.asinh(NAN)))
def testAtan(self):
self.assertRaises(TypeError, math.atan)
self.ftest('atan(-1)', math.atan(-1), -math.pi/4)
self.ftest('atan(0)', math.atan(0), 0)
self.ftest('atan(1)', math.atan(1), math.pi/4)
self.ftest('atan(inf)', math.atan(INF), math.pi/2)
self.ftest('atan(-inf)', math.atan(-INF), -math.pi/2)
self.assert_(math.isnan(math.atan(NAN)))
def testAtanh(self):
self.assertRaises(TypeError, math.atan)
self.ftest('atanh(0)', math.atanh(0), 0)
self.ftest('atanh(0.5)', math.atanh(0.5), 0.54930614433405489)
self.ftest('atanh(-0.5)', math.atanh(-0.5), -0.54930614433405489)
self.assertRaises(ValueError, math.atanh, 1)
self.assertRaises(ValueError, math.atanh, -1)
self.assertRaises(ValueError, math.atanh, INF)
self.assertRaises(ValueError, math.atanh, NINF)
self.assert_(math.isnan(math.atanh(NAN)))
def testAtan2(self):
self.assertRaises(TypeError, math.atan2)
@ -61,6 +136,9 @@ class MathTests(unittest.TestCase):
self.ftest('ceil(-0.5)', math.ceil(-0.5), 0)
self.ftest('ceil(-1.0)', math.ceil(-1.0), -1)
self.ftest('ceil(-1.5)', math.ceil(-1.5), -1)
self.assertEquals(math.ceil(INF), INF)
self.assertEquals(math.ceil(NINF), NINF)
self.assert_(math.isnan(math.ceil(NAN)))
class TestCeil(object):
def __float__(self):
@ -75,17 +153,55 @@ class MathTests(unittest.TestCase):
self.assertRaises(TypeError, math.ceil, t)
self.assertRaises(TypeError, math.ceil, t, 0)
if float.__getformat__("double").startswith("IEEE"):
def testCopysign(self):
self.assertRaises(TypeError, math.copysign)
# copysign should let us distinguish signs of zeros
self.assertEquals(copysign(1., 0.), 1.)
self.assertEquals(copysign(1., -0.), -1.)
self.assertEquals(copysign(INF, 0.), INF)
self.assertEquals(copysign(INF, -0.), NINF)
self.assertEquals(copysign(NINF, 0.), INF)
self.assertEquals(copysign(NINF, -0.), NINF)
# and of infinities
self.assertEquals(copysign(1., INF), 1.)
self.assertEquals(copysign(1., NINF), -1.)
self.assertEquals(copysign(INF, INF), INF)
self.assertEquals(copysign(INF, NINF), NINF)
self.assertEquals(copysign(NINF, INF), INF)
self.assertEquals(copysign(NINF, NINF), NINF)
self.assert_(math.isnan(copysign(NAN, 1.)))
self.assert_(math.isnan(copysign(NAN, INF)))
self.assert_(math.isnan(copysign(NAN, NINF)))
self.assert_(math.isnan(copysign(NAN, NAN)))
# copysign(INF, NAN) may be INF or it may be NINF, since
# we don't know whether the sign bit of NAN is set on any
# given platform.
self.assert_(math.isinf(copysign(INF, NAN)))
# similarly, copysign(2., NAN) could be 2. or -2.
self.assertEquals(abs(copysign(2., NAN)), 2.)
def testCos(self):
self.assertRaises(TypeError, math.cos)
self.ftest('cos(-pi/2)', math.cos(-math.pi/2), 0)
self.ftest('cos(0)', math.cos(0), 1)
self.ftest('cos(pi/2)', math.cos(math.pi/2), 0)
self.ftest('cos(pi)', math.cos(math.pi), -1)
try:
self.assert_(math.isnan(math.cos(INF)))
self.assert_(math.isnan(math.cos(NINF)))
except ValueError:
self.assertRaises(ValueError, math.cos, INF)
self.assertRaises(ValueError, math.cos, NINF)
self.assert_(math.isnan(math.cos(NAN)))
def testCosh(self):
self.assertRaises(TypeError, math.cosh)
self.ftest('cosh(0)', math.cosh(0), 1)
self.ftest('cosh(2)-2*cosh(1)**2', math.cosh(2)-2*math.cosh(1)**2, -1) # Thanks to Lambert
self.assertEquals(math.cosh(INF), INF)
self.assertEquals(math.cosh(NINF), INF)
self.assert_(math.isnan(math.cosh(NAN)))
def testDegrees(self):
self.assertRaises(TypeError, math.degrees)
@ -98,6 +214,9 @@ class MathTests(unittest.TestCase):
self.ftest('exp(-1)', math.exp(-1), 1/math.e)
self.ftest('exp(0)', math.exp(0), 1)
self.ftest('exp(1)', math.exp(1), math.e)
self.assertEquals(math.exp(INF), INF)
self.assertEquals(math.exp(NINF), 0.)
self.assert_(math.isnan(math.exp(NAN)))
def testFabs(self):
self.assertRaises(TypeError, math.fabs)
@ -121,6 +240,9 @@ class MathTests(unittest.TestCase):
# This fails on some platforms - so check it here
self.ftest('floor(1.23e167)', math.floor(1.23e167), 1.23e167)
self.ftest('floor(-1.23e167)', math.floor(-1.23e167), -1.23e167)
self.assertEquals(math.ceil(INF), INF)
self.assertEquals(math.ceil(NINF), NINF)
self.assert_(math.isnan(math.floor(NAN)))
class TestFloor(object):
def __float__(self):
@ -143,6 +265,19 @@ class MathTests(unittest.TestCase):
self.ftest('fmod(-10,1)', math.fmod(-10,1), 0)
self.ftest('fmod(-10,0.5)', math.fmod(-10,0.5), 0)
self.ftest('fmod(-10,1.5)', math.fmod(-10,1.5), -1)
self.assert_(math.isnan(math.fmod(NAN, 1.)))
self.assert_(math.isnan(math.fmod(1., NAN)))
self.assert_(math.isnan(math.fmod(NAN, NAN)))
self.assertRaises(ValueError, math.fmod, 1., 0.)
self.assertRaises(ValueError, math.fmod, INF, 1.)
self.assertRaises(ValueError, math.fmod, NINF, 1.)
self.assertRaises(ValueError, math.fmod, INF, 0.)
self.assertEquals(math.fmod(3.0, INF), 3.0)
self.assertEquals(math.fmod(-3.0, INF), -3.0)
self.assertEquals(math.fmod(3.0, NINF), 3.0)
self.assertEquals(math.fmod(-3.0, NINF), -3.0)
self.assertEquals(math.fmod(0.0, 3.0), 0.0)
self.assertEquals(math.fmod(0.0, NINF), 0.0)
def testFrexp(self):
self.assertRaises(TypeError, math.frexp)
@ -157,10 +292,20 @@ class MathTests(unittest.TestCase):
testfrexp('frexp(1)', math.frexp(1), (0.5, 1))
testfrexp('frexp(2)', math.frexp(2), (0.5, 2))
self.assertEquals(math.frexp(INF)[0], INF)
self.assertEquals(math.frexp(NINF)[0], NINF)
self.assert_(math.isnan(math.frexp(NAN)[0]))
def testHypot(self):
self.assertRaises(TypeError, math.hypot)
self.ftest('hypot(0,0)', math.hypot(0,0), 0)
self.ftest('hypot(3,4)', math.hypot(3,4), 5)
self.assertEqual(math.hypot(NAN, INF), INF)
self.assertEqual(math.hypot(INF, NAN), INF)
self.assertEqual(math.hypot(NAN, NINF), INF)
self.assertEqual(math.hypot(NINF, NAN), INF)
self.assert_(math.isnan(math.hypot(1.0, NAN)))
self.assert_(math.isnan(math.hypot(NAN, -2.0)))
def testLdexp(self):
self.assertRaises(TypeError, math.ldexp)
@ -168,6 +313,13 @@ class MathTests(unittest.TestCase):
self.ftest('ldexp(1,1)', math.ldexp(1,1), 2)
self.ftest('ldexp(1,-1)', math.ldexp(1,-1), 0.5)
self.ftest('ldexp(-1,1)', math.ldexp(-1,1), -2)
self.assertRaises(OverflowError, math.ldexp, 1., 1000000)
self.assertRaises(OverflowError, math.ldexp, -1., 1000000)
self.assertEquals(math.ldexp(1., -1000000), 0.)
self.assertEquals(math.ldexp(-1., -1000000), -0.)
self.assertEquals(math.ldexp(INF, 30), INF)
self.assertEquals(math.ldexp(NINF, -213), NINF)
self.assert_(math.isnan(math.ldexp(NAN, 0)))
def testLog(self):
self.assertRaises(TypeError, math.log)
@ -177,12 +329,31 @@ class MathTests(unittest.TestCase):
self.ftest('log(32,2)', math.log(32,2), 5)
self.ftest('log(10**40, 10)', math.log(10**40, 10), 40)
self.ftest('log(10**40, 10**20)', math.log(10**40, 10**20), 2)
self.assertEquals(math.log(INF), INF)
self.assertRaises(ValueError, math.log, NINF)
self.assert_(math.isnan(math.log(NAN)))
def testLog1p(self):
self.assertRaises(TypeError, math.log1p)
self.ftest('log1p(1/e -1)', math.log1p(1/math.e-1), -1)
self.ftest('log1p(0)', math.log1p(0), 0)
self.ftest('log1p(e-1)', math.log1p(math.e-1), 1)
self.ftest('log1p(1)', math.log1p(1), math.log(2))
self.assertEquals(math.log1p(INF), INF)
self.assertRaises(ValueError, math.log1p, NINF)
self.assert_(math.isnan(math.log1p(NAN)))
n= 2**90
self.assertAlmostEquals(math.log1p(n), 62.383246250395075)
self.assertAlmostEquals(math.log1p(n), math.log1p(float(n)))
def testLog10(self):
self.assertRaises(TypeError, math.log10)
self.ftest('log10(0.1)', math.log10(0.1), -1)
self.ftest('log10(1)', math.log10(1), 0)
self.ftest('log10(10)', math.log10(10), 1)
self.assertEquals(math.log(INF), INF)
self.assertRaises(ValueError, math.log10, NINF)
self.assert_(math.isnan(math.log10(NAN)))
def testModf(self):
self.assertRaises(TypeError, math.modf)
@ -195,12 +366,35 @@ class MathTests(unittest.TestCase):
testmodf('modf(1.5)', math.modf(1.5), (0.5, 1.0))
testmodf('modf(-1.5)', math.modf(-1.5), (-0.5, -1.0))
self.assertEquals(math.modf(INF), (0.0, INF))
self.assertEquals(math.modf(NINF), (-0.0, NINF))
modf_nan = math.modf(NAN)
self.assert_(math.isnan(modf_nan[0]))
self.assert_(math.isnan(modf_nan[1]))
def testPow(self):
self.assertRaises(TypeError, math.pow)
self.ftest('pow(0,1)', math.pow(0,1), 0)
self.ftest('pow(1,0)', math.pow(1,0), 1)
self.ftest('pow(2,1)', math.pow(2,1), 2)
self.ftest('pow(2,-1)', math.pow(2,-1), 0.5)
self.assertEqual(math.pow(INF, 1), INF)
self.assertEqual(math.pow(NINF, 1), NINF)
self.assertEqual((math.pow(1, INF)), 1.)
self.assertEqual((math.pow(1, NINF)), 1.)
self.assert_(math.isnan(math.pow(NAN, 1)))
self.assert_(math.isnan(math.pow(2, NAN)))
self.assert_(math.isnan(math.pow(0, NAN)))
self.assertEqual(math.pow(1, NAN), 1)
self.assertEqual(1**NAN, 1)
self.assertEqual(1**INF, 1)
self.assertEqual(1**NINF, 1)
self.assertEqual(1**0, 1)
self.assertEqual(1.**NAN, 1)
self.assertEqual(1.**INF, 1)
self.assertEqual(1.**NINF, 1)
self.assertEqual(1.**0, 1)
def testRadians(self):
self.assertRaises(TypeError, math.radians)
@ -213,29 +407,52 @@ class MathTests(unittest.TestCase):
self.ftest('sin(0)', math.sin(0), 0)
self.ftest('sin(pi/2)', math.sin(math.pi/2), 1)
self.ftest('sin(-pi/2)', math.sin(-math.pi/2), -1)
try:
self.assert_(math.isnan(math.sin(INF)))
self.assert_(math.isnan(math.sin(NINF)))
except ValueError:
self.assertRaises(ValueError, math.sin, INF)
self.assertRaises(ValueError, math.sin, NINF)
self.assert_(math.isnan(math.sin(NAN)))
def testSinh(self):
self.assertRaises(TypeError, math.sinh)
self.ftest('sinh(0)', math.sinh(0), 0)
self.ftest('sinh(1)**2-cosh(1)**2', math.sinh(1)**2-math.cosh(1)**2, -1)
self.ftest('sinh(1)+sinh(-1)', math.sinh(1)+math.sinh(-1), 0)
self.assertEquals(math.sinh(INF), INF)
self.assertEquals(math.sinh(-INF), -INF)
self.assert_(math.isnan(math.sinh(NAN)))
def testSqrt(self):
self.assertRaises(TypeError, math.sqrt)
self.ftest('sqrt(0)', math.sqrt(0), 0)
self.ftest('sqrt(1)', math.sqrt(1), 1)
self.ftest('sqrt(4)', math.sqrt(4), 2)
self.assertEquals(math.sqrt(INF), INF)
self.assertRaises(ValueError, math.sqrt, NINF)
self.assert_(math.isnan(math.sqrt(NAN)))
def testTan(self):
self.assertRaises(TypeError, math.tan)
self.ftest('tan(0)', math.tan(0), 0)
self.ftest('tan(pi/4)', math.tan(math.pi/4), 1)
self.ftest('tan(-pi/4)', math.tan(-math.pi/4), -1)
try:
self.assert_(math.isnan(math.tan(INF)))
self.assert_(math.isnan(math.tan(NINF)))
except:
self.assertRaises(ValueError, math.tan, INF)
self.assertRaises(ValueError, math.tan, NINF)
self.assert_(math.isnan(math.tan(NAN)))
def testTanh(self):
self.assertRaises(TypeError, math.tanh)
self.ftest('tanh(0)', math.tanh(0), 0)
self.ftest('tanh(1)+tanh(-1)', math.tanh(1)+math.tanh(-1), 0)
self.ftest('tanh(inf)', math.tanh(INF), 1)
self.ftest('tanh(-inf)', math.tanh(NINF), -1)
self.assert_(math.isnan(math.tanh(NAN)))
def test_trunc(self):
self.assertEqual(math.trunc(1), 1)
@ -329,9 +546,27 @@ class MathTests(unittest.TestCase):
else:
self.fail("sqrt(-1) didn't raise ValueError")
def test_testfile(self):
if not float.__getformat__("double").startswith("IEEE"):
return
for id, fn, ar, ai, er, ei, flags in parse_testfile(test_file):
# Skip if either the input or result is complex, or if
# flags is nonempty
if ai != 0. or ei != 0. or flags:
continue
if fn in ['rect', 'polar']:
# no real versions of rect, polar
continue
func = getattr(math, fn)
result = func(ar)
self.ftest("%s:%s(%r)" % (id, fn, ar), result, er)
def test_main():
run_unittest(MathTests)
from doctest import DocFileSuite
suite = unittest.TestSuite()
suite.addTest(unittest.makeSuite(MathTests))
suite.addTest(DocFileSuite("ieee754.txt"))
run_unittest(suite)
if __name__ == '__main__':
test_main()

2
Makefile.pre.in
View File

@ -277,6 +277,7 @@ PYTHON_OBJS= \
Python/peephole.o \
Python/pyarena.o \
Python/pyfpe.o \
Python/pymath.o \
Python/pystate.o \
Python/pythonrun.o \
Python/structmember.o \
@ -627,6 +628,7 @@ PYTHON_HEADERS= \
Include/pydebug.h \
Include/pyerrors.h \
Include/pyfpe.h \
Include/pymath.h \
Include/pygetopt.h \
Include/pymem.h \
Include/pyport.h \

11
Misc/NEWS
View File

@ -18,6 +18,12 @@ Core and builtins
Extensions Modules
------------------
- Added phase(z) -> phi, polar(z) -> r, phi and rect(r, phi) -> z to the cmath
module.
- Four new methods were added to the math and cmath modules:
acosh, asinh, atanh and log1p.
- zlib.decompressobj().flush(value) no longer crashes the interpreter when
passed a value less than or equal to zero.
@ -108,6 +114,11 @@ Build
C API
-----
- Added implementation of copysign, acosh, asinh, atanh and log1p
to the new files Include/pymath.h and Python/pymath.h for
platforms which provide the functions through their libm. The
files also contains several helpers and constants for math.
What's New in Python 2.6 alpha 2?
=================================

1019
Modules/cmathmodule.c
View File

File diff suppressed because it is too large Load Diff

442
Modules/mathmodule.c
View File

@ -1,17 +1,60 @@
/* Math module -- standard C math library functions, pi and e */
/* Here are some comments from Tim Peters, extracted from the
discussion attached to http://bugs.python.org/issue1640. They
describe the general aims of the math module with respect to
special values, IEEE-754 floating-point exceptions, and Python
exceptions.
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.
*/
/*
In general, on an IEEE-754 platform the aim is to follow the C99
standard, including Annex 'F', whenever possible. Where the
standard recommends raising the 'divide-by-zero' or 'invalid'
floating-point exceptions, Python should raise a ValueError. Where
the standard recommends raising 'overflow', Python should raise an
OverflowError. In all other circumstances a value should be
returned.
*/
#include "Python.h"
#include "longintrepr.h" /* just for SHIFT */
#ifndef _MSC_VER
#ifndef __STDC__
extern double fmod (double, double);
extern double frexp (double, int *);
extern double ldexp (double, int);
extern double modf (double, double *);
#endif /* __STDC__ */
#endif /* _MSC_VER */
#ifdef _OSF_SOURCE
/* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
extern double copysign(double, double);
@ -52,28 +95,97 @@ is_error(double x)
return result;
}
/*
math_1 is used to wrap a libm function f that takes a double
arguments and returns a double.
The error reporting follows these rules, which are designed to do
the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
platforms.
- a NaN result from non-NaN inputs causes ValueError to be raised
- an infinite result from finite inputs causes OverflowError to be
raised if can_overflow is 1, or raises ValueError if can_overflow
is 0.
- if the result is finite and errno == EDOM then ValueError is
raised
- if the result is finite and nonzero and errno == ERANGE then
OverflowError is raised
The last rule is used to catch overflow on platforms which follow
C89 but for which HUGE_VAL is not an infinity.
For the majority of one-argument functions these rules are enough
to ensure that Python's functions behave as specified in 'Annex F'
of the C99 standard, with the 'invalid' and 'divide-by-zero'
floating-point exceptions mapping to Python's ValueError and the
'overflow' floating-point exception mapping to OverflowError.
math_1 only works for functions that don't have singularities *and*
the possibility of overflow; fortunately, that covers everything we
care about right now.
*/
static PyObject *
math_1(PyObject *arg, double (*func) (double))
math_1(PyObject *arg, double (*func) (double), int can_overflow)
{
double x = PyFloat_AsDouble(arg);
double x, r;
x = PyFloat_AsDouble(arg);
if (x == -1.0 && PyErr_Occurred())
return NULL;
errno = 0;
PyFPE_START_PROTECT("in math_1", return 0)
x = (*func)(x);
PyFPE_END_PROTECT(x)
Py_SET_ERRNO_ON_MATH_ERROR(x);
if (errno && is_error(x))
PyFPE_START_PROTECT("in math_1", return 0);
r = (*func)(x);
PyFPE_END_PROTECT(r);
if (Py_IS_NAN(r)) {
if (!Py_IS_NAN(x))
errno = EDOM;
else
errno = 0;
}
else if (Py_IS_INFINITY(r)) {
if (Py_IS_FINITE(x))
errno = can_overflow ? ERANGE : EDOM;
else
errno = 0;
}
if (errno && is_error(r))
return NULL;
else
return PyFloat_FromDouble(x);
return PyFloat_FromDouble(r);
}
/*
math_2 is used to wrap a libm function f that takes two double
arguments and returns a double.
The error reporting follows these rules, which are designed to do
the right thing on C89/C99 platforms and IEEE 754/non IEEE 754
platforms.
- a NaN result from non-NaN inputs causes ValueError to be raised
- an infinite result from finite inputs causes OverflowError to be
raised.
- if the result is finite and errno == EDOM then ValueError is
raised
- if the result is finite and nonzero and errno == ERANGE then
OverflowError is raised
The last rule is used to catch overflow on platforms which follow
C89 but for which HUGE_VAL is not an infinity.
For most two-argument functions (copysign, fmod, hypot, atan2)
these rules are enough to ensure that Python's functions behave as
specified in 'Annex F' of the C99 standard, with the 'invalid' and
'divide-by-zero' floating-point exceptions mapping to Python's
ValueError and the 'overflow' floating-point exception mapping to
OverflowError.
*/
static PyObject *
math_2(PyObject *args, double (*func) (double, double), char *funcname)
{
PyObject *ox, *oy;
double x, y;
double x, y, r;
if (! PyArg_UnpackTuple(args, funcname, 2, 2, &ox, &oy))
return NULL;
x = PyFloat_AsDouble(ox);
@ -81,19 +193,30 @@ math_2(PyObject *args, double (*func) (double, double), char *funcname)
if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
return NULL;
errno = 0;
PyFPE_START_PROTECT("in math_2", return 0)
x = (*func)(x, y);
PyFPE_END_PROTECT(x)
Py_SET_ERRNO_ON_MATH_ERROR(x);
if (errno && is_error(x))
PyFPE_START_PROTECT("in math_2", return 0);
r = (*func)(x, y);
PyFPE_END_PROTECT(r);
if (Py_IS_NAN(r)) {
if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
errno = EDOM;
else
errno = 0;
}
else if (Py_IS_INFINITY(r)) {
if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
errno = ERANGE;
else
errno = 0;
}
if (errno && is_error(r))
return NULL;
else
return PyFloat_FromDouble(x);
return PyFloat_FromDouble(r);
}
#define FUNC1(funcname, func, docstring) \
#define FUNC1(funcname, func, can_overflow, docstring) \
static PyObject * math_##funcname(PyObject *self, PyObject *args) { \
return math_1(args, func); \
return math_1(args, func, can_overflow); \
}\
PyDoc_STRVAR(math_##funcname##_doc, docstring);
@ -103,55 +226,49 @@ math_2(PyObject *args, double (*func) (double, double), char *funcname)
}\
PyDoc_STRVAR(math_##funcname##_doc, docstring);
FUNC1(acos, acos,
FUNC1(acos, acos, 0,
"acos(x)\n\nReturn the arc cosine (measured in radians) of x.")
FUNC1(asin, asin,
FUNC1(acosh, acosh, 0,
"acosh(x)\n\nReturn the hyperbolic arc cosine (measured in radians) of x.")
FUNC1(asin, asin, 0,
"asin(x)\n\nReturn the arc sine (measured in radians) of x.")
FUNC1(atan, atan,
FUNC1(asinh, asinh, 0,
"asinh(x)\n\nReturn the hyperbolic arc sine (measured in radians) of x.")
FUNC1(atan, atan, 0,
"atan(x)\n\nReturn the arc tangent (measured in radians) of x.")
FUNC2(atan2, atan2,
"atan2(y, x)\n\nReturn the arc tangent (measured in radians) of y/x.\n"
"Unlike atan(y/x), the signs of both x and y are considered.")
FUNC1(ceil, ceil,
FUNC1(atanh, atanh, 0,
"atanh(x)\n\nReturn the hyperbolic arc tangent (measured in radians) of x.")
FUNC1(ceil, ceil, 0,
"ceil(x)\n\nReturn the ceiling of x as a float.\n"
"This is the smallest integral value >= x.")
FUNC1(cos, cos,
"cos(x)\n\nReturn the cosine of x (measured in radians).")
FUNC1(cosh, cosh,
"cosh(x)\n\nReturn the hyperbolic cosine of x.")
#ifdef MS_WINDOWS
# define copysign _copysign
# define HAVE_COPYSIGN 1
#endif
#ifdef HAVE_COPYSIGN
FUNC2(copysign, copysign,
"copysign(x,y)\n\nReturn x with the sign of y.");
#endif
FUNC1(exp, exp,
"copysign(x,y)\n\nReturn x with the sign of y.")
FUNC1(cos, cos, 0,
"cos(x)\n\nReturn the cosine of x (measured in radians).")
FUNC1(cosh, cosh, 1,
"cosh(x)\n\nReturn the hyperbolic cosine of x.")
FUNC1(exp, exp, 1,
"exp(x)\n\nReturn e raised to the power of x.")
FUNC1(fabs, fabs,
FUNC1(fabs, fabs, 0,
"fabs(x)\n\nReturn the absolute value of the float x.")
FUNC1(floor, floor,
FUNC1(floor, floor, 0,
"floor(x)\n\nReturn the floor of x as a float.\n"
"This is the largest integral value <= x.")
FUNC2(fmod, fmod,
"fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
" x % y may differ.")
FUNC2(hypot, hypot,
"hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).")
FUNC2(pow, pow,
"pow(x,y)\n\nReturn x**y (x to the power of y).")
FUNC1(sin, sin,
FUNC1(log1p, log1p, 1,
"log1p(x)\n\nReturn the natural logarithm of 1+x (base e).\n\
The result is computed in a way which is accurate for x near zero.")
FUNC1(sin, sin, 0,
"sin(x)\n\nReturn the sine of x (measured in radians).")
FUNC1(sinh, sinh,
FUNC1(sinh, sinh, 1,
"sinh(x)\n\nReturn the hyperbolic sine of x.")
FUNC1(sqrt, sqrt,
FUNC1(sqrt, sqrt, 0,
"sqrt(x)\n\nReturn the square root of x.")
FUNC1(tan, tan,
FUNC1(tan, tan, 0,
"tan(x)\n\nReturn the tangent of x (measured in radians).")
FUNC1(tanh, tanh,
FUNC1(tanh, tanh, 0,
"tanh(x)\n\nReturn the hyperbolic tangent of x.")
static PyObject *
@ -172,13 +289,17 @@ math_frexp(PyObject *self, PyObject *arg)
double x = PyFloat_AsDouble(arg);
if (x == -1.0 && PyErr_Occurred())
return NULL;
errno = 0;
x = frexp(x, &i);
Py_SET_ERRNO_ON_MATH_ERROR(x);
if (errno && is_error(x))
return NULL;
else
return Py_BuildValue("(di)", x, i);
/* deal with special cases directly, to sidestep platform
differences */
if (Py_IS_NAN(x) || Py_IS_INFINITY(x) || !x) {
i = 0;
}
else {
PyFPE_START_PROTECT("in math_frexp", return 0);
x = frexp(x, &i);
PyFPE_END_PROTECT(x);
}
return Py_BuildValue("(di)", x, i);
}
PyDoc_STRVAR(math_frexp_doc,
@ -191,19 +312,24 @@ PyDoc_STRVAR(math_frexp_doc,
static PyObject *
math_ldexp(PyObject *self, PyObject *args)
{
double x;
double x, r;
int exp;
if (! PyArg_ParseTuple(args, "di:ldexp", &x, &exp))
return NULL;
errno = 0;
PyFPE_START_PROTECT("ldexp", return 0)
x = ldexp(x, exp);
PyFPE_END_PROTECT(x)
Py_SET_ERRNO_ON_MATH_ERROR(x);
if (errno && is_error(x))
PyFPE_START_PROTECT("in math_ldexp", return 0)
r = ldexp(x, exp);
PyFPE_END_PROTECT(r)
if (Py_IS_FINITE(x) && Py_IS_INFINITY(r))
errno = ERANGE;
/* Windows MSVC8 sets errno = EDOM on ldexp(NaN, i);
we unset it to avoid raising a ValueError here. */
if (errno == EDOM)
errno = 0;
if (errno && is_error(r))
return NULL;
else
return PyFloat_FromDouble(x);
return PyFloat_FromDouble(r);
}
PyDoc_STRVAR(math_ldexp_doc,
@ -216,12 +342,10 @@ math_modf(PyObject *self, PyObject *arg)
if (x == -1.0 && PyErr_Occurred())
return NULL;
errno = 0;
PyFPE_START_PROTECT("in math_modf", return 0);
x = modf(x, &y);
Py_SET_ERRNO_ON_MATH_ERROR(x);
if (errno && is_error(x))
return NULL;
else
return Py_BuildValue("(dd)", x, y);
PyFPE_END_PROTECT(x);
return Py_BuildValue("(dd)", x, y);
}
PyDoc_STRVAR(math_modf_doc,
@ -260,7 +384,7 @@ loghelper(PyObject* arg, double (*func)(double), char *funcname)
}
/* Else let libm handle it by itself. */
return math_1(arg, func);
return math_1(arg, func, 0);
}
static PyObject *
@ -303,6 +427,141 @@ math_log10(PyObject *self, PyObject *arg)
PyDoc_STRVAR(math_log10_doc,
"log10(x) -> the base 10 logarithm of x.");
static PyObject *
math_fmod(PyObject *self, PyObject *args)
{
PyObject *ox, *oy;
double r, x, y;
if (! PyArg_UnpackTuple(args, "fmod", 2, 2, &ox, &oy))
return NULL;
x = PyFloat_AsDouble(ox);
y = PyFloat_AsDouble(oy);
if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
return NULL;
/* fmod(x, +/-Inf) returns x for finite x. */
if (Py_IS_INFINITY(y) && Py_IS_FINITE(x))
return PyFloat_FromDouble(x);
errno = 0;
PyFPE_START_PROTECT("in math_fmod", return 0);
r = fmod(x, y);
PyFPE_END_PROTECT(r);
if (Py_IS_NAN(r)) {
if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
errno = EDOM;
else
errno = 0;
}
if (errno && is_error(r))
return NULL;
else
return PyFloat_FromDouble(r);
}
PyDoc_STRVAR(math_fmod_doc,
"fmod(x,y)\n\nReturn fmod(x, y), according to platform C."
" x % y may differ.");
static PyObject *
math_hypot(PyObject *self, PyObject *args)
{
PyObject *ox, *oy;
double r, x, y;
if (! PyArg_UnpackTuple(args, "hypot", 2, 2, &ox, &oy))
return NULL;
x = PyFloat_AsDouble(ox);
y = PyFloat_AsDouble(oy);
if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
return NULL;
/* hypot(x, +/-Inf) returns Inf, even if x is a NaN. */
if (Py_IS_INFINITY(x))
return PyFloat_FromDouble(fabs(x));
if (Py_IS_INFINITY(y))
return PyFloat_FromDouble(fabs(y));
errno = 0;
PyFPE_START_PROTECT("in math_hypot", return 0);
r = hypot(x, y);
PyFPE_END_PROTECT(r);
if (Py_IS_NAN(r)) {
if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
errno = EDOM;
else
errno = 0;
}
else if (Py_IS_INFINITY(r)) {
if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
errno = ERANGE;
else
errno = 0;
}
if (errno && is_error(r))
return NULL;
else
return PyFloat_FromDouble(r);
}
PyDoc_STRVAR(math_hypot_doc,
"hypot(x,y)\n\nReturn the Euclidean distance, sqrt(x*x + y*y).");
/* pow can't use math_2, but needs its own wrapper: the problem is
that an infinite result can arise either as a result of overflow
(in which case OverflowError should be raised) or as a result of
e.g. 0.**-5. (for which ValueError needs to be raised.)
*/
static PyObject *
math_pow(PyObject *self, PyObject *args)
{
PyObject *ox, *oy;
double r, x, y;
if (! PyArg_UnpackTuple(args, "pow", 2, 2, &ox, &oy))
return NULL;
x = PyFloat_AsDouble(ox);
y = PyFloat_AsDouble(oy);
if ((x == -1.0 || y == -1.0) && PyErr_Occurred())
return NULL;
/* 1**x and x**0 return 1., even if x is a NaN or infinity. */
if (x == 1.0 || y == 0.0)
return PyFloat_FromDouble(1.);
errno = 0;
PyFPE_START_PROTECT("in math_pow", return 0);
r = pow(x, y);
PyFPE_END_PROTECT(r);
if (Py_IS_NAN(r)) {
if (!Py_IS_NAN(x) && !Py_IS_NAN(y))
errno = EDOM;
else
errno = 0;
}
/* an infinite result arises either from:
(A) (+/-0.)**negative,
(B) overflow of x**y with both x and y finite (and x nonzero)
(C) (+/-inf)**positive, or
(D) x**inf with |x| > 1, or x**-inf with |x| < 1.
In case (A) we want ValueError to be raised. In case (B)
OverflowError should be raised. In cases (C) and (D) the infinite
result should be returned.
*/
else if (Py_IS_INFINITY(r)) {
if (x == 0.)
errno = EDOM;
else if (Py_IS_FINITE(x) && Py_IS_FINITE(y))
errno = ERANGE;
else
errno = 0;
}
if (errno && is_error(r))
return NULL;
else
return PyFloat_FromDouble(r);
}
PyDoc_STRVAR(math_pow_doc,
"pow(x,y)\n\nReturn x**y (x to the power of y).");
static const double degToRad = Py_MATH_PI / 180.0;
static const double radToDeg = 180.0 / Py_MATH_PI;
@ -356,16 +615,16 @@ PyDoc_STRVAR(math_isinf_doc,
"isinf(x) -> bool\n\
Checks if float x is infinite (positive or negative)");
static PyMethodDef math_methods[] = {
{"acos", math_acos, METH_O, math_acos_doc},
{"acosh", math_acosh, METH_O, math_acosh_doc},
{"asin", math_asin, METH_O, math_asin_doc},
{"asinh", math_asinh, METH_O, math_asinh_doc},
{"atan", math_atan, METH_O, math_atan_doc},
{"atan2", math_atan2, METH_VARARGS, math_atan2_doc},
{"atanh", math_atanh, METH_O, math_atanh_doc},
{"ceil", math_ceil, METH_O, math_ceil_doc},
#ifdef HAVE_COPYSIGN
{"copysign", math_copysign, METH_VARARGS, math_copysign_doc},
#endif
{"cos", math_cos, METH_O, math_cos_doc},
{"cosh", math_cosh, METH_O, math_cosh_doc},
{"degrees", math_degrees, METH_O, math_degrees_doc},
@ -379,6 +638,7 @@ static PyMethodDef math_methods[] = {
{"isnan", math_isnan, METH_O, math_isnan_doc},
{"ldexp", math_ldexp, METH_VARARGS, math_ldexp_doc},
{"log", math_log, METH_VARARGS, math_log_doc},
{"log1p", math_log1p, METH_O, math_log1p_doc},
{"log10", math_log10, METH_O, math_log10_doc},
{"modf", math_modf, METH_O, math_modf_doc},
{"pow", math_pow, METH_VARARGS, math_pow_doc},
@ -400,27 +660,15 @@ PyDoc_STRVAR(module_doc,
PyMODINIT_FUNC
initmath(void)
{
PyObject *m, *d, *v;
PyObject *m;
m = Py_InitModule3("math", math_methods, module_doc);
if (m == NULL)
goto finally;
d = PyModule_GetDict(m);
if (d == NULL)
goto finally;
if (!(v = PyFloat_FromDouble(Py_MATH_PI)))
goto finally;
if (PyDict_SetItemString(d, "pi", v) < 0)
goto finally;
Py_DECREF(v);
PyModule_AddObject(m, "pi", PyFloat_FromDouble(Py_MATH_PI));
PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));
if (!(v = PyFloat_FromDouble(Py_MATH_E)))
goto finally;
if (PyDict_SetItemString(d, "e", v) < 0)
goto finally;
Py_DECREF(v);
finally:
finally:
return;
}

65
Objects/complexobject.c
View File

@ -187,6 +187,38 @@ c_powi(Py_complex x, long n)
}
double
c_abs(Py_complex z)
{
/* sets errno = ERANGE on overflow; otherwise errno = 0 */
double result;
if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {
/* C99 rules: if either the real or the imaginary part is an
infinity, return infinity, even if the other part is a
NaN. */
if (Py_IS_INFINITY(z.real)) {
result = fabs(z.real);
errno = 0;
return result;
}
if (Py_IS_INFINITY(z.imag)) {
result = fabs(z.imag);
errno = 0;
return result;
}
/* either the real or imaginary part is a NaN,
and neither is infinite. Result should be NaN. */
return Py_NAN;
}
result = hypot(z.real, z.imag);
if (!Py_IS_FINITE(result))
errno = ERANGE;
else
errno = 0;
return result;
}
static PyObject *
complex_subtype_from_c_complex(PyTypeObject *type, Py_complex cval)
{
@ -329,8 +361,7 @@ complex_to_buf(char *buf, int bufsz, PyComplexObject *v, int precision)
if (!Py_IS_FINITE(v->cval.imag)) {
if (Py_IS_NAN(v->cval.imag))
strncpy(buf, "nan*j", 6);
/* else if (copysign(1, v->cval.imag) == 1) */
else if (v->cval.imag > 0)
else if (copysign(1, v->cval.imag) == 1)
strncpy(buf, "inf*j", 6);
else
strncpy(buf, "-inf*j", 7);
@ -492,6 +523,7 @@ static PyObject *
complex_div(PyComplexObject *v, PyComplexObject *w)
{
Py_complex quot;
PyFPE_START_PROTECT("complex_div", return 0)
errno = 0;
quot = c_quot(v->cval,w->cval);
@ -655,9 +687,16 @@ static PyObject *
complex_abs(PyComplexObject *v)
{
double result;
PyFPE_START_PROTECT("complex_abs", return 0)
result = hypot(v->cval.real,v->cval.imag);
result = c_abs(v->cval);
PyFPE_END_PROTECT(result)
if (errno == ERANGE) {
PyErr_SetString(PyExc_OverflowError,
"absolute value too large");
return NULL;
}
return PyFloat_FromDouble(result);
}
@ -786,9 +825,29 @@ complex_getnewargs(PyComplexObject *v)
return Py_BuildValue("(D)", &v->cval);
}
#if 0
static PyObject *
complex_is_finite(PyObject *self)
{
Py_complex c;
c = ((PyComplexObject *)self)->cval;
return PyBool_FromLong((long)(Py_IS_FINITE(c.real) &&
Py_IS_FINITE(c.imag)));
}
PyDoc_STRVAR(complex_is_finite_doc,
"complex.is_finite() -> bool\n"
"\n"
"Returns True if the real and the imaginary part is finite.");
#endif
static PyMethodDef complex_methods[] = {
{"conjugate", (PyCFunction)complex_conjugate, METH_NOARGS,
complex_conjugate_doc},
#if 0
{"is_finite", (PyCFunction)complex_is_finite, METH_NOARGS,
complex_is_finite_doc},
#endif
{"__getnewargs__", (PyCFunction)complex_getnewargs, METH_NOARGS},
{NULL, NULL} /* sentinel */
};

601
Objects/doubledigits.c
View File

@ -1,601 +0,0 @@
/* Free-format floating point printer
*
* Based on "Floating-Point Printer Sample Code", by Robert G. Burger,
* http://www.cs.indiana.edu/~burger/fp/index.html
*/
#include "Python.h"
#if defined(__alpha) || defined(__i386) || defined(_M_IX86) || defined(_M_X64) || defined(_M_IA64)
#define LITTLE_ENDIAN_IEEE_DOUBLE
#elif !(defined(__ppc__) || defined(sparc) || defined(__sgi) || defined(_IBMR2) || defined(hpux))
#error unknown machine type
#endif
#if defined(_M_IX86)
#define UNSIGNED64 unsigned __int64
#elif defined(__alpha)
#define UNSIGNED64 unsigned long
#else
#define UNSIGNED64 unsigned long long
#endif
#ifndef U32
#define U32 unsigned int
#endif
/* exponent stored + 1024, hidden bit to left of decimal point */
#define bias 1023
#define bitstoright 52
#define m1mask 0xf
#define hidden_bit 0x100000
#ifdef LITTLE_ENDIAN_IEEE_DOUBLE
struct dblflt {
unsigned int m4: 16;
unsigned int m3: 16;
unsigned int m2: 16;
unsigned int m1: 4;
unsigned int e: 11;
unsigned int s: 1;
};
#else
/* Big Endian IEEE Double Floats */
struct dblflt {
unsigned int s: 1;
unsigned int e: 11;
unsigned int m1: 4;
unsigned int m2: 16;
unsigned int m3: 16;
unsigned int m4: 16;
};
#endif
#define float_radix 2.147483648e9
typedef UNSIGNED64 Bigit;
#define BIGSIZE 24
#define MIN_E -1074
#define MAX_FIVE 325
#define B_P1 ((Bigit)1 << 52)
typedef struct {
int l;
Bigit d[BIGSIZE];
} Bignum;
static Bignum R, S, MP, MM, five[MAX_FIVE];
static Bignum S2, S3, S4, S5, S6, S7, S8, S9;
static int ruf, k, s_n, use_mp, qr_shift, sl, slr;
static void mul10(Bignum *x);
static void big_short_mul(Bignum *x, Bigit y, Bignum *z);
/*
static void print_big(Bignum *x);
*/
static int estimate(int n);
static void one_shift_left(int y, Bignum *z);
static void short_shift_left(Bigit x, int y, Bignum *z);
static void big_shift_left(Bignum *x, int y, Bignum *z);
static int big_comp(Bignum *x, Bignum *y);
static int sub_big(Bignum *x, Bignum *y, Bignum *z);
static void add_big(Bignum *x, Bignum *y, Bignum *z);
static int add_cmp(void);
static int qr(void);
/*static int _PyFloat_Digits(char *buf, double v, int *signum);*/
/*static void _PyFloat_DigitsInit(void);*/
#define ADD(x, y, z, k) {\
Bigit x_add, z_add;\
x_add = (x);\
if ((k))\
z_add = x_add + (y) + 1, (k) = (z_add <= x_add);\
else\
z_add = x_add + (y), (k) = (z_add < x_add);\
(z) = z_add;\
}
#define SUB(x, y, z, b) {\
Bigit x_sub, y_sub;\
x_sub = (x); y_sub = (y);\
if ((b))\
(z) = x_sub - y_sub - 1, b = (y_sub >= x_sub);\
else\
(z) = x_sub - y_sub, b = (y_sub > x_sub);\
}
#define MUL(x, y, z, k) {\
Bigit x_mul, low, high;\
x_mul = (x);\
low = (x_mul & 0xffffffff) * (y) + (k);\
high = (x_mul >> 32) * (y) + (low >> 32);\
(k) = high >> 32;\
(z) = (low & 0xffffffff) | (high << 32);\
}
#define SLL(x, y, z, k) {\
Bigit x_sll = (x);\
(z) = (x_sll << (y)) | (k);\
(k) = x_sll >> (64 - (y));\
}
static void
mul10(Bignum *x)
{
int i, l;
Bigit *p, k;
l = x->l;
for (i = l, p = &x->d[0], k = 0; i >= 0; i--)
MUL(*p, 10, *p++, k);
if (k != 0)
*p = k, x->l = l+1;
}
static void
big_short_mul(Bignum *x, Bigit y, Bignum *z)
{
int i, xl, zl;
Bigit *xp, *zp, k;
U32 high, low;
xl = x->l;
xp = &x->d[0];
zl = xl;
zp = &z->d[0];
high = y >> 32;
low = y & 0xffffffff;
for (i = xl, k = 0; i >= 0; i--, xp++, zp++) {
Bigit xlow, xhigh, z0, t, c, z1;
xlow = *xp & 0xffffffff;
xhigh = *xp >> 32;
z0 = (xlow * low) + k; /* Cout is (z0 < k) */
t = xhigh * low;
z1 = (xlow * high) + t;
c = (z1 < t);
t = z0 >> 32;
z1 += t;
c += (z1 < t);
*zp = (z1 << 32) | (z0 & 0xffffffff);
k = (xhigh * high) + (c << 32) + (z1 >> 32) + (z0 < k);
}
if (k != 0)
*zp = k, zl++;
z->l = zl;
}
/*
static void
print_big(Bignum *x)
{
int i;
Bigit *p;
printf("#x");
i = x->l;
p = &x->d[i];
for (p = &x->d[i]; i >= 0; i--) {
Bigit b = *p--;
printf("%08x%08x", (int)(b >> 32), (int)(b & 0xffffffff));
}
}
*/
static int
estimate(int n)
{
if (n < 0)
return (int)(n*0.3010299956639812);
else
return 1+(int)(n*0.3010299956639811);
}
static void
one_shift_left(int y, Bignum *z)
{
int n, m, i;
Bigit *zp;
n = y / 64;
m = y % 64;
zp = &z->d[0];
for (i = n; i > 0; i--) *zp++ = 0;
*zp = (Bigit)1 << m;
z->l = n;
}
static void
short_shift_left(Bigit x, int y, Bignum *z)
{
int n, m, i, zl;
Bigit *zp;
n = y / 64;
m = y % 64;
zl = n;
zp = &(z->d[0]);
for (i = n; i > 0; i--) *zp++ = 0;
if (m == 0)
*zp = x;
else {
Bigit high = x >> (64 - m);
*zp = x << m;
if (high != 0)
*++zp = high, zl++;
}
z->l = zl;
}
static void
big_shift_left(Bignum *x, int y, Bignum *z)
{
int n, m, i, xl, zl;
Bigit *xp, *zp, k;
n = y / 64;
m = y % 64;
xl = x->l;
xp = &(x->d[0]);
zl = xl + n;
zp = &(z->d[0]);
for (i = n; i > 0; i--) *zp++ = 0;
if (m == 0)
for (i = xl; i >= 0; i--) *zp++ = *xp++;
else {
for (i = xl, k = 0; i >= 0; i--)
SLL(*xp++, m, *zp++, k);
if (k != 0)
*zp = k, zl++;
}
z->l = zl;
}
static int
big_comp(Bignum *x, Bignum *y)
{
int i, xl, yl;
Bigit *xp, *yp;
xl = x->l;
yl = y->l;
if (xl > yl) return 1;
if (xl < yl) return -1;
xp = &x->d[xl];
yp = &y->d[xl];
for (i = xl; i >= 0; i--, xp--, yp--) {
Bigit a = *xp;
Bigit b = *yp;
if (a > b) return 1;
else if (a < b) return -1;
}
return 0;
}
static int
sub_big(Bignum *x, Bignum *y, Bignum *z)
{
int xl, yl, zl, b, i;
Bigit *xp, *yp, *zp;
xl = x->l;
yl = y->l;
if (yl > xl) return 1;
xp = &x->d[0];
yp = &y->d[0];
zp = &z->d[0];
for (i = yl, b = 0; i >= 0; i--)
SUB(*xp++, *yp++, *zp++, b);
for (i = xl-yl; b && i > 0; i--) {
Bigit x_sub;
x_sub = *xp++;
*zp++ = x_sub - 1;
b = (x_sub == 0);
}
for (; i > 0; i--) *zp++ = *xp++;
if (b) return 1;
zl = xl;
while (*--zp == 0) zl--;
z->l = zl;
return 0;
}
static void
add_big(Bignum *x, Bignum *y, Bignum *z)
{
int xl, yl, k, i;
Bigit *xp, *yp, *zp;
xl = x->l;
yl = y->l;
if (yl > xl) {
int tl;
Bignum *tn;
tl = xl; xl = yl; yl = tl;
tn = x; x = y; y = tn;
}
xp = &x->d[0];
yp = &y->d[0];
zp = &z->d[0];
for (i = yl, k = 0; i >= 0; i--)
ADD(*xp++, *yp++, *zp++, k);
for (i = xl-yl; k && i > 0; i--) {
Bigit z_add;
z_add = *xp++ + 1;
k = (z_add == 0);
*zp++ = z_add;
}
for (; i > 0; i--) *zp++ = *xp++;
if (k)
*zp = 1, z->l = xl+1;
else
z->l = xl;
}
static int
add_cmp()
{
int rl, ml, sl, suml;
static Bignum sum;
rl = R.l;
ml = (use_mp ? MP.l : MM.l);
sl = S.l;
suml = rl >= ml ? rl : ml;
if ((sl > suml+1) || ((sl == suml+1) && (S.d[sl] > 1))) return -1;
if (sl < suml) return 1;
add_big(&R, (use_mp ? &MP : &MM), &sum);
return big_comp(&sum, &S);
}
static int
qr()
{
if (big_comp(&R, &S5) < 0)
if (big_comp(&R, &S2) < 0)
if (big_comp(&R, &S) < 0)
return 0;
else {
sub_big(&R, &S, &R);
return 1;
}
else if (big_comp(&R, &S3) < 0) {
sub_big(&R, &S2, &R);
return 2;
}
else if (big_comp(&R, &S4) < 0) {
sub_big(&R, &S3, &R);
return 3;
}
else {
sub_big(&R, &S4, &R);
return 4;
}
else if (big_comp(&R, &S7) < 0)
if (big_comp(&R, &S6) < 0) {
sub_big(&R, &S5, &R);
return 5;
}
else {
sub_big(&R, &S6, &R);
return 6;
}
else if (big_comp(&R, &S9) < 0)
if (big_comp(&R, &S8) < 0) {
sub_big(&R, &S7, &R);
return 7;
}
else {
sub_big(&R, &S8, &R);
return 8;
}
else {
sub_big(&R, &S9, &R);
return 9;
}
}
#define OUTDIG(d) { *buf++ = (d) + '0'; *buf = 0; return k; }
int
_PyFloat_Digits(char *buf, double v, int *signum)
{
struct dblflt *x;
int sign, e, f_n, m_n, i, d, tc1, tc2;
Bigit f;
/* decompose float into sign, mantissa & exponent */
x = (struct dblflt *)&v;
sign = x->s;
e = x->e;
f = (Bigit)(x->m1 << 16 | x->m2) << 32 | (U32)(x->m3 << 16 | x->m4);
if (e != 0) {
e = e - bias - bitstoright;
f |= (Bigit)hidden_bit << 32;
}
else if (f != 0)
/* denormalized */
e = 1 - bias - bitstoright;
*signum = sign;
if (f == 0) {
*buf++ = '0';
*buf = 0;
return 0;
}
ruf = !(f & 1); /* ruf = (even? f) */
/* Compute the scaling factor estimate, k */
if (e > MIN_E)
k = estimate(e+52);
else {
int n;
Bigit y;
for (n = e+52, y = (Bigit)1 << 52; f < y; n--) y >>= 1;
k = estimate(n);
}
if (e >= 0)
if (f != B_P1)
use_mp = 0, f_n = e+1, s_n = 1, m_n = e;
else
use_mp = 1, f_n = e+2, s_n = 2, m_n = e;
else
if ((e == MIN_E) || (f != B_P1))
use_mp = 0, f_n = 1, s_n = 1-e, m_n = 0;
else
use_mp = 1, f_n = 2, s_n = 2-e, m_n = 0;
/* Scale it! */
if (k == 0) {
short_shift_left(f, f_n, &R);
one_shift_left(s_n, &S);
one_shift_left(m_n, &MM);
if (use_mp) one_shift_left(m_n+1, &MP);
qr_shift = 1;
}
else if (k > 0) {
s_n += k;
if (m_n >= s_n)
f_n -= s_n, m_n -= s_n, s_n = 0;
else
f_n -= m_n, s_n -= m_n, m_n = 0;
short_shift_left(f, f_n, &R);
big_shift_left(&five[k-1], s_n, &S);
one_shift_left(m_n, &MM);
if (use_mp) one_shift_left(m_n+1, &MP);
qr_shift = 0;
}
else {
Bignum *power = &five[-k-1];
s_n += k;
big_short_mul(power, f, &S);
big_shift_left(&S, f_n, &R);
one_shift_left(s_n, &S);
big_shift_left(power, m_n, &MM);
if (use_mp) big_shift_left(power, m_n+1, &MP);
qr_shift = 1;
}
/* fixup */
if (add_cmp() <= -ruf) {
k--;
mul10(&R);
mul10(&MM);
if (use_mp) mul10(&MP);
}
/*
printf("\nk = %d\n", k);
printf("R = "); print_big(&R);
printf("\nS = "); print_big(&S);
printf("\nM- = "); print_big(&MM);
if (use_mp) printf("\nM+ = "), print_big(&MP);
putchar('\n');
fflush(0);
*/
if (qr_shift) {
sl = s_n / 64;
slr = s_n % 64;
}
else {
big_shift_left(&S, 1, &S2);
add_big(&S2, &S, &S3);
big_shift_left(&S2, 1, &S4);
add_big(&S4, &S, &S5);
add_big(&S4, &S2, &S6);
add_big(&S4, &S3, &S7);
big_shift_left(&S4, 1, &S8);
add_big(&S8, &S, &S9);
}
again:
if (qr_shift) { /* Take advantage of the fact that S = (ash 1 s_n) */
if (R.l < sl)
d = 0;
else if (R.l == sl) {
Bigit *p;
p = &R.d[sl];
d = *p >> slr;
*p &= ((Bigit)1 << slr) - 1;
for (i = sl; (i > 0) && (*p == 0); i--) p--;
R.l = i;
}
else {
Bigit *p;
p = &R.d[sl+1];
d = *p << (64 - slr) | *(p-1) >> slr;
p--;
*p &= ((Bigit)1 << slr) - 1;
for (i = sl; (i > 0) && (*p == 0); i--) p--;
R.l = i;
}
}
else /* We need to do quotient-remainder */
d = qr();
tc1 = big_comp(&R, &MM) < ruf;
tc2 = add_cmp() > -ruf;
if (!tc1)
if (!tc2) {
mul10(&R);
mul10(&MM);
if (use_mp) mul10(&MP);
*buf++ = d + '0';
goto again;
}
else
OUTDIG(d+1)
else
if (!tc2)
OUTDIG(d)
else {
big_shift_left(&R, 1, &MM);
if (big_comp(&MM, &S) == -1)
OUTDIG(d)
else
OUTDIG(d+1)
}
}
void
_PyFloat_DigitsInit()
{
int n, i, l;
Bignum *b;
Bigit *xp, *zp, k;
five[0].l = l = 0;
five[0].d[0] = 5;
for (n = MAX_FIVE-1, b = &five[0]; n > 0; n--) {
xp = &b->d[0];
b++;
zp = &b->d[0];
for (i = l, k = 0; i >= 0; i--)
MUL(*xp++, 5, *zp++, k);
if (k != 0)
*zp = k, l++;
b->l = l;
}
/*
for (n = 1, b = &five[0]; n <= MAX_FIVE; n++) {
big_shift_left(b++, n, &R);
print_big(&R);
putchar('\n');
}
fflush(0);
*/
}

206
Objects/floatobject.c
View File

@ -16,10 +16,6 @@
#include "formatter_string.h"
#if !defined(__STDC__)
extern double fmod(double, double);
extern double pow(double, double);
#endif
#ifdef _OSF_SOURCE
/* OSF1 5.1 doesn't make this available with XOPEN_SOURCE_EXTENDED defined */
@ -251,11 +247,11 @@ PyFloat_FromString(PyObject *v, char **pend)
p++;
}
if (PyOS_strnicmp(p, "inf", 4) == 0) {
return PyFloat_FromDouble(sign * Py_HUGE_VAL);
Py_RETURN_INF(sign);
}
#ifdef Py_NAN
if(PyOS_strnicmp(p, "nan", 4) == 0) {
return PyFloat_FromDouble(Py_NAN);
Py_RETURN_NAN;
}
#endif
PyOS_snprintf(buffer, sizeof(buffer),
@ -405,110 +401,6 @@ PyFloat_AsStringEx(char *buf, PyFloatObject *v, int precision)
format_float(buf, 100, v, precision);
}
#ifdef Py_BROKEN_REPR
/* The following function is based on Tcl_PrintDouble,
* from tclUtil.c.
*/
#define is_infinite(d) ( (d) > DBL_MAX || (d) < -DBL_MAX )
#define is_nan(d) ((d) != (d))
static void
format_double_repr(char *dst, double value)
{
char *p, c;
int exp;
int signum;
char buffer[30];
/*
* Handle NaN.
*/
if (is_nan(value)) {
strcpy(dst, "nan");
return;
}
/*
* Handle infinities.
*/
if (is_infinite(value)) {
if (value < 0) {
strcpy(dst, "-inf");
} else {
strcpy(dst, "inf");
}
return;
}
/*
* Ordinary (normal and denormal) values.
*/
exp = _PyFloat_Digits(buffer, value, &signum)+1;
if (signum) {
*dst++ = '-';
}
p = buffer;
if (exp < -3 || exp > 17) {
/*
* E format for numbers < 1e-3 or >= 1e17.
*/
*dst++ = *p++;
c = *p;
if (c != '\0') {
*dst++ = '.';
while (c != '\0') {
*dst++ = c;
c = *++p;
}
}
sprintf(dst, "e%+d", exp-1);
} else {
/*
* F format for others.
*/
if (exp <= 0) {
*dst++ = '0';
}
c = *p;
while (exp-- > 0) {
if (c != '\0') {
*dst++ = c;
c = *++p;
} else {
*dst++ = '0';
}
}
*dst++ = '.';
if (c == '\0') {
*dst++ = '0';
} else {
while (++exp < 0) {
*dst++ = '0';
}
while (c != '\0') {
*dst++ = c;
c = *++p;
}
}
*dst++ = '\0';
}
}
static void
format_float_repr(char *buf, PyFloatObject *v)
{
assert(PyFloat_Check(v));
format_double_repr(buf, PyFloat_AS_DOUBLE(v));
}
#endif /* Py_BROKEN_REPR */
/* Macro and helper that convert PyObject obj to a C double and store
the value in dbl; this replaces the functionality of the coercion
slot function. If conversion to double raises an exception, obj is
@ -593,13 +485,8 @@ float_print(PyFloatObject *v, FILE *fp, int flags)
static PyObject *
float_repr(PyFloatObject *v)
{
#ifdef Py_BROKEN_REPR
char buf[30];
format_float_repr(buf, v);
#else
char buf[100];
format_float(buf, sizeof(buf), v, PREC_REPR);
#endif
return PyString_FromString(buf);
}
@ -889,10 +776,13 @@ float_div(PyObject *v, PyObject *w)
double a,b;
CONVERT_TO_DOUBLE(v, a);
CONVERT_TO_DOUBLE(w, b);
#ifdef Py_NAN
if (b == 0.0) {
PyErr_SetString(PyExc_ZeroDivisionError, "float division");
PyErr_SetString(PyExc_ZeroDivisionError,
"float division");
return NULL;
}
#endif
PyFPE_START_PROTECT("divide", return 0)
a = a / b;
PyFPE_END_PROTECT(a)
@ -908,10 +798,13 @@ float_classic_div(PyObject *v, PyObject *w)
if (Py_DivisionWarningFlag >= 2 &&
PyErr_Warn(PyExc_DeprecationWarning, "classic float division") < 0)
return NULL;
#ifdef Py_NAN
if (b == 0.0) {
PyErr_SetString(PyExc_ZeroDivisionError, "float division");
PyErr_SetString(PyExc_ZeroDivisionError,
"float division");
return NULL;
}
#endif
PyFPE_START_PROTECT("divide", return 0)
a = a / b;
PyFPE_END_PROTECT(a)
@ -923,12 +816,15 @@ float_rem(PyObject *v, PyObject *w)
{
double vx, wx;
double mod;
CONVERT_TO_DOUBLE(v, vx);
CONVERT_TO_DOUBLE(w, wx);
CONVERT_TO_DOUBLE(v, vx);
CONVERT_TO_DOUBLE(w, wx);
#ifdef Py_NAN
if (wx == 0.0) {
PyErr_SetString(PyExc_ZeroDivisionError, "float modulo");
PyErr_SetString(PyExc_ZeroDivisionError,
"float modulo");
return NULL;
}
#endif
PyFPE_START_PROTECT("modulo", return 0)
mod = fmod(vx, wx);
/* note: checking mod*wx < 0 is incorrect -- underflows to
@ -1032,6 +928,9 @@ float_pow(PyObject *v, PyObject *w, PyObject *z)
}
return PyFloat_FromDouble(0.0);
}
if (iv == 1.0) { /* 1**w is 1, even 1**inf and 1**nan */
return PyFloat_FromDouble(1.0);
}
if (iv < 0.0) {
/* Whether this is an error is a mess, and bumps into libm
* bugs so we have to figure it out ourselves.
@ -1122,6 +1021,57 @@ float_coerce(PyObject **pv, PyObject **pw)
return 1; /* Can't do it */
}
static PyObject *
float_is_integer(PyObject *v)
{
double x = PyFloat_AsDouble(v);
PyObject *o;
if (x == -1.0 && PyErr_Occurred())
return NULL;
if (!Py_IS_FINITE(x))
Py_RETURN_FALSE;
PyFPE_START_PROTECT("is_integer", return NULL)
o = (floor(x) == x) ? Py_True : Py_False;
PyFPE_END_PROTECT(x)
if (errno != 0) {
PyErr_SetFromErrno(errno == ERANGE ? PyExc_OverflowError :
PyExc_ValueError);
return NULL;
}
Py_INCREF(o);
return o;
}
#if 0
static PyObject *
float_is_inf(PyObject *v)
{
double x = PyFloat_AsDouble(v);
if (x == -1.0 && PyErr_Occurred())
return NULL;
return PyBool_FromLong((long)Py_IS_INFINITY(x));
}
static PyObject *
float_is_nan(PyObject *v)
{
double x = PyFloat_AsDouble(v);
if (x == -1.0 && PyErr_Occurred())
return NULL;
return PyBool_FromLong((long)Py_IS_NAN(x));
}
static PyObject *
float_is_finite(PyObject *v)
{
double x = PyFloat_AsDouble(v);
if (x == -1.0 && PyErr_Occurred())
return NULL;
return PyBool_FromLong((long)Py_IS_FINITE(x));
}
#endif
static PyObject *
float_trunc(PyObject *v)
{
@ -1480,12 +1430,22 @@ PyDoc_STRVAR(float__format__doc,
static PyMethodDef float_methods[] = {
{"conjugate", (PyCFunction)float_float, METH_NOARGS,
{"conjugate", (PyCFunction)float_float, METH_NOARGS,
"Returns self, the complex conjugate of any float."},
{"__trunc__", (PyCFunction)float_trunc, METH_NOARGS,
"Returns the Integral closest to x between 0 and x."},
{"as_integer_ratio", (PyCFunction)float_as_integer_ratio, METH_NOARGS,
float_as_integer_ratio_doc},
{"is_integer", (PyCFunction)float_is_integer, METH_NOARGS,
"Returns True if the float is an integer."},
#if 0
{"is_inf", (PyCFunction)float_is_inf, METH_NOARGS,
"Returns True if the float is positive or negative infinite."},
{"is_finite", (PyCFunction)float_is_finite, METH_NOARGS,
"Returns True if the float is finite, neither infinite nor NaN."},
{"is_nan", (PyCFunction)float_is_nan, METH_NOARGS,
"Returns True if the float is not a number (NaN)."},
#endif
{"__getnewargs__", (PyCFunction)float_getnewargs, METH_NOARGS},
{"__getformat__", (PyCFunction)float_getformat,
METH_O|METH_CLASS, float_getformat_doc},
@ -1646,10 +1606,6 @@ _PyFloat_Init(void)
double_format = detected_double_format;
float_format = detected_float_format;
#ifdef Py_BROKEN_REPR
/* Initialize floating point repr */
_PyFloat_DigitsInit();
#endif
/* Init float info */
if (FloatInfoType.tp_name == 0)
PyStructSequence_InitType(&FloatInfoType, &floatinfo_desc);

12
Objects/intobject.c
View File

@ -1143,9 +1143,21 @@ int__format__(PyObject *self, PyObject *args)
return NULL;
}
#if 0
static PyObject *
int_is_finite(PyObject *v)
{
Py_RETURN_TRUE;
}
#endif
static PyMethodDef int_methods[] = {
{"conjugate", (PyCFunction)int_int, METH_NOARGS,
"Returns self, the complex conjugate of any int."},
#if 0
{"is_finite", (PyCFunction)int_is_finite, METH_NOARGS,
"Returns always True."},
#endif
{"__trunc__", (PyCFunction)int_int, METH_NOARGS,
"Truncating an Integral returns itself."},
{"__getnewargs__", (PyCFunction)int_getnewargs, METH_NOARGS},

12
Objects/longobject.c
View File

@ -3441,9 +3441,21 @@ long__format__(PyObject *self, PyObject *args)
return NULL;
}
#if 0
static PyObject *
long_is_finite(PyObject *v)
{
Py_RETURN_TRUE;
}
#endif
static PyMethodDef long_methods[] = {
{"conjugate", (PyCFunction)long_long, METH_NOARGS,
"Returns self, the complex conjugate of any long."},
#if 0
{"is_finite", (PyCFunction)long_is_finite, METH_NOARGS,
"Returns always True."},
#endif
{"__trunc__", (PyCFunction)long_long, METH_NOARGS,
"Truncating an Integral returns itself."},
{"__getnewargs__", (PyCFunction)long_getnewargs, METH_NOARGS},

5
PC/pyconfig.h
View File

@ -211,12 +211,13 @@ typedef _W64 int ssize_t;
#endif /* MS_WIN32 && !MS_WIN64 */
typedef int pid_t;
#define hypot _hypot
#include <float.h>
#define Py_IS_NAN _isnan
#define Py_IS_INFINITY(X) (!_finite(X) && !_isnan(X))
#define Py_IS_FINITE(X) _finite(X)
#define copysign _copysign
#define hypot _hypot
#endif /* _MSC_VER */
@ -396,7 +397,7 @@ Py_NO_ENABLE_SHARED to find out. Also support MS_NO_COREDLL for b/w compat */
/* Fairly standard from here! */
/* Define to 1 if you have the `copysign' function. */
/* #define HAVE_COPYSIGN 1*/
#define HAVE_COPYSIGN 1
/* Define to 1 if you have the `isinf' function. */
#define HAVE_ISINF 1

25
Python/hypot.c
View File

@ -1,25 +0,0 @@
/* hypot() replacement */
#include "Python.h"
#ifndef HAVE_HYPOT
double hypot(double x, double y)
{
double yx;
x = fabs(x);
y = fabs(y);
if (x < y) {
double temp = x;
x = y;
y = temp;
}
if (x == 0.)
return 0.;
else {
yx = y/x;
return x*sqrt(1.+yx*yx);
}
}
#endif /* HAVE_HYPOT */

232
Python/pymath.c Normal file
View File

@ -0,0 +1,232 @@
#include "Python.h"
#ifndef HAVE_HYPOT
double hypot(double x, double y)
{
double yx;
x = fabs(x);
y = fabs(y);
if (x < y) {
double temp = x;
x = y;
y = temp;
}
if (x == 0.)
return 0.;
else {
yx = y/x;
return x*sqrt(1.+yx*yx);
}
}
#endif /* HAVE_HYPOT */
#ifndef HAVE_COPYSIGN
static double
copysign(double x, double y)
{
/* use atan2 to distinguish -0. from 0. */
if (y > 0. || (y == 0. && atan2(y, -1.) > 0.)) {
return fabs(x);
} else {
return -fabs(x);
}
}
#endif /* HAVE_COPYSIGN */
#ifndef HAVE_LOG1P
double
log1p(double x)
{
/* For x small, we use the following approach. Let y be the nearest
float to 1+x, then
1+x = y * (1 - (y-1-x)/y)
so log(1+x) = log(y) + log(1-(y-1-x)/y). Since (y-1-x)/y is tiny,
the second term is well approximated by (y-1-x)/y. If abs(x) >=
DBL_EPSILON/2 or the rounding-mode is some form of round-to-nearest
then y-1-x will be exactly representable, and is computed exactly
by (y-1)-x.
If abs(x) < DBL_EPSILON/2 and the rounding mode is not known to be
round-to-nearest then this method is slightly dangerous: 1+x could
be rounded up to 1+DBL_EPSILON instead of down to 1, and in that
case y-1-x will not be exactly representable any more and the
result can be off by many ulps. But this is easily fixed: for a
floating-point number |x| < DBL_EPSILON/2., the closest
floating-point number to log(1+x) is exactly x.
*/
double y;
if (fabs(x) < DBL_EPSILON/2.) {
return x;
} else if (-0.5 <= x && x <= 1.) {
/* WARNING: it's possible than an overeager compiler
will incorrectly optimize the following two lines
to the equivalent of "return log(1.+x)". If this
happens, then results from log1p will be inaccurate
for small x. */
y = 1.+x;
return log(y)-((y-1.)-x)/y;
} else {
/* NaNs and infinities should end up here */
return log(1.+x);
}
}
#endif /* HAVE_LOG1P */
/*
* ====================================================
* 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.
* ====================================================
*/
static const double ln2 = 6.93147180559945286227E-01;
static const double two_pow_m28 = 3.7252902984619141E-09; /* 2**-28 */
static const double two_pow_p28 = 268435456.0; /* 2**28 */
static const double zero = 0.0;
/* asinh(x)
* Method :
* Based on
* asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
* we have
* asinh(x) := x if 1+x*x=1,
* := sign(x)*(log(x)+ln2)) for large |x|, else
* := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
* := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
*/
#ifndef HAVE_ASINH
double
asinh(double x)
{
double w;
double absx = fabs(x);
if (Py_IS_NAN(x) || Py_IS_INFINITY(x)) {
return x+x;
}
if (absx < two_pow_m28) { /* |x| < 2**-28 */
return x; /* return x inexact except 0 */
}
if (absx > two_pow_p28) { /* |x| > 2**28 */
w = log(absx)+ln2;
}
else if (absx > 2.0) { /* 2 < |x| < 2**28 */
w = log(2.0*absx + 1.0 / (sqrt(x*x + 1.0) + absx));
}
else { /* 2**-28 <= |x| < 2= */
double t = x*x;
w = log1p(absx + t / (1.0 + sqrt(1.0 + t)));
}
return copysign(w, x);
}
#endif /* HAVE_ASINH */
/* acosh(x)
* Method :
* Based on
* acosh(x) = log [ x + sqrt(x*x-1) ]
* we have
* acosh(x) := log(x)+ln2, if x is large; else
* acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
* acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
*
* Special cases:
* acosh(x) is NaN with signal if x<1.
* acosh(NaN) is NaN without signal.
*/
#ifndef HAVE_ACOSH
double
acosh(double x)
{
if (Py_IS_NAN(x)) {
return x+x;
}
if (x < 1.) { /* x < 1; return a signaling NaN */
errno = EDOM;
#ifdef Py_NAN
return Py_NAN;
#else
return (x-x)/(x-x);
#endif
}
else if (x >= two_pow_p28) { /* x > 2**28 */
if (Py_IS_INFINITY(x)) {
return x+x;
} else {
return log(x)+ln2; /* acosh(huge)=log(2x) */
}
}
else if (x == 1.) {
return 0.0; /* acosh(1) = 0 */
}
else if (x > 2.) { /* 2 < x < 2**28 */
double t = x*x;
return log(2.0*x - 1.0 / (x + sqrt(t - 1.0)));
}
else { /* 1 < x <= 2 */
double t = x - 1.0;
return log1p(t + sqrt(2.0*t + t*t));
}
}
#endif /* HAVE_ACOSH */
/* atanh(x)
* Method :
* 1.Reduced x to positive by atanh(-x) = -atanh(x)
* 2.For x>=0.5
* 1 2x x
* atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
* 2 1 - x 1 - x
*
* For x<0.5
* atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
*
* Special cases:
* atanh(x) is NaN if |x| >= 1 with signal;
* atanh(NaN) is that NaN with no signal;
*
*/
#ifndef HAVE_ATANH
double
atanh(double x)
{
double absx;
double t;
if (Py_IS_NAN(x)) {
return x+x;
}
absx = fabs(x);
if (absx >= 1.) { /* |x| >= 1 */
errno = EDOM;
#ifdef Py_NAN
return Py_NAN;
#else
return x/zero;
#endif
}
if (absx < two_pow_m28) { /* |x| < 2**-28 */
return x;
}
if (absx < 0.5) { /* |x| < 0.5 */
t = absx+absx;
t = 0.5 * log1p(t + t*absx / (1.0 - absx));
}
else { /* 0.5 <= |x| <= 1.0 */
t = 0.5 * log1p((absx + absx) / (1.0 - absx));
}
return copysign(t, x);
}
#endif /* HAVE_ATANH */

2
configure vendored
View File

@ -1,5 +1,5 @@
#! /bin/sh
# From configure.in Revision: 61992 .
# From configure.in Revision: 62145 .
# Guess values for system-dependent variables and create Makefiles.
# Generated by GNU Autoconf 2.61 for python 2.6.
#