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 #pragma ident "%Z%%M% %I% %E% SMI" 28 29 /* 30 * _D_cplx_div_ix(b, w) returns (I * b) / w with infinities handled 31 * according to C99. 32 * 33 * If b and w are both finite and w is nonzero, _D_cplx_div_ix(b, w) 34 * delivers the complex quotient q according to the usual formula: 35 * let c = Re(w), and d = Im(w); then q = x + I * y where x = (b * d) 36 * / r and y = (b * c) / r with r = c * c + d * d. This implementa- 37 * tion computes intermediate results in extended precision to avoid 38 * premature underflow or overflow. 39 * 40 * If b is neither NaN nor zero and w is zero, or if b is infinite 41 * and w is finite and nonzero, _D_cplx_div_ix delivers an infinite 42 * result. If b is finite and w is infinite, _D_cplx_div_ix delivers 43 * a zero result. 44 * 45 * If b and w are both zero or both infinite, or if either b or w is 46 * NaN, _D_cplx_div_ix delivers NaN + I * NaN. C99 doesn't specify 47 * these cases. 48 * 49 * This implementation can raise spurious invalid operation, inexact, 50 * and division-by-zero exceptions. C99 allows this. 51 * 52 * Warning: Do not attempt to "optimize" this code by removing multi- 53 * plications by zero. 54 */ 55 56 #if !defined(i386) && !defined(__i386) && !defined(__amd64) 57 #error This code is for x86 only 58 #endif 59 60 /* 61 * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise 62 */ 63 static int 64 testinf(double x) 65 { 66 union { 67 int i[2]; 68 double d; 69 } xx; 70 71 xx.d = x; 72 return (((((xx.i[1] << 1) - 0xffe00000) | xx.i[0]) == 0)? 73 (1 | (xx.i[1] >> 31)) : 0); 74 } 75 76 double _Complex 77 _D_cplx_div_ix(double b, double _Complex w) 78 { 79 double _Complex v; 80 union { 81 int i[2]; 82 double d; 83 } cc, dd; 84 double c, d; 85 long double r, x, y; 86 int i, j; 87 88 /* 89 * The following is equivalent to 90 * 91 * c = creal(w); d = cimag(w); 92 */ 93 /* LINTED alignment */ 94 c = ((double *)&w)[0]; 95 /* LINTED alignment */ 96 d = ((double *)&w)[1]; 97 98 r = (long double)c * c + (long double)d * d; 99 100 if (r == 0.0f) { 101 /* w is zero; multiply b by 1/Re(w) - I * Im(w) */ 102 c = 1.0f / c; 103 j = testinf(b); 104 if (j) { /* b is infinite */ 105 b = j; 106 } 107 /* LINTED alignment */ 108 ((double *)&v)[0] = (b == 0.0f)? b * c : b * d; 109 /* LINTED alignment */ 110 ((double *)&v)[1] = b * c; 111 return (v); 112 } 113 114 r = (long double)b / r; 115 x = (long double)d * r; 116 y = (long double)c * r; 117 118 if (x != x || y != y) { 119 /* 120 * x or y is NaN, so b and w can't both be finite and 121 * nonzero. Since we handled the case w = 0 above, the 122 * only case to check here is when w is infinite. 123 */ 124 i = testinf(c); 125 j = testinf(d); 126 if (i | j) { /* w is infinite */ 127 cc.d = c; 128 dd.d = d; 129 c = (cc.i[1] < 0)? -0.0f : 0.0f; 130 d = (dd.i[1] < 0)? -0.0f : 0.0f; 131 x = (long double)d * b; 132 y = (long double)c * b; 133 } 134 } 135 136 /* 137 * The following is equivalent to 138 * 139 * return x + I * y; 140 */ 141 /* LINTED alignment */ 142 ((double *)&v)[0] = (double)x; 143 /* LINTED alignment */ 144 ((double *)&v)[1] = (double)y; 145 return (v); 146 } 147