usr/src/lib/libc/i386/fp/_X_cplx_mul.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_mul(z, w) returns z * w with infinities handled according
  31  * to C99.
  32  *
  33  * If z and w are both finite, _X_cplx_mul(z, w) delivers the complex
  34  * product according to the usual formula: let a = Re(z), b = Im(z),
  35  * c = Re(w), and d = Im(w); then _X_cplx_mul(z, w) delivers x + I * y
  36  * where x = a * c - b * d and y = a * d + b * c.  Note that if both
  37  * ac and bd overflow, then at least one of ad or bc must also over-
  38  * flow, and vice versa, so that if one component of the product is
  39  * NaN, the other is infinite.  (Such a value is considered infinite
  40  * according to C99.)
  41  *
  42  * If one of z or w is infinite and the other is either finite nonzero
  43  * or infinite, _X_cplx_mul delivers an infinite result.  If one factor
  44  * is infinite and the other is zero, _X_cplx_mul delivers NaN + I * NaN.
  45  * C99 doesn't specify the latter case.
  46  *
  47  * C99 also doesn't specify what should happen if either z or w is a
  48  * complex NaN (i.e., neither finite nor infinite).  This implementation
  49  * delivers NaN + I * NaN in this case.
  50  *
  51  * This implementation can raise spurious underflow, overflow, invalid
  52  * operation, and inexact exceptions.  C99 allows this.
  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 testinfl(long double x)
  71 {
  72         union {
  73                 int             i[3];
  74                 long double     e;
  75         } xx;
  76
  77         xx.e = x;
  78         if ((xx.i[2] & 0x7fff) != 0x7fff || ((xx.i[1] << 1) | xx.i[0]) != 0)
  79                 return (0);
  80         return (1 | ((xx.i[2] << 16) >> 31));
  81 }
  82
  83 long double _Complex
  84 _X_cplx_mul(long double _Complex z, long double _Complex w)
  85 {
  86         long double _Complex    v;
  87         long double             a, b, c, d, x, y;
  88         int                     recalc, i, j;
  89
  90         /*
  91          * The following is equivalent to
  92          *
  93          *  a = creall(z); b = cimagl(z);
  94          *  c = creall(w); d = cimagl(w);
  95          */
  96         a = ((long double *)&z)[0];
  97         b = ((long double *)&z)[1];
  98         c = ((long double *)&w)[0];
  99         d = ((long double *)&w)[1];
 100
 101         x = a * c - b * d;
 102         y = a * d + b * c;
 103
 104         if (x != x && y != y) {
 105                 /*
 106                  * Both x and y are NaN, so z and w can't both be finite.
 107                  * If at least one of z or w is a complex NaN, and neither
 108                  * is infinite, then we might as well deliver NaN + I * NaN.
 109                  * So the only cases to check are when one of z or w is
 110                  * infinite.
 111                  */
 112                 recalc = 0;
 113                 i = testinfl(a);
 114                 j = testinfl(b);
 115                 if (i | j) { /* z is infinite */
 116                         /* "factor out" infinity */
 117                         a = i;
 118                         b = j;
 119                         recalc = 1;
 120                 }
 121                 i = testinfl(c);
 122                 j = testinfl(d);
 123                 if (i | j) { /* w is infinite */
 124                         /* "factor out" infinity */
 125                         c = i;
 126                         d = j;
 127                         recalc = 1;
 128                 }
 129                 if (recalc) {
 130                         x = inf.f * (a * c - b * d);
 131                         y = inf.f * (a * d + b * c);
 132                 }
 133         }
 134
 135         /*
 136          * The following is equivalent to
 137          *
 138          *  return x + I * y;
 139          */
 140         ((long double *)&v)[0] = x;
 141         ((long double *)&v)[1] = y;
 142         return (v);
 143 }