~ubuntu-branches/ubuntu/wily/musl/wily

« back to all changes in this revision

Viewing changes to src/math/fma.c

  • Committer: Package Import Robot
  • Author(s): Kevin Bortis
  • Date: 2013-09-27 23:47:18 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20130927234718-a96bgcnvzx5buf60
Tags: 0.9.14-1
* Import upstream version 0.9.14
* Only build on fully supported architectures
* Point to new homepage in control file (Closes: #724277)
* Revorked debian/rules
* Solved possible problem with postrm script (Closes: #724247)

Show diffs side-by-side

added added

removed removed

Lines of Context:
2
2
#include "libm.h"
3
3
 
4
4
#if LDBL_MANT_DIG==64 && LDBL_MAX_EXP==16384
5
 
union ld80 {
6
 
        long double x;
7
 
        struct {
8
 
                uint64_t m;
9
 
                uint16_t e : 15;
10
 
                uint16_t s : 1;
11
 
                uint16_t pad;
12
 
        } bits;
13
 
};
14
 
 
15
5
/* exact add, assumes exponent_x >= exponent_y */
16
6
static void add(long double *hi, long double *lo, long double x, long double y)
17
7
{
45
35
*/
46
36
static long double adjust(long double hi, long double lo)
47
37
{
48
 
        union ld80 uhi, ulo;
 
38
        union ldshape uhi, ulo;
49
39
 
50
40
        if (lo == 0)
51
41
                return hi;
52
 
        uhi.x = hi;
53
 
        if (uhi.bits.m & 0x3ff)
 
42
        uhi.f = hi;
 
43
        if (uhi.i.m & 0x3ff)
54
44
                return hi;
55
 
        ulo.x = lo;
56
 
        if (uhi.bits.s == ulo.bits.s)
57
 
                uhi.bits.m++;
 
45
        ulo.f = lo;
 
46
        if ((uhi.i.se & 0x8000) == (ulo.i.se & 0x8000))
 
47
                uhi.i.m++;
58
48
        else {
59
 
                uhi.bits.m--;
60
49
                /* handle underflow and take care of ld80 implicit msb */
61
 
                if (uhi.bits.m == (uint64_t)-1/2) {
62
 
                        uhi.bits.m *= 2;
63
 
                        uhi.bits.e--;
 
50
                if (uhi.i.m << 1 == 0) {
 
51
                        uhi.i.m = 0;
 
52
                        uhi.i.se--;
64
53
                }
 
54
                uhi.i.m--;
65
55
        }
66
 
        return uhi.x;
 
56
        return uhi.f;
67
57
}
68
58
 
69
59
/* adjusted add so the result is correct when rounded to double (or less) precision */
82
72
 
83
73
static int getexp(long double x)
84
74
{
85
 
        union ld80 u;
86
 
        u.x = x;
87
 
        return u.bits.e;
 
75
        union ldshape u;
 
76
        u.f = x;
 
77
        return u.i.se & 0x7fff;
88
78
}
89
79
 
90
80
double fma(double x, double y, double z)
242
232
static inline double add_adjusted(double a, double b)
243
233
{
244
234
        struct dd sum;
245
 
        uint64_t hibits, lobits;
 
235
        union {double f; uint64_t i;} uhi, ulo;
246
236
 
247
237
        sum = dd_add(a, b);
248
238
        if (sum.lo != 0) {
249
 
                EXTRACT_WORD64(hibits, sum.hi);
250
 
                if ((hibits & 1) == 0) {
 
239
                uhi.f = sum.hi;
 
240
                if ((uhi.i & 1) == 0) {
251
241
                        /* hibits += (int)copysign(1.0, sum.hi * sum.lo) */
252
 
                        EXTRACT_WORD64(lobits, sum.lo);
253
 
                        hibits += 1 - ((hibits ^ lobits) >> 62);
254
 
                        INSERT_WORD64(sum.hi, hibits);
 
242
                        ulo.f = sum.lo;
 
243
                        uhi.i += 1 - ((uhi.i ^ ulo.i) >> 62);
 
244
                        sum.hi = uhi.f;
255
245
                }
256
246
        }
257
247
        return (sum.hi);
265
255
static inline double add_and_denormalize(double a, double b, int scale)
266
256
{
267
257
        struct dd sum;
268
 
        uint64_t hibits, lobits;
 
258
        union {double f; uint64_t i;} uhi, ulo;
269
259
        int bits_lost;
270
260
 
271
261
        sum = dd_add(a, b);
281
271
         * break the ties manually.
282
272
         */
283
273
        if (sum.lo != 0) {
284
 
                EXTRACT_WORD64(hibits, sum.hi);
285
 
                bits_lost = -((int)(hibits >> 52) & 0x7ff) - scale + 1;
286
 
                if (bits_lost != 1 ^ (int)(hibits & 1)) {
 
274
                uhi.f = sum.hi;
 
275
                bits_lost = -((int)(uhi.i >> 52) & 0x7ff) - scale + 1;
 
276
                if (bits_lost != 1 ^ (int)(uhi.i & 1)) {
287
277
                        /* hibits += (int)copysign(1.0, sum.hi * sum.lo) */
288
 
                        EXTRACT_WORD64(lobits, sum.lo);
289
 
                        hibits += 1 - (((hibits ^ lobits) >> 62) & 2);
290
 
                        INSERT_WORD64(sum.hi, hibits);
 
278
                        ulo.f = sum.lo;
 
279
                        uhi.i += 1 - (((uhi.i ^ ulo.i) >> 62) & 2);
 
280
                        sum.hi = uhi.f;
291
281
                }
292
282
        }
293
283
        return scalbn(sum.hi, scale);