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