usr/src/lib/libc/i386/fp/_X_cplx_div.c

   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  * _X_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, _X_cplx_div delivers
  34  * the complex quotient q according to the usual formula: let a =
  35  * Re(z), b = Im(z), c = Re(w), and d = Im(w); then q = x + I * y
  36  * where x = (a * c + b * d) / r and y = (b * c - a * d) / r with
  37  * r = c * c + d * d.  This implementation scales to avoid premature
  38  * underflow or overflow.
  39  *
  40  * If z is neither NaN nor zero and w is zero, or if z is infinite
  41  * and w is finite and nonzero, _X_cplx_div delivers an infinite
  42  * result.  If z is finite and w is infinite, _X_cplx_div delivers
  43  * a zero result.
  44  *
  45  * If z and w are both zero or both infinite, or if either z or w is
  46  * a complex NaN, _X_cplx_div delivers NaN + I * NaN.  C99 doesn't
  47  * specify these cases.
  48  *
  49  * This implementation can raise spurious underflow, overflow, in-
  50  * valid operation, inexact, and division-by-zero exceptions.  C99
  51  * allows this.
  52  */
  53
  54 #if !defined(i386) && !defined(__i386) && !defined(__amd64)
  55 #error This code is for x86 only
  56 #endif
  57
  58 static union {
  59         int     i;
  60         float   f;
  61 } inf = {
  62         0x7f800000
  63 };
  64
  65 /*
  66  * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise
  67  */
  68 static int
  69 testinfl(long double x)
  70 {
  71         union {
  72                 int             i[3];
  73                 long double     e;
  74         } xx;
  75
  76         xx.e = x;
  77         if ((xx.i[2] & 0x7fff) != 0x7fff || ((xx.i[1] << 1) | xx.i[0]) != 0)
  78                 return (0);
  79         return (1 | ((xx.i[2] << 16) >> 31));
  80 }
  81
  82 long double _Complex
  83 _X_cplx_div(long double _Complex z, long double _Complex w)
  84 {
  85         long double _Complex    v;
  86         union {
  87                 int             i[3];
  88                 long double     e;
  89         } aa, bb, cc, dd, ss;
  90         long double     a, b, c, d, r;
  91         int             ea, eb, ec, ed, ez, ew, es, i, j;
  92
  93         /*
  94          * The following is equivalent to
  95          *
  96          *  a = creall(*z); b = cimagl(*z);
  97          *  c = creall(*w); d = cimagl(*w);
  98          */
  99         a = ((long double *)&z)[0];
 100         b = ((long double *)&z)[1];
 101         c = ((long double *)&w)[0];
 102         d = ((long double *)&w)[1];
 103
 104         /* extract exponents to estimate |z| and |w| */
 105         aa.e = a;
 106         bb.e = b;
 107         ea = aa.i[2] & 0x7fff;
 108         eb = bb.i[2] & 0x7fff;
 109         ez = (ea > eb)? ea : eb;
 110
 111         cc.e = c;
 112         dd.e = d;
 113         ec = cc.i[2] & 0x7fff;
 114         ed = dd.i[2] & 0x7fff;
 115         ew = (ec > ed)? ec : ed;
 116
 117         /* check for special cases */
 118         if (ew >= 0x7fff) { /* w is inf or nan */
 119                 r = 0.0f;
 120                 i = testinfl(c);
 121                 j = testinfl(d);
 122                 if (i | j) { /* w is infinite */
 123                         /*
 124                          * "factor out" infinity, being careful to preserve
 125                          * signs of finite values
 126                          */
 127                         c = i? i : (((cc.i[2] << 16) < 0)? -0.0f : 0.0f);
 128                         d = j? j : (((dd.i[2] << 16) < 0)? -0.0f : 0.0f);
 129                         if (ez >= 0x7ffe) {
 130                                 /* scale to avoid overflow below */
 131                                 c *= 0.5f;
 132                                 d *= 0.5f;
 133                         }
 134                 }
 135                 ((long double *)&v)[0] = (a * c + b * d) * r;
 136                 ((long double *)&v)[1] = (b * c - a * d) * r;
 137                 return (v);
 138         }
 139
 140         if (ew == 0 && (cc.i[1] | cc.i[0] | dd.i[1] | dd.i[0]) == 0) {
 141                 /* w is zero; multiply z by 1/Re(w) - I * Im(w) */
 142                 c = 1.0f / c;
 143                 i = testinfl(a);
 144                 j = testinfl(b);
 145                 if (i | j) { /* z is infinite */
 146                         a = i;
 147                         b = j;
 148                 }
 149                 ((long double *)&v)[0] = a * c + b * d;
 150                 ((long double *)&v)[1] = b * c - a * d;
 151                 return (v);
 152         }
 153
 154         if (ez >= 0x7fff) { /* z is inf or nan */
 155                 i = testinfl(a);
 156                 j = testinfl(b);
 157                 if (i | j) { /* z is infinite */
 158                         a = i;
 159                         b = j;
 160                         r = inf.f;
 161                 }
 162                 ((long double *)&v)[0] = a * c + b * d;
 163                 ((long double *)&v)[1] = b * c - a * d;
 164                 return (v);
 165         }
 166
 167         /*
 168          * Scale c and d to compute 1/|w|^2 and the real and imaginary
 169          * parts of the quotient.
 170          */
 171         es = ((ew >> 2) - ew) + 0x6ffd;
 172         if (ez < 0x0086) { /* |z| < 2^-16249 */
 173                 if (((ew - 0x3efe) | (0x4083 - ew)) >= 0)
 174                         es = ((0x4083 - ew) >> 1) + 0x3fff;
 175         }
 176         ss.i[2] = es;
 177         ss.i[1] = 0x80000000;
 178         ss.i[0] = 0;
 179
 180         c *= ss.e;
 181         d *= ss.e;
 182         r = 1.0f / (c * c + d * d);
 183
 184         c *= ss.e;
 185         d *= ss.e;
 186
 187         ((long double *)&v)[0] = (a * c + b * d) * r;
 188         ((long double *)&v)[1] = (b * c - a * d) * r;
 189         return (v);
 190 }