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