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 } 139