blob: 4e08722bccac54f83874f6376f5ca502f6038bdc [file] [log] [blame]
Guido van Rossum71aa32f1996-01-12 01:34:57 +00001/* Complex math module */
2
3/* much code borrowed from mathmodule.c */
4
Roger E. Masse24070ca1996-12-09 22:59:53 +00005#include "Python.h"
Guido van Rossum71aa32f1996-01-12 01:34:57 +00006
Guido van Rossum71aa32f1996-01-12 01:34:57 +00007#include "mymath.h"
8
9#ifdef i860
10/* Cray APP has bogus definition of HUGE_VAL in <math.h> */
11#undef HUGE_VAL
12#endif
13
14#ifdef HUGE_VAL
15#define CHECK(x) if (errno != 0) ; \
16 else if (-HUGE_VAL <= (x) && (x) <= HUGE_VAL) ; \
17 else errno = ERANGE
18#else
19#define CHECK(x) /* Don't know how to check */
20#endif
21
22#ifndef M_PI
23#define M_PI (3.141592653589793239)
24#endif
25
26/* First, the C functions that do the real work */
27
28/* constants */
Guido van Rossum9e720e31996-07-21 02:31:35 +000029static Py_complex c_1 = {1., 0.};
30static Py_complex c_half = {0.5, 0.};
31static Py_complex c_i = {0., 1.};
32static Py_complex c_i2 = {0., 0.5};
Guido van Rossuma376cc51996-12-05 23:43:35 +000033#if 0
Guido van Rossum9e720e31996-07-21 02:31:35 +000034static Py_complex c_mi = {0., -1.};
35static Py_complex c_pi2 = {M_PI/2., 0.};
Guido van Rossuma376cc51996-12-05 23:43:35 +000036#endif
Guido van Rossum71aa32f1996-01-12 01:34:57 +000037
38/* forward declarations */
Guido van Rossum9e720e31996-07-21 02:31:35 +000039staticforward Py_complex c_log();
40staticforward Py_complex c_prodi();
41staticforward Py_complex c_sqrt();
Guido van Rossum71aa32f1996-01-12 01:34:57 +000042
43
Guido van Rossum9e720e31996-07-21 02:31:35 +000044static Py_complex c_acos(x)
45 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +000046{
47 return c_neg(c_prodi(c_log(c_sum(x,c_prod(c_i,
48 c_sqrt(c_diff(c_1,c_prod(x,x))))))));
49}
50
Guido van Rossum9e720e31996-07-21 02:31:35 +000051static Py_complex c_acosh(x)
52 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +000053{
54 return c_log(c_sum(x,c_prod(c_i,
55 c_sqrt(c_diff(c_1,c_prod(x,x))))));
56}
57
Guido van Rossum9e720e31996-07-21 02:31:35 +000058static Py_complex c_asin(x)
59 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +000060{
61 return c_neg(c_prodi(c_log(c_sum(c_prod(c_i,x),
62 c_sqrt(c_diff(c_1,c_prod(x,x)))))));
63}
64
Guido van Rossum9e720e31996-07-21 02:31:35 +000065static Py_complex c_asinh(x)
66 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +000067{
68 return c_neg(c_log(c_diff(c_sqrt(c_sum(c_1,c_prod(x,x))),x)));
69}
70
Guido van Rossum9e720e31996-07-21 02:31:35 +000071static Py_complex c_atan(x)
72 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +000073{
74 return c_prod(c_i2,c_log(c_quot(c_sum(c_i,x),c_diff(c_i,x))));
75}
76
Guido van Rossum9e720e31996-07-21 02:31:35 +000077static Py_complex c_atanh(x)
78 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +000079{
80 return c_prod(c_half,c_log(c_quot(c_sum(c_1,x),c_diff(c_1,x))));
81}
82
Guido van Rossum9e720e31996-07-21 02:31:35 +000083static Py_complex c_cos(x)
84 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +000085{
Guido van Rossum9e720e31996-07-21 02:31:35 +000086 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +000087 r.real = cos(x.real)*cosh(x.imag);
88 r.imag = -sin(x.real)*sinh(x.imag);
89 return r;
90}
91
Guido van Rossum9e720e31996-07-21 02:31:35 +000092static Py_complex c_cosh(x)
93 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +000094{
Guido van Rossum9e720e31996-07-21 02:31:35 +000095 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +000096 r.real = cos(x.imag)*cosh(x.real);
97 r.imag = sin(x.imag)*sinh(x.real);
98 return r;
99}
100
Guido van Rossum9e720e31996-07-21 02:31:35 +0000101static Py_complex c_exp(x)
102 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000103{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000104 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000105 double l = exp(x.real);
106 r.real = l*cos(x.imag);
107 r.imag = l*sin(x.imag);
108 return r;
109}
110
Guido van Rossum9e720e31996-07-21 02:31:35 +0000111static Py_complex c_log(x)
112 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000113{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000114 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000115 double l = hypot(x.real,x.imag);
116 r.imag = atan2(x.imag, x.real);
117 r.real = log(l);
118 return r;
119}
120
Guido van Rossum9e720e31996-07-21 02:31:35 +0000121static Py_complex c_log10(x)
122 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000123{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000124 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000125 double l = hypot(x.real,x.imag);
126 r.imag = atan2(x.imag, x.real)/log(10.);
127 r.real = log10(l);
128 return r;
129}
130
Guido van Rossum9e720e31996-07-21 02:31:35 +0000131static Py_complex c_prodi(x)
132 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000133{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000134 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000135 r.real = -x.imag;
136 r.imag = x.real;
137 return r;
138}
139
Guido van Rossum9e720e31996-07-21 02:31:35 +0000140static Py_complex c_sin(x)
141 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000142{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000143 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000144 r.real = sin(x.real)*cosh(x.imag);
145 r.imag = cos(x.real)*sinh(x.imag);
146 return r;
147}
148
Guido van Rossum9e720e31996-07-21 02:31:35 +0000149static Py_complex c_sinh(x)
150 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000151{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000152 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000153 r.real = cos(x.imag)*sinh(x.real);
154 r.imag = sin(x.imag)*cosh(x.real);
155 return r;
156}
157
Guido van Rossum9e720e31996-07-21 02:31:35 +0000158static Py_complex c_sqrt(x)
159 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000160{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000161 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000162 double s,d;
163 if (x.real == 0. && x.imag == 0.)
164 r = x;
165 else {
166 s = sqrt(0.5*(fabs(x.real) + hypot(x.real,x.imag)));
167 d = 0.5*x.imag/s;
168 if (x.real > 0.) {
169 r.real = s;
170 r.imag = d;
171 }
172 else if (x.imag >= 0.) {
173 r.real = d;
174 r.imag = s;
175 }
176 else {
177 r.real = -d;
178 r.imag = -s;
179 }
180 }
181 return r;
182}
183
Guido van Rossum9e720e31996-07-21 02:31:35 +0000184static Py_complex c_tan(x)
185 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000186{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000187 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000188 double sr,cr,shi,chi;
189 double rs,is,rc,ic;
190 double d;
191 sr = sin(x.real);
192 cr = cos(x.real);
193 shi = sinh(x.imag);
194 chi = cosh(x.imag);
195 rs = sr*chi;
196 is = cr*shi;
197 rc = cr*chi;
198 ic = -sr*shi;
199 d = rc*rc + ic*ic;
200 r.real = (rs*rc+is*ic)/d;
201 r.imag = (is*rc-rs*ic)/d;
202 return r;
203}
204
Guido van Rossum9e720e31996-07-21 02:31:35 +0000205static Py_complex c_tanh(x)
206 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000207{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000208 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000209 double si,ci,shr,chr;
210 double rs,is,rc,ic;
211 double d;
212 si = sin(x.imag);
213 ci = cos(x.imag);
214 shr = sinh(x.real);
215 chr = cosh(x.real);
216 rs = ci*shr;
217 is = si*chr;
218 rc = ci*chr;
219 ic = si*shr;
220 d = rc*rc + ic*ic;
221 r.real = (rs*rc+is*ic)/d;
222 r.imag = (is*rc-rs*ic)/d;
223 return r;
224}
225
226
227/* And now the glue to make them available from Python: */
228
Roger E. Masse24070ca1996-12-09 22:59:53 +0000229static PyObject *
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000230math_error()
231{
232 if (errno == EDOM)
Roger E. Masse24070ca1996-12-09 22:59:53 +0000233 PyErr_SetString(PyExc_ValueError, "math domain error");
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000234 else if (errno == ERANGE)
Roger E. Masse24070ca1996-12-09 22:59:53 +0000235 PyErr_SetString(PyExc_OverflowError, "math range error");
236 else /* Unexpected math error */
237 PyErr_SetFromErrno(PyExc_ValueError);
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000238 return NULL;
239}
240
Roger E. Masse24070ca1996-12-09 22:59:53 +0000241static PyObject *
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000242math_1(args, func)
Roger E. Masse24070ca1996-12-09 22:59:53 +0000243 PyObject *args;
244 Py_complex (*func) Py_FPROTO((Py_complex));
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000245{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000246 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000247 if (!PyArg_ParseTuple(args, "D", &x))
248 return NULL;
249 errno = 0;
Guido van Rossum52fa3a61997-02-14 22:59:58 +0000250 PyFPE_START_PROTECT("complex function", return 0)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000251 x = (*func)(x);
Guido van Rossum45b83911997-03-14 04:32:50 +0000252 PyFPE_END_PROTECT(x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000253 CHECK(x.real);
254 CHECK(x.imag);
255 if (errno != 0)
256 return math_error();
257 else
Roger E. Masse24070ca1996-12-09 22:59:53 +0000258 return PyComplex_FromCComplex(x);
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000259}
260
261#define FUNC1(stubname, func) \
Roger E. Masse24070ca1996-12-09 22:59:53 +0000262 static PyObject * stubname(self, args) PyObject *self, *args; { \
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000263 return math_1(args, func); \
264 }
265
266FUNC1(cmath_acos, c_acos)
267FUNC1(cmath_acosh, c_acosh)
268FUNC1(cmath_asin, c_asin)
269FUNC1(cmath_asinh, c_asinh)
270FUNC1(cmath_atan, c_atan)
271FUNC1(cmath_atanh, c_atanh)
272FUNC1(cmath_cos, c_cos)
273FUNC1(cmath_cosh, c_cosh)
274FUNC1(cmath_exp, c_exp)
275FUNC1(cmath_log, c_log)
276FUNC1(cmath_log10, c_log10)
277FUNC1(cmath_sin, c_sin)
278FUNC1(cmath_sinh, c_sinh)
279FUNC1(cmath_sqrt, c_sqrt)
280FUNC1(cmath_tan, c_tan)
281FUNC1(cmath_tanh, c_tanh)
282
283
Roger E. Masse24070ca1996-12-09 22:59:53 +0000284static PyMethodDef cmath_methods[] = {
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000285 {"acos", cmath_acos, 1},
286 {"acosh", cmath_acosh, 1},
287 {"asin", cmath_asin, 1},
288 {"asinh", cmath_asinh, 1},
289 {"atan", cmath_atan, 1},
290 {"atanh", cmath_atanh, 1},
291 {"cos", cmath_cos, 1},
292 {"cosh", cmath_cosh, 1},
293 {"exp", cmath_exp, 1},
294 {"log", cmath_log, 1},
295 {"log10", cmath_log10, 1},
296 {"sin", cmath_sin, 1},
297 {"sinh", cmath_sinh, 1},
298 {"sqrt", cmath_sqrt, 1},
299 {"tan", cmath_tan, 1},
300 {"tanh", cmath_tanh, 1},
301 {NULL, NULL} /* sentinel */
302};
303
304void
305initcmath()
306{
Roger E. Masse24070ca1996-12-09 22:59:53 +0000307 PyObject *m, *d, *v;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000308
309 m = Py_InitModule("cmath", cmath_methods);
Roger E. Masse24070ca1996-12-09 22:59:53 +0000310 d = PyModule_GetDict(m);
311 PyDict_SetItemString(d, "pi",
312 v = PyFloat_FromDouble(atan(1.0) * 4.0));
313 Py_DECREF(v);
314 PyDict_SetItemString(d, "e", v = PyFloat_FromDouble(exp(1.0)));
315 Py_DECREF(v);
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000316}