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 2003 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 * _F_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, _F_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 double 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, _F_cplx_div delivers an infinite 43 * result. If z is finite and w is infinite, _F_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, _F_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(sparc) && !defined(__sparc) 58 #error This code is for SPARC only 59 #endif 60 61 static union { 62 int i[2]; 63 double d; 64 } inf = { 65 0x7ff00000, 0 66 }; 67 68 /* 69 * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise 70 */ 71 static int 72 testinff(float x) 73 { 74 union { 75 int i; 76 float f; 77 } xx; 78 79 xx.f = x; 80 return ((((xx.i << 1) - 0xff000000) == 0)? (1 | (xx.i >> 31)) : 0); 81 } 82 83 float _Complex 84 _F_cplx_div(float _Complex z, float _Complex w) 85 { 86 float _Complex v; 87 union { 88 int i; 89 float f; 90 } cc, dd; 91 float a, b, c, d; 92 double r, x, y; 93 int i, j, recalc; 94 95 /* 96 * The following is equivalent to 97 * 98 * a = crealf(z); b = cimagf(z); 99 * c = crealf(w); d = cimagf(w); 100 */ 101 a = ((float *)&z)[0]; 102 b = ((float *)&z)[1]; 103 c = ((float *)&w)[0]; 104 d = ((float *)&w)[1]; 105 106 r = (double)c * c + (double)d * d; 107 108 if (r == 0.0) { 109 /* w is zero; multiply z by 1/Re(w) - I * Im(w) */ 110 c = 1.0f / c; 111 i = testinff(a); 112 j = testinff(b); 113 if (i | j) { /* z is infinite */ 114 a = i; 115 b = j; 116 } 117 ((float *)&v)[0] = a * c + b * d; 118 ((float *)&v)[1] = b * c - a * d; 119 return (v); 120 } 121 122 r = 1.0 / r; 123 x = ((double)a * c + (double)b * d) * r; 124 y = ((double)b * c - (double)a * d) * r; 125 126 if (x != x && y != y) { 127 /* 128 * Both x and y are NaN, so z and w can't both be finite 129 * and nonzero. Since we handled the case w = 0 above, 130 * the only cases to check here are when one of z or w 131 * is infinite. 132 */ 133 r = 1.0; 134 recalc = 0; 135 i = testinff(a); 136 j = testinff(b); 137 if (i | j) { /* z is infinite */ 138 /* "factor out" infinity */ 139 a = i; 140 b = j; 141 r = inf.d; 142 recalc = 1; 143 } 144 i = testinff(c); 145 j = testinff(d); 146 if (i | j) { /* w is infinite */ 147 /* 148 * "factor out" infinity, being careful to preserve 149 * signs of finite values 150 */ 151 cc.f = c; 152 dd.f = d; 153 c = i? i : ((cc.i < 0)? -0.0f : 0.0f); 154 d = j? j : ((dd.i < 0)? -0.0f : 0.0f); 155 r *= 0.0; 156 recalc = 1; 157 } 158 if (recalc) { 159 x = ((double)a * c + (double)b * d) * r; 160 y = ((double)b * c - (double)a * d) * r; 161 } 162 } 163 164 /* 165 * The following is equivalent to 166 * 167 * return x + I * y; 168 */ 169 ((float *)&v)[0] = (float)x; 170 ((float *)&v)[1] = (float)y; 171 return (v); 172 } 173