1*2f2c0062Sguenther /* $OpenBSD: s_cpow.c,v 1.7 2016/09/12 19:47:02 guenther Exp $ */
27b36286aSmartynas /*
37b36286aSmartynas * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
47b36286aSmartynas *
57b36286aSmartynas * Permission to use, copy, modify, and distribute this software for any
67b36286aSmartynas * purpose with or without fee is hereby granted, provided that the above
77b36286aSmartynas * copyright notice and this permission notice appear in all copies.
87b36286aSmartynas *
97b36286aSmartynas * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
107b36286aSmartynas * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
117b36286aSmartynas * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
127b36286aSmartynas * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
137b36286aSmartynas * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
147b36286aSmartynas * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
157b36286aSmartynas * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
167b36286aSmartynas */
177b36286aSmartynas
187b36286aSmartynas /* cpow
197b36286aSmartynas *
207b36286aSmartynas * Complex power function
217b36286aSmartynas *
227b36286aSmartynas *
237b36286aSmartynas *
247b36286aSmartynas * SYNOPSIS:
257b36286aSmartynas *
267b36286aSmartynas * double complex cpow();
277b36286aSmartynas * double complex a, z, w;
287b36286aSmartynas *
297b36286aSmartynas * w = cpow (a, z);
307b36286aSmartynas *
317b36286aSmartynas *
327b36286aSmartynas *
337b36286aSmartynas * DESCRIPTION:
347b36286aSmartynas *
357b36286aSmartynas * Raises complex A to the complex Zth power.
367b36286aSmartynas * Definition is per AMS55 # 4.2.8,
377b36286aSmartynas * analytically equivalent to cpow(a,z) = cexp(z clog(a)).
387b36286aSmartynas *
397b36286aSmartynas * ACCURACY:
407b36286aSmartynas *
417b36286aSmartynas * Relative error:
427b36286aSmartynas * arithmetic domain # trials peak rms
437b36286aSmartynas * IEEE -10,+10 30000 9.4e-15 1.5e-15
447b36286aSmartynas *
457b36286aSmartynas */
467b36286aSmartynas
477b36286aSmartynas #include <complex.h>
48de3697aaSmartynas #include <float.h>
497b36286aSmartynas #include <math.h>
507b36286aSmartynas
517b36286aSmartynas double complex
cpow(double complex a,double complex z)527b36286aSmartynas cpow(double complex a, double complex z)
537b36286aSmartynas {
547b36286aSmartynas double complex w;
557b36286aSmartynas double x, y, r, theta, absa, arga;
567b36286aSmartynas
577b36286aSmartynas x = creal (z);
587b36286aSmartynas y = cimag (z);
597b36286aSmartynas absa = cabs (a);
607b36286aSmartynas if (absa == 0.0) {
617b36286aSmartynas return (0.0 + 0.0 * I);
627b36286aSmartynas }
637b36286aSmartynas arga = carg (a);
647b36286aSmartynas r = pow (absa, x);
657b36286aSmartynas theta = x * arga;
667b36286aSmartynas if (y != 0.0) {
677b36286aSmartynas r = r * exp (-y * arga);
687b36286aSmartynas theta = theta + y * log (absa);
697b36286aSmartynas }
707b36286aSmartynas w = r * cos (theta) + (r * sin (theta)) * I;
717b36286aSmartynas return (w);
727b36286aSmartynas }
73*2f2c0062Sguenther DEF_STD(cpow);
74*2f2c0062Sguenther LDBL_MAYBE_UNUSED_CLONE(cpow);
75