usr/src/lib/libc/i386/fp/_F_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  * _F_cplx_mul(z, w) returns z * w with infinities handled according
  31  * to C99.
  32  *
  33  * If z and w are both finite, _F_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 _F_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, _F_cplx_mul delivers an infinite result.  If one factor
  42  * is infinite and the other is zero, _F_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 testinff(float x)
  69 {
  70         union {
  71                 int     i;
  72                 float   f;
  73         } xx;
  74
  75         xx.f = x;
  76         return ((((xx.i << 1) - 0xff000000) == 0)? (1 | (xx.i >> 31)) : 0);
  77 }
  78
  79 float _Complex
  80 _F_cplx_mul(float _Complex z, float _Complex w)
  81 {
  82         float _Complex  v;
  83         float           a, b, c, d;
  84         long double     x, y;
  85         int             recalc, i, j;
  86
  87         /*
  88          * The following is equivalent to
  89          *
  90          *  a = crealf(z); b = cimagf(z);
  91          *  c = crealf(w); d = cimagf(w);
  92          */
  93         a = ((float *)&z)[0];
  94         b = ((float *)&z)[1];
  95         c = ((float *)&w)[0];
  96         d = ((float *)&w)[1];
  97
  98         x = (long double)a * c - (long double)b * d;
  99         y = (long double)a * d + (long double)b * c;
 100
 101         if (x != x && y != y) {
 102                 /*
 103                  * Both x and y are NaN, so z and w can't both be finite.
 104                  * If at least one of z or w is a complex NaN, and neither
 105                  * is infinite, then we might as well deliver NaN + I * NaN.
 106                  * So the only cases to check are when one of z or w is
 107                  * infinite.
 108                  */
 109                 recalc = 0;
 110                 i = testinff(a);
 111                 j = testinff(b);
 112                 if (i | j) { /* z is infinite */
 113                         /* "factor out" infinity */
 114                         a = i;
 115                         b = j;
 116                         recalc = 1;
 117                 }
 118                 i = testinff(c);
 119                 j = testinff(d);
 120                 if (i | j) { /* w is infinite */
 121                         /* "factor out" infinity */
 122                         c = i;
 123                         d = j;
 124                         recalc = 1;
 125                 }
 126                 if (recalc) {
 127                         x = inf.f * ((long double)a * c - (long double)b * d);
 128                         y = inf.f * ((long double)a * d + (long double)b * c);
 129                 }
 130         }
 131
 132         /*
 133          * The following is equivalent to
 134          *
 135          *  return x + I * y;
 136          */
 137         ((float *)&v)[0] = (float)x;
 138         ((float *)&v)[1] = (float)y;
 139         return (v);
 140 }