← Back to branch summary

~ubuntu-branches/ubuntu/trusty/eglibc/trusty

« back to all changes in this revision

Viewing changes to sysdeps/ieee754/dbl-64/s_fma.c

  • Committer: Package Import Robot
  • Author(s): Adam Conrad
  • Date: 2013-01-10 18:39:35 UTC
  • mfrom: (1.5.2) (4.4.24 experimental)
  • Revision ID: package-import@ubuntu.com-20130110183935-afsgfxkmg7wk5eaj
Tags: 2.17-0ubuntu1
https://launchpad.net/bugs/1085342
https://launchpad.net/bugs/613662
* Merge with Debian, bringing in a new upstream and many small fixes:
  - patches/any/cvs-malloc-deadlock.diff: Dropped, merged upstream.
  - patches/ubuntu/lddebug-scopes.diff: Rebase for upstream changes.
  - patches/ubuntu/local-CVE-2012-3406.diff: Rebased against upstream.
  - patches/ubuntu/no-asm-mtune-i686.diff: Fixed in recent binutils.
* This upstream merge fixes a nasty hang in pulseaudio (LP: #1085342)
* Bump MIN_KERNEL_SUPPORTED to 2.6.32 on ARM, now that we no longer
  have to support shonky 2.6.31 kernels on imx51 babbage builders.
* Drop patches/ubuntu/local-disable-nscd-host-caching.diff, as these
  issues were apparently resolved upstream a while ago (LP: #613662)
* Fix the compiled-in bug URL to point to launchpad.net, not Debian.

Show diffs side-by-side

added added

removed removed

Lines of Context:
sysdeps/ieee754/dbl-64/s_fma.c
22
22
#include <fenv.h>
23
23
#include <ieee754.h>
24
24
#include <math_private.h>
 
25
#include <tininess.h>
25
26
 
26
27
/* This implementation uses rounding to odd to avoid problems with
27
28
   double rounding.  See a paper by Boldo and Melquiond:
49
50
          && u.ieee.exponent != 0x7ff
50
51
          && v.ieee.exponent != 0x7ff)
51
52
        return (z + x) + y;
52
 
      /* If x or y or z is Inf/NaN, or if fma will certainly overflow,
53
 
         or if x * y is less than half of DBL_DENORM_MIN,
54
 
         compute as x * y + z.  */
 
53
      /* If z is zero and x are y are nonzero, compute the result
 
54
         as x * y to avoid the wrong sign of a zero result if x * y
 
55
         underflows to 0.  */
 
56
      if (z == 0 && x != 0 && y != 0)
 
57
        return x * y;
 
58
      /* If x or y or z is Inf/NaN, or if x * y is zero, compute as
 
59
         x * y + z.  */
55
60
      if (u.ieee.exponent == 0x7ff
56
61
          || v.ieee.exponent == 0x7ff
57
62
          || w.ieee.exponent == 0x7ff
58
 
          || u.ieee.exponent + v.ieee.exponent
59
 
             > 0x7ff + IEEE754_DOUBLE_BIAS
60
 
          || u.ieee.exponent + v.ieee.exponent
61
 
             < IEEE754_DOUBLE_BIAS - DBL_MANT_DIG - 2)
 
63
          || x == 0
 
64
          || y == 0)
62
65
        return x * y + z;
 
66
      /* If fma will certainly overflow, compute as x * y.  */
 
67
      if (u.ieee.exponent + v.ieee.exponent > 0x7ff + IEEE754_DOUBLE_BIAS)
 
68
        return x * y;
 
69
      /* If x * y is less than 1/4 of DBL_DENORM_MIN, neither the
 
70
         result nor whether there is underflow depends on its exact
 
71
         value, only on its sign.  */
 
72
      if (u.ieee.exponent + v.ieee.exponent
 
73
          < IEEE754_DOUBLE_BIAS - DBL_MANT_DIG - 2)
 
74
        {
 
75
          int neg = u.ieee.negative ^ v.ieee.negative;
 
76
          double tiny = neg ? -0x1p-1074 : 0x1p-1074;
 
77
          if (w.ieee.exponent >= 3)
 
78
            return tiny + z;
 
79
          /* Scaling up, adding TINY and scaling down produces the
 
80
             correct result, because in round-to-nearest mode adding
 
81
             TINY has no effect and in other modes double rounding is
 
82
             harmless.  But it may not produce required underflow
 
83
             exceptions.  */
 
84
          v.d = z * 0x1p54 + tiny;
 
85
          if (TININESS_AFTER_ROUNDING
 
86
              ? v.ieee.exponent < 55
 
87
              : (w.ieee.exponent == 0
 
88
                 || (w.ieee.exponent == 1
 
89
                     && w.ieee.negative != neg
 
90
                     && w.ieee.mantissa1 == 0
 
91
                     && w.ieee.mantissa0 == 0)))
 
92
            {
 
93
              volatile double force_underflow = x * y;
 
94
              (void) force_underflow;
 
95
            }
 
96
          return v.d * 0x1p-54;
 
97
        }
63
98
      if (u.ieee.exponent + v.ieee.exponent
64
99
          >= 0x7ff + IEEE754_DOUBLE_BIAS - DBL_MANT_DIG)
65
100
        {
79
114
        {
80
115
          /* Similarly.
81
116
             If z exponent is very large and x and y exponents are
82
 
             very small, it doesn't matter if we don't adjust it.  */
83
 
          if (u.ieee.exponent > v.ieee.exponent)
 
117
             very small, adjust them up to avoid spurious underflows,
 
118
             rather than down.  */
 
119
          if (u.ieee.exponent + v.ieee.exponent
 
120
              <= IEEE754_DOUBLE_BIAS + DBL_MANT_DIG)
 
121
            {
 
122
              if (u.ieee.exponent > v.ieee.exponent)
 
123
                u.ieee.exponent += 2 * DBL_MANT_DIG + 2;
 
124
              else
 
125
                v.ieee.exponent += 2 * DBL_MANT_DIG + 2;
 
126
            }
 
127
          else if (u.ieee.exponent > v.ieee.exponent)
84
128
            {
85
129
              if (u.ieee.exponent > DBL_MANT_DIG)
86
130
                u.ieee.exponent -= DBL_MANT_DIG;
110
154
                  <= IEEE754_DOUBLE_BIAS + DBL_MANT_DIG) */
111
155
        {
112
156
          if (u.ieee.exponent > v.ieee.exponent)
113
 
            u.ieee.exponent += 2 * DBL_MANT_DIG;
 
157
            u.ieee.exponent += 2 * DBL_MANT_DIG + 2;
114
158
          else
115
 
            v.ieee.exponent += 2 * DBL_MANT_DIG;
116
 
          if (w.ieee.exponent <= 4 * DBL_MANT_DIG + 4)
 
159
            v.ieee.exponent += 2 * DBL_MANT_DIG + 2;
 
160
          if (w.ieee.exponent <= 4 * DBL_MANT_DIG + 6)
117
161
            {
118
162
              if (w.ieee.exponent)
119
 
                w.ieee.exponent += 2 * DBL_MANT_DIG;
 
163
                w.ieee.exponent += 2 * DBL_MANT_DIG + 2;
120
164
              else
121
 
                w.d *= 0x1p106;
 
165
                w.d *= 0x1p108;
122
166
              adjust = -1;
123
167
            }
124
168
          /* Otherwise x * y should just affect inexact
128
172
      y = v.d;
129
173
      z = w.d;
130
174
    }
 
175
 
 
176
  /* Ensure correct sign of exact 0 + 0.  */
 
177
  if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0))
 
178
    return x * y + z;
 
179
 
 
180
  fenv_t env;
 
181
  libc_feholdexcept_setround (&env, FE_TONEAREST);
 
182
 
131
183
  /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
132
184
#define C ((1 << (DBL_MANT_DIG + 1) / 2) + 1)
133
185
  double x1 = x * C;
146
198
  t1 = m1 - t1;
147
199
  t2 = z - t2;
148
200
  double a2 = t1 + t2;
149
 
 
150
 
  fenv_t env;
151
 
  libc_feholdexcept_setround (&env, FE_TOWARDZERO);
 
201
  feclearexcept (FE_INEXACT);
 
202
 
 
203
  /* If the result is an exact zero, ensure it has the correct
 
204
     sign.  */
 
205
  if (a1 == 0 && m2 == 0)
 
206
    {
 
207
      libc_feupdateenv (&env);
 
208
      /* Ensure that round-to-nearest value of z + m1 is not
 
209
         reused.  */
 
210
      asm volatile ("" : "=m" (z) : "m" (z));
 
211
      return z + m1;
 
212
    }
 
213
 
 
214
  libc_fesetround (FE_TOWARDZERO);
152
215
 
153
216
  /* Perform m2 + a2 addition with round to odd.  */
154
217
  u.d = a2 + m2;
184
247
      /* If a1 + u.d is exact, the only rounding happens during
185
248
         scaling down.  */
186
249
      if (j == 0)
187
 
        return v.d * 0x1p-106;
 
250
        return v.d * 0x1p-108;
188
251
      /* If result rounded to zero is not subnormal, no double
189
252
         rounding will occur.  */
190
 
      if (v.ieee.exponent > 106)
191
 
        return (a1 + u.d) * 0x1p-106;
192
 
      /* If v.d * 0x1p-106 with round to zero is a subnormal above
193
 
         or equal to DBL_MIN / 2, then v.d * 0x1p-106 shifts mantissa
 
253
      if (v.ieee.exponent > 108)
 
254
        return (a1 + u.d) * 0x1p-108;
 
255
      /* If v.d * 0x1p-108 with round to zero is a subnormal above
 
256
         or equal to DBL_MIN / 2, then v.d * 0x1p-108 shifts mantissa
194
257
         down just by 1 bit, which means v.ieee.mantissa1 |= j would
195
258
         change the round bit, not sticky or guard bit.
196
 
         v.d * 0x1p-106 never normalizes by shifting up,
 
259
         v.d * 0x1p-108 never normalizes by shifting up,
197
260
         so round bit plus sticky bit should be already enough
198
261
         for proper rounding.  */
199
 
      if (v.ieee.exponent == 106)
 
262
      if (v.ieee.exponent == 108)
200
263
        {
 
264
          /* If the exponent would be in the normal range when
 
265
             rounding to normal precision with unbounded exponent
 
266
             range, the exact result is known and spurious underflows
 
267
             must be avoided on systems detecting tininess after
 
268
             rounding.  */
 
269
          if (TININESS_AFTER_ROUNDING)
 
270
            {
 
271
              w.d = a1 + u.d;
 
272
              if (w.ieee.exponent == 109)
 
273
                return w.d * 0x1p-108;
 
274
            }
201
275
          /* v.ieee.mantissa1 & 2 is LSB bit of the result before rounding,
202
276
             v.ieee.mantissa1 & 1 is the round bit and j is our sticky
203
 
             bit.  In round-to-nearest 001 rounds down like 00,
204
 
             011 rounds up, even though 01 rounds down (thus we need
205
 
             to adjust), 101 rounds down like 10 and 111 rounds up
206
 
             like 11.  */
207
 
          if ((v.ieee.mantissa1 & 3) == 1)
208
 
            {
209
 
              v.d *= 0x1p-106;
210
 
              if (v.ieee.negative)
211
 
                return v.d - 0x1p-1074 /* __DBL_DENORM_MIN__ */;
212
 
              else
213
 
                return v.d + 0x1p-1074 /* __DBL_DENORM_MIN__ */;
214
 
            }
215
 
          else
216
 
            return v.d * 0x1p-106;
 
277
             bit.  */
 
278
          w.d = 0.0;
 
279
          w.ieee.mantissa1 = ((v.ieee.mantissa1 & 3) << 1) | j;
 
280
          w.ieee.negative = v.ieee.negative;
 
281
          v.ieee.mantissa1 &= ~3U;
 
282
          v.d *= 0x1p-108;
 
283
          w.d *= 0x1p-2;
 
284
          return v.d + w.d;
217
285
        }
218
286
      v.ieee.mantissa1 |= j;
219
 
      return v.d * 0x1p-106;
 
287
      return v.d * 0x1p-108;
220
288
    }
221
289
}
222
290
#ifndef __fma