بهترین مفسر پایتون برای برنامه نویسی — راهنمای کاربردی

1#include "Python.h"
2#include 
3#include "_math.h"
4
5
6#if !defined(HAVE_ACOSH) || !defined(HAVE_ASINH)
7static const double ln2 = 6.93147180559945286227E-01;
8static const double two_pow_p28 = 268435456.0; /* 2**28 */
9#endif
10#if !defined(HAVE_ASINH) || !defined(HAVE_ATANH)
11static const double two_pow_m28 = 3.7252902984619141E-09; /* 2**-28 */
12#endif
13#if !defined(HAVE_ATANH) && !defined(Py_NAN)
14static const double zero = 0.0;
15#endif
16
17double
18_Py_acosh(double x)
19{
20    if (Py_IS_NAN(x)) {
21        return x+x;
22    }
23    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 */
24        if (Py_IS_INFINITY(x)) {
25            return x+x;
26        }
27        else {
28            return log(x) + ln2;          /* acosh(huge)=log(2x) */
29        }
30    }
31    else if (x == 1.) {
32        return 0.0;                     /* acosh(1) = 0 */
33    }
34    else if (x > 2.) {                  /* 2 < x < 2**28 */
35        double t = x * x;
36        return log(2.0 * x - 1.0 / (x + sqrt(t - 1.0)));
37    }
38    else {                              /* 1 < x <= 2 */ double t = x - 1.0; return m_log1p(t + sqrt(2.0 * t + t * t)); } } #endif /* HAVE_ACOSH */ #ifndef HAVE_ASINH /* 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
39 *               := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
40 */
41
42double
43_Py_asinh(double x)
44{
45    double w;
46    double absx = fabs(x);
47
48    if (Py_IS_NAN(x) || Py_IS_INFINITY(x)) {
49        return x+x;
50    }
51    if (absx < two_pow_m28) {           /* |x| < 2**-28 */ return x; /* return x inexact except 0 */ } if (absx > two_pow_p28) {           /* |x| > 2**28 */
52        w = log(absx) + ln2;
53    }
54    else if (absx > 2.0) {              /* 2 < |x| < 2**28 */
55        w = log(2.0 * absx + 1.0 / (sqrt(x * x + 1.0) + absx));
56    }
57    else {                              /* 2**-28 <= |x| < 2= */ double t = x*x; w = m_log1p(absx + t / (1.0 + sqrt(1.0 + t))); } return copysign(w, x); } #endif /* HAVE_ASINH */ #ifndef HAVE_ATANH /* atanh(x) * Method : * 1.Reduced x to positive by atanh(-x) = -atanh(x) * 2.For x>=0.5
58 *                  1              2x                          x
59 *      atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * -------)
60 *                  2             1 - x                      1 - x
61 *
62 *      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;
63 *      atanh(NaN) is that NaN with no signal;
64 *
65 */
66
67double
68_Py_atanh(double x)
69{
70    double absx;
71    double t;
72
73    if (Py_IS_NAN(x)) {
74        return x+x;
75    }
76    absx = fabs(x);
77    if (absx >= 1.) {                   /* |x| >= 1 */
78        errno = EDOM;
79#ifdef Py_NAN
80        return Py_NAN;
81#else
82        return x / zero;
83#endif
84    }
85    if (absx < two_pow_m28) {           /* |x| < 2**-28 */
86        return x;
87    }
88    if (absx < 0.5) {                   /* |x| < 0.5 */
89        t = absx+absx;
90        t = 0.5 * m_log1p(t + t*absx / (1.0 - absx));
91    }
92    else {                              /* 0.5 <= |x| <= 1.0 */ t = 0.5 * m_log1p((absx + absx) / (1.0 - absx)); } return copysign(t, x); } #endif /* HAVE_ATANH */ #ifndef HAVE_EXPM1 /* Mathematically, expm1(x) = exp(x) - 1. The expm1 function is designed to avoid the significant loss of precision that arises from direct evaluation of the expression exp(x) - 1, for x near 0. */ double _Py_expm1(double x) { /* For abs(x) >= log(2), it's safe to evaluate exp(x) - 1 directly; this
93       also works fine for infinities and nans.
94       For smaller x, we can use a method due to Kahan that achieves close to
95       full accuracy.
96    
97    double y;
98    if (fabs(x) < DBL_EPSILON / 2.) {
99        return x;
100    }
101    else if (-0.5 <= x && x <= 1.) {
102        /* WARNING: it's possible that an overeager compiler
103           will incorrectly optimize the following two lines
104           to the equivalent of "return log(1.+x)". If this
105           happens, then results from log1p will be inaccurate
106           for small x. */
107        y = 1.+x;
108        return log(y) - ((y - 1.) - x) / y;
109    }
110    else {
111        /* NaNs and infinities should end up here */
112        return log(1.+x);
113    }
114}
115

1# Package initialisation
2from pypy.interpreter.mixedmodule import MixedModule
3
4class Module(MixedModule):
5    appleveldefs = {
6       'factorial' : 'app_math.factorial'
7    }
8
9    interpleveldefs = {
10       'e'              : 'interp_math.get(space).w_e', 
11       'pi'             : 'interp_math.get(space).w_pi', 
12       'pow'            : 'interp_math.pow',
13       'cosh'           : 'interp_math.cosh',
14       'copysign'       : 'interp_math.copysign',
15       'ldexp'          : 'interp_math.ldexp',
16       'hypot'          : 'interp_math.hypot',
17       'tan'            : 'interp_math.tan',
18       'asin'           : 'interp_math.asin',
19       'fabs'           : 'interp_math.fabs',
20       'floor'          : 'interp_math.floor',
21       'sqrt'           : 'interp_math.sqrt',
22       'frexp'          : 'interp_math.frexp',
23       'degrees'        : 'interp_math.degrees',
24       'log'            : 'interp_math.log',
25       'log10'          : 'interp_math.log10',
26       'fmod'           : 'interp_math.fmod',
27       'atan'           : 'interp_math.atan',
28       'ceil'           : 'interp_math.ceil',
29       'sinh'           : 'interp_math.sinh',
30       'cos'            : 'interp_math.cos',
31       'tanh'           : 'interp_math.tanh',
32       'radians'        : 'interp_math.radians',
33       'sin'            : 'interp_math.sin',
34       'atan2'          : 'interp_math.atan2',
35       'modf'           : 'interp_math.modf',
36       'exp'            : 'interp_math.exp',
37       'expm1'          : 'interp_math.expm1',
38       'acos'           : 'interp_math.acos',
39       'isinf'          : 'interp_math.isinf',
40       'isnan'          : 'interp_math.isnan',
41       'trunc'          : 'interp_math.trunc',
42       'fsum'           : 'interp_math.fsum',
43       'asinh'          : 'interp_math.asinh',
44       'acosh'          : 'interp_math.acosh',
45       'atanh'          : 'interp_math.atanh',
46       'log1p'          : 'interp_math.log1p',
47       'expm1'          : 'interp_math.expm1',
48       'erf'            : 'interp_math.erf',
49       'erfc'           : 'interp_math.erfc',
50       'gamma'          : 'interp_math.gamma',
51       'lgamma'         : 'interp_math.lgamma',
52}
53
54import sys, math
55
56INFINITY = 1e200 * 1e200
57NAN = abs(INFINITY / INFINITY)
58
59def isinf(x):
60    return x == INFINITY or x == -INFINITY
61
62def isnan(v):
63    return v != v
64
65def positiveinf(x):
66    return isinf(x) and x > 0.0
67
68def negativeinf(x):
69    return isinf(x) and x < 0.0
70
71def finite(x):
72    return not isinf(x) and not isnan(x)
73
74unary_math_functions = ['acos', 'asin', 'atan',
75                        'ceil', 'cos', 'cosh', 'exp', 'fabs', 'floor',
76                        'sin', 'sinh', 'sqrt', 'tan', 'tanh', 'log', 'log10',
77                        'acosh', 'asinh', 'atanh', 'log1p', 'expm1']
78binary_math_functions = ['atan2', 'fmod', 'hypot', 'pow']
79
80class MathTests:
81
82    REGCASES = []
83    for name in unary_math_functions:
84        try:
85            input, output = (0.3,), getattr(math, name)(0.3)
86        except AttributeError:
87            # cannot test this function
88            continue
89        except ValueError:
90            input, output = (1.3,), getattr(math, name)(1.3)
91        REGCASES.append((name, input, output))
92
93    IRREGCASES = [
94        ('atan2', (0.31, 0.123), math.atan2(0.31, 0.123)),
95        ('fmod',  (0.31, 0.123), math.fmod(0.31, 0.123)),
96        ('hypot', (0.31, 0.123), math.hypot(0.31, 0.123)),
97        ('pow',   (0.31, 0.123), math.pow(0.31, 0.123)),
98        ('pow',   (-0.31, 0.123), ValueError),
99        ('pow',   (-0.5, 2.0), 0.25),
100        ('pow',   (-0.5, 1.0), -0.5),
101        ('pow',   (-0.5, 0.0), 1.0),
102        ('pow',   (-0.5, -1.0), -2.0),
103        ('ldexp', (3.375, 2), 13.5),
104        ('ldexp', (1.0, -10000), 0.0),   # underflow
105        ('frexp', (-1.25,), lambda x: x == (-0.625, 1)),
106        ('modf',  (4.25,), lambda x: x == (0.25, 4.0)),
107        ('modf',  (-4.25,), lambda x: x == (-0.25, -4.0)),
108        ('copysign', (1.5, 0.0), 1.5),
109        ('copysign', (1.5, -0.0), -1.5),
110        ('copysign', (1.5, INFINITY), 1.5),
111        ('copysign', (1.5, -INFINITY), -1.5),
112        ('copysign', (1.5, NAN), 1.5),
113        ('copysign', (1.75, -NAN), -1.75),      # special case for -NAN here
114        ]
115
116    OVFCASES = [
117        ('cosh', (9999.9,), OverflowError),
118        ('sinh', (9999.9,), OverflowError),
119        ('exp', (9999.9,), OverflowError),
120        ('pow', (10.0, 40000.0), OverflowError),
121        ('ldexp', (10.0, 40000), OverflowError),
122        ('log', (0.0,), ValueError), #cpython does it this way
123        ('log1p', (-1.0,), OverflowError),
124        ('log', (-1.,), ValueError),
125        ('log10', (0.0,), ValueError), #cpython does it this way
126        ]
127
128    INFCASES = [
129        ('acos', (INFINITY,), ValueError),
130        ('acos', (-INFINITY,), ValueError),
131        ('asin', (INFINITY,), ValueError),
132        ('asin', (-INFINITY,), ValueError),
133        ('atan', (INFINITY,), math.pi / 2),
134        ('atan', (-INFINITY,), -math.pi / 2),
135        ('atanh', (INFINITY,), ValueError),
136        ('atanh', (-INFINITY,), ValueError),
137        ('ceil', (INFINITY,), positiveinf),
138        ('ceil', (-INFINITY,), negativeinf),
139        ('cos', (INFINITY,), ValueError),
140        ('cos', (-INFINITY,), ValueError),
141        ('cosh', (INFINITY,), positiveinf),
142        ('cosh', (-INFINITY,), positiveinf),
143        ('exp', (INFINITY,), positiveinf),
144        ('exp', (-INFINITY,), 0.0),
145        ('fabs', (INFINITY,), positiveinf),
146        ('fabs', (-INFINITY,), positiveinf),
147        ('floor', (INFINITY,), positiveinf),
148        ('floor', (-INFINITY,), negativeinf),
149        ('sin', (INFINITY,), ValueError),
150        ('sin', (-INFINITY,), ValueError),
151        ('sinh', (INFINITY,), positiveinf),
152        ('sinh', (-INFINITY,), negativeinf),
153        ('sqrt', (INFINITY,), positiveinf),
154        ('sqrt', (-INFINITY,), ValueError),
155        ('tan', (INFINITY,), ValueError),
156        ('tan', (-INFINITY,), ValueError),
157        ('tanh', (INFINITY,), 1.0),
158        ('tanh', (-INFINITY,), -1.0),
159        ('log', (INFINITY,), positiveinf),
160        ('log', (-INFINITY,), ValueError),
161        ('log10', (INFINITY,), positiveinf),
162        ('log10', (-INFINITY,), ValueError),
163        ('frexp', (INFINITY,), lambda x: isinf(x[0])),
164        ('ldexp', (INFINITY, 3), positiveinf),
165        ('ldexp', (-INFINITY, 3), negativeinf),
166        ('modf',  (INFINITY,), lambda x: positiveinf(x[1])),
167        ('modf',  (-INFINITY,), lambda x: negativeinf(x[1])),
168        ('pow', (INFINITY, 0.0), 1.0),
169        ('pow', (INFINITY, 0.001), positiveinf),
170        ('pow', (INFINITY, -0.001), 0.0),
171        ('pow', (-INFINITY, 0.0), 1.0),
172        ('pow', (-INFINITY, 0.001), positiveinf),
173        ('pow', (-INFINITY, -0.001), 0.0),
174        ('pow', (-INFINITY, 3.0), negativeinf),
175        ('pow', (-INFINITY, 6.0), positiveinf),
176        ('pow', (-INFINITY, -13.0), -0.0),
177        ('pow', (-INFINITY, -128.0), 0.0),
178        ('pow', (1.001, INFINITY), positiveinf),
179        ('pow', (1.0,   INFINITY), 1.0),
180        ('pow', (0.999, INFINITY), 0.0),
181        ('pow', (-0.999,INFINITY), 0.0),
182        #('pow', (-1.0, INFINITY), 1.0, but strange, could also be -1.0),
183        ('pow', (-1.001,INFINITY), positiveinf),
184        ('pow', (1.001, -INFINITY), 0.0),
185        ('pow', (1.0,   -INFINITY), 1.0),
186        #('pow', (0.999, -INFINITY), positiveinf, but get OverflowError),
187        #('pow', (INFINITY, INFINITY), positiveinf, but get OverflowError),
188        ('pow', (INFINITY, -INFINITY), 0.0),
189        ('pow', (-INFINITY, INFINITY), positiveinf),
190        ]
191
192    IRREGERRCASES = [
193        #
194        ('atan2', (INFINITY, -2.3), math.pi / 2),
195        ('atan2', (INFINITY, 0.0), math.pi / 2),
196        ('atan2', (INFINITY, 3.0), math.pi / 2),
197        #('atan2', (INFINITY, INFINITY), ?strange),
198        ('atan2', (2.1, INFINITY), 0.0),
199        ('atan2', (0.0, INFINITY), 0.0),
200        ('atan2', (-0.1, INFINITY), -0.0),
201        ('atan2', (-INFINITY, 0.4), -math.pi / 2),
202        ('atan2', (2.1, -INFINITY), math.pi),
203        ('atan2', (0.0, -INFINITY), math.pi),
204        ('atan2', (-0.1, -INFINITY), -math.pi),
205        #
206        ('fmod', (INFINITY, 1.0), ValueError),
207        ('fmod', (1.0, INFINITY), 1.0),
208        ('fmod', (1.0, -INFINITY), 1.0),
209        ('fmod', (INFINITY, INFINITY), ValueError),
210        #
211        ('hypot', (-INFINITY, 1.0), positiveinf),
212        ('hypot', (-2.3, -INFINITY), positiveinf),
213        ]
214
215    binary_math_functions = ['atan2', 'fmod', 'hypot', 'pow']
216
217    NANCASES1 = [
218        (name, (NAN,), isnan) for name in unary_math_functions]
219    NANCASES2 = [
220        (name, (NAN, 0.1), isnan) for name in binary_math_functions]
221    NANCASES3 = [
222        (name, (-0.2, NAN), isnan) for name in binary_math_functions]
223    NANCASES4 = [
224        (name, (NAN, -INFINITY), isnan) for name in binary_math_functions
225                                        if name != 'hypot']
226    NANCASES5 = [
227        (name, (INFINITY, NAN), isnan) for name in binary_math_functions
228                                       if name != 'hypot']
229    NANCASES6 = [
230        ('frexp', (NAN,), lambda x: isnan(x[0])),
231        ('ldexp', (NAN, 2), isnan),
232        ('hypot', (NAN, INFINITY), positiveinf),
233        ('hypot', (NAN, -INFINITY), positiveinf),
234        ('hypot', (INFINITY, NAN), positiveinf),
235        ('hypot', (-INFINITY, NAN), positiveinf),
236        ('modf', (NAN,), lambda x: (isnan(x[0]) and isnan(x[1]))),
237        ]
238
239    # The list of test cases.  Note that various tests import this,
240    # notably in rpython/lltypesystem/module and in translator/c/test.
241    TESTCASES = (REGCASES + IRREGCASES + OVFCASES + INFCASES + IRREGERRCASES
242                 + NANCASES1 + NANCASES2 + NANCASES3 + NANCASES4 + NANCASES5
243                 + NANCASES6)
244
245
246def get_tester(expected):
247    if type(expected) is type(Exception):
248        def tester(value):
249            return False
250    elif callable(expected):
251        def tester(value):
252            ok = expected(value)
253            assert isinstance(ok, bool)
254            return ok
255    else:
256        assert finite(expected), "badly written test"
257        def tester(got):
258            gotsign = expectedsign = 1
259            if got < 0.0: gotsign = -gotsign
260            if expected < 0.0: expectedsign = -expectedsign
261            return finite(got) and (got == expected and
262                                    gotsign == expectedsign)
263    return tester
264

1package org.python.modules;
2
3import org.python.core.Py;
4import org.python.core.PyComplex;
5import org.python.core.PyException;
6import org.python.core.PyFloat;
7import org.python.core.PyInstance;
8import org.python.core.PyObject;
9import org.python.core.PyTuple;
10
11public class cmath {
12
13    public static final PyFloat pi = new PyFloat(Math.PI);
14    public static final PyFloat e = new PyFloat(Math.E);
15
16    /** 2<sup>-½</sup> (Ref: Abramowitz & Stegun [1972], p2). */
17    private static final double ROOT_HALF = 0.70710678118654752440;
18    /** ln({@link Double#MAX_VALUE}) or a little less */
19    private static final double NEARLY_LN_DBL_MAX = 709.4361393;
20    
21    private static final double ATLEAST_27LN2 = 18.72;
22    private static final double HALF_E2 = 0.5 * Math.E * Math.E;
23
24    /** log<sub>10</sub>e (Ref: Abramowitz & Stegun [1972], p3). */
25    private static final double LOG10E = 0.43429448190325182765;
26
27    private static PyComplex complexFromPyObject(PyObject obj) {
28        // If op is already of type PyComplex_Type, return its value
29        if (obj instanceof PyComplex) {
30            return (PyComplex)obj;
31        }
32
33        // If not, use op's __complex__ method, if it exists
34        PyObject newObj = null;
35        if (obj instanceof PyInstance) {
36            // this can go away in python 3000
37            if (obj.__findattr__("__complex__") != null) {
38                newObj = obj.invoke("__complex__");
39            }
40            // else try __float__
41        } else {
42            PyObject complexFunc = obj.getType().lookup("__complex__");
43            if (complexFunc != null) {
44                newObj = complexFunc.__call__(obj);
45            }
46        }
47
48        if (newObj != null) {
49            if (!(newObj instanceof PyComplex)) {
50                throw Py.TypeError("__complex__ should return a complex object");
51            }
52            return (PyComplex)newObj;
53        }
54
55        // If neither of the above works, interpret op as a float giving the real part of
56        // the result, and fill in the imaginary part as 0
57        return new PyComplex(obj.asDouble(), 0);
58    }
59
60    public static PyComplex acos(PyObject w) {
61        return _acos(complexFromPyObject(w));
62    }
63
64    
65    private static PyComplex _acos(PyComplex w) {
66
67        // Let z = x + iy and w = u + iv.
68        double x, y, u = w.real, v = w.imag;
69
70        if (Math.abs(u) > 0x1p27 || Math.abs(v) > 0x1p27) {
71            /*
72             * w is large: approximate 2cos(z) by exp(i(x+iy)) or exp(-i(x+iy)), whichever
73             * dominates. Hence, z = x+iy = i ln(2(u+iv)) or -i ln(2(u+iv))
74             */
75            x = Math.atan2(Math.abs(v), u);
76            y = Math.copySign(logHypot(u, v) + math.LN2, -v);
77
78        } else if (Double.isNaN(v)) {
79            // Special cases
80            x = (u == 0.) ? Math.PI / 2. : v;
81            y = v;
82
83        } else {
84            // Normal case, without risk of overflow.
85            PyComplex a = sqrt(new PyComplex(1. - u, -v)); // a = sqrt(1-w) = sqrt(2) sin(z/2)
86            PyComplex b = sqrt(new PyComplex(1 + u, v));   // b = sqrt(1+w) = sqrt(2) cos(z/2)
87            // Arguments here are sin(x/2)cosh(y/2), cos(x/2)cosh(y/2) giving tan(x/2)
88            x = 2. * Math.atan2(a.real, b.real);
89            // 2 (cos(x/2)**2+sin(x/2)**2) sinh(y/2)cosh(y/2) = sinh y
90            y = math.asinh(a.imag * b.real - a.real * b.imag);
91        }
92
93        // If that generated a nan, and there wasn't one in the argument, raise a domain error.
94        return exceptNaN(new PyComplex(x, y), w);
95    }
96
97    
98    public static PyComplex acosh(PyObject w) {
99        return _acosh(complexFromPyObject(w));
100    }
101
102    
103    private static PyComplex _acosh(PyComplex w) {
104
105        // Let z = x + iy and w = u + iv.
106        double x, y, u = w.real, v = w.imag;
107
108        if (Math.abs(u) > 0x1p27 || Math.abs(v) > 0x1p27) {
109            /*
110             * w is large: approximate 2cosh(z) by exp(x+iy) or exp(-x-iy), whichever dominates.
111             * Hence, z = x+iy = ln(2(u+iv)) or -ln(2(u+iv))
112             */
113            x = logHypot(u, v) + math.LN2;
114            y = Math.atan2(v, u);
115
116        } else if (v == 0. && !Double.isNaN(u)) {
117            /*
118             * We're on the real axis (and maybe the branch cut). u = cosh x cos y. In all cases,
119             * the sign of y follows v.
120             */
121            if (u >= 1.) {
122                // As real library, cos y = 1, u = cosh x.
123                x = math.acosh(u);
124                y = v;
125            } else if (u < -1.) {
126                // Left part of cut: cos y = -1, u = -cosh x
127                x = math.acosh(-u);
128                y = Math.copySign(Math.PI, v);
129            } else {
130                // -1 <= u <= 1: cosh x = 1, u = cos y.
131                x = 0.;
132                y = Math.copySign(Math.acos(u), v);
133            }
134
135        } else {
136            // Normal case, without risk of overflow.
137            PyComplex a = sqrt(new PyComplex(u - 1., v)); // a = sqrt(w-1) = sqrt(2) sinh(z/2)
138            PyComplex b = sqrt(new PyComplex(u + 1., v)); // b = sqrt(w+1) = sqrt(2) cosh(z/2)
139            // 2 sinh(x/2)cosh(x/2) (cos(y/2)**2+sin(y/2)**2) = sinh x
140            x = math.asinh(a.real * b.real + a.imag * b.imag);
141            // Arguments here are cosh(x/2)sin(y/2) and cosh(x/2)cos(y/2) giving tan y/2
142            y = 2. * Math.atan2(a.imag, b.real);
143        }
144
145        // If that generated a nan, and there wasn't one in the argument, raise a domain error.
146        return exceptNaN(new PyComplex(x, y), w);
147    }
148
149    public static PyComplex asin(PyObject w) {
150        return asinOrAsinh(complexFromPyObject(w), false);
151    }
152
153    public static PyComplex asinh(PyObject w) {
154        return asinOrAsinh(complexFromPyObject(w), true);
155    }
156
157    private static PyComplex asinOrAsinh(PyComplex w, boolean h) {
158        double u, v, x, y;
159        PyComplex z;
160
161        if (h) {
162            // We compute z = asinh(w). Let z = x + iy and w = u + iv.
163            u = w.real;
164            v = w.imag;
165            // Then the function body computes x + iy = asinh(w).
166        } else {
167            // We compute w = asin(z). Unusually, let w = u - iv, so u + iv = iw.
168            v = w.real;
169            u = -w.imag;
170            // Then as before, the function body computes asinh(u+iv) = asinh(iw) = i asin(w),
171            // but we finally return z = y - ix = -i asinh(iw) = asin(w).
172        }
173
174        if (Double.isNaN(u)) {
175            // Special case for nan in real part. Default clause deals naturally with v=nan.
176            if (v == 0.) {
177                x = u;
178                y = v;
179            } else if (Double.isInfinite(v)) {
180                x = Double.POSITIVE_INFINITY;
181                y = u;
182            } else { // Any other value of v -> nan+nanj
183                x = y = u;
184            }
185
186        } else if (Math.abs(u) > 0x1p27 || Math.abs(v) > 0x1p27) {
187            /*
188             * w is large: approximate 2sinh(z) by exp(x+iy) or -exp(-x-iy), whichever dominates.
189             * Hence, z = x+iy = ln(2(u+iv)) or -ln(-2(u+iv))
190             */
191            x = logHypot(u, v) + math.LN2;
192            if (Math.copySign(1., u) > 0.) {
193                y = Math.atan2(v, u);
194            } else {
195                // Adjust for sign, choosing the angle so that -pi/2 < y < pi/2
196                x = -x;
197                y = Math.atan2(v, -u);
198            }
199
200        } else {
201            // Normal case, without risk of overflow.
202            PyComplex a = sqrt(new PyComplex(1. + v, -u)); // a = sqrt(1-iw)
203            PyComplex b = sqrt(new PyComplex(1. - v, u));  // b = sqrt(1+iw)
204            // Combine the parts so as that terms in y cancel, leaving us with sinh x:
205            x = math.asinh(a.real * b.imag - a.imag * b.real);
206            // The arguments are v = cosh x sin y, and cosh x cos y
207            y = Math.atan2(v, a.real * b.real - a.imag * b.imag);
208        }
209
210        // Compose the result w according to whether we're computing asin(w) or asinh(w).
211        if (h) {
212            z = new PyComplex(x, y);    // z = x + iy = asinh(u+iv).
213        } else {
214            z = new PyComplex(y, -x);   // z = y - ix = -i asinh(v-iu) = asin(w)
215        }
216
217        // If that generated a nan, and there wasn't one in the argument, raise a domain error.
218        return exceptNaN(z, w);
219    }
220
221
222    public static PyComplex atan(PyObject w) {
223        return atanOrAtanh(complexFromPyObject(w), false);
224    }
225
226    
227    public static PyComplex atanh(PyObject w) {
228        return atanOrAtanh(complexFromPyObject(w), true);
229    }
230
231    private static PyComplex atanOrAtanh(PyComplex w, boolean h) {
232        double u, v, x, y;
233        PyComplex z;
234
235        if (h) {
236            // We compute z = atanh(w). Let z = x + iy and w = u + iv.
237            u = w.real;
238            v = w.imag;
239            // Then the function body computes x + iy = atanh(w).
240        } else {
241            // We compute w = atan(z). Unusually, let w = u - iv, so u + iv = iw.
242            v = w.real;
243            u = -w.imag;
244            // Then as before, the function body computes atanh(u+iv) = atanh(iw) = i atan(w),
245            // but we finally return z = y - ix = -i atanh(iw) = atan(w).
246        }
247
248        double absu = Math.abs(u), absv = Math.abs(v);
249
250        if (absu >= 0x1p511 || absv >= 0x1p511) {
251            // w is large: approximate atanh(w) by 1/w + i pi/2. 1/w = conjg(w)/|w|**2.
252            if (Double.isInfinite(absu) || Double.isInfinite(absv)) {
253                x = Math.copySign(0., u);
254            } else {
255                // w is also too big to square, carry a 2**-N scaling factor.
256                int N = 520;
257                double uu = Math.scalb(u, -N), vv = Math.scalb(v, -N);
258                double mod2w = uu * uu + vv * vv;
259                x = Math.scalb(uu / mod2w, -N);
260            }
261            // We don't need the imaginary part of 1/z. Just pi/2 with the sign of v. (If not nan.)
262            if (Double.isNaN(v)) {
263                y = v;
264            } else {
265                y = Math.copySign(Math.PI / 2., v);
266            }
267
268        } else if (absu < 0x1p-53) {
269            // u is small enough that u**2 may be neglected relative to 1.
270            if (absv > 0x1p-27) {
271                // v is not small, but is not near overflow either.
272                double v2 = v * v;
273                double d = 1. + v2;
274                x = Math.copySign(Math.log1p(4. * absu / d), u) * 0.25;
275                y = Math.atan2(2. * v, 1. - v2) * 0.5;
276            } else {
277                // v is also small enough that v**2 may be neglected (or is nan). So z = w.
278                x = u;
279                y = v;
280            }
281
282        } else if (absu == 1. && absv < 0x1p-27) {
283            // w is close to +1 or -1: needs a different expression, good as v->0
284            x = Math.copySign(Math.log(absv) - math.LN2, u) * 0.5;
285            if (v == 0.) {
286                y = Double.NaN;
287            } else {
288                y = Math.copySign(Math.atan2(2., absv), v) * 0.5;
289            }
290
291        } else {
292            /*
293             * Normal case, without risk of overflow. The basic expression is z =
294             * 0.5*ln((1+w)/(1-w)), which for positive u we rearrange as 0.5*ln(1+2w/(1-w)) and for
295             * negative u as -0.5*ln(1-2w/(1+w)). By use of absu, we reduce the difference between
296             * the expressions fo u>=0 and u<0 to a sign transfer.
297             */
298            double lmu = (1. - absu), lpu = (1. + absu), v2 = v * v;
299            double d = lmu * lmu + v2;
300            x = Math.copySign(Math.log1p(4. * absu / d), u) * 0.25;
301            y = Math.atan2(2. * v, lmu * lpu - v2) * 0.5;
302        }
303
304        // Compose the result w according to whether we're computing atan(w) or atanh(w).
305        if (h) {
306            z = new PyComplex(x, y);    // z = x + iy = atanh(u+iv).
307        } else {
308            z = new PyComplex(y, -x);   // z = y - ix = -i atanh(v-iu) = atan(w)
309        }
310
311        // If that generated a nan, and there wasn't one in the argument, raise a domain error.
312        return exceptNaN(z, w);
313    }

1using System;
2using System.Collections;
3using System.Collections.Generic;
4using System.Diagnostics;
5using System.Linq;
6
7using Microsoft.Scripting;
8using Microsoft.Scripting.Runtime;
9using Microsoft.Scripting.Utils;
10
11using IronPython.Runtime;
12using IronPython.Runtime.Operations;
13using IronPython.Runtime.Types;
14
15using System.Numerics;
16
17[assembly: PythonModule("math", typeof(IronPython.Modules.PythonMath))]
18namespace IronPython.Modules {
19    public static partial class PythonMath {
20        public const string __doc__ = "Provides common mathematical functions.";
21
22        public const double pi = Math.PI;
23        public const double e = Math.E;
24
25        private const double degreesToRadians = Math.PI / 180.0;
26        private const int Bias = 0x3FE;
27
28        public static double degrees(double radians) {
29            return Check(radians, radians / degreesToRadians);
30        }
31
32        public static double radians(double degrees) {
33            return Check(degrees, degrees * degreesToRadians);
34        }
35
36        public static double fmod(double v, double w) {
37            return Check(v, w, v % w);
38        }
39
40        private static double sum(List<double> partials) {
41            // sum the partials the same was as CPython does
42            var n = partials.Count;
43            var hi = 0.0;
44
45            if (n == 0) return hi;
46
47            var lo = 0.0;
48
49            // sum exact
50            while (n > 0) {
51                var x = hi;
52                var y = partials[--n];
53                hi = x + y;
54                lo = y - (hi - x);
55                if (lo != 0.0)
56                    break;
57            }
58
59            if (n == 0) return hi;
60
61            // half-even rounding
62            if (lo < 0.0 && partials[n - 1] < 0.0 || lo > 0.0 && partials[n - 1] > 0.0) {
63                var y = lo * 2.0;
64                var x = hi + y;
65                var yr = x - hi;
66                if (y == yr)
67                    hi = x;
68            }
69            return hi;
70        }
71
72        public static double fsum(IEnumerable e) {
73            // msum from https://code.activestate.com/recipes/393090/
74            var partials = new List<double>();
75            foreach (var v in e.Cast<object>().Select(o => Converter.ConvertToDouble(o))) {
76                var x = v;
77                var i = 0;
78                for (var j = 0; j < partials.Count; j++) {
79                    var y = partials[j];
80                    if (Math.Abs(x) < Math.Abs(y)) {
81                        var t = x;
82                        x = y;
83                        y = t;
84                    }
85                    var hi = x + y;
86                    var lo = y - (hi - x);
87                    if (lo != 0) {
88                        partials[i++] = lo;
89                    }
90                    x = hi;
91                }
92                partials.RemoveRange(i, partials.Count - i);
93                partials.Add(x);
94            }
95
96            return sum(partials);
97        }
98
99        public static PythonTuple frexp(double v) {
100            if (Double.IsInfinity(v) || Double.IsNaN(v)) {
101                return PythonTuple.MakeTuple(v, 0.0);
102            }
103            int exponent = 0;
104            double mantissa = 0;
105
106            if (v == 0) {
107                mantissa = 0;
108                exponent = 0;
109            } else {
110                byte[] vb = BitConverter.GetBytes(v);
111                if (BitConverter.IsLittleEndian) {
112                    DecomposeLe(vb, out mantissa, out exponent);
113                } else {
114                    throw new NotImplementedException();
115                }
116            }
117
118            return PythonTuple.MakeTuple(mantissa, exponent);
119        }
120
121        public static PythonTuple modf(double v) {
122            if (double.IsInfinity(v)) {
123                return PythonTuple.MakeTuple(0.0, v);
124            }
125            double w = v % 1.0;
126            v -= w;
127            return PythonTuple.MakeTuple(w, v);
128        }
129
130        public static double ldexp(double v, BigInteger w) {
131            if (v == 0.0 || double.IsInfinity(v)) {
132                return v;
133            }
134            return Check(v, v * Math.Pow(2.0, (double)w));
135        }
136
137        public static double hypot(double v, double w) {
138            if (double.IsInfinity(v) || double.IsInfinity(w)) {
139                return double.PositiveInfinity;
140            }
141            return Check(v, w, MathUtils.Hypot(v, w));
142        }
143
144        public static double pow(double v, double exp) {
145            if (v == 1.0 || exp == 0.0) {
146                return 1.0;
147            } else if (double.IsNaN(v) || double.IsNaN(exp)) {
148                return double.NaN;
149            } else if (v == 0.0) {
150                if (exp > 0.0) {
151                    return 0.0;
152                }
153                throw PythonOps.ValueError("math domain error");
154            } else if (double.IsPositiveInfinity(exp)) {
155                if (v > 1.0 || v < -1.0) {
156                    return double.PositiveInfinity;
157                } else if (v == -1.0) {
158                    return 1.0;
159                } else {
160                    return 0.0;
161                }
162            } else if (double.IsNegativeInfinity(exp)) {
163                if (v > 1.0 || v < -1.0) {
164                    return 0.0;
165                } else if (v == -1.0) {
166                    return 1.0;
167                } else {
168                    return double.PositiveInfinity;
169                }
170            }
171            return Check(v, exp, Math.Pow(v, exp));
172        }
173
174        public static double log(double v0) {
175            if (v0 <= 0.0) {
176                throw PythonOps.ValueError("math domain error");
177            }
178            return Check(v0, Math.Log(v0));
179        }
180
181        public static double log(double v0, double v1) {
182            if (v0 <= 0.0 || v1 == 0.0) {
183                throw PythonOps.ValueError("math domain error");
184            } else if (v1 == 1.0) {
185                throw PythonOps.ZeroDivisionError("float division");
186            } else if (v1 == Double.PositiveInfinity) {
187                return 0.0;
188            }
189            return Check(Math.Log(v0, v1));
190        }
191
192        public static double log(BigInteger value) {
193            if (value.Sign <= 0) {
194                throw PythonOps.ValueError("math domain error");
195            }
196            return value.Log();
197        }
198
199        public static double log(object value) {
200            // CPython tries float first, then double, so we need
201            // an explicit overload which properly matches the order here
202            double val;
203            if (Converter.TryConvertToDouble(value, out val)) {
204                return log(val);
205            } else {
206                return log(Converter.ConvertToBigInteger(value));
207            }
208        }
209
210        public static double log(BigInteger value, double newBase) {
211            if (newBase <= 0.0 || value <= 0) {
212                throw PythonOps.ValueError("math domain error");
213            } else if (newBase == 1.0) {
214                throw PythonOps.ZeroDivisionError("float division");
215            } else if (newBase == Double.PositiveInfinity) {
216                return 0.0;
217            }
218            return Check(value.Log(newBase));
219        }
220
221        public static double log(object value, double newBase) {
222            // CPython tries float first, then double, so we need
223            // an explicit overload which properly matches the order here
224            double val;
225            if (Converter.TryConvertToDouble(value, out val)) {
226                return log(val, newBase);
227            } else {
228                return log(Converter.ConvertToBigInteger(value), newBase);
229            }
230        }
231
232        public static double log10(double v0) {
233            if (v0 <= 0.0) {
234                throw PythonOps.ValueError("math domain error");
235            }
236            return Check(v0, Math.Log10(v0));
237        }