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