Print this page
11210 libm should be cstyle(1ONBLD) clean
   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 }
   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 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 
  26 /*
  27  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __catanf = catanf
  32 
  33 #include "libm.h"
  34 #include "complex_wrapper.h"
  35 
  36 #if defined(__i386) && !defined(__amd64)
  37 extern int __swapRP(int);
  38 #endif
  39 
  40 static const float pi_2 = 1.570796326794896558e+00F,

  41         zero = 0.0F,
  42         half = 0.5F,
  43         two = 2.0F,
  44         one = 1.0F;
  45 
  46 fcomplex
  47 catanf(fcomplex z)
  48 {
  49         fcomplex ans;
  50         float x, y, ax, ay, t;
  51         double dx, dy, dt;
  52         int hx, hy, ix, iy;
  53 
  54         x = F_RE(z);
  55         y = F_IM(z);
  56         ax = fabsf(x);
  57         ay = fabsf(y);
  58         hx = THE_WORD(x);
  59         hy = THE_WORD(y);
  60         ix = hx & 0x7fffffff;
  61         iy = hy & 0x7fffffff;
  62 
  63         if (ix >= 0x7f800000) {              /* x is inf or NaN */
  64                 if (ix == 0x7f800000) {
  65                         F_RE(ans) = pi_2;
  66                         F_IM(ans) = zero;
  67                 } else {
  68                         F_RE(ans) = x * x;
  69 
  70                         if (iy == 0 || iy == 0x7f800000)
  71                                 F_IM(ans) = zero;
  72                         else
  73                                 F_IM(ans) = (fabsf(y) - ay) / (fabsf(y) - ay);
  74                 }
  75         } else if (iy >= 0x7f800000) {       /* y is inf or NaN */
  76                 if (iy == 0x7f800000) {
  77                         F_RE(ans) = pi_2;
  78                         F_IM(ans) = zero;
  79                 } else {
  80                         F_RE(ans) = (fabsf(x) - ax) / (fabsf(x) - ax);
  81                         F_IM(ans) = y * y;
  82                 }
  83         } else if (ix == 0) {
  84                 /* BEGIN CSTYLED */
  85                 /*
  86                  * x = 0
  87                  *      1                            1
  88                  * A = --- * atan2(2x, 1-x*x-y*y) = --- atan2(0,1-|y|)
  89                  *      2                            2
  90                  *
  91                  *     1     [ (y+1)*(y+1) ]   1          2      1         2y
  92                  * B = - log [ ----------- ] = - log (1+ ---) or - log(1+ ----)
  93                  *     4     [ (y-1)*(y-1) ]   2         y-1     2         1-y
  94                  */
  95                 /* END CSTYLED */
  96                 t = one - ay;
  97 
  98                 if (iy == 0x3f800000) {
  99                         /* y=1: catan(0,1)=(0,+inf) with 1/0 signal */
 100                         F_IM(ans) = ay / ax;
 101                         F_RE(ans) = zero;
 102                 } else if (iy > 0x3f800000) {        /* y>1 */
 103                         F_IM(ans) = half * log1pf(two / (-t));
 104                         F_RE(ans) = pi_2;
 105                 } else {                        /* y<1 */
 106                         F_IM(ans) = half * log1pf((ay + ay) / t);
 107                         F_RE(ans) = zero;
 108                 }
 109         } else {
 110                 /* BEGIN CSTYLED */
 111                 /*
 112                  * use double precision x,y
 113                  *      1
 114                  * A = --- * atan2(2x, 1-x*x-y*y)
 115                  *      2
 116                  *
 117                  *     1     [ x*x+(y+1)*(y+1) ]   1               4y
 118                  * B = - log [ --------------- ] = - log (1+ -----------------)
 119                  *     4     [ x*x+(y-1)*(y-1) ]   4         x*x + (y-1)*(y-1)
 120                  */
 121                 /* END CSTYLED */
 122 
 123 #if defined(__i386) && !defined(__amd64)
 124                 int rp = __swapRP(fp_extended);
 125 #endif
 126                 dx = (double)ax;
 127                 dy = (double)ay;
 128                 F_RE(ans) = (float)(0.5 * atan2(dx + dx, 1.0 - dx * dx - dy *
 129                     dy));
 130                 dt = dy - 1.0;
 131                 F_IM(ans) = (float)(0.25 * log1p(4.0 * dy / (dx * dx + dt *
 132                     dt)));
 133 
 134 #if defined(__i386) && !defined(__amd64)
 135                 if (rp != fp_extended)
 136                         (void) __swapRP(rp);
 137 #endif
 138         }
 139 
 140         if (hx < 0)
 141                 F_RE(ans) = -F_RE(ans);
 142 
 143         if (hy < 0)
 144                 F_IM(ans) = -F_IM(ans);
 145 
 146         return (ans);
 147 }