usr/src/lib/libm/common/complex/catanf.c

   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 /*
  22  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  23  */
  24 /*
  25  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  26  * Use is subject to license terms.
  27  */
  28
  29 #pragma weak __catanf = catanf
  30
  31 #include "libm.h"
  32 #include "complex_wrapper.h"
  33
  34 #if defined(__i386) && !defined(__amd64)
  35 extern int __swapRP(int);
  36 #endif
  37
  38 static const float
  39         pi_2 = 1.570796326794896558e+00F,
  40         zero = 0.0F,
  41         half = 0.5F,
  42         two = 2.0F,
  43         one = 1.0F;
  44
  45 fcomplex
  46 catanf(fcomplex z) {
  47         fcomplex        ans;
  48         float           x, y, ax, ay, t;
  49         double          dx, dy, dt;
  50         int             hx, hy, ix, iy;
  51
  52         x = F_RE(z);
  53         y = F_IM(z);
  54         ax = fabsf(x);
  55         ay = fabsf(y);
  56         hx = THE_WORD(x);
  57         hy = THE_WORD(y);
  58         ix = hx & 0x7fffffff;
  59         iy = hy & 0x7fffffff;
  60
  61         if (ix >= 0x7f800000) {         /* x is inf or NaN */
  62                 if (ix == 0x7f800000) {
  63                         F_RE(ans) = pi_2;
  64                         F_IM(ans) = zero;
  65                 } else {
  66                         F_RE(ans) = x * x;
  67                         if (iy == 0 || iy == 0x7f800000)
  68                                 F_IM(ans) = zero;
  69                         else
  70                                 F_IM(ans) = (fabsf(y) - ay) / (fabsf(y) - ay);
  71                 }
  72         } else if (iy >= 0x7f800000) {  /* y is inf or NaN */
  73                 if (iy == 0x7f800000) {
  74                         F_RE(ans) = pi_2;
  75                         F_IM(ans) = zero;
  76                 } else {
  77                         F_RE(ans) = (fabsf(x) - ax) / (fabsf(x) - ax);
  78                         F_IM(ans) = y * y;
  79                 }
  80         } else if (ix == 0) {
  81                 /* INDENT OFF */
  82                 /*
  83                  * x = 0
  84                  *      1                            1
  85                  * A = --- * atan2(2x, 1-x*x-y*y) = --- atan2(0,1-|y|)
  86                  *      2                            2
  87                  *
  88                  *     1     [ (y+1)*(y+1) ]   1          2      1         2y
  89                  * B = - log [ ----------- ] = - log (1+ ---) or - log(1+ ----)
  90                  *     4     [ (y-1)*(y-1) ]   2         y-1     2         1-y
  91                  */
  92                 /* INDENT ON */
  93                 t = one - ay;
  94                 if (iy == 0x3f800000) {
  95                         /* y=1: catan(0,1)=(0,+inf) with 1/0 signal */
  96                         F_IM(ans) = ay / ax;
  97                         F_RE(ans) = zero;
  98                 } else if (iy > 0x3f800000) {   /* y>1 */
  99                         F_IM(ans) = half * log1pf(two / (-t));
 100                         F_RE(ans) = pi_2;
 101                 } else {                /* y<1 */
 102                         F_IM(ans) = half * log1pf((ay + ay) / t);
 103                         F_RE(ans) = zero;
 104                 }
 105         } else {
 106                 /* INDENT OFF */
 107                 /*
 108                  * use double precision x,y
 109                  *      1
 110                  * A = --- * atan2(2x, 1-x*x-y*y)
 111                  *      2
 112                  *
 113                  *     1     [ x*x+(y+1)*(y+1) ]   1               4y
 114                  * B = - log [ --------------- ] = - log (1+ -----------------)
 115                  *     4     [ x*x+(y-1)*(y-1) ]   4         x*x + (y-1)*(y-1)
 116                  */
 117                 /* INDENT ON */
 118 #if defined(__i386) && !defined(__amd64)
 119                 int     rp = __swapRP(fp_extended);
 120 #endif
 121                 dx = (double)ax;
 122                 dy = (double)ay;
 123                 F_RE(ans) = (float)(0.5 * atan2(dx + dx,
 124                     1.0 - dx * dx - dy * dy));
 125                 dt = dy - 1.0;
 126                 F_IM(ans) = (float)(0.25 * log1p(4.0 * dy /
 127                     (dx * dx + dt * dt)));
 128 #if defined(__i386) && !defined(__amd64)
 129                 if (rp != fp_extended)
 130                         (void) __swapRP(rp);
 131 #endif
 132         }
 133         if (hx < 0)
 134                 F_RE(ans) = -F_RE(ans);
 135         if (hy < 0)
 136                 F_IM(ans) = -F_IM(ans);
 137         return (ans);
 138 }