cmathmodule.c 7.84 KB
Newer Older
Guido van Rossum's avatar
Guido van Rossum committed
1 2 3 4
/* Complex math module */

/* much code borrowed from mathmodule.c */

Roger E. Masse's avatar
Roger E. Masse committed
5
#include "Python.h"
Guido van Rossum's avatar
Guido van Rossum committed
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

#ifdef i860
/* Cray APP has bogus definition of HUGE_VAL in <math.h> */
#undef HUGE_VAL
#endif

#ifdef HUGE_VAL
#define CHECK(x) if (errno != 0) ; \
	else if (-HUGE_VAL <= (x) && (x) <= HUGE_VAL) ; \
	else errno = ERANGE
#else
#define CHECK(x) /* Don't know how to check */
#endif

#ifndef M_PI
#define M_PI (3.141592653589793239)
#endif

/* First, the C functions that do the real work */

/* constants */
27
static Py_complex c_one = {1., 0.};
Guido van Rossum's avatar
Guido van Rossum committed
28 29
static Py_complex c_half = {0.5, 0.};
static Py_complex c_i = {0., 1.};
30
static Py_complex c_halfi = {0., 0.5};
Guido van Rossum's avatar
Guido van Rossum committed
31 32

/* forward declarations */
33 34 35
staticforward Py_complex c_log(Py_complex);
staticforward Py_complex c_prodi(Py_complex);
staticforward Py_complex c_sqrt(Py_complex);
Guido van Rossum's avatar
Guido van Rossum committed
36 37


38 39
static Py_complex
c_acos(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
40 41
{
	return c_neg(c_prodi(c_log(c_sum(x,c_prod(c_i,
42
		    c_sqrt(c_diff(c_one,c_prod(x,x))))))));
Guido van Rossum's avatar
Guido van Rossum committed
43 44
}

45 46 47 48
static char c_acos_doc[] =
"acos(x)\n"
"\n"
"Return the arc cosine of x.";
49 50


51 52
static Py_complex
c_acosh(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
53
{
54 55
	Py_complex z;
	z = c_sqrt(c_half);
56 57
	z = c_log(c_prod(z, c_sum(c_sqrt(c_sum(x,c_one)),
				  c_sqrt(c_diff(x,c_one)))));
58
	return c_sum(z, z);
Guido van Rossum's avatar
Guido van Rossum committed
59 60
}

61 62 63 64
static char c_acosh_doc[] =
"acosh(x)\n"
"\n"
"Return the hyperbolic arccosine of x.";
65 66


67 68
static Py_complex
c_asin(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
69
{
70 71 72
	/* -i * log[(sqrt(1-x**2) + i*x] */
	const Py_complex squared = c_prod(x, x);
	const Py_complex sqrt_1_minus_x_sq = c_sqrt(c_diff(c_one, squared));
73 74 75
        return c_neg(c_prodi(c_log(
        		c_sum(sqrt_1_minus_x_sq, c_prodi(x))
		    )       )     );
Guido van Rossum's avatar
Guido van Rossum committed
76 77
}

78 79 80 81
static char c_asin_doc[] =
"asin(x)\n"
"\n"
"Return the arc sine of x.";
82 83


84 85
static Py_complex
c_asinh(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
86
{
87
	Py_complex z;
88 89 90 91
	z = c_sqrt(c_half);
	z = c_log(c_prod(z, c_sum(c_sqrt(c_sum(x, c_i)),
				  c_sqrt(c_diff(x, c_i)))));
	return c_sum(z, z);
Guido van Rossum's avatar
Guido van Rossum committed
92 93
}

94 95 96 97
static char c_asinh_doc[] =
"asinh(x)\n"
"\n"
"Return the hyperbolic arc sine of x.";
98 99


100 101
static Py_complex
c_atan(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
102
{
103
	return c_prod(c_halfi,c_log(c_quot(c_sum(c_i,x),c_diff(c_i,x))));
Guido van Rossum's avatar
Guido van Rossum committed
104 105
}

106 107 108 109
static char c_atan_doc[] =
"atan(x)\n"
"\n"
"Return the arc tangent of x.";
110 111


112 113
static Py_complex
c_atanh(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
114
{
115
	return c_prod(c_half,c_log(c_quot(c_sum(c_one,x),c_diff(c_one,x))));
Guido van Rossum's avatar
Guido van Rossum committed
116 117
}

118 119 120 121
static char c_atanh_doc[] =
"atanh(x)\n"
"\n"
"Return the hyperbolic arc tangent of x.";
122 123


124 125
static Py_complex
c_cos(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
126
{
Guido van Rossum's avatar
Guido van Rossum committed
127
	Py_complex r;
Guido van Rossum's avatar
Guido van Rossum committed
128 129 130 131 132
	r.real = cos(x.real)*cosh(x.imag);
	r.imag = -sin(x.real)*sinh(x.imag);
	return r;
}

133 134 135 136
static char c_cos_doc[] =
"cos(x)\n"
"n"
"Return the cosine of x.";
137 138


139 140
static Py_complex
c_cosh(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
141
{
Guido van Rossum's avatar
Guido van Rossum committed
142
	Py_complex r;
Guido van Rossum's avatar
Guido van Rossum committed
143 144 145 146 147
	r.real = cos(x.imag)*cosh(x.real);
	r.imag = sin(x.imag)*sinh(x.real);
	return r;
}

148 149 150 151
static char c_cosh_doc[] =
"cosh(x)\n"
"n"
"Return the hyperbolic cosine of x.";
152 153


154 155
static Py_complex
c_exp(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
156
{
Guido van Rossum's avatar
Guido van Rossum committed
157
	Py_complex r;
Guido van Rossum's avatar
Guido van Rossum committed
158 159 160 161 162 163
	double l = exp(x.real);
	r.real = l*cos(x.imag);
	r.imag = l*sin(x.imag);
	return r;
}

164 165 166 167
static char c_exp_doc[] =
"exp(x)\n"
"\n"
"Return the exponential value e**x.";
168 169


170 171
static Py_complex
c_log(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
172
{
Guido van Rossum's avatar
Guido van Rossum committed
173
	Py_complex r;
Guido van Rossum's avatar
Guido van Rossum committed
174 175 176 177 178 179
	double l = hypot(x.real,x.imag);
	r.imag = atan2(x.imag, x.real);
	r.real = log(l);
	return r;
}

180 181 182 183
static char c_log_doc[] =
"log(x)\n"
"\n"
"Return the natural logarithm of x.";
184 185


186 187
static Py_complex
c_log10(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
188
{
Guido van Rossum's avatar
Guido van Rossum committed
189
	Py_complex r;
Guido van Rossum's avatar
Guido van Rossum committed
190 191 192 193 194 195
	double l = hypot(x.real,x.imag);
	r.imag = atan2(x.imag, x.real)/log(10.);
	r.real = log10(l);
	return r;
}

196 197 198 199
static char c_log10_doc[] =
"log10(x)\n"
"\n"
"Return the base-10 logarithm of x.";
200 201 202


/* internal function not available from Python */
203 204
static Py_complex
c_prodi(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
205
{
Guido van Rossum's avatar
Guido van Rossum committed
206
	Py_complex r;
Guido van Rossum's avatar
Guido van Rossum committed
207 208 209 210 211
	r.real = -x.imag;
	r.imag = x.real;
	return r;
}

212

213 214
static Py_complex
c_sin(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
215
{
Guido van Rossum's avatar
Guido van Rossum committed
216
	Py_complex r;
217 218
	r.real = sin(x.real) * cosh(x.imag);
	r.imag = cos(x.real) * sinh(x.imag);
Guido van Rossum's avatar
Guido van Rossum committed
219 220 221
	return r;
}

222 223 224 225
static char c_sin_doc[] =
"sin(x)\n"
"\n"
"Return the sine of x.";
226 227


228 229
static Py_complex
c_sinh(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
230
{
Guido van Rossum's avatar
Guido van Rossum committed
231
	Py_complex r;
232 233
	r.real = cos(x.imag) * sinh(x.real);
	r.imag = sin(x.imag) * cosh(x.real);
Guido van Rossum's avatar
Guido van Rossum committed
234 235 236
	return r;
}

237 238 239 240
static char c_sinh_doc[] =
"sinh(x)\n"
"\n"
"Return the hyperbolic sine of x.";
241 242


243 244
static Py_complex
c_sqrt(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
245
{
Guido van Rossum's avatar
Guido van Rossum committed
246
	Py_complex r;
Guido van Rossum's avatar
Guido van Rossum committed
247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268
	double s,d;
	if (x.real == 0. && x.imag == 0.)
		r = x;
	else {
		s = sqrt(0.5*(fabs(x.real) + hypot(x.real,x.imag)));
		d = 0.5*x.imag/s;
		if (x.real > 0.) {
			r.real = s;
			r.imag = d;
		}
		else if (x.imag >= 0.) {
			r.real = d;
			r.imag = s;
		}
		else {
			r.real = -d;
			r.imag = -s;
		}
	}
	return r;
}

269 270 271 272
static char c_sqrt_doc[] =
"sqrt(x)\n"
"\n"
"Return the square root of x.";
273 274


275 276
static Py_complex
c_tan(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
277
{
Guido van Rossum's avatar
Guido van Rossum committed
278
	Py_complex r;
Guido van Rossum's avatar
Guido van Rossum committed
279 280 281 282 283 284 285
	double sr,cr,shi,chi;
	double rs,is,rc,ic;
	double d;
	sr = sin(x.real);
	cr = cos(x.real);
	shi = sinh(x.imag);
	chi = cosh(x.imag);
286 287 288 289 290 291 292
	rs = sr * chi;
	is = cr * shi;
	rc = cr * chi;
	ic = -sr * shi;
	d = rc*rc + ic * ic;
	r.real = (rs*rc + is*ic) / d;
	r.imag = (is*rc - rs*ic) / d;
Guido van Rossum's avatar
Guido van Rossum committed
293 294 295
	return r;
}

296 297 298 299
static char c_tan_doc[] =
"tan(x)\n"
"\n"
"Return the tangent of x.";
300 301


302 303
static Py_complex
c_tanh(Py_complex x)
Guido van Rossum's avatar
Guido van Rossum committed
304
{
Guido van Rossum's avatar
Guido van Rossum committed
305
	Py_complex r;
Guido van Rossum's avatar
Guido van Rossum committed
306 307 308 309 310 311 312
	double si,ci,shr,chr;
	double rs,is,rc,ic;
	double d;
	si = sin(x.imag);
	ci = cos(x.imag);
	shr = sinh(x.real);
	chr = cosh(x.real);
313 314 315 316
	rs = ci * shr;
	is = si * chr;
	rc = ci * chr;
	ic = si * shr;
Guido van Rossum's avatar
Guido van Rossum committed
317
	d = rc*rc + ic*ic;
318 319
	r.real = (rs*rc + is*ic) / d;
	r.imag = (is*rc - rs*ic) / d;
Guido van Rossum's avatar
Guido van Rossum committed
320 321 322
	return r;
}

323 324 325 326
static char c_tanh_doc[] =
"tanh(x)\n"
"\n"
"Return the hyperbolic tangent of x.";
327

Guido van Rossum's avatar
Guido van Rossum committed
328 329 330

/* And now the glue to make them available from Python: */

Roger E. Masse's avatar
Roger E. Masse committed
331
static PyObject *
332
math_error(void)
Guido van Rossum's avatar
Guido van Rossum committed
333 334
{
	if (errno == EDOM)
Roger E. Masse's avatar
Roger E. Masse committed
335
		PyErr_SetString(PyExc_ValueError, "math domain error");
Guido van Rossum's avatar
Guido van Rossum committed
336
	else if (errno == ERANGE)
Roger E. Masse's avatar
Roger E. Masse committed
337 338
		PyErr_SetString(PyExc_OverflowError, "math range error");
	else    /* Unexpected math error */
339
		PyErr_SetFromErrno(PyExc_ValueError);
Guido van Rossum's avatar
Guido van Rossum committed
340 341 342
	return NULL;
}

Roger E. Masse's avatar
Roger E. Masse committed
343
static PyObject *
Peter Schneider-Kamp's avatar
Peter Schneider-Kamp committed
344
math_1(PyObject *args, Py_complex (*func)(Py_complex))
Guido van Rossum's avatar
Guido van Rossum committed
345
{
Guido van Rossum's avatar
Guido van Rossum committed
346
	Py_complex x;
Guido van Rossum's avatar
Guido van Rossum committed
347 348 349
	if (!PyArg_ParseTuple(args, "D", &x))
		return NULL;
	errno = 0;
350
	PyFPE_START_PROTECT("complex function", return 0)
Guido van Rossum's avatar
Guido van Rossum committed
351
	x = (*func)(x);
352
	PyFPE_END_PROTECT(x)
Guido van Rossum's avatar
Guido van Rossum committed
353 354 355 356 357
	CHECK(x.real);
	CHECK(x.imag);
	if (errno != 0)
		return math_error();
	else
Roger E. Masse's avatar
Roger E. Masse committed
358
		return PyComplex_FromCComplex(x);
Guido van Rossum's avatar
Guido van Rossum committed
359 360 361
}

#define FUNC1(stubname, func) \
Peter Schneider-Kamp's avatar
Peter Schneider-Kamp committed
362
	static PyObject * stubname(PyObject *self, PyObject *args) { \
Guido van Rossum's avatar
Guido van Rossum committed
363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383
		return math_1(args, func); \
	}

FUNC1(cmath_acos, c_acos)
FUNC1(cmath_acosh, c_acosh)
FUNC1(cmath_asin, c_asin)
FUNC1(cmath_asinh, c_asinh)
FUNC1(cmath_atan, c_atan)
FUNC1(cmath_atanh, c_atanh)
FUNC1(cmath_cos, c_cos)
FUNC1(cmath_cosh, c_cosh)
FUNC1(cmath_exp, c_exp)
FUNC1(cmath_log, c_log)
FUNC1(cmath_log10, c_log10)
FUNC1(cmath_sin, c_sin)
FUNC1(cmath_sinh, c_sinh)
FUNC1(cmath_sqrt, c_sqrt)
FUNC1(cmath_tan, c_tan)
FUNC1(cmath_tanh, c_tanh)


384 385 386
static char module_doc[] =
"This module is always available. It provides access to mathematical\n"
"functions for complex numbers.";
387

Roger E. Masse's avatar
Roger E. Masse committed
388
static PyMethodDef cmath_methods[] = {
389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404
	{"acos",   cmath_acos,  METH_VARARGS, c_acos_doc},
	{"acosh",  cmath_acosh, METH_VARARGS, c_acosh_doc},
	{"asin",   cmath_asin,  METH_VARARGS, c_asin_doc},
	{"asinh",  cmath_asinh, METH_VARARGS, c_asinh_doc},
	{"atan",   cmath_atan,  METH_VARARGS, c_atan_doc},
	{"atanh",  cmath_atanh, METH_VARARGS, c_atanh_doc},
	{"cos",    cmath_cos,   METH_VARARGS, c_cos_doc},
	{"cosh",   cmath_cosh,  METH_VARARGS, c_cosh_doc},
	{"exp",    cmath_exp,   METH_VARARGS, c_exp_doc},
	{"log",    cmath_log,   METH_VARARGS, c_log_doc},
	{"log10",  cmath_log10, METH_VARARGS, c_log10_doc},
	{"sin",    cmath_sin,   METH_VARARGS, c_sin_doc},
	{"sinh",   cmath_sinh,  METH_VARARGS, c_sinh_doc},
	{"sqrt",   cmath_sqrt,  METH_VARARGS, c_sqrt_doc},
	{"tan",    cmath_tan,   METH_VARARGS, c_tan_doc},
	{"tanh",   cmath_tanh,  METH_VARARGS, c_tanh_doc},
Guido van Rossum's avatar
Guido van Rossum committed
405 406 407
	{NULL,		NULL}		/* sentinel */
};

408
DL_EXPORT(void)
409
initcmath(void)
Guido van Rossum's avatar
Guido van Rossum committed
410
{
Roger E. Masse's avatar
Roger E. Masse committed
411
	PyObject *m, *d, *v;
412

413
	m = Py_InitModule3("cmath", cmath_methods, module_doc);
Roger E. Masse's avatar
Roger E. Masse committed
414 415 416 417 418 419
	d = PyModule_GetDict(m);
	PyDict_SetItemString(d, "pi",
			     v = PyFloat_FromDouble(atan(1.0) * 4.0));
	Py_DECREF(v);
	PyDict_SetItemString(d, "e", v = PyFloat_FromDouble(exp(1.0)));
	Py_DECREF(v);
Guido van Rossum's avatar
Guido van Rossum committed
420
}