1 /* Complex math module */

     2 

     3 /* much code borrowed from mathmodule.c */

     4 

     5 #include "Python.h"

     6 /* we need DBL_MAX, DBL_MIN, DBL_EPSILON, DBL_MANT_DIG and FLT_RADIX from

     7    float.h.  We assume that FLT_RADIX is either 2 or 16. */

     8 #include <float.h>

     9 

    10 #if (FLT_RADIX != 2 && FLT_RADIX != 16)

    11 #error "Modules/cmathmodule.c expects FLT_RADIX to be 2 or 16"

    12 #endif

    13 

    14 #ifndef M_LN2

    15 #define M_LN2 (0.6931471805599453094) /* natural log of 2 */

    16 #endif

    17 

    18 #ifndef M_LN10

    19 #define M_LN10 (2.302585092994045684) /* natural log of 10 */

    20 #endif

    21 

    22 /*

    23    CM_LARGE_DOUBLE is used to avoid spurious overflow in the sqrt, log,

    24    inverse trig and inverse hyperbolic trig functions.  Its log is used in the

    25    evaluation of exp, cos, cosh, sin, sinh, tan, and tanh to avoid unecessary

    26    overflow.

    27  */

    28 

    29 #define CM_LARGE_DOUBLE (DBL_MAX/4.)

    30 #define CM_SQRT_LARGE_DOUBLE (sqrt(CM_LARGE_DOUBLE))

    31 #define CM_LOG_LARGE_DOUBLE (log(CM_LARGE_DOUBLE))

    32 #define CM_SQRT_DBL_MIN (sqrt(DBL_MIN))

    33 

    34 /* 

    35    CM_SCALE_UP is an odd integer chosen such that multiplication by

    36    2**CM_SCALE_UP is sufficient to turn a subnormal into a normal.

    37    CM_SCALE_DOWN is (-(CM_SCALE_UP+1)/2).  These scalings are used to compute

    38    square roots accurately when the real and imaginary parts of the argument

    39    are subnormal.

    40 */

    41 

    42 #if FLT_RADIX==2

    43 #define CM_SCALE_UP (2*(DBL_MANT_DIG/2) + 1)

    44 #elif FLT_RADIX==16

    45 #define CM_SCALE_UP (4*DBL_MANT_DIG+1)

    46 #endif

    47 #define CM_SCALE_DOWN (-(CM_SCALE_UP+1)/2)

    48 

    49 /* forward declarations */

    50 static Py_complex c_asinh(Py_complex);

    51 static Py_complex c_atanh(Py_complex);

    52 static Py_complex c_cosh(Py_complex);

    53 static Py_complex c_sinh(Py_complex);

    54 static Py_complex c_sqrt(Py_complex);

    55 static Py_complex c_tanh(Py_complex);

    56 static PyObject * math_error(void);

    57 

    58 /* Code to deal with special values (infinities, NaNs, etc.). */

    59 

    60 /* special_type takes a double and returns an integer code indicating

    61    the type of the double as follows:

    62 */

    63 

    64 enum special_types {

    65 	ST_NINF,	/* 0, negative infinity */

    66 	ST_NEG,		/* 1, negative finite number (nonzero) */

    67 	ST_NZERO,	/* 2, -0. */

    68 	ST_PZERO,	/* 3, +0. */

    69 	ST_POS,		/* 4, positive finite number (nonzero) */

    70 	ST_PINF,	/* 5, positive infinity */

    71 	ST_NAN,		/* 6, Not a Number */

    72 };

    73 

    74 static enum special_types

    75 special_type(double d)

    76 {

    77 	if (Py_IS_FINITE(d)) {

    78 		if (d != 0) {

    79 			if (copysign(1., d) == 1.)

    80 				return ST_POS;

    81 			else

    82 				return ST_NEG;

    83 		}

    84 		else {

    85 			if (copysign(1., d) == 1.)

    86 				return ST_PZERO;

    87 			else

    88 				return ST_NZERO;

    89 		}

    90 	}

    91 	if (Py_IS_NAN(d))

    92 		return ST_NAN;

    93 	if (copysign(1., d) == 1.)

    94 		return ST_PINF;

    95 	else

    96 		return ST_NINF;

    97 }

    98 

    99 #define SPECIAL_VALUE(z, table)						\

   100 	if (!Py_IS_FINITE((z).real) || !Py_IS_FINITE((z).imag)) {	\

   101 		errno = 0;                                              \

   102 		return table[special_type((z).real)]	                \

   103 			    [special_type((z).imag)];			\

   104 	}

   105 

   106 #define P Py_MATH_PI

   107 #define P14 0.25*Py_MATH_PI

   108 #define P12 0.5*Py_MATH_PI

   109 #define P34 0.75*Py_MATH_PI

   110 #define INF Py_HUGE_VAL

   111 #define N Py_NAN

   112 #define U -9.5426319407711027e33 /* unlikely value, used as placeholder */

   113 

   114 /* First, the C functions that do the real work.  Each of the c_*

   115    functions computes and returns the C99 Annex G recommended result

   116    and also sets errno as follows: errno = 0 if no floating-point

   117    exception is associated with the result; errno = EDOM if C99 Annex

   118    G recommends raising divide-by-zero or invalid for this result; and

   119    errno = ERANGE where the overflow floating-point signal should be

   120    raised.

   121 */

   122 

   123 static Py_complex acos_special_values[7][7];

   124 

   125 static Py_complex

   126 c_acos(Py_complex z)

   127 {

   128 	Py_complex s1, s2, r;

   129 

   130 	SPECIAL_VALUE(z, acos_special_values);

   131 

   132 	if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {

   133 		/* avoid unnecessary overflow for large arguments */

   134 		r.real = atan2(fabs(z.imag), z.real);

   135 		/* split into cases to make sure that the branch cut has the

   136 		   correct continuity on systems with unsigned zeros */

   137 		if (z.real < 0.) {

   138 			r.imag = -copysign(log(hypot(z.real/2., z.imag/2.)) +

   139 					   M_LN2*2., z.imag);

   140 		} else {

   141 			r.imag = copysign(log(hypot(z.real/2., z.imag/2.)) +

   142 					  M_LN2*2., -z.imag);

   143 		}

   144 	} else {

   145 		s1.real = 1.-z.real;

   146 		s1.imag = -z.imag;

   147 		s1 = c_sqrt(s1);

   148 		s2.real = 1.+z.real;

   149 		s2.imag = z.imag;

   150 		s2 = c_sqrt(s2);

   151 		r.real = 2.*atan2(s1.real, s2.real);

   152 		r.imag = asinh(s2.real*s1.imag - s2.imag*s1.real);

   153 	}

   154 	errno = 0;

   155 	return r;

   156 }

   157 

   158 PyDoc_STRVAR(c_acos_doc,

   159 "acos(x)\n"

   160 "\n"

   161 "Return the arc cosine of x.");

   162 

   163 

   164 static Py_complex acosh_special_values[7][7];

   165 

   166 static Py_complex

   167 c_acosh(Py_complex z)

   168 {

   169 	Py_complex s1, s2, r;

   170 

   171 	SPECIAL_VALUE(z, acosh_special_values);

   172 

   173 	if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {

   174 		/* avoid unnecessary overflow for large arguments */

   175 		r.real = log(hypot(z.real/2., z.imag/2.)) + M_LN2*2.;

   176 		r.imag = atan2(z.imag, z.real);

   177 	} else {

   178 		s1.real = z.real - 1.;

   179 		s1.imag = z.imag;

   180 		s1 = c_sqrt(s1);

   181 		s2.real = z.real + 1.;

   182 		s2.imag = z.imag;

   183 		s2 = c_sqrt(s2);

   184 		r.real = asinh(s1.real*s2.real + s1.imag*s2.imag);

   185 		r.imag = 2.*atan2(s1.imag, s2.real);

   186 	}

   187 	errno = 0;

   188 	return r;

   189 }

   190 

   191 PyDoc_STRVAR(c_acosh_doc,

   192 "acosh(x)\n"

   193 "\n"

   194 "Return the hyperbolic arccosine of x.");

   195 

   196 

   197 static Py_complex

   198 c_asin(Py_complex z)

   199 {

   200 	/* asin(z) = -i asinh(iz) */

   201 	Py_complex s, r;

   202 	s.real = -z.imag;

   203 	s.imag = z.real;

   204 	s = c_asinh(s);

   205 	r.real = s.imag;

   206 	r.imag = -s.real;

   207 	return r;

   208 }

   209 

   210 PyDoc_STRVAR(c_asin_doc,

   211 "asin(x)\n"

   212 "\n"

   213 "Return the arc sine of x.");

   214 

   215 

   216 static Py_complex asinh_special_values[7][7];

   217 

   218 static Py_complex

   219 c_asinh(Py_complex z)

   220 {

   221 	Py_complex s1, s2, r;

   222 

   223 	SPECIAL_VALUE(z, asinh_special_values);

   224 

   225 	if (fabs(z.real) > CM_LARGE_DOUBLE || fabs(z.imag) > CM_LARGE_DOUBLE) {

   226 		if (z.imag >= 0.) {

   227 			r.real = copysign(log(hypot(z.real/2., z.imag/2.)) +

   228 					  M_LN2*2., z.real);

   229 		} else {

   230 			r.real = -copysign(log(hypot(z.real/2., z.imag/2.)) +

   231 					   M_LN2*2., -z.real);

   232 		}

   233 		r.imag = atan2(z.imag, fabs(z.real));

   234 	} else {

   235 		s1.real = 1.+z.imag;

   236 		s1.imag = -z.real;

   237 		s1 = c_sqrt(s1);

   238 		s2.real = 1.-z.imag;

   239 		s2.imag = z.real;

   240 		s2 = c_sqrt(s2);

   241 		r.real = asinh(s1.real*s2.imag-s2.real*s1.imag);

   242 		r.imag = atan2(z.imag, s1.real*s2.real-s1.imag*s2.imag);

   243 	}

   244 	errno = 0;

   245 	return r;

   246 }

   247 

   248 PyDoc_STRVAR(c_asinh_doc,

   249 "asinh(x)\n"

   250 "\n"

   251 "Return the hyperbolic arc sine of x.");

   252 

   253 

   254 static Py_complex

   255 c_atan(Py_complex z)

   256 {

   257 	/* atan(z) = -i atanh(iz) */

   258 	Py_complex s, r;

   259 	s.real = -z.imag;

   260 	s.imag = z.real;

   261 	s = c_atanh(s);

   262 	r.real = s.imag;

   263 	r.imag = -s.real;

   264 	return r;

   265 }

   266 

   267 /* Windows screws up atan2 for inf and nan, and alpha Tru64 5.1 doesn't follow

   268    C99 for atan2(0., 0.). */

   269 static double

   270 c_atan2(Py_complex z)

   271 {

   272 	if (Py_IS_NAN(z.real) || Py_IS_NAN(z.imag))

   273 		return Py_NAN;

   274 	if (Py_IS_INFINITY(z.imag)) {

   275 		if (Py_IS_INFINITY(z.real)) {

   276 			if (copysign(1., z.real) == 1.)

   277 				/* atan2(+-inf, +inf) == +-pi/4 */

   278 				return copysign(0.25*Py_MATH_PI, z.imag);

   279 			else

   280 				/* atan2(+-inf, -inf) == +-pi*3/4 */

   281 				return copysign(0.75*Py_MATH_PI, z.imag);

   282 		}

   283 		/* atan2(+-inf, x) == +-pi/2 for finite x */

   284 		return copysign(0.5*Py_MATH_PI, z.imag);

   285 	}

   286 	if (Py_IS_INFINITY(z.real) || z.imag == 0.) {

   287 		if (copysign(1., z.real) == 1.)

   288 			/* atan2(+-y, +inf) = atan2(+-0, +x) = +-0. */

   289 			return copysign(0., z.imag);

   290 		else

   291 			/* atan2(+-y, -inf) = atan2(+-0., -x) = +-pi. */

   292 			return copysign(Py_MATH_PI, z.imag);

   293 	}

   294 	return atan2(z.imag, z.real);

   295 }

   296 

   297 PyDoc_STRVAR(c_atan_doc,

   298 "atan(x)\n"

   299 "\n"

   300 "Return the arc tangent of x.");

   301 

   302 

   303 static Py_complex atanh_special_values[7][7];

   304 

   305 static Py_complex

   306 c_atanh(Py_complex z)

   307 {

   308 	Py_complex r;

   309 	double ay, h;

   310 

   311 	SPECIAL_VALUE(z, atanh_special_values);

   312 

   313 	/* Reduce to case where z.real >= 0., using atanh(z) = -atanh(-z). */

   314 	if (z.real < 0.) {

   315 		return c_neg(c_atanh(c_neg(z)));

   316 	}

   317 

   318 	ay = fabs(z.imag);

   319 	if (z.real > CM_SQRT_LARGE_DOUBLE || ay > CM_SQRT_LARGE_DOUBLE) {

   320 		/*

   321 		   if abs(z) is large then we use the approximation

   322 		   atanh(z) ~ 1/z +/- i*pi/2 (+/- depending on the sign

   323 		   of z.imag)

   324 		*/

   325 		h = hypot(z.real/2., z.imag/2.);  /* safe from overflow */

   326 		r.real = z.real/4./h/h;

   327 		/* the two negations in the next line cancel each other out

   328 		   except when working with unsigned zeros: they're there to

   329 		   ensure that the branch cut has the correct continuity on

   330 		   systems that don't support signed zeros */

   331 		r.imag = -copysign(Py_MATH_PI/2., -z.imag);

   332 		errno = 0;

   333 	} else if (z.real == 1. && ay < CM_SQRT_DBL_MIN) {

   334 		/* C99 standard says:  atanh(1+/-0.) should be inf +/- 0i */

   335 		if (ay == 0.) {

   336 			r.real = INF;

   337 			r.imag = z.imag;

   338 			errno = EDOM;

   339 		} else {

   340 			r.real = -log(sqrt(ay)/sqrt(hypot(ay, 2.)));

   341 			r.imag = copysign(atan2(2., -ay)/2, z.imag);

   342 			errno = 0;

   343 		}

   344 	} else {

   345 		r.real = log1p(4.*z.real/((1-z.real)*(1-z.real) + ay*ay))/4.;

   346 		r.imag = -atan2(-2.*z.imag, (1-z.real)*(1+z.real) - ay*ay)/2.;

   347 		errno = 0;

   348 	}

   349 	return r;

   350 }

   351 

   352 PyDoc_STRVAR(c_atanh_doc,

   353 "atanh(x)\n"

   354 "\n"

   355 "Return the hyperbolic arc tangent of x.");

   356 

   357 

   358 static Py_complex

   359 c_cos(Py_complex z)

   360 {

   361 	/* cos(z) = cosh(iz) */

   362 	Py_complex r;

   363 	r.real = -z.imag;

   364 	r.imag = z.real;

   365 	r = c_cosh(r);

   366 	return r;

   367 }

   368 

   369 PyDoc_STRVAR(c_cos_doc,

   370 "cos(x)\n"

   371 "n"

   372 "Return the cosine of x.");

   373 

   374 

   375 /* cosh(infinity + i*y) needs to be dealt with specially */

   376 static Py_complex cosh_special_values[7][7];

   377 

   378 static Py_complex

   379 c_cosh(Py_complex z)

   380 {

   381 	Py_complex r;

   382 	double x_minus_one;

   383 

   384 	/* special treatment for cosh(+/-inf + iy) if y is not a NaN */

   385 	if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {

   386 		if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag) &&

   387 		    (z.imag != 0.)) {

   388 			if (z.real > 0) {

   389 				r.real = copysign(INF, cos(z.imag));

   390 				r.imag = copysign(INF, sin(z.imag));

   391 			}

   392 			else {

   393 				r.real = copysign(INF, cos(z.imag));

   394 				r.imag = -copysign(INF, sin(z.imag));

   395 			}

   396 		}

   397 		else {

   398 			r = cosh_special_values[special_type(z.real)]

   399 				               [special_type(z.imag)];

   400 		}

   401 		/* need to set errno = EDOM if y is +/- infinity and x is not

   402 		   a NaN */

   403 		if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))

   404 			errno = EDOM;

   405 		else

   406 			errno = 0;

   407 		return r;

   408 	}

   409 

   410 	if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {

   411 		/* deal correctly with cases where cosh(z.real) overflows but

   412 		   cosh(z) does not. */

   413 		x_minus_one = z.real - copysign(1., z.real);

   414 		r.real = cos(z.imag) * cosh(x_minus_one) * Py_MATH_E;

   415 		r.imag = sin(z.imag) * sinh(x_minus_one) * Py_MATH_E;

   416 	} else {

   417 		r.real = cos(z.imag) * cosh(z.real);

   418 		r.imag = sin(z.imag) * sinh(z.real);

   419 	}

   420 	/* detect overflow, and set errno accordingly */

   421 	if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))

   422 		errno = ERANGE;

   423 	else

   424 		errno = 0;

   425 	return r;

   426 }

   427 

   428 PyDoc_STRVAR(c_cosh_doc,

   429 "cosh(x)\n"

   430 "n"

   431 "Return the hyperbolic cosine of x.");

   432 

   433 

   434 /* exp(infinity + i*y) and exp(-infinity + i*y) need special treatment for

   435    finite y */

   436 static Py_complex exp_special_values[7][7];

   437 

   438 static Py_complex

   439 c_exp(Py_complex z)

   440 {

   441 	Py_complex r;

   442 	double l;

   443 

   444 	if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {

   445 		if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)

   446 		    && (z.imag != 0.)) {

   447 			if (z.real > 0) {

   448 				r.real = copysign(INF, cos(z.imag));

   449 				r.imag = copysign(INF, sin(z.imag));

   450 			}

   451 			else {

   452 				r.real = copysign(0., cos(z.imag));

   453 				r.imag = copysign(0., sin(z.imag));

   454 			}

   455 		}

   456 		else {

   457 			r = exp_special_values[special_type(z.real)]

   458 				              [special_type(z.imag)];

   459 		}

   460 		/* need to set errno = EDOM if y is +/- infinity and x is not

   461 		   a NaN and not -infinity */

   462 		if (Py_IS_INFINITY(z.imag) &&

   463 		    (Py_IS_FINITE(z.real) ||

   464 		     (Py_IS_INFINITY(z.real) && z.real > 0)))

   465 			errno = EDOM;

   466 		else

   467 			errno = 0;

   468 		return r;

   469 	}

   470 

   471 	if (z.real > CM_LOG_LARGE_DOUBLE) {

   472 		l = exp(z.real-1.);

   473 		r.real = l*cos(z.imag)*Py_MATH_E;

   474 		r.imag = l*sin(z.imag)*Py_MATH_E;

   475 	} else {

   476 		l = exp(z.real);

   477 		r.real = l*cos(z.imag);

   478 		r.imag = l*sin(z.imag);

   479 	}

   480 	/* detect overflow, and set errno accordingly */

   481 	if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))

   482 		errno = ERANGE;

   483 	else

   484 		errno = 0;

   485 	return r;

   486 }

   487 

   488 PyDoc_STRVAR(c_exp_doc,

   489 "exp(x)\n"

   490 "\n"

   491 "Return the exponential value e**x.");

   492 

   493 

   494 static Py_complex log_special_values[7][7];

   495 

   496 static Py_complex

   497 c_log(Py_complex z)

   498 {

   499 	/*

   500 	   The usual formula for the real part is log(hypot(z.real, z.imag)).

   501 	   There are four situations where this formula is potentially

   502 	   problematic:

   503 

   504 	   (1) the absolute value of z is subnormal.  Then hypot is subnormal,

   505 	   so has fewer than the usual number of bits of accuracy, hence may

   506 	   have large relative error.  This then gives a large absolute error

   507 	   in the log.  This can be solved by rescaling z by a suitable power

   508 	   of 2.

   509 

   510 	   (2) the absolute value of z is greater than DBL_MAX (e.g. when both

   511 	   z.real and z.imag are within a factor of 1/sqrt(2) of DBL_MAX)

   512 	   Again, rescaling solves this.

   513 

   514 	   (3) the absolute value of z is close to 1.  In this case it's

   515 	   difficult to achieve good accuracy, at least in part because a

   516 	   change of 1ulp in the real or imaginary part of z can result in a

   517 	   change of billions of ulps in the correctly rounded answer.

   518 

   519 	   (4) z = 0.  The simplest thing to do here is to call the

   520 	   floating-point log with an argument of 0, and let its behaviour

   521 	   (returning -infinity, signaling a floating-point exception, setting

   522 	   errno, or whatever) determine that of c_log.  So the usual formula

   523 	   is fine here.

   524 

   525 	 */

   526 

   527 	Py_complex r;

   528 	double ax, ay, am, an, h;

   529 

   530 	SPECIAL_VALUE(z, log_special_values);

   531 

   532 	ax = fabs(z.real);

   533 	ay = fabs(z.imag);

   534 

   535 	if (ax > CM_LARGE_DOUBLE || ay > CM_LARGE_DOUBLE) {

   536 		r.real = log(hypot(ax/2., ay/2.)) + M_LN2;

   537 	} else if (ax < DBL_MIN && ay < DBL_MIN) {

   538 		if (ax > 0. || ay > 0.) {

   539 			/* catch cases where hypot(ax, ay) is subnormal */

   540 			r.real = log(hypot(ldexp(ax, DBL_MANT_DIG),

   541 				 ldexp(ay, DBL_MANT_DIG))) - DBL_MANT_DIG*M_LN2;

   542 		}

   543 		else {

   544 			/* log(+/-0. +/- 0i) */

   545 			r.real = -INF;

   546 			r.imag = atan2(z.imag, z.real);

   547 			errno = EDOM;

   548 			return r;

   549 		}

   550 	} else {

   551 		h = hypot(ax, ay);

   552 		if (0.71 <= h && h <= 1.73) {

   553 			am = ax > ay ? ax : ay;  /* max(ax, ay) */

   554 			an = ax > ay ? ay : ax;  /* min(ax, ay) */

   555 			r.real = log1p((am-1)*(am+1)+an*an)/2.;

   556 		} else {

   557 			r.real = log(h);

   558 		}

   559 	}

   560 	r.imag = atan2(z.imag, z.real);

   561 	errno = 0;

   562 	return r;

   563 }

   564 

   565 

   566 static Py_complex

   567 c_log10(Py_complex z)

   568 {

   569 	Py_complex r;

   570 	int errno_save;

   571 

   572 	r = c_log(z);

   573 	errno_save = errno; /* just in case the divisions affect errno */

   574 	r.real = r.real / M_LN10;

   575 	r.imag = r.imag / M_LN10;

   576 	errno = errno_save;

   577 	return r;

   578 }

   579 

   580 PyDoc_STRVAR(c_log10_doc,

   581 "log10(x)\n"

   582 "\n"

   583 "Return the base-10 logarithm of x.");

   584 

   585 

   586 static Py_complex

   587 c_sin(Py_complex z)

   588 {

   589 	/* sin(z) = -i sin(iz) */

   590 	Py_complex s, r;

   591 	s.real = -z.imag;

   592 	s.imag = z.real;

   593 	s = c_sinh(s);

   594 	r.real = s.imag;

   595 	r.imag = -s.real;

   596 	return r;

   597 }

   598 

   599 PyDoc_STRVAR(c_sin_doc,

   600 "sin(x)\n"

   601 "\n"

   602 "Return the sine of x.");

   603 

   604 

   605 /* sinh(infinity + i*y) needs to be dealt with specially */

   606 static Py_complex sinh_special_values[7][7];

   607 

   608 static Py_complex

   609 c_sinh(Py_complex z)

   610 {

   611 	Py_complex r;

   612 	double x_minus_one;

   613 

   614 	/* special treatment for sinh(+/-inf + iy) if y is finite and

   615 	   nonzero */

   616 	if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {

   617 		if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)

   618 		    && (z.imag != 0.)) {

   619 			if (z.real > 0) {

   620 				r.real = copysign(INF, cos(z.imag));

   621 				r.imag = copysign(INF, sin(z.imag));

   622 			}

   623 			else {

   624 				r.real = -copysign(INF, cos(z.imag));

   625 				r.imag = copysign(INF, sin(z.imag));

   626 			}

   627 		}

   628 		else {

   629 			r = sinh_special_values[special_type(z.real)]

   630 				               [special_type(z.imag)];

   631 		}

   632 		/* need to set errno = EDOM if y is +/- infinity and x is not

   633 		   a NaN */

   634 		if (Py_IS_INFINITY(z.imag) && !Py_IS_NAN(z.real))

   635 			errno = EDOM;

   636 		else

   637 			errno = 0;

   638 		return r;

   639 	}

   640 

   641 	if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {

   642 		x_minus_one = z.real - copysign(1., z.real);

   643 		r.real = cos(z.imag) * sinh(x_minus_one) * Py_MATH_E;

   644 		r.imag = sin(z.imag) * cosh(x_minus_one) * Py_MATH_E;

   645 	} else {

   646 		r.real = cos(z.imag) * sinh(z.real);

   647 		r.imag = sin(z.imag) * cosh(z.real);

   648 	}

   649 	/* detect overflow, and set errno accordingly */

   650 	if (Py_IS_INFINITY(r.real) || Py_IS_INFINITY(r.imag))

   651 		errno = ERANGE;

   652 	else

   653 		errno = 0;

   654 	return r;

   655 }

   656 

   657 PyDoc_STRVAR(c_sinh_doc,

   658 "sinh(x)\n"

   659 "\n"

   660 "Return the hyperbolic sine of x.");

   661 

   662 

   663 static Py_complex sqrt_special_values[7][7];

   664 

   665 static Py_complex

   666 c_sqrt(Py_complex z)

   667 {

   668 	/*

   669 	   Method: use symmetries to reduce to the case when x = z.real and y

   670 	   = z.imag are nonnegative.  Then the real part of the result is

   671 	   given by

   672 

   673 	     s = sqrt((x + hypot(x, y))/2)

   674 

   675 	   and the imaginary part is

   676 

   677 	     d = (y/2)/s

   678 

   679 	   If either x or y is very large then there's a risk of overflow in

   680 	   computation of the expression x + hypot(x, y).  We can avoid this

   681 	   by rewriting the formula for s as:

   682 

   683 	     s = 2*sqrt(x/8 + hypot(x/8, y/8))

   684 

   685 	   This costs us two extra multiplications/divisions, but avoids the

   686 	   overhead of checking for x and y large.

   687 

   688 	   If both x and y are subnormal then hypot(x, y) may also be

   689 	   subnormal, so will lack full precision.  We solve this by rescaling

   690 	   x and y by a sufficiently large power of 2 to ensure that x and y

   691 	   are normal.

   692 	*/

   693 

   694 

   695 	Py_complex r;

   696 	double s,d;

   697 	double ax, ay;

   698 

   699 	SPECIAL_VALUE(z, sqrt_special_values);

   700 

   701 	if (z.real == 0. && z.imag == 0.) {

   702 		r.real = 0.;

   703 		r.imag = z.imag;

   704 		return r;

   705 	}

   706 

   707 	ax = fabs(z.real);

   708 	ay = fabs(z.imag);

   709 

   710 	if (ax < DBL_MIN && ay < DBL_MIN && (ax > 0. || ay > 0.)) {

   711 		/* here we catch cases where hypot(ax, ay) is subnormal */

   712 		ax = ldexp(ax, CM_SCALE_UP);

   713 		s = ldexp(sqrt(ax + hypot(ax, ldexp(ay, CM_SCALE_UP))),

   714 			  CM_SCALE_DOWN);

   715 	} else {

   716 		ax /= 8.;

   717 		s = 2.*sqrt(ax + hypot(ax, ay/8.));

   718 	}

   719 	d = ay/(2.*s);

   720 

   721 	if (z.real >= 0.) {

   722 		r.real = s;

   723 		r.imag = copysign(d, z.imag);

   724 	} else {

   725 		r.real = d;

   726 		r.imag = copysign(s, z.imag);

   727 	}

   728 	errno = 0;

   729 	return r;

   730 }

   731 

   732 PyDoc_STRVAR(c_sqrt_doc,

   733 "sqrt(x)\n"

   734 "\n"

   735 "Return the square root of x.");

   736 

   737 

   738 static Py_complex

   739 c_tan(Py_complex z)

   740 {

   741 	/* tan(z) = -i tanh(iz) */

   742 	Py_complex s, r;

   743 	s.real = -z.imag;

   744 	s.imag = z.real;

   745 	s = c_tanh(s);

   746 	r.real = s.imag;

   747 	r.imag = -s.real;

   748 	return r;

   749 }

   750 

   751 PyDoc_STRVAR(c_tan_doc,

   752 "tan(x)\n"

   753 "\n"

   754 "Return the tangent of x.");

   755 

   756 

   757 /* tanh(infinity + i*y) needs to be dealt with specially */

   758 static Py_complex tanh_special_values[7][7];

   759 

   760 static Py_complex

   761 c_tanh(Py_complex z)

   762 {

   763 	/* Formula:

   764 

   765 	   tanh(x+iy) = (tanh(x)(1+tan(y)^2) + i tan(y)(1-tanh(x))^2) /

   766 	   (1+tan(y)^2 tanh(x)^2)

   767 

   768 	   To avoid excessive roundoff error, 1-tanh(x)^2 is better computed

   769 	   as 1/cosh(x)^2.  When abs(x) is large, we approximate 1-tanh(x)^2

   770 	   by 4 exp(-2*x) instead, to avoid possible overflow in the

   771 	   computation of cosh(x).

   772 

   773 	*/

   774 

   775 	Py_complex r;

   776 	double tx, ty, cx, txty, denom;

   777 

   778 	/* special treatment for tanh(+/-inf + iy) if y is finite and

   779 	   nonzero */

   780 	if (!Py_IS_FINITE(z.real) || !Py_IS_FINITE(z.imag)) {

   781 		if (Py_IS_INFINITY(z.real) && Py_IS_FINITE(z.imag)

   782 		    && (z.imag != 0.)) {

   783 			if (z.real > 0) {

   784 				r.real = 1.0;

   785 				r.imag = copysign(0.,

   786 						  2.*sin(z.imag)*cos(z.imag));

   787 			}

   788 			else {

   789 				r.real = -1.0;

   790 				r.imag = copysign(0.,

   791 						  2.*sin(z.imag)*cos(z.imag));

   792 			}

   793 		}

   794 		else {

   795 			r = tanh_special_values[special_type(z.real)]

   796 				               [special_type(z.imag)];

   797 		}

   798 		/* need to set errno = EDOM if z.imag is +/-infinity and

   799 		   z.real is finite */

   800 		if (Py_IS_INFINITY(z.imag) && Py_IS_FINITE(z.real))

   801 			errno = EDOM;

   802 		else

   803 			errno = 0;

   804 		return r;

   805 	}

   806 

   807 	/* danger of overflow in 2.*z.imag !*/

   808 	if (fabs(z.real) > CM_LOG_LARGE_DOUBLE) {

   809 		r.real = copysign(1., z.real);

   810 		r.imag = 4.*sin(z.imag)*cos(z.imag)*exp(-2.*fabs(z.real));

   811 	} else {

   812 		tx = tanh(z.real);

   813 		ty = tan(z.imag);

   814 		cx = 1./cosh(z.real);

   815 		txty = tx*ty;

   816 		denom = 1. + txty*txty;

   817 		r.real = tx*(1.+ty*ty)/denom;

   818 		r.imag = ((ty/denom)*cx)*cx;

   819 	}

   820 	errno = 0;

   821 	return r;

   822 }

   823 

   824 PyDoc_STRVAR(c_tanh_doc,

   825 "tanh(x)\n"

   826 "\n"

   827 "Return the hyperbolic tangent of x.");

   828 

   829 

   830 static PyObject *

   831 cmath_log(PyObject *self, PyObject *args)

   832 {

   833 	Py_complex x;

   834 	Py_complex y;

   835 

   836 	if (!PyArg_ParseTuple(args, "D|D", &x, &y))

   837 		return NULL;

   838 

   839 	errno = 0;

   840 	PyFPE_START_PROTECT("complex function", return 0)

   841 	x = c_log(x);

   842 	if (PyTuple_GET_SIZE(args) == 2) {

   843 		y = c_log(y);

   844 		x = c_quot(x, y);

   845 	}

   846 	PyFPE_END_PROTECT(x)

   847 	if (errno != 0)

   848 		return math_error();

   849 	return PyComplex_FromCComplex(x);

   850 }

   851 

   852 PyDoc_STRVAR(cmath_log_doc,

   853 "log(x[, base]) -> the logarithm of x to the given base.\n\

   854 If the base not specified, returns the natural logarithm (base e) of x.");

   855 

   856 

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

   858 

   859 static PyObject *

   860 math_error(void)

   861 {

   862 	if (errno == EDOM)

   863 		PyErr_SetString(PyExc_ValueError, "math domain error");

   864 	else if (errno == ERANGE)

   865 		PyErr_SetString(PyExc_OverflowError, "math range error");

   866 	else    /* Unexpected math error */

   867 		PyErr_SetFromErrno(PyExc_ValueError);

   868 	return NULL;

   869 }

   870 

   871 static PyObject *

   872 math_1(PyObject *args, Py_complex (*func)(Py_complex))

   873 {

   874 	Py_complex x,r ;

   875 	if (!PyArg_ParseTuple(args, "D", &x))

   876 		return NULL;

   877 	errno = 0;

   878 	PyFPE_START_PROTECT("complex function", return 0);

   879 	r = (*func)(x);

   880 	PyFPE_END_PROTECT(r);

   881 	if (errno == EDOM) {

   882 		PyErr_SetString(PyExc_ValueError, "math domain error");

   883 		return NULL;

   884 	}

   885 	else if (errno == ERANGE) {

   886 		PyErr_SetString(PyExc_OverflowError, "math range error");

   887 		return NULL;

   888 	}

   889 	else {

   890 		return PyComplex_FromCComplex(r);

   891 	}

   892 }

   893 

   894 #define FUNC1(stubname, func) \

   895 	static PyObject * stubname(PyObject *self, PyObject *args) { \

   896 		return math_1(args, func); \

   897 	}

   898 

   899 FUNC1(cmath_acos, c_acos)

   900 FUNC1(cmath_acosh, c_acosh)

   901 FUNC1(cmath_asin, c_asin)

   902 FUNC1(cmath_asinh, c_asinh)

   903 FUNC1(cmath_atan, c_atan)

   904 FUNC1(cmath_atanh, c_atanh)

   905 FUNC1(cmath_cos, c_cos)

   906 FUNC1(cmath_cosh, c_cosh)

   907 FUNC1(cmath_exp, c_exp)

   908 FUNC1(cmath_log10, c_log10)

   909 FUNC1(cmath_sin, c_sin)

   910 FUNC1(cmath_sinh, c_sinh)

   911 FUNC1(cmath_sqrt, c_sqrt)

   912 FUNC1(cmath_tan, c_tan)

   913 FUNC1(cmath_tanh, c_tanh)

   914 

   915 static PyObject *

   916 cmath_phase(PyObject *self, PyObject *args)

   917 {

   918 	Py_complex z;

   919 	double phi;

   920 	if (!PyArg_ParseTuple(args, "D:phase", &z))

   921 		return NULL;

   922 	errno = 0;

   923 	PyFPE_START_PROTECT("arg function", return 0)

   924 	phi = c_atan2(z);

   925 	PyFPE_END_PROTECT(phi)

   926 	if (errno != 0)

   927 		return math_error();

   928 	else

   929 		return PyFloat_FromDouble(phi);

   930 }

   931 

   932 PyDoc_STRVAR(cmath_phase_doc,

   933 "phase(z) -> float\n\n\

   934 Return argument, also known as the phase angle, of a complex.");

   935 

   936 static PyObject *

   937 cmath_polar(PyObject *self, PyObject *args)

   938 {

   939 	Py_complex z;

   940 	double r, phi;

   941 	if (!PyArg_ParseTuple(args, "D:polar", &z))

   942 		return NULL;

   943 	PyFPE_START_PROTECT("polar function", return 0)

   944 	phi = c_atan2(z); /* should not cause any exception */

   945 	r = c_abs(z); /* sets errno to ERANGE on overflow;  otherwise 0 */

   946 	PyFPE_END_PROTECT(r)

   947 	if (errno != 0)

   948 		return math_error();

   949 	else

   950 		return Py_BuildValue("dd", r, phi);

   951 }

   952 

   953 PyDoc_STRVAR(cmath_polar_doc,

   954 "polar(z) -> r: float, phi: float\n\n\

   955 Convert a complex from rectangular coordinates to polar coordinates. r is\n\

   956 the distance from 0 and phi the phase angle.");

   957 

   958 /*

   959   rect() isn't covered by the C99 standard, but it's not too hard to

   960   figure out 'spirit of C99' rules for special value handing:

   961 

   962     rect(x, t) should behave like exp(log(x) + it) for positive-signed x

   963     rect(x, t) should behave like -exp(log(-x) + it) for negative-signed x

   964     rect(nan, t) should behave like exp(nan + it), except that rect(nan, 0)

   965       gives nan +- i0 with the sign of the imaginary part unspecified.

   966 

   967 */

   968 

   969 static Py_complex rect_special_values[7][7];

   970 

   971 static PyObject *

   972 cmath_rect(PyObject *self, PyObject *args)

   973 {

   974 	Py_complex z;

   975 	double r, phi;

   976 	if (!PyArg_ParseTuple(args, "dd:rect", &r, &phi))

   977 		return NULL;

   978 	errno = 0;

   979 	PyFPE_START_PROTECT("rect function", return 0)

   980 

   981 	/* deal with special values */

   982 	if (!Py_IS_FINITE(r) || !Py_IS_FINITE(phi)) {

   983 		/* if r is +/-infinity and phi is finite but nonzero then

   984 		   result is (+-INF +-INF i), but we need to compute cos(phi)

   985 		   and sin(phi) to figure out the signs. */

   986 		if (Py_IS_INFINITY(r) && (Py_IS_FINITE(phi)

   987 					  && (phi != 0.))) {

   988 			if (r > 0) {

   989 				z.real = copysign(INF, cos(phi));

   990 				z.imag = copysign(INF, sin(phi));

   991 			}

   992 			else {

   993 				z.real = -copysign(INF, cos(phi));

   994 				z.imag = -copysign(INF, sin(phi));

   995 			}

   996 		}

   997 		else {

   998 			z = rect_special_values[special_type(r)]

   999 				               [special_type(phi)];

  1000 		}

  1001 		/* need to set errno = EDOM if r is a nonzero number and phi

  1002 		   is infinite */

  1003 		if (r != 0. && !Py_IS_NAN(r) && Py_IS_INFINITY(phi))

  1004 			errno = EDOM;

  1005 		else

  1006 			errno = 0;

  1007 	}

  1008 	else {

  1009 		z.real = r * cos(phi);

  1010 		z.imag = r * sin(phi);

  1011 		errno = 0;

  1012 	}

  1013 

  1014 	PyFPE_END_PROTECT(z)

  1015 	if (errno != 0)

  1016 		return math_error();

  1017 	else

  1018 		return PyComplex_FromCComplex(z);

  1019 }

  1020 

  1021 PyDoc_STRVAR(cmath_rect_doc,

  1022 "rect(r, phi) -> z: complex\n\n\

  1023 Convert from polar coordinates to rectangular coordinates.");

  1024 

  1025 static PyObject *

  1026 cmath_isnan(PyObject *self, PyObject *args)

  1027 {

  1028 	Py_complex z;

  1029 	if (!PyArg_ParseTuple(args, "D:isnan", &z))

  1030 		return NULL;

  1031 	return PyBool_FromLong(Py_IS_NAN(z.real) || Py_IS_NAN(z.imag));

  1032 }

  1033 

  1034 PyDoc_STRVAR(cmath_isnan_doc,

  1035 "isnan(z) -> bool\n\

  1036 Checks if the real or imaginary part of z not a number (NaN)");

  1037 

  1038 static PyObject *

  1039 cmath_isinf(PyObject *self, PyObject *args)

  1040 {

  1041 	Py_complex z;

  1042 	if (!PyArg_ParseTuple(args, "D:isnan", &z))

  1043 		return NULL;

  1044 	return PyBool_FromLong(Py_IS_INFINITY(z.real) ||

  1045 			       Py_IS_INFINITY(z.imag));

  1046 }

  1047 

  1048 PyDoc_STRVAR(cmath_isinf_doc,

  1049 "isinf(z) -> bool\n\

  1050 Checks if the real or imaginary part of z is infinite.");

  1051 

  1052 

  1053 PyDoc_STRVAR(module_doc,

  1054 "This module is always available. It provides access to mathematical\n"

  1055 "functions for complex numbers.");

  1056 

  1057 static PyMethodDef cmath_methods[] = {

  1058 	{"acos",   cmath_acos,  METH_VARARGS, c_acos_doc},

  1059 	{"acosh",  cmath_acosh, METH_VARARGS, c_acosh_doc},

  1060 	{"asin",   cmath_asin,  METH_VARARGS, c_asin_doc},

  1061 	{"asinh",  cmath_asinh, METH_VARARGS, c_asinh_doc},

  1062 	{"atan",   cmath_atan,  METH_VARARGS, c_atan_doc},

  1063 	{"atanh",  cmath_atanh, METH_VARARGS, c_atanh_doc},

  1064 	{"cos",    cmath_cos,   METH_VARARGS, c_cos_doc},

  1065 	{"cosh",   cmath_cosh,  METH_VARARGS, c_cosh_doc},

  1066 	{"exp",    cmath_exp,   METH_VARARGS, c_exp_doc},

  1067 	{"isinf",  cmath_isinf, METH_VARARGS, cmath_isinf_doc},

  1068 	{"isnan",  cmath_isnan, METH_VARARGS, cmath_isnan_doc},

  1069 	{"log",    cmath_log,   METH_VARARGS, cmath_log_doc},

  1070 	{"log10",  cmath_log10, METH_VARARGS, c_log10_doc},

  1071 	{"phase",  cmath_phase, METH_VARARGS, cmath_phase_doc},

  1072 	{"polar",  cmath_polar, METH_VARARGS, cmath_polar_doc},

  1073 	{"rect",   cmath_rect,  METH_VARARGS, cmath_rect_doc},

  1074 	{"sin",    cmath_sin,   METH_VARARGS, c_sin_doc},

  1075 	{"sinh",   cmath_sinh,  METH_VARARGS, c_sinh_doc},

  1076 	{"sqrt",   cmath_sqrt,  METH_VARARGS, c_sqrt_doc},

  1077 	{"tan",    cmath_tan,   METH_VARARGS, c_tan_doc},

  1078 	{"tanh",   cmath_tanh,  METH_VARARGS, c_tanh_doc},

  1079 	{NULL,		NULL}		/* sentinel */

  1080 };

  1081 

  1082 PyMODINIT_FUNC

  1083 initcmath(void)

  1084 {

  1085 	PyObject *m;

  1086 

  1087 	m = Py_InitModule3("cmath", cmath_methods, module_doc);

  1088 	if (m == NULL)

  1089 		return;

  1090 

  1091 	PyModule_AddObject(m, "pi",

  1092                            PyFloat_FromDouble(Py_MATH_PI));

  1093 	PyModule_AddObject(m, "e", PyFloat_FromDouble(Py_MATH_E));

  1094 

  1095 	/* initialize special value tables */

  1096 

  1097 #define INIT_SPECIAL_VALUES(NAME, BODY) { Py_complex* p = (Py_complex*)NAME; BODY }

  1098 #define C(REAL, IMAG) p->real = REAL; p->imag = IMAG; ++p;

  1099 

  1100 	INIT_SPECIAL_VALUES(acos_special_values, {

  1101 	  C(P34,INF) C(P,INF)  C(P,INF)  C(P,-INF)  C(P,-INF)  C(P34,-INF) C(N,INF)

  1102 	  C(P12,INF) C(U,U)    C(U,U)    C(U,U)     C(U,U)     C(P12,-INF) C(N,N)

  1103 	  C(P12,INF) C(U,U)    C(P12,0.) C(P12,-0.) C(U,U)     C(P12,-INF) C(P12,N)

  1104 	  C(P12,INF) C(U,U)    C(P12,0.) C(P12,-0.) C(U,U)     C(P12,-INF) C(P12,N)

  1105 	  C(P12,INF) C(U,U)    C(U,U)    C(U,U)     C(U,U)     C(P12,-INF) C(N,N)

  1106 	  C(P14,INF) C(0.,INF) C(0.,INF) C(0.,-INF) C(0.,-INF) C(P14,-INF) C(N,INF)

  1107 	  C(N,INF)   C(N,N)    C(N,N)    C(N,N)     C(N,N)     C(N,-INF)   C(N,N)

  1108 	})

  1109 

  1110 	INIT_SPECIAL_VALUES(acosh_special_values, {

  1111 	  C(INF,-P34) C(INF,-P)  C(INF,-P)  C(INF,P)  C(INF,P)  C(INF,P34) C(INF,N)

  1112 	  C(INF,-P12) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,P12) C(N,N)

  1113 	  C(INF,-P12) C(U,U)     C(0.,-P12) C(0.,P12) C(U,U)    C(INF,P12) C(N,N)

  1114 	  C(INF,-P12) C(U,U)     C(0.,-P12) C(0.,P12) C(U,U)    C(INF,P12) C(N,N)

  1115 	  C(INF,-P12) C(U,U)     C(U,U)     C(U,U)    C(U,U)    C(INF,P12) C(N,N)

  1116 	  C(INF,-P14) C(INF,-0.) C(INF,-0.) C(INF,0.) C(INF,0.) C(INF,P14) C(INF,N)

  1117 	  C(INF,N)    C(N,N)     C(N,N)     C(N,N)    C(N,N)    C(INF,N)   C(N,N)

  1118 	})

  1119 

  1120 	INIT_SPECIAL_VALUES(asinh_special_values, {

  1121 	  C(-INF,-P14) C(-INF,-0.) C(-INF,-0.) C(-INF,0.) C(-INF,0.) C(-INF,P14) C(-INF,N)

  1122 	  C(-INF,-P12) C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(-INF,P12) C(N,N)

  1123 	  C(-INF,-P12) C(U,U)      C(-0.,-0.)  C(-0.,0.)  C(U,U)     C(-INF,P12) C(N,N)

  1124 	  C(INF,-P12)  C(U,U)      C(0.,-0.)   C(0.,0.)   C(U,U)     C(INF,P12)  C(N,N)

  1125 	  C(INF,-P12)  C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(INF,P12)  C(N,N)

  1126 	  C(INF,-P14)  C(INF,-0.)  C(INF,-0.)  C(INF,0.)  C(INF,0.)  C(INF,P14)  C(INF,N)

  1127 	  C(INF,N)     C(N,N)      C(N,-0.)    C(N,0.)    C(N,N)     C(INF,N)    C(N,N)

  1128 	})

  1129 

  1130 	INIT_SPECIAL_VALUES(atanh_special_values, {

  1131 	  C(-0.,-P12) C(-0.,-P12) C(-0.,-P12) C(-0.,P12) C(-0.,P12) C(-0.,P12) C(-0.,N)

  1132 	  C(-0.,-P12) C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(-0.,P12) C(N,N)

  1133 	  C(-0.,-P12) C(U,U)      C(-0.,-0.)  C(-0.,0.)  C(U,U)     C(-0.,P12) C(-0.,N)

  1134 	  C(0.,-P12)  C(U,U)      C(0.,-0.)   C(0.,0.)   C(U,U)     C(0.,P12)  C(0.,N)

  1135 	  C(0.,-P12)  C(U,U)      C(U,U)      C(U,U)     C(U,U)     C(0.,P12)  C(N,N)