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 clogf = __clogf 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 fcomplex 39 clogf(fcomplex z) { 40 fcomplex ans; 41 float x, y, ax, ay; 42 double dx, dy; 43 int ix, iy, hx, hy; 44 45 x = F_RE(z); 46 y = F_IM(z); 47 hx = THE_WORD(x); 48 hy = THE_WORD(y); 49 ix = hx & 0x7fffffff; 50 iy = hy & 0x7fffffff; 51 ay = fabsf(y); 52 ax = fabsf(x); 53 F_IM(ans) = atan2f(y, x); 54 if (ix >= 0x7f800000 || iy >= 0x7f800000) { 55 /* x or y is Inf or NaN */ 56 if (iy == 0x7f800000) 57 F_RE(ans) = ay; 58 else if (ix == 0x7f800000) 59 F_RE(ans) = ax; 60 else 61 F_RE(ans) = ax + ay; 62 } else { 63 #if defined(__i386) && !defined(__amd64) 64 int rp = __swapRP(fp_extended); 65 #endif 66 dx = (double)ax; 67 dy = (double)ay; 68 if (ix == 0x3f800000) 69 F_RE(ans) = (float)(0.5 * log1p(dy * dy)); 70 else if (iy == 0x3f800000) 71 F_RE(ans) = (float)(0.5 * log1p(dx * dx)); 72 else if ((ix | iy) == 0) 73 F_RE(ans) = -1.0f / ax; 74 else 75 F_RE(ans) = (float)(0.5 * log(dx * dx + dy * dy)); 76 #if defined(__i386) && !defined(__amd64) 77 if (rp != fp_extended) 78 (void) __swapRP(rp); 79 #endif 80 } 81 return (ans); 82 } 83