compiler/mlib/s_fmal.c

   1 /*-
   2  * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
   3  * All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   8  * 1. Redistributions of source code must retain the above copyright
   9  *    notice, this list of conditions and the following disclaimer.
  10  * 2. Redistributions in binary form must reproduce the above copyright
  11  *    notice, this list of conditions and the following disclaimer in the
  12  *    documentation and/or other materials provided with the distribution.
  13  *
  14  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  15  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  16  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  17  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  18  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  19  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  20  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  21  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  22  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  23  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  24  * SUCH DAMAGE.
  25  */
  26
  27 #include <sys/cdefs.h>
  28 __FBSDID("$FreeBSD: src/lib/msun/src/s_fmal.c,v 1.3 2007/01/07 07:54:21 das Exp $");
  29
  30 #include <fenv.h>
  31 #include <float.h>
  32 #include <math.h>
  33
  34 /*
  35  * Fused multiply-add: Compute x * y + z with a single rounding error.
  36  *
  37  * We use scaling to avoid overflow/underflow, along with the
  38  * canonical precision-doubling technique adapted from:
  39  *
  40  *      Dekker, T.  A Floating-Point Technique for Extending the
  41  *      Available Precision.  Numer. Math. 18, 224-242 (1971).
  42  */
  43 long double
  44 fmal(long double x, long double y, long double z)
  45 {
  46 #if LDBL_MANT_DIG == 64
  47         static const long double split = 0x1p32L + 1.0;
  48 #elif LDBL_MANT_DIG == 113
  49         static const long double split = 0x1p57L + 1.0;
  50 #endif
  51         long double xs, ys, zs;
  52         long double c, cc, hx, hy, p, q, tx, ty;
  53         long double r, rr, s;
  54         int oround;
  55         int ex, ey, ez;
  56         int spread;
  57
  58         if (z == 0.0)
  59                 return (x * y);
  60         if (x == 0.0 || y == 0.0)
  61                 return (x * y + z);
  62
  63         /* Results of frexp() are undefined for these cases. */
  64         if (!isfinite(x) || !isfinite(y) || !isfinite(z))
  65                 return (x * y + z);
  66
  67         xs = frexpl(x, &ex);
  68         ys = frexpl(y, &ey);
  69         zs = frexpl(z, &ez);
  70         oround = fegetround();
  71         spread = ex + ey - ez;
  72
  73         /*
  74          * If x * y and z are many orders of magnitude apart, the scaling
  75          * will overflow, so we handle these cases specially.  Rounding
  76          * modes other than FE_TONEAREST are painful.
  77          */
  78         if (spread > LDBL_MANT_DIG * 2) {
  79                 fenv_t env;
  80                 feraiseexcept(FE_INEXACT);
  81                 switch(oround) {
  82                 case FE_TONEAREST:
  83                         return (x * y);
  84                 case FE_TOWARDZERO:
  85                         if (x > 0.0 ^ y < 0.0 ^ z < 0.0)
  86                                 return (x * y);
  87                         feholdexcept(&env);
  88                         r = x * y;
  89                         if (!fetestexcept(FE_INEXACT))
  90                                 r = nextafterl(r, 0);
  91                         feupdateenv(&env);
  92                         return (r);
  93                 case FE_DOWNWARD:
  94                         if (z > 0.0)
  95                                 return (x * y);
  96                         feholdexcept(&env);
  97                         r = x * y;
  98                         if (!fetestexcept(FE_INEXACT))
  99                                 r = nextafterl(r, -INFINITY);
 100                         feupdateenv(&env);
 101                         return (r);
 102                 default:        /* FE_UPWARD */
 103                         if (z < 0.0)
 104                                 return (x * y);
 105                         feholdexcept(&env);
 106                         r = x * y;
 107                         if (!fetestexcept(FE_INEXACT))
 108                                 r = nextafterl(r, INFINITY);
 109                         feupdateenv(&env);
 110                         return (r);
 111                 }
 112         }
 113         if (spread < -LDBL_MANT_DIG) {
 114                 feraiseexcept(FE_INEXACT);
 115                 if (!isnormal(z))
 116                         feraiseexcept(FE_UNDERFLOW);
 117                 switch (oround) {
 118                 case FE_TONEAREST:
 119                         return (z);
 120                 case FE_TOWARDZERO:
 121                         if (x > 0.0 ^ y < 0.0 ^ z < 0.0)
 122                                 return (z);
 123                         else
 124                                 return (nextafterl(z, 0));
 125                 case FE_DOWNWARD:
 126                         if (x > 0.0 ^ y < 0.0)
 127                                 return (z);
 128                         else
 129                                 return (nextafterl(z, -INFINITY));
 130                 default:        /* FE_UPWARD */
 131                         if (x > 0.0 ^ y < 0.0)
 132                                 return (nextafterl(z, INFINITY));
 133                         else
 134                                 return (z);
 135                 }
 136         }
 137
 138         /*
 139          * Use Dekker's algorithm to perform the multiplication and
 140          * subsequent addition in twice the machine precision.
 141          * Arrange so that x * y = c + cc, and x * y + z = r + rr.
 142          */
 143         fesetround(FE_TONEAREST);
 144
 145         p = xs * split;
 146         hx = xs - p;
 147         hx += p;
 148         tx = xs - hx;
 149
 150         p = ys * split;
 151         hy = ys - p;
 152         hy += p;
 153         ty = ys - hy;
 154
 155         p = hx * hy;
 156         q = hx * ty + tx * hy;
 157         c = p + q;
 158         cc = p - c + q + tx * ty;
 159
 160         zs = ldexpl(zs, -spread);
 161         r = c + zs;
 162         s = r - c;
 163         rr = (c - (r - s)) + (zs - s) + cc;
 164
 165         spread = ex + ey;
 166         if (spread + ilogbl(r) > -16383) {
 167                 fesetround(oround);
 168                 r = r + rr;
 169         } else {
 170                 /*
 171                  * The result is subnormal, so we round before scaling to
 172                  * avoid double rounding.
 173                  */
 174                 p = ldexpl(copysignl(0x1p-16382L, r), -spread);
 175                 c = r + p;
 176                 s = c - r;
 177                 cc = (r - (c - s)) + (p - s) + rr;
 178                 fesetround(oround);
 179                 r = (c + cc) - p;
 180         }
 181         return (ldexpl(r, spread));
 182 }