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 /* cpowl
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 * long double complex cpowl();
26 6813d08fSMatt Macy * long double complex a, z, w;
27 6813d08fSMatt Macy *
28 6813d08fSMatt Macy * w = cpowl (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 <complex.h>
47 6813d08fSMatt Macy #include <math.h>
48 *5a4c3b83SDimitry Andric #include "math_private.h"
49 6813d08fSMatt Macy
50 6813d08fSMatt Macy long double complex
cpowl(long double complex a,long double complex z)51 6813d08fSMatt Macy cpowl(long double complex a, long double complex z)
52 6813d08fSMatt Macy {
53 6813d08fSMatt Macy long double complex w;
54 6813d08fSMatt Macy long double x, y, r, theta, absa, arga;
55 6813d08fSMatt Macy
56 6813d08fSMatt Macy x = creall(z);
57 6813d08fSMatt Macy y = cimagl(z);
58 6813d08fSMatt Macy absa = cabsl(a);
59 6813d08fSMatt Macy if (absa == 0.0L) {
60 *5a4c3b83SDimitry Andric return (CMPLXL(0.0L, 0.0L));
61 6813d08fSMatt Macy }
62 6813d08fSMatt Macy arga = cargl(a);
63 6813d08fSMatt Macy r = powl(absa, x);
64 6813d08fSMatt Macy theta = x * arga;
65 6813d08fSMatt Macy if (y != 0.0L) {
66 6813d08fSMatt Macy r = r * expl(-y * arga);
67 6813d08fSMatt Macy theta = theta + y * logl(absa);
68 6813d08fSMatt Macy }
69 *5a4c3b83SDimitry Andric w = CMPLXL(r * cosl(theta), r * sinl(theta));
70 6813d08fSMatt Macy return (w);
71 6813d08fSMatt Macy }
72