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