blob: adf76b8a331dd5e3b228c283d91dd3ad3e1d97e3 [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#ifndef M_PI
8#define M_PI (3.141592653589793239)
9#endif
10
11/* First, the C functions that do the real work */
12
13/* constants */
Tim Peters8cb82bd2001-02-20 20:36:38 +000014static Py_complex c_one = {1., 0.};
Guido van Rossum9e720e31996-07-21 02:31:35 +000015static Py_complex c_half = {0.5, 0.};
16static Py_complex c_i = {0., 1.};
Tim Peters8cb82bd2001-02-20 20:36:38 +000017static Py_complex c_halfi = {0., 0.5};
Guido van Rossum71aa32f1996-01-12 01:34:57 +000018
19/* forward declarations */
Thomas Wouterseb61b6e2000-07-24 11:17:40 +000020staticforward Py_complex c_log(Py_complex);
21staticforward Py_complex c_prodi(Py_complex);
22staticforward Py_complex c_sqrt(Py_complex);
Guido van Rossum71aa32f1996-01-12 01:34:57 +000023
24
Tim Peters14e26402001-02-20 20:15:19 +000025static Py_complex
26c_acos(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +000027{
28 return c_neg(c_prodi(c_log(c_sum(x,c_prod(c_i,
Tim Peters8cb82bd2001-02-20 20:36:38 +000029 c_sqrt(c_diff(c_one,c_prod(x,x))))))));
Guido van Rossum71aa32f1996-01-12 01:34:57 +000030}
31
Tim Peters14e26402001-02-20 20:15:19 +000032static char c_acos_doc[] =
33"acos(x)\n"
34"\n"
35"Return the arc cosine of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +000036
37
Tim Peters14e26402001-02-20 20:15:19 +000038static Py_complex
39c_acosh(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +000040{
Guido van Rossumf385c5e2000-06-30 02:29:22 +000041 Py_complex z;
42 z = c_sqrt(c_half);
Tim Peters8cb82bd2001-02-20 20:36:38 +000043 z = c_log(c_prod(z, c_sum(c_sqrt(c_sum(x,c_one)),
44 c_sqrt(c_diff(x,c_one)))));
Guido van Rossumf385c5e2000-06-30 02:29:22 +000045 return c_sum(z, z);
Guido van Rossum71aa32f1996-01-12 01:34:57 +000046}
47
Tim Peters14e26402001-02-20 20:15:19 +000048static char c_acosh_doc[] =
49"acosh(x)\n"
50"\n"
51"Return the hyperbolic arccosine of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +000052
53
Tim Peters14e26402001-02-20 20:15:19 +000054static Py_complex
55c_asin(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +000056{
Tim Petersedf22102001-02-21 03:22:39 +000057 /* -i * log[(sqrt(1-x**2) + i*x] */
58 const Py_complex squared = c_prod(x, x);
59 const Py_complex sqrt_1_minus_x_sq = c_sqrt(c_diff(c_one, squared));
Tim Peters1c8e1f02001-02-22 19:51:56 +000060 return c_neg(c_prodi(c_log(
61 c_sum(sqrt_1_minus_x_sq, c_prodi(x))
62 ) ) );
Guido van Rossum71aa32f1996-01-12 01:34:57 +000063}
64
Tim Peters14e26402001-02-20 20:15:19 +000065static char c_asin_doc[] =
66"asin(x)\n"
67"\n"
68"Return the arc sine of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +000069
70
Tim Peters14e26402001-02-20 20:15:19 +000071static Py_complex
72c_asinh(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +000073{
Guido van Rossum11a50711999-01-14 19:11:11 +000074 Py_complex z;
Tim Petersedf22102001-02-21 03:22:39 +000075 z = c_sqrt(c_half);
76 z = c_log(c_prod(z, c_sum(c_sqrt(c_sum(x, c_i)),
77 c_sqrt(c_diff(x, c_i)))));
78 return c_sum(z, z);
Guido van Rossum71aa32f1996-01-12 01:34:57 +000079}
80
Tim Peters14e26402001-02-20 20:15:19 +000081static char c_asinh_doc[] =
82"asinh(x)\n"
83"\n"
84"Return the hyperbolic arc sine of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +000085
86
Tim Peters14e26402001-02-20 20:15:19 +000087static Py_complex
88c_atan(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +000089{
Tim Peters8cb82bd2001-02-20 20:36:38 +000090 return c_prod(c_halfi,c_log(c_quot(c_sum(c_i,x),c_diff(c_i,x))));
Guido van Rossum71aa32f1996-01-12 01:34:57 +000091}
92
Tim Peters14e26402001-02-20 20:15:19 +000093static char c_atan_doc[] =
94"atan(x)\n"
95"\n"
96"Return the arc tangent of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +000097
98
Tim Peters14e26402001-02-20 20:15:19 +000099static Py_complex
100c_atanh(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000101{
Tim Peters8cb82bd2001-02-20 20:36:38 +0000102 return c_prod(c_half,c_log(c_quot(c_sum(c_one,x),c_diff(c_one,x))));
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000103}
104
Tim Peters14e26402001-02-20 20:15:19 +0000105static char c_atanh_doc[] =
106"atanh(x)\n"
107"\n"
108"Return the hyperbolic arc tangent of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000109
110
Tim Peters14e26402001-02-20 20:15:19 +0000111static Py_complex
112c_cos(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 r.real = cos(x.real)*cosh(x.imag);
116 r.imag = -sin(x.real)*sinh(x.imag);
117 return r;
118}
119
Tim Peters14e26402001-02-20 20:15:19 +0000120static char c_cos_doc[] =
121"cos(x)\n"
122"n"
123"Return the cosine of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000124
125
Tim Peters14e26402001-02-20 20:15:19 +0000126static Py_complex
127c_cosh(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000128{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000129 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000130 r.real = cos(x.imag)*cosh(x.real);
131 r.imag = sin(x.imag)*sinh(x.real);
132 return r;
133}
134
Tim Peters14e26402001-02-20 20:15:19 +0000135static char c_cosh_doc[] =
136"cosh(x)\n"
137"n"
138"Return the hyperbolic cosine of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000139
140
Tim Peters14e26402001-02-20 20:15:19 +0000141static Py_complex
142c_exp(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000143{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000144 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000145 double l = exp(x.real);
146 r.real = l*cos(x.imag);
147 r.imag = l*sin(x.imag);
148 return r;
149}
150
Tim Peters14e26402001-02-20 20:15:19 +0000151static char c_exp_doc[] =
152"exp(x)\n"
153"\n"
154"Return the exponential value e**x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000155
156
Tim Peters14e26402001-02-20 20:15:19 +0000157static Py_complex
158c_log(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000159{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000160 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000161 double l = hypot(x.real,x.imag);
Martin v. Löwis387c5472001-09-06 08:16:17 +0000162 r.imag = atan2(x.imag, x.real);
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000163 r.real = log(l);
164 return r;
165}
166
Tim Peters14e26402001-02-20 20:15:19 +0000167static char c_log_doc[] =
168"log(x)\n"
169"\n"
170"Return the natural logarithm of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000171
172
Tim Peters14e26402001-02-20 20:15:19 +0000173static Py_complex
174c_log10(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000175{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000176 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000177 double l = hypot(x.real,x.imag);
Martin v. Löwis387c5472001-09-06 08:16:17 +0000178 r.imag = atan2(x.imag, x.real)/log(10.);
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000179 r.real = log10(l);
180 return r;
181}
182
Tim Peters14e26402001-02-20 20:15:19 +0000183static char c_log10_doc[] =
184"log10(x)\n"
185"\n"
186"Return the base-10 logarithm of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000187
188
189/* internal function not available from Python */
Tim Peters14e26402001-02-20 20:15:19 +0000190static Py_complex
191c_prodi(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000192{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000193 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000194 r.real = -x.imag;
195 r.imag = x.real;
196 return r;
197}
198
Guido van Rossumc6e22901998-12-04 19:26:43 +0000199
Tim Peters14e26402001-02-20 20:15:19 +0000200static Py_complex
201c_sin(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000202{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000203 Py_complex r;
Tim Peters14e26402001-02-20 20:15:19 +0000204 r.real = sin(x.real) * cosh(x.imag);
205 r.imag = cos(x.real) * sinh(x.imag);
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000206 return r;
207}
208
Tim Peters14e26402001-02-20 20:15:19 +0000209static char c_sin_doc[] =
210"sin(x)\n"
211"\n"
212"Return the sine of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000213
214
Tim Peters14e26402001-02-20 20:15:19 +0000215static Py_complex
216c_sinh(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000217{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000218 Py_complex r;
Tim Peters14e26402001-02-20 20:15:19 +0000219 r.real = cos(x.imag) * sinh(x.real);
220 r.imag = sin(x.imag) * cosh(x.real);
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000221 return r;
222}
223
Tim Peters14e26402001-02-20 20:15:19 +0000224static char c_sinh_doc[] =
225"sinh(x)\n"
226"\n"
227"Return the hyperbolic sine of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000228
229
Tim Peters14e26402001-02-20 20:15:19 +0000230static Py_complex
231c_sqrt(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000232{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000233 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000234 double s,d;
235 if (x.real == 0. && x.imag == 0.)
236 r = x;
237 else {
238 s = sqrt(0.5*(fabs(x.real) + hypot(x.real,x.imag)));
239 d = 0.5*x.imag/s;
240 if (x.real > 0.) {
241 r.real = s;
242 r.imag = d;
243 }
244 else if (x.imag >= 0.) {
245 r.real = d;
246 r.imag = s;
247 }
248 else {
249 r.real = -d;
250 r.imag = -s;
251 }
252 }
253 return r;
254}
255
Tim Peters14e26402001-02-20 20:15:19 +0000256static char c_sqrt_doc[] =
257"sqrt(x)\n"
258"\n"
259"Return the square root of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000260
261
Tim Peters14e26402001-02-20 20:15:19 +0000262static Py_complex
263c_tan(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000264{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000265 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000266 double sr,cr,shi,chi;
267 double rs,is,rc,ic;
268 double d;
269 sr = sin(x.real);
270 cr = cos(x.real);
271 shi = sinh(x.imag);
272 chi = cosh(x.imag);
Tim Peters14e26402001-02-20 20:15:19 +0000273 rs = sr * chi;
274 is = cr * shi;
275 rc = cr * chi;
276 ic = -sr * shi;
277 d = rc*rc + ic * ic;
278 r.real = (rs*rc + is*ic) / d;
279 r.imag = (is*rc - rs*ic) / d;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000280 return r;
281}
282
Tim Peters14e26402001-02-20 20:15:19 +0000283static char c_tan_doc[] =
284"tan(x)\n"
285"\n"
286"Return the tangent of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000287
288
Tim Peters14e26402001-02-20 20:15:19 +0000289static Py_complex
290c_tanh(Py_complex x)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000291{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000292 Py_complex r;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000293 double si,ci,shr,chr;
294 double rs,is,rc,ic;
295 double d;
296 si = sin(x.imag);
297 ci = cos(x.imag);
298 shr = sinh(x.real);
299 chr = cosh(x.real);
Tim Peters14e26402001-02-20 20:15:19 +0000300 rs = ci * shr;
301 is = si * chr;
302 rc = ci * chr;
303 ic = si * shr;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000304 d = rc*rc + ic*ic;
Tim Peters14e26402001-02-20 20:15:19 +0000305 r.real = (rs*rc + is*ic) / d;
306 r.imag = (is*rc - rs*ic) / d;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000307 return r;
308}
309
Tim Peters14e26402001-02-20 20:15:19 +0000310static char c_tanh_doc[] =
311"tanh(x)\n"
312"\n"
313"Return the hyperbolic tangent of x.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000314
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000315
316/* And now the glue to make them available from Python: */
317
Roger E. Masse24070ca1996-12-09 22:59:53 +0000318static PyObject *
Thomas Woutersf3f33dc2000-07-21 06:00:07 +0000319math_error(void)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000320{
321 if (errno == EDOM)
Roger E. Masse24070ca1996-12-09 22:59:53 +0000322 PyErr_SetString(PyExc_ValueError, "math domain error");
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000323 else if (errno == ERANGE)
Roger E. Masse24070ca1996-12-09 22:59:53 +0000324 PyErr_SetString(PyExc_OverflowError, "math range error");
325 else /* Unexpected math error */
Tim Peters14e26402001-02-20 20:15:19 +0000326 PyErr_SetFromErrno(PyExc_ValueError);
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000327 return NULL;
328}
329
Roger E. Masse24070ca1996-12-09 22:59:53 +0000330static PyObject *
Peter Schneider-Kampf1ca8982000-07-10 09:31:34 +0000331math_1(PyObject *args, Py_complex (*func)(Py_complex))
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000332{
Guido van Rossum9e720e31996-07-21 02:31:35 +0000333 Py_complex x;
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000334 if (!PyArg_ParseTuple(args, "D", &x))
335 return NULL;
336 errno = 0;
Guido van Rossum52fa3a61997-02-14 22:59:58 +0000337 PyFPE_START_PROTECT("complex function", return 0)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000338 x = (*func)(x);
Guido van Rossum45b83911997-03-14 04:32:50 +0000339 PyFPE_END_PROTECT(x)
Tim Petersa40c7932001-09-05 22:36:56 +0000340 Py_SET_ERANGE_IF_OVERFLOW(x.real);
341 Py_SET_ERANGE_IF_OVERFLOW(x.imag);
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000342 if (errno != 0)
343 return math_error();
344 else
Roger E. Masse24070ca1996-12-09 22:59:53 +0000345 return PyComplex_FromCComplex(x);
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000346}
347
348#define FUNC1(stubname, func) \
Peter Schneider-Kampf1ca8982000-07-10 09:31:34 +0000349 static PyObject * stubname(PyObject *self, PyObject *args) { \
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000350 return math_1(args, func); \
351 }
352
353FUNC1(cmath_acos, c_acos)
354FUNC1(cmath_acosh, c_acosh)
355FUNC1(cmath_asin, c_asin)
356FUNC1(cmath_asinh, c_asinh)
357FUNC1(cmath_atan, c_atan)
358FUNC1(cmath_atanh, c_atanh)
359FUNC1(cmath_cos, c_cos)
360FUNC1(cmath_cosh, c_cosh)
361FUNC1(cmath_exp, c_exp)
362FUNC1(cmath_log, c_log)
363FUNC1(cmath_log10, c_log10)
364FUNC1(cmath_sin, c_sin)
365FUNC1(cmath_sinh, c_sinh)
366FUNC1(cmath_sqrt, c_sqrt)
367FUNC1(cmath_tan, c_tan)
368FUNC1(cmath_tanh, c_tanh)
369
370
Tim Peters14e26402001-02-20 20:15:19 +0000371static char module_doc[] =
372"This module is always available. It provides access to mathematical\n"
373"functions for complex numbers.";
Guido van Rossumc6e22901998-12-04 19:26:43 +0000374
Roger E. Masse24070ca1996-12-09 22:59:53 +0000375static PyMethodDef cmath_methods[] = {
Tim Peters14e26402001-02-20 20:15:19 +0000376 {"acos", cmath_acos, METH_VARARGS, c_acos_doc},
377 {"acosh", cmath_acosh, METH_VARARGS, c_acosh_doc},
378 {"asin", cmath_asin, METH_VARARGS, c_asin_doc},
379 {"asinh", cmath_asinh, METH_VARARGS, c_asinh_doc},
380 {"atan", cmath_atan, METH_VARARGS, c_atan_doc},
381 {"atanh", cmath_atanh, METH_VARARGS, c_atanh_doc},
382 {"cos", cmath_cos, METH_VARARGS, c_cos_doc},
383 {"cosh", cmath_cosh, METH_VARARGS, c_cosh_doc},
384 {"exp", cmath_exp, METH_VARARGS, c_exp_doc},
385 {"log", cmath_log, METH_VARARGS, c_log_doc},
386 {"log10", cmath_log10, METH_VARARGS, c_log10_doc},
387 {"sin", cmath_sin, METH_VARARGS, c_sin_doc},
388 {"sinh", cmath_sinh, METH_VARARGS, c_sinh_doc},
389 {"sqrt", cmath_sqrt, METH_VARARGS, c_sqrt_doc},
390 {"tan", cmath_tan, METH_VARARGS, c_tan_doc},
391 {"tanh", cmath_tanh, METH_VARARGS, c_tanh_doc},
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000392 {NULL, NULL} /* sentinel */
393};
394
Guido van Rossum3886bb61998-12-04 18:50:17 +0000395DL_EXPORT(void)
Thomas Woutersf3f33dc2000-07-21 06:00:07 +0000396initcmath(void)
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000397{
Roger E. Masse24070ca1996-12-09 22:59:53 +0000398 PyObject *m, *d, *v;
Tim Peters14e26402001-02-20 20:15:19 +0000399
Guido van Rossumc6e22901998-12-04 19:26:43 +0000400 m = Py_InitModule3("cmath", cmath_methods, module_doc);
Roger E. Masse24070ca1996-12-09 22:59:53 +0000401 d = PyModule_GetDict(m);
402 PyDict_SetItemString(d, "pi",
403 v = PyFloat_FromDouble(atan(1.0) * 4.0));
404 Py_DECREF(v);
405 PyDict_SetItemString(d, "e", v = PyFloat_FromDouble(exp(1.0)));
406 Py_DECREF(v);
Guido van Rossum71aa32f1996-01-12 01:34:57 +0000407}