compiler/stdc/math/s_casinl.c

   1 /*      $OpenBSD: s_casinl.c,v 1.3 2011/07/20 21:02:51 martynas Exp $   */
   2
   3 /*
   4  * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
   5  *
   6  * Permission to use, copy, modify, and distribute this software for any
   7  * purpose with or without fee is hereby granted, provided that the above
   8  * copyright notice and this permission notice appear in all copies.
   9  *
  10  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  11  * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
  12  * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
  13  * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  14  * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
  15  * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
  16  * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  17  */
  18
  19 /*                                                      casinl()
  20  *
  21  *      Complex circular arc sine
  22  *
  23  *
  24  *
  25  * SYNOPSIS:
  26  *
  27  * long double complex casinl();
  28  * long double complex z, w;
  29  *
  30  * w = casinl( z );
  31  *
  32  *
  33  *
  34  * DESCRIPTION:
  35  *
  36  * Inverse complex sine:
  37  *
  38  *                               2
  39  * w = -i clog( iz + csqrt( 1 - z ) ).
  40  *
  41  *
  42  * ACCURACY:
  43  *
  44  *                      Relative error:
  45  * arithmetic   domain     # trials      peak         rms
  46  *    DEC       -10,+10     10100       2.1e-15     3.4e-16
  47  *    IEEE      -10,+10     30000       2.2e-14     2.7e-15
  48  * Larger relative error can be observed for z near zero.
  49  * Also tested by csin(casin(z)) = z.
  50  */
  51
  52 #include <float.h>
  53 #include <complex.h>
  54 #include "math.h"
  55
  56 #if     LDBL_MANT_DIG == 64
  57 static const long double MACHEPL= 5.42101086242752217003726400434970855712890625E-20L;
  58 #elif   LDBL_MANT_DIG == 113
  59 static const long double MACHEPL = 9.629649721936179265279889712924636592690508e-35L;
  60 #endif
  61
  62 static const long double PIO2L = 1.570796326794896619231321691639751442098585L;
  63
  64 long double complex
  65 casinl(long double complex z)
  66 {
  67         long double complex w;
  68         long double x, y, b;
  69         long double complex ca, ct, zz, z2;
  70
  71         x = creall(z);
  72         y = cimagl(z);
  73
  74         if (y == 0.0L) {
  75                 if (fabsl(x) > 1.0L) {
  76                         w = PIO2L + 0.0L * I;
  77                         /*mtherr( "casinl", DOMAIN );*/
  78                 }
  79                 else {
  80                         w = asinl(x) + 0.0L * I;
  81                 }
  82                 return (w);
  83         }
  84
  85         /* Power series expansion */
  86         b = cabsl(z);
  87         if (b < 0.125L) {
  88                 long double complex sum;
  89                 long double n, cn;
  90
  91                 z2 = (x - y) * (x + y) + (2.0L * x * y) * I;
  92                 cn = 1.0L;
  93                 n = 1.0L;
  94                 ca = x + y * I;
  95                 sum = x + y * I;
  96                 do {
  97                         ct = z2 * ca;
  98                         ca = ct;
  99
 100                         cn *= n;
 101                         n += 1.0L;
 102                         cn /= n;
 103                         n += 1.0L;
 104                         b = cn/n;
 105
 106                         ct *= b;
 107                         sum += ct;
 108                         b = cabsl(ct);
 109                 }
 110
 111                 while (b > MACHEPL);
 112                 w = sum;
 113                 return w;
 114         }
 115
 116         ca = x + y * I;
 117         ct = ca * I;    /* iz */
 118         /* sqrt(1 - z*z) */
 119         /* cmul(&ca, &ca, &zz) */
 120         /* x * x  -  y * y */
 121         zz = (x - y) * (x + y) + (2.0L * x * y) * I;
 122         zz = 1.0L - creall(zz) - cimagl(zz) * I;
 123         z2 = csqrtl(zz);
 124
 125         zz = ct + z2;
 126         zz = clogl(zz);
 127         /* multiply by 1/i = -i */
 128         w = zz * (-1.0L * I);
 129         return (w);
 130 }