svn commit: r251836 - user/peterj
Peter Jeremy
peterj at FreeBSD.org
Mon Jun 17 08:53:56 UTC 2013
Author: peterj
Date: Mon Jun 17 08:53:55 2013
New Revision: 251836
URL: http://svnweb.freebsd.org/changeset/base/251836
Log:
Initial work-in-progress commit of a cpow(3) implementation.
This implements a subset of the desired special cases and then punts
to the fallback: cpow(z, w) => cexp(w * clog(z))
Added:
user/peterj/cpow.c (contents, props changed)
Modified:
user/peterj/README
Modified: user/peterj/README
==============================================================================
--- user/peterj/README Mon Jun 17 08:49:08 2013 (r251835)
+++ user/peterj/README Mon Jun 17 08:53:55 2013 (r251836)
@@ -1,6 +1,6 @@
My odds & ends
+cpow.c Complex pow() implementation. This is a work-in-progress.
+
ctest.c test libm exception conditions listed in WG14/N1256 G.6
against a set of float, double and long double functions
-
-$FreeBSD$
Added: user/peterj/cpow.c
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ user/peterj/cpow.c Mon Jun 17 08:53:55 2013 (r251836)
@@ -0,0 +1,240 @@
+/*-
+ * Copyright (c) 2012 Peter Jeremy <peterj at FreeBSD.ORG>
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <complex.h>
+#include <math.h>
+#include "math_private.h"
+
+/*
+ * Raise a complex number to a complex power.
+ *
+ * Algorithmetically:
+ * cpow(z, w) = z ** w = cexp(w * clog(z))
+ * cpow(rz + I*iz, rw + I*iw) = cexp((rw + I*iw) * clog(rz + I*iz))
+ * Let: rad = hypot(rz, iz)
+ * th = atan2(iz, rz)
+ * clog(rz + I*iz) = log(hypot(rz, iz)) + I*atan2(iz, rz)
+ * = log(rad) + I*th
+ * then:
+ * cpow(rz + I*iz, rw + I*iw) = cexp((rw + I*iw) * (log(rad) + I*th))
+ * = cexp((rw*log(rad) - iw*th) + I*(rw*th + iw*log(rad)))
+ * = exp(rw*log(rad)) / exp(iw*th) * cis(rw*th + iw*log(rad))
+ * = pow(rad, rw) / exp(iw*th) * cis(rw*th + iw*log(rad))
+ *
+ * cexp(x + I*y) = exp(x) * cis(y)
+ * clog(x + I*y) = log(hypot(x, y)) + I*atan2(y, x)
+ *
+ * Special cases:
+ * iw = 0:
+ * = pow(rad, rw) * cis(rw*th)
+ * rw = 0
+ *
+ * There are no special exceptions defined for cpow(), instead implementations
+ * raise exceptions as appropriate for the calculations. WG14/N1256 explicitly
+ * allows implementations to return cexp(y * clog(x)) but this loses precision
+ * roughly equal to the number of bits in the floating point type's exponent.
+ *
+ * Special cases & exceptions:
+ * 'int' integer
+ * 'rat' rational
+ * 'odd' odd integer
+ * 'ANY' includes NaN
+ * z.r z.i w.r w.i
+ * ±0 0 odd<0 0 ±Inf + I*0, Div0 [F.9.4.4]
+ * ±0 0 <0,!odd 0 +Inf + I*0, Div0 [F.9.4.4]
+ * ±0 0 odd>0 0 ±0 + I*0 [F.9.4.4]
+ * ±0 0 >0,!odd 0 +0 + I*0 [F.9.4.4]
+ * -1 0 ±Inf 0 +1 + I*0 [F.9.4.4] #
+ * +1 0 ANY ANY +1 + I*0 [F.9.4.4] #
+ * ANY ANY ±0 0 +1 + I*0 [F.9.4.4]+
+ * <0 0 !int 0 NaN + I*0, Invalid [F.9.4.4]*
+ * |z|<1 -Inf 0 +Inf + I*0 [F.9.4.4]+ #
+ * |z|>1 -Inf 0 +0 + I*0 [F.9.4.4]+ #
+ * |z|<1 +Inf 0 +0 + I*0 [F.9.4.4]+ #
+ * |z|>1 +Inf 0 +Inf + I*0 [F.9.4.4]+ #
+ * -Inf 0 odd<0 0 -0 + I*0 [F.9.4.4]
+ * -Inf 0 <0,!odd 0 +0 + I*0 [F.9.4.4]
+ * -Inf 0 odd>0 0 -Inf + I*0 [F.9.4.4]
+ * -Inf 0 >0,!odd 0 +Inf + I*0 [F.9.4.4]
+ * +Inf 0 <0 0 +0 + I*0 [F.9.4.4]
+ * +Inf 0 >0 0 +Inf + I*0 [F.9.4.4]
+ * ? ? x Inf NaN + I*NaN, Invalid [G.6.3.1]
+ * ? ? x NaN NaN + I*NaN, optInvalid [G.6.3.1]
+ *
+ * Special cases (evaluated in order):
+ * 0. +1 ** (anything) is 1
+ * 1. (anything) ** 0 is 1
+ * 2. (anything) ** 1 is itself
+ * 3. (anything) ** NAN is NAN
+ * 4. NAN ** (anything except 0) is NAN
+ * 5. (|x| > 1) ** +INF is +INF
+ * 6. (|x| > 1) ** -INF is +0
+ * 7. (|x| < 1) ** +INF is +0
+ * 8. (|x| < 1) ** -INF is +INF
+ * 9. (|x| == 1) ** +-INF is NAN
+#* 10. +0 ** (+anything except 0, NAN) is +0
+#* 11. -0 ** (+anything except 0, NAN, odd integer) is +0
+#* 12. +0 ** (-anything except 0, NAN) is +INF
+#* 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
+#* 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
+#* 15. +INF ** (+anything except 0,NAN) is +INF
+#* 16. +INF ** (-anything except 0,NAN) is +0
+#* 17. -INF ** (anything) = -0 ** (-anything)
+#* 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
+#* 19. (-anything except 0 and inf) ** (non-integer) is NAN
+^ '#' means it's not implemented yet
+ *
+ */
+double complex
+cpow(double complex z, double complex w)
+{
+ double iw, iz, rw, rz; /* Imaginary & real parts of parameters */
+ int32_t hiw, hiz, hrw, hrz; /* high words of imag & real parts */
+ int32_t liw, liz, lrw, lrz; /* low words of imag & real parts */
+ int32_t miw, miz, mrw, mrz; /* masked high words */
+ double rad, th; /* z in polar form */
+ double t1;
+ int32_t hrad, lrad; /* breakdown of rad */
+ int32_t j, k, wisint;
+
+ rw = creal(w);
+ iw = cimag(w);
+ EXTRACT_WORDS(hrw, lrw, rw);
+ EXTRACT_WORDS(hiw, liw, iw);
+ mrw = 0x7fffffff & hrw;
+ miw = 0x7fffffff & hiw;
+
+ rz = creal(z);
+ iz = cimag(z);
+ EXTRACT_WORDS(hrz, lrz, rz);
+ EXTRACT_WORDS(hiz, liz, iz);
+ mrz = 0x7fffffff & hrz;
+ miz = 0x7fffffff & hiz;
+
+ /* +1 ±I*0 ** ANY is +1 +I*0 */
+ if (((hrz - 0x3ff00000) | lrz | miz | liz) == 0)
+ return(cpack(1.0, 0.0));
+
+ if ((miw | liw) == 0) { /* w is real */
+ /* (anything) ** 0 is +1 */
+ if ((mrw | lrw) == 0)
+ return (cpack(1.0, 0.0));
+
+ /* (anything) ** +1 is itself */
+ if (((hrw - 0x3ff00000) | lrw) == 0)
+ return (z);
+
+ /* Use csqrt() for (anything) ** +0.5 */
+ if (((hrw - 0x3fe00000) | lrw) == 0)
+ return (csqrt(z));
+
+ /* Use pow() for (positive real) ** real */
+ if (hrz >= 0 && (miz | liz) == 0)
+ return (cpack(pow(rz, rw), 0.0));
+
+ rad = cabs(z); /*### can be Inf even though z is finite */
+ EXTRACT_WORDS(hrad, lrad, rad);
+
+ /* NAN ** (anything except 0) is NAN */
+ if (hrad > 0x7ff00000 || (hrad == 0x7ff00000 && lrad != 0))
+ return (cpack(rad, rad));
+
+ if (((mrw - 0x7ff00000) | lrw) == 0) { /* Inf */
+ /* (Anything != +1 on the unit circle) ** Inf is NaN */
+ if (((hrad - 0x3ff00000) | lrad) == 0) /* 1 */
+ return (cpack(NAN, NAN));
+
+ /*
+ * (Anything inside the unit circle) ** -Inf is +Inf.
+ * (Anything outside the unit circle) ** +Inf is +Inf.
+ */
+ if ((hrad < 0x3ff00000) == (hrw < 0))
+ return (cpack(INFINITY, 0.0));
+
+ /*
+ * Anything inside the unit circle ** +Inf is +0.
+ * Anything outside the unit circle ** -Inf is +0.
+ */
+ return (cpack(0.0, 0.0));
+ }
+
+ /* (anything) ** NaN is NaN */
+ if (mrw >= 0x7ff00000)
+ return (cpack(rz + rw, iz + rw));
+
+ /*
+ * At this point, the exponent is a real, finite number
+ * and the base is a, possibly infinite, number, excluding
+ * positive reals.
+ */
+
+ /* Determine if the exponent is an odd integer:
+ * wisint = 0 ... w is not an integer
+ * wisint = 1 ... w is an odd int
+ * wisint = 2 ... w is an even int
+ */
+ wisint = 0;
+ if (mrw >= 0x43400000) /* must be even integer */
+ wisint = 2;
+ else if (mrw >= 0x3ff00000) {
+ k = (mrw >> 20) - 0x3ff; /* exponent */
+ if (k > 20) {
+ j = lrw >> (52 - k);
+ if ((j << (52 - k)) == lrw)
+ wisint = 2 - (j & 1);
+ } else if (lrw == 0) {
+ j = mrw >> (20 - k);
+ if ((j << (20 - k)) == mrw)
+ wisint = 2 - (j & 1);
+ }
+ }
+
+ /*###### Special case w is a small integer */
+
+ /*######*/
+ rad = pow(rad, rw);
+ th = atan2(iz, rz) * rw;
+ return (cpack(rad * cos(th), rad * sin(th)));
+ }
+
+ /* w is not pure real */
+
+ rad = cabs(z); /*### can be Inf even though z is finite */
+ EXTRACT_WORDS(hrad, lrad, rad);
+
+ /* NAN ** (anything except 0) is NAN */
+ if (hrad > 0x7ff00000 || (hrad == 0x7ff00000 && lrad == 0))
+ return (cpack(rad, rad));
+
+ /*######*/
+ th = atan2(iz, rz);
+ t1 = rw * th + iw * log(rad);
+ rad = pow(rad, rw) / exp(iw * th);
+ return (cpack(rad * cos(t1), rad * sin(t1)));
+}
More information about the svn-src-user
mailing list