16813d08fSMatt Macy /*- 26813d08fSMatt Macy * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 36813d08fSMatt Macy * 46813d08fSMatt Macy * Permission to use, copy, modify, and distribute this software for any 56813d08fSMatt Macy * purpose with or without fee is hereby granted, provided that the above 66813d08fSMatt Macy * copyright notice and this permission notice appear in all copies. 76813d08fSMatt Macy * 86813d08fSMatt Macy * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 96813d08fSMatt Macy * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 106813d08fSMatt Macy * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 116813d08fSMatt Macy * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 126813d08fSMatt Macy * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 136813d08fSMatt Macy * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 146813d08fSMatt Macy * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 156813d08fSMatt Macy */ 166813d08fSMatt Macy 176813d08fSMatt Macy /* cpowf 186813d08fSMatt Macy * 196813d08fSMatt Macy * Complex power function 206813d08fSMatt Macy * 216813d08fSMatt Macy * 226813d08fSMatt Macy * 236813d08fSMatt Macy * SYNOPSIS: 246813d08fSMatt Macy * 256813d08fSMatt Macy * float complex cpowf(); 266813d08fSMatt Macy * float complex a, z, w; 276813d08fSMatt Macy * 286813d08fSMatt Macy * w = cpowf (a, z); 296813d08fSMatt Macy * 306813d08fSMatt Macy * 316813d08fSMatt Macy * 326813d08fSMatt Macy * DESCRIPTION: 336813d08fSMatt Macy * 346813d08fSMatt Macy * Raises complex A to the complex Zth power. 356813d08fSMatt Macy * Definition is per AMS55 # 4.2.8, 366813d08fSMatt Macy * analytically equivalent to cpow(a,z) = cexp(z clog(a)). 376813d08fSMatt Macy * 386813d08fSMatt Macy * ACCURACY: 396813d08fSMatt Macy * 406813d08fSMatt Macy * Relative error: 416813d08fSMatt Macy * arithmetic domain # trials peak rms 426813d08fSMatt Macy * IEEE -10,+10 30000 9.4e-15 1.5e-15 436813d08fSMatt Macy * 446813d08fSMatt Macy */ 456813d08fSMatt Macy 466813d08fSMatt Macy #include <sys/cdefs.h> 476813d08fSMatt Macy __FBSDID("$FreeBSD$"); 486813d08fSMatt Macy 496813d08fSMatt Macy #include <complex.h> 506813d08fSMatt Macy #include <math.h> 51*5a4c3b83SDimitry Andric #include "math_private.h" 526813d08fSMatt Macy 536813d08fSMatt Macy float complex 546813d08fSMatt Macy cpowf(float complex a, float complex z) 556813d08fSMatt Macy { 566813d08fSMatt Macy float complex w; 576813d08fSMatt Macy float x, y, r, theta, absa, arga; 586813d08fSMatt Macy 596813d08fSMatt Macy x = crealf(z); 606813d08fSMatt Macy y = cimagf(z); 616813d08fSMatt Macy absa = cabsf (a); 626813d08fSMatt Macy if (absa == 0.0f) { 63*5a4c3b83SDimitry Andric return (CMPLXF(0.0f, 0.0f)); 646813d08fSMatt Macy } 656813d08fSMatt Macy arga = cargf (a); 666813d08fSMatt Macy r = powf (absa, x); 676813d08fSMatt Macy theta = x * arga; 686813d08fSMatt Macy if (y != 0.0f) { 696813d08fSMatt Macy r = r * expf (-y * arga); 706813d08fSMatt Macy theta = theta + y * logf (absa); 716813d08fSMatt Macy } 72*5a4c3b83SDimitry Andric w = CMPLXF(r * cosf (theta), r * sinf (theta)); 736813d08fSMatt Macy return (w); 746813d08fSMatt Macy } 75