+2015-06-13 Keith Marshall <keithmarshall@users.sourceforge.net>
+
+ Prefer our hypot() implementation in complex maths operations.
+
+ * mingwex/complex/cabs_generic.c: New file; it replaces...
+ * mingwex/complex/cabs.c mingwex/complex/cabsf.c: ...all of...
+ * mingwex/complex/cabsl.c: ...these.
+
+ * mingwex/complex/catan_generic.c: New file; it replaces...
+ * mingwex/complex/catan.c mingwex/complex/catanf.c: ...all of...
+ * mingwex/complex/catanl.c: ...these.
+
+ * mingwex/complex/clog_generic.c: New file; it replaces...
+ * mingwex/complex/clog.c mingwex/complex/clogf.c: ...all of...
+ * mingwex/complex/clogl.c: ...these.
+
+ * mingwex/complex/cpow_generic.c: New file; it replaces...
+ * mingwex/complex/cpow.c mingwex/complex/cpowf.c: ...all of...
+ * mingwex/complex/cpowl.c: ...these.
+
+ * mingwex/complex/csqrt_generic.c: Do not use...
+ (_hypot): ...this Microsoft form of function reference; use...
+ (hypot): ...this ANSI standard form instead.
+
2015-06-10 Keith Marshall <keithmarshall@users.sourceforge.net>
Correct C++ compilation anomaly with hypotf() in cmath header.
+++ /dev/null
-#include <math.h>
-#include <complex.h>
-
-double cabs (double complex Z)
-{
- return _hypot ( __real__ Z, __imag__ Z);
-}
--- /dev/null
+/*
+ * cabs_generic.c
+ *
+ * $Id$
+ *
+ * Compute the modulus of a complex number; this provides a generic
+ * implementation for the cabs(), cabsf(), and cabsl() functions.
+ *
+ * Written by Keith Marshall <keithmarshall@users.sourceforge.net>
+ * This is an adaptation of original contributions by Danny Smith
+ * Copyright (C) 2003, 2015, MinGW.org Project
+ *
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice, this permission notice, and the following
+ * disclaimer shall be included in all copies or substantial portions of
+ * the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OF OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *
+ *
+ * This is a generic implementation for all of the csqrt(), csqrtl(),
+ * and csqrth() functions; each is to be compiled separately, i.e.
+ *
+ * gcc -D FUNCTION=cabs -o cabs.o cabs_generic.c
+ * gcc -D FUNCTION=cabsl -o cabsl.o cabs_generic.c
+ * gcc -D FUNCTION=cabsf -o cabsf.o cabs_generic.c
+ *
+ */
+#include <math.h>
+#include <complex.h>
+
+#ifndef FUNCTION
+/* If user neglected to specify it, the default compilation is for
+ * the cabs() function.
+ */
+# define FUNCTION cabs
+#endif
+
+#define argtype_cabs double
+#define argtype_cabsl long double
+#define argtype_cabsf float
+
+#define PASTE(PREFIX,SUFFIX) PREFIX##SUFFIX
+#define mapfunc(NAME) mapname(mapfunc_,FUNCTION)(NAME)
+#define mapname(PREFIX,SUFFIX) PASTE(PREFIX,SUFFIX)
+
+#define ARGTYPE(FUNCTION) PASTE(argtype_,FUNCTION)
+
+#define mapfunc_cabs(NAME) NAME
+#define mapfunc_cabsl(NAME) PASTE(NAME,l)
+#define mapfunc_cabsf(NAME) PASTE(NAME,f)
+
+/* Define the generic function implementation.
+ */
+ARGTYPE(FUNCTION) FUNCTION( ARGTYPE(FUNCTION) complex Z )
+{
+ return mapfunc(hypot)( __real__ Z, __imag__ Z );
+}
+
+/* $RCSfile$: end of file */
+++ /dev/null
-#include <math.h>
-#include <complex.h>
-
-float cabsf (float complex Z)
-{
- return (float) _hypot ( __real__ Z, __imag__ Z);
-}
+++ /dev/null
-#include <math.h>
-#include <complex.h>
-
-long double cabsl (long double complex Z)
-{
- return hypotl ( __real__ Z, __imag__ Z);
-}
+++ /dev/null
-/* catan.c */
-
-/*
- Contributed by Danny Smith
- 2003-10-17
-
- FIXME: This needs some serious numerical analysis.
-*/
-
-#include <math.h>
-#include <complex.h>
-#include <errno.h>
-
-/* catan (z) = -I/2 * clog ((I + z) / (I - z)) */
-
-double complex
-catan (double complex Z)
-{
- double complex Res;
- double complex Tmp;
- double x = __real__ Z;
- double y = __imag__ Z;
-
- if ( x == 0.0 && (1.0 - fabs (y)) == 0.0)
- {
- errno = ERANGE;
- __real__ Res = HUGE_VAL;
- __imag__ Res = HUGE_VAL;
- }
- else if (isinf (_hypot (x, y)))
- {
- __real__ Res = (x > 0 ? M_PI_2 : -M_PI_2);
- __imag__ Res = 0.0;
- }
- else
- {
- __real__ Tmp = - x;
- __imag__ Tmp = 1.0 - y;
-
- __real__ Res = x;
- __imag__ Res = y + 1.0;
-
- Tmp = clog (Res/Tmp);
- __real__ Res = - 0.5 * __imag__ Tmp;
- __imag__ Res = 0.5 * __real__ Tmp;
- }
-
- return Res;
-}
--- /dev/null
+/*
+ * catan_generic.c
+ *
+ * $Id$
+ *
+ * Compute the complex arctan corresponding to a complex tangent value;
+ * a generic implementation for catan(), catanf(), and catanl().
+ *
+ * Written by Keith Marshall <keithmarshall@users.sourceforge.net>
+ * This is an adaptation of an original contribution by Danny Smith
+ * Copyright (C) 2003-2005, 2015, MinGW.org Project
+ *
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice, this permission notice, and the following
+ * disclaimer shall be included in all copies or substantial portions of
+ * the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OF OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *
+ *
+ * This is a generic implementation for all of the catan(), catanl(),
+ * and catanh() functions; each is to be compiled separately, i.e.
+ *
+ * gcc -D FUNCTION=catan -o catan.o catan_generic.c
+ * gcc -D FUNCTION=catanl -o catanl.o catan_generic.c
+ * gcc -D FUNCTION=catanf -o catanf.o catan_generic.c
+ *
+ */
+#include <math.h>
+#include <complex.h>
+
+#ifndef FUNCTION
+/* If user neglected to specify it, the default compilation is for
+ * the catan() function.
+ */
+# define FUNCTION catan
+#endif
+
+#define argtype_catan double
+#define argtype_catanl long double
+#define argtype_catanf float
+
+#define PASTE(PREFIX,SUFFIX) PREFIX##SUFFIX
+#define mapname(PREFIX,SUFFIX) PASTE(PREFIX,SUFFIX)
+
+#define ARGTYPE(FUNCTION) PASTE(argtype_,FUNCTION)
+#define ARGCAST(VALUE) mapname(argcast_,FUNCTION)(VALUE)
+
+#define argcast_catan(VALUE) VALUE
+#define argcast_catanl(VALUE) PASTE(VALUE,L)
+#define argcast_catanf(VALUE) PASTE(VALUE,F)
+
+#define mapfunc(NAME) mapname(mapfunc_,FUNCTION)(NAME)
+
+#define mapfunc_catan(NAME) NAME
+#define mapfunc_catanl(NAME) PASTE(NAME,l)
+#define mapfunc_catanf(NAME) PASTE(NAME,f)
+
+/* Define the generic function implementation.
+ */
+ARGTYPE(FUNCTION) complex FUNCTION( ARGTYPE(FUNCTION) complex Z )
+{
+ ARGTYPE(FUNCTION) x = __real__ Z; /* real part of Z saved as x */
+ ARGTYPE(FUNCTION) y = __imag__ Z; /* imaginary part of Z saved as iy */
+
+ /* Having thus saved the original value of Z as separate real and
+ * imaginary parts, we may reuse Z to accumulate the arctan value,
+ * noting that...
+ */
+ if( isfinite( mapfunc(hypot)( x, y )) )
+ {
+ /* ...for any complex number with finite modulus, its argtan may
+ * be computed as:
+ *
+ * catan(Z) = (clog(1.0 + iZ) - clog(1.0 - iZ)) / 2.0i
+ *
+ * For computational convenience, we may introduce an additional
+ * intermediate accumulator variable...
+ */
+ ARGTYPE(FUNCTION) complex tmp;
+
+ /* ...in which we initially compute the (1.0 - iZ) term...
+ */
+ __real__ tmp = ARGCAST(1.0) + y;
+ __imag__ tmp = -x;
+
+ /* ...while computing (1.0 + iZ) directly in Z itself.
+ */
+ __real__ Z = ARGCAST(1.0) - y;
+ __imag__ Z = x;
+
+ /* We may then compute and combine the complex logarithms
+ * of this pair of sub-expressions, while we simultaneously
+ * perform the REAL division by 2.0...
+ */
+ tmp = ARGCAST(0.5) * (mapfunc(clog)( Z ) - mapfunc(clog)( tmp ));
+
+ /* ...and finally, the division by i has the effect of
+ * exchanging the real and imaginary parts of the result,
+ * with a change of sign in the resultant imaginary part.
+ */
+ __real__ Z = __imag__ tmp;
+ __imag__ Z = -(__real__ tmp);
+ }
+
+ else
+ { /* Conversely, if the modulus of Z is in any way undefined,
+ * we return a real value, with modulus of PI / 2.0, and with
+ * the same sign as the original real part of Z.
+ */
+ __real__ Z = mapfunc(copysign)( ARGCAST(M_PI_2), x );
+ __imag__ Z = ARGCAST(0.0);
+ }
+
+ /* In either case, Z now represents the arctan of its original
+ * value, which we are required to return.
+ */
+ return Z;
+}
+
+/* $RCSfile$: end of file */
+++ /dev/null
-/* catanf.c */
-
-/*
- Contributed by Danny Smith
- 2004-12-24
-
- FIXME: This needs some serious numerical analysis.
-*/
-
-#include <math.h>
-#include <complex.h>
-#include <errno.h>
-
-/* catan (z) = -I/2 * clog ((I + z) / (I - z)) */
-
-float complex
-catanf (float complex Z)
-{
- float complex Res;
- float complex Tmp;
- float x = __real__ Z;
- float y = __imag__ Z;
-
- if ( x == 0.0f && (1.0f - fabsf (y)) == 0.0f)
- {
- errno = ERANGE;
- __real__ Res = HUGE_VALF;
- __imag__ Res = HUGE_VALF;
- }
- else if (isinf (hypotf (x, y)))
- {
- __real__ Res = (x > 0 ? M_PI_2 : -M_PI_2);
- __imag__ Res = 0.0f;
- }
- else
- {
- __real__ Tmp = - x;
- __imag__ Tmp = 1.0f - y;
-
- __real__ Res = x;
- __imag__ Res = y + 1.0f;
-
- Tmp = clogf (Res/Tmp);
- __real__ Res = - 0.5f * __imag__ Tmp;
- __imag__ Res = 0.5f * __real__ Tmp;
- }
-
- return Res;
-}
+++ /dev/null
-/* catanl.c */
-
-/*
- Contributed by Danny Smith
- 2005-01-04
-
- FIXME: This needs some serious numerical analysis.
-*/
-
-#include <math.h>
-#include <complex.h>
-#include <errno.h>
-
-/* catan (z) = -I/2 * clog ((I + z) / (I - z)) */
-
-#ifndef _M_PI_2L
-#define _M_PI_2L 1.5707963267948966192313L
-#endif
-
-long double complex
-catanl (long double complex Z)
-{
- long double complex Res;
- long double complex Tmp;
- long double x = __real__ Z;
- long double y = __imag__ Z;
-
- if ( x == 0.0L && (1.0L - fabsl (y)) == 0.0L)
- {
- errno = ERANGE;
- __real__ Res = HUGE_VALL;
- __imag__ Res = HUGE_VALL;
- }
- else if (isinf (hypotl (x, y)))
- {
- __real__ Res = (x > 0 ? _M_PI_2L : -_M_PI_2L);
- __imag__ Res = 0.0L;
- }
- else
- {
- __real__ Tmp = - x;
- __imag__ Tmp = 1.0L - y;
-
- __real__ Res = x;
- __imag__ Res = y + 1.0L;
-
- Tmp = clogl (Res/Tmp);
- __real__ Res = - 0.5L * __imag__ Tmp;
- __imag__ Res = 0.5L * __real__ Tmp;
- }
-
- return Res;
-}
+++ /dev/null
-/*
- clog.c
- Contributed by Danny Smith
- 2003-10-20
-*/
-
-/* clog (x + I * y) = log (hypot (x, y)) + I * atan2 (y, x) */
-
-#include <math.h>
-#include <complex.h>
-
-double complex clog (double complex Z)
-{
- double complex Res;
- __real__ Res = log (_hypot (__real__ Z, __imag__ Z));
- __imag__ Res = carg (Z);
- return Res;
-}
-
--- /dev/null
+/*
+ * clog_generic.c
+ *
+ * $Id$
+ *
+ * Compute the logarithm of a complex number; this provides a generic
+ * implementation for the clog(), clogf(), and clogl() functions.
+ *
+ * Written by Keith Marshall <keithmarshall@users.sourceforge.net>
+ * This is an adaptation of original contributions by Danny Smith
+ * Copyright (C) 2003-2005, 2015, MinGW.org Project
+ *
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice, this permission notice, and the following
+ * disclaimer shall be included in all copies or substantial portions of
+ * the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OF OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *
+ *
+ * This is a generic implementation for all of the clog(), clogl(),
+ * and clogf() functions; each is to be compiled separately, i.e.
+ *
+ * gcc -D FUNCTION=clog -o clog.o clog_generic.c
+ * gcc -D FUNCTION=clogl -o clogl.o clog_generic.c
+ * gcc -D FUNCTION=clogf -o clogf.o clog_generic.c
+ *
+ */
+#include <math.h>
+#include <complex.h>
+
+#ifndef FUNCTION
+/* If user neglected to specify it, the default compilation is for
+ * the clog() function.
+ */
+# define FUNCTION clog
+#endif
+
+#define argtype_clog double
+#define argtype_clogl long double
+#define argtype_clogf float
+
+#define PASTE(PREFIX,SUFFIX) PREFIX##SUFFIX
+#define mapname(PREFIX,SUFFIX) PASTE(PREFIX,SUFFIX)
+
+#define ARGTYPE(FUNCTION) PASTE(argtype_,FUNCTION)
+#define ARGCAST(VALUE) mapname(argcast_,FUNCTION)(VALUE)
+
+#define argcast_clog(VALUE) VALUE
+#define argcast_clogl(VALUE) PASTE(VALUE,L)
+#define argcast_clogf(VALUE) PASTE(VALUE,F)
+
+#define mapfunc(NAME) mapname(mapfunc_,FUNCTION)(NAME)
+
+#define mapfunc_clog(NAME) NAME
+#define mapfunc_clogl(NAME) PASTE(NAME,l)
+#define mapfunc_clogf(NAME) PASTE(NAME,f)
+
+/* Prefer the fast version of the log() function.
+ */
+#include "../math/fastmath.h"
+#define log __fast_log
+
+/* Define the generic function implementation.
+ */
+ARGTYPE(FUNCTION) complex FUNCTION( ARGTYPE(FUNCTION) complex Z )
+{
+ /* Compute the value to return; real part is the logarithm of the
+ * modulus of Z, (equivalent to the length of the hypotenuse of the
+ * triangle formed when representing Z by cartesian co-ordinates in
+ * the complex number plane); imaginary part is its phase angle,
+ * (a.k.a. its argument, in the polar representation).
+ *
+ * Note that the applicable POSIX and ISO-C standards do not require
+ * any specific error return, and none is provided; (in particular,
+ * whereas the modulus of Z is ALWAYS >= 0.0, and thus there is no
+ * scope for a log() domain error to occur, modulus of Z == 0.0
+ * causes a log() range error for which errno is NOT set).
+ */
+ ARGTYPE(FUNCTION) complex result;
+ __real__ result = mapfunc(log)( mapfunc(hypot)( __real__ Z, __imag__ Z ) );
+ __imag__ result = mapfunc(carg)( Z );
+ return result;
+}
+
+/* $RCSfile$: end of file */
+++ /dev/null
-/*
- clogf.c
- Contributed by Danny Smith
- 2004-12-24
-*/
-
-/* clog (x + I * y) = log (hypot (x, y)) + I * atan2 (y, x) */
-
-#include <math.h>
-#include <complex.h>
-
-float complex clogf (float complex Z)
-{
- float complex Res;
- __real__ Res = logf (_hypot (__real__ Z, __imag__ Z));
- __imag__ Res = cargf (Z);
- return Res;
-}
-
+++ /dev/null
-/*
- clogl.c
- Contributed by Danny Smith
- 2005-01-04
-*/
-
-/* clog (x + I * y) = log (hypot (x, y)) + I * atan2 (y, x) */
-
-#include <math.h>
-#include <complex.h>
-
-long double complex clogl (long double complex Z)
-{
- long double complex Res;
- __real__ Res = logl (hypotl (__real__ Z, __imag__ Z));
- __imag__ Res = cargl (Z);
- return Res;
-}
-
+++ /dev/null
-/* cpow.c */
-/*
- Contributed by Danny Smith
- 2003-10-20
-*/
-
-/* cpow(X, Y) = cexp(X * clog(Y)) */
-
-#include <math.h>
-#include <complex.h>
-
-/* Use dll version of pow */
-extern double (*_imp__pow) (double, double);
-#define pow (*_imp__pow)
-
-double complex cpow (double complex X, double complex Y)
-{
- double complex Res;
- double i;
- double r = hypot (__real__ X, __imag__ X);
- if (r == 0.0)
- {
- __real__ Res = __imag__ Res = 0.0;
- }
- else
- {
- double rho;
- double theta;
- i = carg (X);
- theta = i * __real__ Y;
-
- if (__imag__ Y == 0.0)
- /* This gives slightly more accurate results in these cases. */
- rho = pow (r, __real__ Y);
- else
- {
- r = log (r);
- /* rearrangement of cexp(X * clog(Y)) */
- theta += r * __imag__ Y;
- rho = exp (r * __real__ Y - i * __imag__ Y);
- }
-
- __real__ Res = rho * cos (theta);
- __imag__ Res = rho * sin (theta);
- }
- return Res;
-}
-
--- /dev/null
+/*
+ * cpow_generic.c
+ *
+ * $Id$
+ *
+ * Compute the result of raising a complex number, Z, to a complex
+ * power, N; this provides a generic implementation for each of the
+ * cpow(), cpowf(), and cpowl() functions.
+ *
+ * Written by Keith Marshall <keithmarshall@users.sourceforge.net>
+ * This is an adaptation of an original contribution by Danny Smith
+ * Copyright (C) 2003, 2014, MinGW.org Project
+ *
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice, this permission notice, and the following
+ * disclaimer shall be included in all copies or substantial portions of
+ * the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OF OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *
+ *
+ * This is a generic implementation for all of the cpow(), cpowl(),
+ * and cpowh() functions; each is to be compiled separately, i.e.
+ *
+ * gcc -D FUNCTION=cpow -o cpow.o cpow_generic.c
+ * gcc -D FUNCTION=cpowl -o cpowl.o cpow_generic.c
+ * gcc -D FUNCTION=cpowf -o cpowf.o cpow_generic.c
+ *
+ */
+#include <math.h>
+#include <complex.h>
+
+#ifndef FUNCTION
+/* If user neglected to specify it, the default compilation is for
+ * the cpow() function.
+ */
+# define FUNCTION cpow
+#endif
+
+#define argtype_cpow double
+#define argtype_cpowl long double
+#define argtype_cpowf float
+
+#define PASTE(PREFIX,SUFFIX) PREFIX##SUFFIX
+#define mapname(PREFIX,SUFFIX) PASTE(PREFIX,SUFFIX)
+
+#define ARGTYPE(FUNCTION) PASTE(argtype_,FUNCTION)
+#define ARGCAST(VALUE) mapname(argcast_,FUNCTION)(VALUE)
+
+#define argcast_cpow(VALUE) VALUE
+#define argcast_cpowl(VALUE) PASTE(VALUE,L)
+#define argcast_cpowf(VALUE) PASTE(VALUE,F)
+
+#define mapfunc(NAME) mapname(mapfunc_,FUNCTION)(NAME)
+
+#define mapfunc_cpow(NAME) NAME
+#define mapfunc_cpowl(NAME) PASTE(NAME,l)
+#define mapfunc_cpowf(NAME) PASTE(NAME,f)
+
+/* Danny Smith's original implementation explicitly referenced the
+ * DLL entry point for the pow() function, in the case of the cpow()
+ * function only; the following variation on mapfunc() reproduces
+ * this implementation model.
+ */
+extern double (*_imp__pow)( double, double );
+#define powfunc mapname(powfunc_,FUNCTION)
+#define powfunc_cpow (*_imp__pow)
+#define powfunc_cpowf powf
+#define powfunc_cpowl powl
+
+/* Define the generic function implementation.
+ */
+ARGTYPE(FUNCTION) complex FUNCTION
+( ARGTYPE(FUNCTION) complex Z, ARGTYPE(FUNCTION) complex N )
+{
+ ARGTYPE(FUNCTION) radius;
+ if( (radius = mapfunc(hypot)( __real__ Z, __imag__ Z )) == ARGCAST(0.0) )
+ /*
+ * Modulus of Z is zero, hence result is also zero; force it.
+ */
+ __real__ Z = __imag__ Z = ARGCAST(0.0);
+
+ else
+ { /* Compute the complex power function result, in terms of its
+ * polar (r, theta) representation in the complex plane.
+ */
+ ARGTYPE(FUNCTION) r = mapfunc(carg)( Z );
+ ARGTYPE(FUNCTION) theta = r * __real__ N;
+
+ if( __imag__ N == ARGCAST(0.0) )
+ /*
+ * For entirely real N, the power function tends to yield
+ * more accuracy than exponential expansion.
+ */
+ r = powfunc( radius, __real__ N );
+
+ else
+ {
+ theta += (radius = mapfunc(log)( radius )) * __imag__ N;
+ r = mapfunc(exp)( radius * __real__ N - r * __imag__ N );
+ }
+
+ /* Convert polar to cartesian representation for return.
+ */
+ __real__ Z = r * mapfunc(cos)( theta );
+ __imag__ Z = r * mapfunc(sin)( theta );
+ }
+
+ /* Regardless of how we got here, the original base value has
+ * now been raised to the specified power; return it.
+ */
+ return Z;
+}
+
+/* $RCSfile$: end of file */
+++ /dev/null
-/* cpowf.c */
-/*
- Contributed by Danny Smith
- 2004-12-24
-*/
-
-/* cpow(X, Y) = cexp(X * clog(Y)) */
-
-#include <math.h>
-#include <complex.h>
-
-float complex cpowf (float complex X, float complex Y)
-{
- float complex Res;
- float i;
- float r = _hypot (__real__ X, __imag__ X);
- if (r == 0.0f)
- {
- __real__ Res = __imag__ Res = 0.0;
- }
- else
- {
- float rho;
- float theta;
- i = cargf (X);
- theta = i * __real__ Y;
-
- if (__imag__ Y == 0.0f)
- /* This gives slightly more accurate results in these cases. */
- rho = powf (r, __real__ Y);
- else
- {
- r = logf (r);
- /* rearrangement of cexp(X * clog(Y)) */
- theta += r * __imag__ Y;
- rho = expf (r * __real__ Y - i * __imag__ Y);
- }
-
- __real__ Res = rho * cosf (theta);
- __imag__ Res = rho * sinf (theta);
- }
- return Res;
-}
+++ /dev/null
-/* cpowl.c */
-/*
- Contributed by Danny Smith
- 2005-01-04
-*/
-
-/* cpow(X, Y) = cexp(X * clog(Y)) */
-
-#include <math.h>
-#include <complex.h>
-
-long double complex cpowl (long double complex X, long double complex Y)
-{
- long double complex Res;
- long double i;
- long double r = hypotl (__real__ X, __imag__ X);
- if (r == 0.0L)
- {
- __real__ Res = __imag__ Res = 0.0L;
- }
- else
- {
- long double rho;
- long double theta;
- i = cargl (X);
- theta = i * __real__ Y;
-
- if (__imag__ Y == 0.0L)
- /* This gives slightly more accurate results in these cases. */
- rho = powl (r, __real__ Y);
- else
- {
- r = logl (r);
- /* rearrangement of cexp(X * clog(Y)) */
- theta += r * __imag__ Y;
- rho = expl (r * __real__ Y - i * __imag__ Y);
- }
-
- __real__ Res = rho * cosl (theta);
- __imag__ Res = rho * sinl (theta);
- }
- return Res;
-}
*
* $Id$
*
- * Compute the complex arcsin corresponding to a complex sine value;
- * a generic implementation for casin(), casinf(), and casinl().
+ * Compute the principal square root of a complex number; this provides
+ * a generic implementation for csqrt(), csqrtf(), and csqrtl().
*
* Written by Keith Marshall <keithmarshall@users.sourceforge.net>
* This is an adaptation of an original contribution by Danny Smith
- * Copyright (C) 2003, 2014, MinGW.org Project
+ * Copyright (C) 2003, 2014, 2015, MinGW.org Project
*
*
* Permission is hereby granted, free of charge, to any person obtaining a
#define log1p __fast_log1p
#define sqrt __fast_sqrt
-/* To compute the modulus of complex numbers, we use the hypot()
- * function, representing it generically by reference to MSVCRT.DLL's
- * double _hypot(); for the float and long double variants, we use
- * our own implementations, which have no leading underscore.
- */
-#define _hypotf hypotf
-#define _hypotl hypotl
-
/* Define the generic function implementation.
*/
ARGTYPE(FUNCTION) complex FUNCTION( ARGTYPE(FUNCTION) complex Z )
}
else
{ /* Finally, this represents the general complex number case; we
- * may apply (12) and (13) directly, (noting that _hypot(x, y) is
- * functionally equivalent to sqrt(y^2 + x^2), with the advantage
- * that it protects against overflow when x or y is large).
+ * may apply (12) and (13) directly, (noting that hypot(x, y) is
+ * functionally equivalent to sqrt(y^2 + x^2), with the benefit
+ * of protection against overflow when x or y is large).
*/
- __real__ Z = x = mapfunc(sqrt)(ARGCAST(0.5) * (mapfunc(_hypot)(x, y) + x));
+ __real__ Z = x = mapfunc(sqrt)(ARGCAST(0.5) * (mapfunc(hypot)(x, y) + x));
__imag__ Z = ARGCAST(0.5) * y / x;
}
return Z;
OSDN Git Service
