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 * On SPARC V8, _Q_cplx_div_rx(v, a, w) sets *v = *a / *w with in-
31 * finities handling according to C99.
32 *
33 * On SPARC V9, _Q_cplx_div_rx(a, w) returns *a / *w with infinities
34 * handled according to C99.
35 *
36 * If a and w are both finite and w is nonzero, _Q_cplx_div_rx de-
37 * livers the complex quotient q according to the usual formula: let
38 * c = Re(w), and d = Im(w); then q = x + I * y where x = (a * c) / r
39 * and y = (-a * d) / r with r = c * c + d * d.  This implementation
40 * scales to avoid premature underflow or overflow.
41 *
42 * If a is neither NaN nor zero and w is zero, or if a is infinite
43 * and w is finite and nonzero, _Q_cplx_div_rx delivers an infinite
44 * result.  If a is finite and w is infinite, _Q_cplx_div_rx delivers
45 * a zero result.
46 *
47 * If a and w are both zero or both infinite, or if either a or w is
48 * NaN, _Q_cplx_div_rx delivers NaN + I * NaN.  C99 doesn't specify
49 * these cases.
50 *
51 * This implementation can raise spurious underflow, overflow, in-
52 * valid operation, inexact, and division-by-zero exceptions.  C99
53 * allows this.
54 */
55
56#if !defined(sparc) && !defined(__sparc)
57#error This code is for SPARC only
58#endif
59
60extern void _Q_scl(long double *, int);
61extern void _Q_scle(long double *, int);
62
63/*
64 * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise
65 */
66static int
67testinfl(long double x)
68{
69	union {
70		int		i[4];
71		long double	q;
72	} xx;
73
74	xx.q = x;
75	return (((((xx.i[0] << 1) - 0xfffe0000) | xx.i[1] | xx.i[2] | xx.i[3])
76		== 0)? (1 | (xx.i[0] >> 31)) : 0);
77}
78
79#ifdef __sparcv9
80long double _Complex
81_Q_cplx_div_rx(const long double *pa, const long double _Complex *w)
82{
83	long double _Complex	v;
84#else
85void
86_Q_cplx_div_rx(long double _Complex *v, const long double *pa,
87	const long double _Complex *w)
88{
89#endif
90	union {
91		int		i[4];
92		long double	q;
93	} aa, cc, dd;
94	long double	a, c, d, sc, sd, r;
95	int		ha, hc, hd, hw, i, j;
96
97	a = *pa;
98
99	/*
100	 * The following is equivalent to
101	 *
102	 *  c = creall(*w); d = cimagl(*w);
103	 */
104	c = ((long double *)w)[0];
105	d = ((long double *)w)[1];
106
107	/* extract high-order words to estimate |a| and |w| */
108	aa.q = a;
109	ha = aa.i[0] & ~0x80000000;
110
111	cc.q = c;
112	dd.q = d;
113	hc = cc.i[0] & ~0x80000000;
114	hd = dd.i[0] & ~0x80000000;
115	hw = (hc > hd)? hc : hd;
116
117	/* check for special cases */
118	if (hw >= 0x7fff0000) { /* w is inf or nan */
119		i = testinfl(c);
120		j = testinfl(d);
121		if (i | j) { /* w is infinite */
122			c = (cc.i[0] < 0)? -0.0l : 0.0l;
123			d = (dd.i[0] < 0)? -0.0l : 0.0l;
124		} else /* w is nan */
125			a += c + d;
126		c *= a;
127		d *= -a;
128		goto done;
129	}
130
131	if (hw == 0 && (cc.i[1] | cc.i[2] | cc.i[3] |
132		dd.i[1] | dd.i[2] | dd.i[3]) == 0) {
133		/* w is zero; multiply a by 1/Re(w) - I * Im(w) */
134		c = 1.0l / c;
135		i = testinfl(a);
136		if (i) { /* a is infinite */
137			a = i;
138		}
139		c *= a;
140		d = (a == 0.0l)? c : -a * d;
141		goto done;
142	}
143
144	if (ha >= 0x7fff0000) { /* a is inf or nan */
145		c *= a;
146		d *= -a;
147		goto done;
148	}
149
150	/*
151	 * Compute the real and imaginary parts of the quotient,
152	 * scaling to avoid overflow or underflow.
153	 */
154	hw = (hw - 0x3fff0000) >> 16;
155	sc = c;
156	sd = d;
157	_Q_scl(&sc, -hw);
158	_Q_scl(&sd, -hw);
159	r = sc * sc + sd * sd;
160
161	ha = (ha - 0x3fff0000) >> 16;
162	_Q_scl(&a, -ha);
163	a /= r;
164	ha -= (hw + hw);
165
166	hc = (hc - 0x3fff0000) >> 16;
167	_Q_scl(&c, -hc);
168	c *= a;
169	hc += ha;
170
171	hd = (hd - 0x3fff0000) >> 16;
172	_Q_scl(&d, -hd);
173	d *= -a;
174	hd += ha;
175
176	/* compensate for scaling */
177	_Q_scle(&c, hc);
178	_Q_scle(&d, hd);
179
180done:
181#ifdef __sparcv9
182	((long double *)&v)[0] = c;
183	((long double *)&v)[1] = d;
184	return (v);
185#else
186	((long double *)v)[0] = c;
187	((long double *)v)[1] = d;
188#endif
189}
190