mingwrt/mingwex/complex/casin_generic.c

   1 /*
   2  * casin_generic.c
   3  *
   4  * $Id$
   5  *
   6  * Compute the complex arcsin corresponding to a complex sine value;
   7  * a generic implementation for casin(), casinf(), and casinl().
   8  *
   9  * Written by Keith Marshall <keith@users.osdn.me>
  10  * This is an adaptation of an original contribution by Danny Smith
  11  * Copyright (C) 2003, 2014, 2022, MinGW.OSDN Project
  12  *
  13  *
  14  * Permission is hereby granted, free of charge, to any person obtaining a
  15  * copy of this software and associated documentation files (the "Software"),
  16  * to deal in the Software without restriction, including without limitation
  17  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
  18  * and/or sell copies of the Software, and to permit persons to whom the
  19  * Software is furnished to do so, subject to the following conditions:
  20  *
  21  * The above copyright notice, this permission notice, and the following
  22  * disclaimer shall be included in all copies or substantial portions of
  23  * the Software.
  24  *
  25  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  26  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  27  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
  28  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  29  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  30  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OF OR OTHER
  31  * DEALINGS IN THE SOFTWARE.
  32  *
  33  *
  34  * This is a generic implementation for all of the casin(), casinl(),
  35  * and casinf() functions; each is to be compiled separately, i.e.
  36  *
  37  *   gcc -D FUNCTION=casin  -o casin.o  casin_generic.c
  38  *   gcc -D FUNCTION=casinl -o casinl.o casin_generic.c
  39  *   gcc -D FUNCTION=casinf -o casinf.o casin_generic.c
  40  *
  41  */
  42 #include <math.h>
  43 #include <complex.h>
  44
  45 #ifndef FUNCTION
  46 /* If user neglected to specify it, the default compilation is for
  47  * the casin() function.
  48  */
  49 # define FUNCTION casin
  50 #endif
  51
  52 #define argtype_casin  double
  53 #define argtype_casinl long double
  54 #define argtype_casinf float
  55
  56 #define PASTE(PREFIX,SUFFIX)   PREFIX##SUFFIX
  57 #define mapname(PREFIX,SUFFIX) PASTE(PREFIX,SUFFIX)
  58
  59 #define ARGTYPE(FUNCTION) PASTE(argtype_,FUNCTION)
  60 #define ARGCAST(VALUE)    mapname(argcast_,FUNCTION)(VALUE)
  61
  62 #define argcast_casin(VALUE)  VALUE
  63 #define argcast_casinl(VALUE) PASTE(VALUE,L)
  64 #define argcast_casinf(VALUE) PASTE(VALUE,F)
  65
  66 #define mapfunc(NAME) mapname(mapfunc_,FUNCTION)(NAME)
  67
  68 #define mapfunc_casin(NAME)  NAME
  69 #define mapfunc_casinl(NAME) PASTE(NAME,l)
  70 #define mapfunc_casinf(NAME) PASTE(NAME,f)
  71
  72 /* Define the generic function implementation.
  73  */
  74 ARGTYPE(FUNCTION) complex FUNCTION( ARGTYPE(FUNCTION) complex Z )
  75 {
  76   ARGTYPE(FUNCTION) x = __real__ Z;     /* real part of Z saved as x */
  77   ARGTYPE(FUNCTION) y = __imag__ Z;     /* imaginary part of Z saved as iy */
  78
  79   /* Having thus saved the original value of Z as separate real and
  80    * imaginary parts, we may reuse Z to compute the arcsin value...
  81    */
  82   if( (y == ARGCAST(0.0)) && (ARGCAST(1.0) >= mapfunc(fabs)(x)) )
  83   {
  84     /* ...but when Z is entirely real, and lies within the valid
  85      * domain, (-1.0 .. +1.0), for real sine values, we may simply
  86      * optimize this call, using the real arcsin function...
  87      */
  88     __real__ Z = mapfunc(asin)(x);
  89   }
  90   else
  91   { /* Z is complex, (i.e. it has a non-zero imaginary part), or it
  92      * is entirely real but lies outside the valid domain for a real
  93      * sine value; we must evaluate the complex arcsin expression:
  94      *
  95      *   casin(Z) = -i * clog(i * Z + csqrt(1.0 - Z * Z))
  96      *
  97      * Here, we use a local accumulator for the evaluation:
  98      */
  99     ARGTYPE(FUNCTION) complex ZZ;
 100
 101     /* First, note that we may express 1.0 - Z * Z as equivalent to
 102      * 1.0 - (x + iy) * (x + iy), which expands to:
 103      *
 104      *   1.0 - (x * x + i * y * i * y + 2.0 * i * x * y)
 105      *
 106      * and, noting that i * i is equivalent to -1, this simplifies
 107      * to become:
 108      *
 109      *   1.0 - (x * x - y * y + 2.0 * i * x * y)
 110      *
 111      * Separating real and imaginary parts yields:
 112      *
 113      *   (1.0 - x * x - y * y) - (i * 2.0 * x * y)
 114      *
 115      * and factorization of the real part then yields:
 116      *
 117      *   (1.0 - (x - y) * (x + y)) - (i * 2.0 * x * y)
 118      *
 119      * hence...
 120      */
 121     __real__ ZZ = ARGCAST(1.0) - (x - y) * (x + y);
 122     __imag__ ZZ = ARGCAST(-2.0) * x * y;
 123
 124     /* ...and the complex square root term may be computed as...
 125      */
 126     ZZ = mapfunc(csqrt)(ZZ);
 127
 128     /* ...to which we then add i * (x + iy), (i.e. i times the
 129      * original value of Z); here, we note that i * iy is equal to
 130      * -1 * y, and hence adding it is equivalent to subtracting y
 131      *  from the real part of our recomputed Z, while adding i * x
 132      *  is equivalent to adding x to the imaginary part...
 133      */
 134     __real__ ZZ -= y;
 135     __imag__ ZZ += x;
 136
 137     /* ...thus obtaining the intermediate expression value from
 138      * which we must compute the complex logarithm...
 139      */
 140     ZZ = mapfunc(clog)(ZZ);
 141
 142     /* ...and ultimately, the arcsin value to be returned is the
 143      * value of this complex logarithm multiplied by -1 * i; (the
 144      * imaginary part of the intermediate result, multiplied by i
 145      * becomes the real part of the result, with inverted sign,
 146      * which is then inverted again on multiplication by -1...
 147      */
 148     __real__ Z = __imag__ ZZ;
 149
 150     /* ...while the real part of the intermediate simply becomes
 151      * the imaginary part of the result, with sign preserved, but
 152      * then inverted by multiplication by -1).
 153      */
 154     __imag__ Z = ARGCAST(-1.0) * __real__ ZZ;
 155   }
 156   return Z;
 157 }
 158
 159 /* $RCSfile$: end of file */