usr/src/lib/libc/i386/fp/_D_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  * _D_cplx_mul(z, w) returns z * w with infinities handled according
  31  * to C99.
  32  *
  33  * If z and w are both finite, _D_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 _D_cplx_mul(z, w) delivers x + I * y
  36  * where x = a * c - b * d and y = a * d + b * c.  This implementation
  37  * uses extended precision to form these expressions, so none of the
  38  * intermediate products can overflow.
  39  *
  40  * If one of z or w is infinite and the other is either finite nonzero
  41  * or infinite, _D_cplx_mul delivers an infinite result.  If one factor
  42  * is infinite and the other is zero, _D_cplx_mul delivers NaN + I * NaN.
  43  * C99 doesn't specify the latter case.
  44  *
  45  * C99 also doesn't specify what should happen if either z or w is a
  46  * complex NaN (i.e., neither finite nor infinite).  This implementation
  47  * delivers NaN + I * NaN in this case.
  48  *
  49  * This implementation can raise spurious invalid operation and inexact
  50  * exceptions.  C99 allows this.
  51  */
  52
  53 #if !defined(i386) && !defined(__i386) && !defined(__amd64)
  54 #error This code is for x86 only
  55 #endif
  56
  57 static union {
  58         int     i;
  59         float   f;
  60 } inf = {
  61         0x7f800000
  62 };
  63
  64 /*
  65  * Return +1 if x is +Inf, -1 if x is -Inf, and 0 otherwise
  66  */
  67 static int
  68 testinf(double x)
  69 {
  70         union {
  71                 int     i[2];
  72                 double  d;
  73         } xx;
  74
  75         xx.d = x;
  76         return (((((xx.i[1] << 1) - 0xffe00000) | xx.i[0]) == 0)?
  77                 (1 | (xx.i[1] >> 31)) : 0);
  78 }
  79
  80 double _Complex
  81 _D_cplx_mul(double _Complex z, double _Complex w)
  82 {
  83         double _Complex v;
  84         double          a, b, c, d;
  85         long double     x, y;
  86         int             recalc, i, j;
  87
  88         /*
  89          * The following is equivalent to
  90          *
  91          *  a = creal(z); b = cimag(z);
  92          *  c = creal(w); d = cimag(w);
  93          */
  94         /* LINTED alignment */
  95         a = ((double *)&z)[0];
  96         /* LINTED alignment */
  97         b = ((double *)&z)[1];
  98         /* LINTED alignment */
  99         c = ((double *)&w)[0];
 100         /* LINTED alignment */
 101         d = ((double *)&w)[1];
 102
 103         x = (long double)a * c - (long double)b * d;
 104         y = (long double)a * d + (long double)b * c;
 105
 106         if (x != x && y != y) {
 107                 /*
 108                  * Both x and y are NaN, so z and w can't both be finite.
 109                  * If at least one of z or w is a complex NaN, and neither
 110                  * is infinite, then we might as well deliver NaN + I * NaN.
 111                  * So the only cases to check are when one of z or w is
 112                  * infinite.
 113                  */
 114                 recalc = 0;
 115                 i = testinf(a);
 116                 j = testinf(b);
 117                 if (i | j) { /* z is infinite */
 118                         /* "factor out" infinity */
 119                         a = i;
 120                         b = j;
 121                         recalc = 1;
 122                 }
 123                 i = testinf(c);
 124                 j = testinf(d);
 125                 if (i | j) { /* w is infinite */
 126                         /* "factor out" infinity */
 127                         c = i;
 128                         d = j;
 129                         recalc = 1;
 130                 }
 131                 if (recalc) {
 132                         x = inf.f * ((long double)a * c - (long double)b * d);
 133                         y = inf.f * ((long double)a * d + (long double)b * c);
 134                 }
 135         }
 136
 137         /*
 138          * The following is equivalent to
 139          *
 140          *  return x + I * y;
 141          */
 142         /* LINTED alignment */
 143         ((double *)&v)[0] = (double)x;
 144         /* LINTED alignment */
 145         ((double *)&v)[1] = (double)y;
 146         return (v);
 147 }