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, Version 1.0 only 6 * (the "License"). You may not use this file except in compliance 7 * with the License. 8 * 9 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 10 * or http://www.opensolaris.org/os/licensing. 11 * See the License for the specific language governing permissions 12 * and limitations under the License. 13 * 14 * When distributing Covered Code, include this CDDL HEADER in each 15 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 16 * If applicable, add the following below this CDDL HEADER, with the 17 * fields enclosed by brackets "[]" replaced with your own identifying 18 * information: Portions Copyright [yyyy] [name of copyright owner] 19 * 20 * CDDL HEADER END 21 */ 22 /* 23 * Copyright 2004 Sun Microsystems, Inc. All rights reserved. 24 * Use is subject to license terms. 25 */ 26 27 /* 28 * _F_cplx_div_rx(a, w) returns a / w with infinities handled according 29 * to C99. 30 * 31 * If a and w are both finite and w is nonzero, _F_cplx_div_rx(a, w) 32 * delivers the complex quotient q according to the usual formula: 33 * let c = Re(w), and d = Im(w); then q = x + I * y where x = (a * c) 34 * / r and y = (-a * d) / r with r = c * c + d * d. This implementa- 35 * tion computes intermediate results in extended precision to avoid 36 * premature underflow or overflow. 37 * 38 * If a is neither NaN nor zero and w is zero, or if a is infinite 39 * and w is finite and nonzero, _F_cplx_div_rx delivers an infinite 40 * result. If a is finite and w is infinite, _F_cplx_div_rx delivers 41 * a zero result. 42 * 43 * If a and w are both zero or both infinite, or if either a or w is 44 * NaN, _F_cplx_div_rx delivers NaN + I * NaN. C99 doesn't specify 45 * these cases. 46 * 47 * This implementation can raise spurious invalid operation, inexact, 48 * and division-by-zero exceptions. C99 allows this. 49 * 50 * Warning: Do not attempt to "optimize" this code by removing multi- 51 * plications by zero. 52 */ 53 54 #if !defined(i386) && !defined(__i386) && !defined(__amd64) 55 #error This code is for x86 only 56 #endif 57 58 /* 59 * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise 60 */ 61 static int 62 testinff(float x) 63 { 64 union { 65 int i; 66 float f; 67 } xx; 68 69 xx.f = x; 70 return ((((xx.i << 1) - 0xff000000) == 0)? (1 | (xx.i >> 31)) : 0); 71 } 72 73 float _Complex 74 _F_cplx_div_rx(float a, float _Complex w) 75 { 76 float _Complex v; 77 union { 78 int i; 79 float f; 80 } cc, dd; 81 float c, d; 82 long double r, x, y; 83 int i, j; 84 85 /* 86 * The following is equivalent to 87 * 88 * c = crealf(w); d = cimagf(w); 89 */ 90 c = ((float *)&w)[0]; 91 d = ((float *)&w)[1]; 92 93 r = (long double)c * c + (long double)d * d; 94 95 if (r == 0.0f) { 96 /* w is zero; multiply a by 1/Re(w) - I * Im(w) */ 97 c = 1.0f / c; 98 i = testinff(a); 99 if (i) { /* a is infinite */ 100 a = i; 101 } 102 ((float *)&v)[0] = a * c; 103 ((float *)&v)[1] = (a == 0.0f)? a * c : -a * d; 104 return (v); 105 } 106 107 r = (long double)a / r; 108 x = (long double)c * r; 109 y = (long double)-d * r; 110 111 if (x != x || y != y) { 112 /* 113 * x or y is NaN, so a and w can't both be finite and 114 * nonzero. Since we handled the case w = 0 above, the 115 * only case to check here is when w is infinite. 116 */ 117 i = testinff(c); 118 j = testinff(d); 119 if (i | j) { /* w is infinite */ 120 cc.f = c; 121 dd.f = d; 122 c = (cc.i < 0)? -0.0f : 0.0f; 123 d = (dd.i < 0)? -0.0f : 0.0f; 124 x = (long double)c * a; 125 y = (long double)-d * a; 126 } 127 } 128 129 /* 130 * The following is equivalent to 131 * 132 * return x + I * y; 133 */ 134 ((float *)&v)[0] = (float)x; 135 ((float *)&v)[1] = (float)y; 136 return (v); 137 } 138