mingwrt/mingwex/complex/cpow_generic.c

   1 /*
   2  * cpow_generic.c
   3  *
   4  * $Id$
   5  *
   6  * Compute the result of raising a complex number, Z, to a complex
   7  * power, N; this provides a generic implementation for each of the
   8  * cpow(), cpowf(), and cpowl() functions.
   9  *
  10  * Written by Keith Marshall <keithmarshall@users.sourceforge.net>
  11  * This is an adaptation of an original contribution by Danny Smith
  12  * Copyright (C) 2003, 2014, MinGW.org Project
  13  *
  14  *
  15  * Permission is hereby granted, free of charge, to any person obtaining a
  16  * copy of this software and associated documentation files (the "Software"),
  17  * to deal in the Software without restriction, including without limitation
  18  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  19  * and/or sell copies of the Software, and to permit persons to whom the
  20  * Software is furnished to do so, subject to the following conditions:
  21  *
  22  * The above copyright notice, this permission notice, and the following
  23  * disclaimer shall be included in all copies or substantial portions of
  24  * the Software.
  25  *
  26  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  27  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  28  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  29  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  30  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  31  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OF OR OTHER
  32  * DEALINGS IN THE SOFTWARE.
  33  *
  34  *
  35  * This is a generic implementation for all of the cpow(), cpowl(),
  36  * and cpowh() functions; each is to be compiled separately, i.e.
  37  *
  38  *   gcc -D FUNCTION=cpow  -o cpow.o  cpow_generic.c
  39  *   gcc -D FUNCTION=cpowl -o cpowl.o cpow_generic.c
  40  *   gcc -D FUNCTION=cpowf -o cpowf.o cpow_generic.c
  41  *
  42  */
  43 #include <math.h>
  44 #include <complex.h>
  45
  46 #ifndef FUNCTION
  47 /* If user neglected to specify it, the default compilation is for
  48  * the cpow() function.
  49  */
  50 # define FUNCTION cpow
  51 #endif
  52
  53 #define argtype_cpow  double
  54 #define argtype_cpowl long double
  55 #define argtype_cpowf float
  56
  57 #define PASTE(PREFIX,SUFFIX)   PREFIX##SUFFIX
  58 #define mapname(PREFIX,SUFFIX) PASTE(PREFIX,SUFFIX)
  59
  60 #define ARGTYPE(FUNCTION) PASTE(argtype_,FUNCTION)
  61 #define ARGCAST(VALUE)    mapname(argcast_,FUNCTION)(VALUE)
  62
  63 #define argcast_cpow(VALUE)  VALUE
  64 #define argcast_cpowl(VALUE) PASTE(VALUE,L)
  65 #define argcast_cpowf(VALUE) PASTE(VALUE,F)
  66
  67 #define mapfunc(NAME) mapname(mapfunc_,FUNCTION)(NAME)
  68
  69 #define mapfunc_cpow(NAME)  NAME
  70 #define mapfunc_cpowl(NAME) PASTE(NAME,l)
  71 #define mapfunc_cpowf(NAME) PASTE(NAME,f)
  72
  73 /* Danny Smith's original implementation explicitly referenced the
  74  * DLL entry point for the pow() function, in the case of the cpow()
  75  * function only; the following variation on mapfunc() reproduces
  76  * this implementation model.
  77  */
  78 extern double (*_imp__pow)( double, double );
  79 #define powfunc mapname(powfunc_,FUNCTION)
  80 #define powfunc_cpow  (*_imp__pow)
  81 #define powfunc_cpowf  powf
  82 #define powfunc_cpowl  powl
  83
  84 /* Define the generic function implementation.
  85  */
  86 ARGTYPE(FUNCTION) complex FUNCTION
  87 ( ARGTYPE(FUNCTION) complex Z, ARGTYPE(FUNCTION) complex N )
  88 {
  89   ARGTYPE(FUNCTION) radius;
  90   if( (radius = mapfunc(hypot)( __real__ Z, __imag__ Z )) == ARGCAST(0.0) )
  91     /*
  92      * Modulus of Z is zero, hence result is also zero; force it.
  93      */
  94     __real__ Z = __imag__ Z = ARGCAST(0.0);
  95
  96   else
  97   { /* Compute the complex power function result, in terms of its
  98      * polar (r, theta) representation in the complex plane.
  99      */
 100     ARGTYPE(FUNCTION) r = mapfunc(carg)( Z );
 101     ARGTYPE(FUNCTION) theta = r * __real__ N;
 102
 103     if( __imag__ N == ARGCAST(0.0) )
 104       /*
 105        * For entirely real N, the power function tends to yield
 106        * more accuracy than exponential expansion.
 107        */
 108       r = powfunc( radius, __real__ N );
 109
 110     else
 111     {
 112       theta += (radius = mapfunc(log)( radius )) * __imag__ N;
 113       r = mapfunc(exp)( radius * __real__ N - r * __imag__ N );
 114     }
 115
 116     /* Convert polar to cartesian representation for return.
 117      */
 118     __real__ Z = r * mapfunc(cos)( theta );
 119     __imag__ Z = r * mapfunc(sin)( theta );
 120   }
 121
 122   /* Regardless of how we got here, the original base value has
 123    * now been raised to the specified power; return it.
 124    */
 125   return Z;
 126 }
 127
 128 /* $RCSfile$: end of file */