← Back to branch summary

~vorlon/ubuntu/natty/eglibc/multiarch

« back to all changes in this revision

Viewing changes to sysdeps/ieee754/ldbl-96/s_fma.c

  • Committer: Steve Langasek
  • Date: 2011-02-18 21:18:44 UTC
  • mfrom: (103.1.7 eglibc)
  • Revision ID: steve.langasek@linaro.org-20110218211844-lodmi8b1qhyq3f3x
Tags: 2.13~pre1-0ubuntu1+multiarch.1
merge from natty

Show diffs side-by-side

added added

removed removed

Lines of Context:
sysdeps/ieee754/ldbl-96/s_fma.c
 
1
/* Compute x * y + z as ternary operation.
 
2
   Copyright (C) 2010 Free Software Foundation, Inc.
 
3
   This file is part of the GNU C Library.
 
4
   Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.
 
5
 
 
6
   The GNU C Library is free software; you can redistribute it and/or
 
7
   modify it under the terms of the GNU Lesser General Public
 
8
   License as published by the Free Software Foundation; either
 
9
   version 2.1 of the License, or (at your option) any later version.
 
10
 
 
11
   The GNU C Library is distributed in the hope that it will be useful,
 
12
   but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
14
   Lesser General Public License for more details.
 
15
 
 
16
   You should have received a copy of the GNU Lesser General Public
 
17
   License along with the GNU C Library; if not, write to the Free
 
18
   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
 
19
   02111-1307 USA.  */
 
20
 
 
21
#include <float.h>
 
22
#include <math.h>
 
23
#include <fenv.h>
 
24
#include <ieee754.h>
 
25
 
 
26
/* This implementation uses rounding to odd to avoid problems with
 
27
   double rounding.  See a paper by Boldo and Melquiond:
 
28
   http://www.lri.fr/~melquion/doc/08-tc.pdf  */
 
29
 
 
30
double
 
31
__fma (double x, double y, double z)
 
32
{
 
33
  if (__builtin_expect (isinf (z), 0))
 
34
    {
 
35
      /* If z is Inf, but x and y are finite, the result should be
 
36
         z rather than NaN.  */
 
37
      if (finite (x) && finite (y))
 
38
        return (z + x) + y;
 
39
      return (x * y) + z;
 
40
    }
 
41
 
 
42
  /* Multiplication m1 + m2 = x * y using Dekker's algorithm.  */
 
43
#define C ((1ULL << (LDBL_MANT_DIG + 1) / 2) + 1)
 
44
  long double x1 = (long double) x * C;
 
45
  long double y1 = (long double) y * C;
 
46
  long double m1 = (long double) x * y;
 
47
  x1 = (x - x1) + x1;
 
48
  y1 = (y - y1) + y1;
 
49
  long double x2 = x - x1;
 
50
  long double y2 = y - y1;
 
51
  long double m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;
 
52
 
 
53
  /* Addition a1 + a2 = z + m1 using Knuth's algorithm.  */
 
54
  long double a1 = z + m1;
 
55
  long double t1 = a1 - z;
 
56
  long double t2 = a1 - t1;
 
57
  t1 = m1 - t1;
 
58
  t2 = z - t2;
 
59
  long double a2 = t1 + t2;
 
60
 
 
61
  fenv_t env;
 
62
  feholdexcept (&env);
 
63
  fesetround (FE_TOWARDZERO);
 
64
  /* Perform m2 + a2 addition with round to odd.  */
 
65
  a2 = a2 + m2;
 
66
 
 
67
  /* Add that to a1 again using rounding to odd.  */
 
68
  union ieee854_long_double u;
 
69
  u.d = a1 + a2;
 
70
  if ((u.ieee.mantissa1 & 1) == 0 && u.ieee.exponent != 0x7fff)
 
71
    u.ieee.mantissa1 |= fetestexcept (FE_INEXACT) != 0;
 
72
  feupdateenv (&env);
 
73
 
 
74
  /* Add finally round to double precision.  */
 
75
  return u.d;
 
76
}
 
77
#ifndef __fma
 
78
weak_alias (__fma, fma)
 
79
#endif