Date: Mon, 17 Jun 2013 08:53:55 +0000 (UTC) From: Peter Jeremy <peterj@FreeBSD.org> To: src-committers@freebsd.org, svn-src-user@freebsd.org Subject: svn commit: r251836 - user/peterj Message-ID: <201306170853.r5H8rtWp003891@svn.freebsd.org>
next in thread | raw e-mail | index | archive | help
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@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))); +}
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?201306170853.r5H8rtWp003891>