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