Home | History | Annotate | Download | only in src 
      1 /*-
      2  * Copyright (c) 2008 Stephen L. Moshier <steve (at) moshier.net>
      3  *
      4  * Permission to use, copy, modify, and distribute this software for any
      5  * purpose with or without fee is hereby granted, provided that the above
      6  * copyright notice and this permission notice appear in all copies.
      7  *
      8  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
      9  * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
     10  * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
     11  * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
     12  * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
     13  * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
     14  * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
     15  */
     16 
     17 /*							cpow
     18  *
     19  *	Complex power function
     20  *
     21  *
     22  *
     23  * SYNOPSIS:
     24  *
     25  * double complex cpow();
     26  * double complex a, z, w;
     27  *
     28  * w = cpow (a, z);
     29  *
     30  *
     31  *
     32  * DESCRIPTION:
     33  *
     34  * Raises complex A to the complex Zth power.
     35  * Definition is per AMS55 # 4.2.8,
     36  * analytically equivalent to cpow(a,z) = cexp(z clog(a)).
     37  *
     38  * ACCURACY:
     39  *
     40  *                      Relative error:
     41  * arithmetic   domain     # trials      peak         rms
     42  *    IEEE      -10,+10     30000       9.4e-15     1.5e-15
     43  *
     44  */
     45 
     46 #include <sys/cdefs.h>
     47 __FBSDID("$FreeBSD: head/lib/msun/src/s_cpow.c 336563 2018-07-20 18:27:30Z dim $");
     48 
     49 #include <complex.h>
     50 #include <float.h>
     51 #include <math.h>
     52 #include "math_private.h"
     53 
     54 double complex
     55 cpow(double complex a, double complex z)
     56 {
     57 	double complex w;
     58 	double x, y, r, theta, absa, arga;
     59 
     60 	x = creal (z);
     61 	y = cimag (z);
     62 	absa = cabs (a);
     63 	if (absa == 0.0) {
     64 		return (CMPLX(0.0, 0.0));
     65 	}
     66 	arga = carg (a);
     67 	r = pow (absa, x);
     68 	theta = x * arga;
     69 	if (y != 0.0) {
     70 		r = r * exp (-y * arga);
     71 		theta = theta + y * log (absa);
     72 	}
     73 	w = CMPLX(r * cos (theta),  r * sin (theta));
     74 	return (w);
     75 }
     76