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 * _D_cplx_div_ix(b, w) returns (I * b) / w with infinities handled 29 * according to C99. 30 * 31 * If b and w are both finite and w is nonzero, _D_cplx_div_ix(b, 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 = (b * d) 34 * / r and y = (b * c) / 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 b is neither NaN nor zero and w is zero, or if b is infinite 39 * and w is finite and nonzero, _D_cplx_div_ix delivers an infinite 40 * result. If b is finite and w is infinite, _D_cplx_div_ix delivers 41 * a zero result. 42 * 43 * If b and w are both zero or both infinite, or if either b or w is 44 * NaN, _D_cplx_div_ix 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 testinf(double x) 63 { 64 union { 65 int i[2]; 66 double d; 67 } xx; 68 69 xx.d = x; 70 return (((((xx.i[1] << 1) - 0xffe00000) | xx.i[0]) == 0)? 71 (1 | (xx.i[1] >> 31)) : 0); 72 } 73 74 double _Complex 75 _D_cplx_div_ix(double b, double _Complex w) 76 { 77 double _Complex v; 78 union { 79 int i[2]; 80 double d; 81 } cc, dd; 82 double c, d; 83 long double r, x, y; 84 int i, j; 85 86 /* 87 * The following is equivalent to 88 * 89 * c = creal(w); d = cimag(w); 90 */ 91 /* LINTED alignment */ 92 c = ((double *)&w)[0]; 93 /* LINTED alignment */ 94 d = ((double *)&w)[1]; 95 96 r = (long double)c * c + (long double)d * d; 97 98 if (r == 0.0f) { 99 /* w is zero; multiply b by 1/Re(w) - I * Im(w) */ 100 c = 1.0f / c; 101 j = testinf(b); 102 if (j) { /* b is infinite */ 103 b = j; 104 } 105 /* LINTED alignment */ 106 ((double *)&v)[0] = (b == 0.0f)? b * c : b * d; 107 /* LINTED alignment */ 108 ((double *)&v)[1] = b * c; 109 return (v); 110 } 111 112 r = (long double)b / r; 113 x = (long double)d * r; 114 y = (long double)c * r; 115 116 if (x != x || y != y) { 117 /* 118 * x or y is NaN, so b and w can't both be finite and 119 * nonzero. Since we handled the case w = 0 above, the 120 * only case to check here is when w is infinite. 121 */ 122 i = testinf(c); 123 j = testinf(d); 124 if (i | j) { /* w is infinite */ 125 cc.d = c; 126 dd.d = d; 127 c = (cc.i[1] < 0)? -0.0f : 0.0f; 128 d = (dd.i[1] < 0)? -0.0f : 0.0f; 129 x = (long double)d * b; 130 y = (long double)c * b; 131 } 132 } 133 134 /* 135 * The following is equivalent to 136 * 137 * return x + I * y; 138 */ 139 /* LINTED alignment */ 140 ((double *)&v)[0] = (double)x; 141 /* LINTED alignment */ 142 ((double *)&v)[1] = (double)y; 143 return (v); 144 } 145