Mercurial > hg > CbC > CbC_gcc

annotate libquadmath/math/fmaq.c @ 158:494b0b89df80 default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
...
author Shinji KONO <kono@ie.u-ryukyu.ac.jp>
date Mon, 25 May 2020 18:13:55 +0900
parents 1830386684a0
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
1 /* Compute x * y + z as ternary operation.
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
2 Copyright (C) 2010-2018 Free Software Foundation, Inc.
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
3 This file is part of the GNU C Library.
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
4 Contributed by Jakub Jelinek <jakub@redhat.com>, 2010.
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
5
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
6 The GNU C Library is free software; you can redistribute it and/or
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
7 modify it under the terms of the GNU Lesser General Public
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
8 License as published by the Free Software Foundation; either
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
9 version 2.1 of the License, or (at your option) any later version.
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
10
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
11 The GNU C Library is distributed in the hope that it will be useful,
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
14 Lesser General Public License for more details.
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
15
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
16 You should have received a copy of the GNU Lesser General Public
111
kono
parents: 68
diff changeset
17 License along with the GNU C Library; if not, see
kono
parents: 68
diff changeset
18 <http://www.gnu.org/licenses/>. */
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
19
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
20 #include "quadmath-imp.h"
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
21
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
22 /* This implementation uses rounding to odd to avoid problems with
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
23 double rounding. See a paper by Boldo and Melquiond:
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
24 http://www.lri.fr/~melquion/doc/08-tc.pdf */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
25
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
26 __float128
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
27 fmaq (__float128 x, __float128 y, __float128 z)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
28 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
29 ieee854_float128 u, v, w;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
30 int adjust = 0;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
31 u.value = x;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
32 v.value = y;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
33 w.value = z;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
34 if (__builtin_expect (u.ieee.exponent + v.ieee.exponent
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
35 >= 0x7fff + IEEE854_FLOAT128_BIAS
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
36 - FLT128_MANT_DIG, 0)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
37 || __builtin_expect (u.ieee.exponent >= 0x7fff - FLT128_MANT_DIG, 0)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
38 || __builtin_expect (v.ieee.exponent >= 0x7fff - FLT128_MANT_DIG, 0)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
39 || __builtin_expect (w.ieee.exponent >= 0x7fff - FLT128_MANT_DIG, 0)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
40 || __builtin_expect (u.ieee.exponent + v.ieee.exponent
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
41 <= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG, 0))
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
42 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
43 /* If z is Inf, but x and y are finite, the result should be
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
44 z rather than NaN. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
45 if (w.ieee.exponent == 0x7fff
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
46 && u.ieee.exponent != 0x7fff
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
47 && v.ieee.exponent != 0x7fff)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
48 return (z + x) + y;
111
kono
parents: 68
diff changeset
49 /* If z is zero and x are y are nonzero, compute the result
kono
parents: 68
diff changeset
50 as x * y to avoid the wrong sign of a zero result if x * y
kono
parents: 68
diff changeset
51 underflows to 0. */
kono
parents: 68
diff changeset
52 if (z == 0 && x != 0 && y != 0)
kono
parents: 68
diff changeset
53 return x * y;
kono
parents: 68
diff changeset
54 /* If x or y or z is Inf/NaN, or if x * y is zero, compute as
kono
parents: 68
diff changeset
55 x * y + z. */
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
56 if (u.ieee.exponent == 0x7fff
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
57 || v.ieee.exponent == 0x7fff
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
58 || w.ieee.exponent == 0x7fff
111
kono
parents: 68
diff changeset
59 || x == 0
kono
parents: 68
diff changeset
60 || y == 0)
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
61 return x * y + z;
111
kono
parents: 68
diff changeset
62 /* If fma will certainly overflow, compute as x * y. */
kono
parents: 68
diff changeset
63 if (u.ieee.exponent + v.ieee.exponent
kono
parents: 68
diff changeset
64 > 0x7fff + IEEE854_FLOAT128_BIAS)
kono
parents: 68
diff changeset
65 return x * y;
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
66 /* If x * y is less than 1/4 of FLT128_TRUE_MIN, neither the
111
kono
parents: 68
diff changeset
67 result nor whether there is underflow depends on its exact
kono
parents: 68
diff changeset
68 value, only on its sign. */
kono
parents: 68
diff changeset
69 if (u.ieee.exponent + v.ieee.exponent
kono
parents: 68
diff changeset
70 < IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG - 2)
kono
parents: 68
diff changeset
71 {
kono
parents: 68
diff changeset
72 int neg = u.ieee.negative ^ v.ieee.negative;
kono
parents: 68
diff changeset
73 __float128 tiny = neg ? -0x1p-16494Q : 0x1p-16494Q;
kono
parents: 68
diff changeset
74 if (w.ieee.exponent >= 3)
kono
parents: 68
diff changeset
75 return tiny + z;
kono
parents: 68
diff changeset
76 /* Scaling up, adding TINY and scaling down produces the
kono
parents: 68
diff changeset
77 correct result, because in round-to-nearest mode adding
kono
parents: 68
diff changeset
78 TINY has no effect and in other modes double rounding is
kono
parents: 68
diff changeset
79 harmless. But it may not produce required underflow
kono
parents: 68
diff changeset
80 exceptions. */
kono
parents: 68
diff changeset
81 v.value = z * 0x1p114Q + tiny;
kono
parents: 68
diff changeset
82 if (TININESS_AFTER_ROUNDING
kono
parents: 68
diff changeset
83 ? v.ieee.exponent < 115
kono
parents: 68
diff changeset
84 : (w.ieee.exponent == 0
kono
parents: 68
diff changeset
85 || (w.ieee.exponent == 1
kono
parents: 68
diff changeset
86 && w.ieee.negative != neg
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
87 && w.ieee.mantissa3 == 0
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
88 && w.ieee.mantissa2 == 0
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
89 && w.ieee.mantissa1 == 0
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
90 && w.ieee.mantissa0 == 0)))
111
kono
parents: 68
diff changeset
91 {
kono
parents: 68
diff changeset
92 __float128 force_underflow = x * y;
kono
parents: 68
diff changeset
93 math_force_eval (force_underflow);
kono
parents: 68
diff changeset
94 }
kono
parents: 68
diff changeset
95 return v.value * 0x1p-114Q;
kono
parents: 68
diff changeset
96 }
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
97 if (u.ieee.exponent + v.ieee.exponent
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
98 >= 0x7fff + IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
99 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
100 /* Compute 1p-113 times smaller result and multiply
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
101 at the end. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
102 if (u.ieee.exponent > v.ieee.exponent)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
103 u.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
104 else
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
105 v.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
106 /* If x + y exponent is very large and z exponent is very small,
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
107 it doesn't matter if we don't adjust it. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
108 if (w.ieee.exponent > FLT128_MANT_DIG)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
109 w.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
110 adjust = 1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
111 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
112 else if (w.ieee.exponent >= 0x7fff - FLT128_MANT_DIG)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
113 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
114 /* Similarly.
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
115 If z exponent is very large and x and y exponents are
111
kono
parents: 68
diff changeset
116 very small, adjust them up to avoid spurious underflows,
kono
parents: 68
diff changeset
117 rather than down. */
kono
parents: 68
diff changeset
118 if (u.ieee.exponent + v.ieee.exponent
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
119 <= IEEE854_FLOAT128_BIAS + 2 * FLT128_MANT_DIG)
111
kono
parents: 68
diff changeset
120 {
kono
parents: 68
diff changeset
121 if (u.ieee.exponent > v.ieee.exponent)
kono
parents: 68
diff changeset
122 u.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
kono
parents: 68
diff changeset
123 else
kono
parents: 68
diff changeset
124 v.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
kono
parents: 68
diff changeset
125 }
kono
parents: 68
diff changeset
126 else if (u.ieee.exponent > v.ieee.exponent)
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
127 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
128 if (u.ieee.exponent > FLT128_MANT_DIG)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
129 u.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
130 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
131 else if (v.ieee.exponent > FLT128_MANT_DIG)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
132 v.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
133 w.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
134 adjust = 1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
135 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
136 else if (u.ieee.exponent >= 0x7fff - FLT128_MANT_DIG)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
137 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
138 u.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
139 if (v.ieee.exponent)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
140 v.ieee.exponent += FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
141 else
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
142 v.value *= 0x1p113Q;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
143 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
144 else if (v.ieee.exponent >= 0x7fff - FLT128_MANT_DIG)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
145 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
146 v.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
147 if (u.ieee.exponent)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
148 u.ieee.exponent += FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
149 else
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
150 u.value *= 0x1p113Q;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
151 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
152 else /* if (u.ieee.exponent + v.ieee.exponent
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
153 <= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG) */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
154 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
155 if (u.ieee.exponent > v.ieee.exponent)
111
kono
parents: 68
diff changeset
156 u.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
157 else
111
kono
parents: 68
diff changeset
158 v.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
kono
parents: 68
diff changeset
159 if (w.ieee.exponent <= 4 * FLT128_MANT_DIG + 6)
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
160 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
161 if (w.ieee.exponent)
111
kono
parents: 68
diff changeset
162 w.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
163 else
111
kono
parents: 68
diff changeset
164 w.value *= 0x1p228Q;
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
165 adjust = -1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
166 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
167 /* Otherwise x * y should just affect inexact
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
168 and nothing else. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
169 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
170 x = u.value;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
171 y = v.value;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
172 z = w.value;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
173 }
111
kono
parents: 68
diff changeset
174
kono
parents: 68
diff changeset
175 /* Ensure correct sign of exact 0 + 0. */
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
176 if (__glibc_unlikely ((x == 0 || y == 0) && z == 0))
111
kono
parents: 68
diff changeset
177 {
kono
parents: 68
diff changeset
178 x = math_opt_barrier (x);
kono
parents: 68
diff changeset
179 return x * y + z;
kono
parents: 68
diff changeset
180 }
kono
parents: 68
diff changeset
181
kono
parents: 68
diff changeset
182 fenv_t env;
kono
parents: 68
diff changeset
183 feholdexcept (&env);
kono
parents: 68
diff changeset
184 fesetround (FE_TONEAREST);
kono
parents: 68
diff changeset
185
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
186 /* Multiplication m1 + m2 = x * y using Dekker's algorithm. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
187 #define C ((1LL << (FLT128_MANT_DIG + 1) / 2) + 1)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
188 __float128 x1 = x * C;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
189 __float128 y1 = y * C;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
190 __float128 m1 = x * y;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
191 x1 = (x - x1) + x1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
192 y1 = (y - y1) + y1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
193 __float128 x2 = x - x1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
194 __float128 y2 = y - y1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
195 __float128 m2 = (((x1 * y1 - m1) + x1 * y2) + x2 * y1) + x2 * y2;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
196
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
197 /* Addition a1 + a2 = z + m1 using Knuth's algorithm. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
198 __float128 a1 = z + m1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
199 __float128 t1 = a1 - z;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
200 __float128 t2 = a1 - t1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
201 t1 = m1 - t1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
202 t2 = z - t2;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
203 __float128 a2 = t1 + t2;
111
kono
parents: 68
diff changeset
204 /* Ensure the arithmetic is not scheduled after feclearexcept call. */
kono
parents: 68
diff changeset
205 math_force_eval (m2);
kono
parents: 68
diff changeset
206 math_force_eval (a2);
kono
parents: 68
diff changeset
207 feclearexcept (FE_INEXACT);
kono
parents: 68
diff changeset
208
kono
parents: 68
diff changeset
209 /* If the result is an exact zero, ensure it has the correct sign. */
kono
parents: 68
diff changeset
210 if (a1 == 0 && m2 == 0)
kono
parents: 68
diff changeset
211 {
kono
parents: 68
diff changeset
212 feupdateenv (&env);
kono
parents: 68
diff changeset
213 /* Ensure that round-to-nearest value of z + m1 is not reused. */
kono
parents: 68
diff changeset
214 z = math_opt_barrier (z);
kono
parents: 68
diff changeset
215 return z + m1;
kono
parents: 68
diff changeset
216 }
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
217
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
218 fesetround (FE_TOWARDZERO);
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
219 /* Perform m2 + a2 addition with round to odd. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
220 u.value = a2 + m2;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
221
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
222 if (__glibc_likely (adjust == 0))
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
223 {
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
224 if ((u.ieee.mantissa3 & 1) == 0 && u.ieee.exponent != 0x7fff)
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
225 u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0;
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
226 feupdateenv (&env);
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
227 /* Result is a1 + u.value. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
228 return a1 + u.value;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
229 }
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
230 else if (__glibc_likely (adjust > 0))
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
231 {
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
232 if ((u.ieee.mantissa3 & 1) == 0 && u.ieee.exponent != 0x7fff)
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
233 u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0;
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
234 feupdateenv (&env);
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
235 /* Result is a1 + u.value, scaled up. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
236 return (a1 + u.value) * 0x1p113Q;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
237 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
238 else
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
239 {
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
240 if ((u.ieee.mantissa3 & 1) == 0)
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
241 u.ieee.mantissa3 |= fetestexcept (FE_INEXACT) != 0;
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
242 v.value = a1 + u.value;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
243 /* Ensure the addition is not scheduled after fetestexcept call. */
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
244 math_force_eval (v.value);
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
245 int j = fetestexcept (FE_INEXACT) != 0;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
246 feupdateenv (&env);
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
247 /* Ensure the following computations are performed in default rounding
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
248 mode instead of just reusing the round to zero computation. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
249 asm volatile ("" : "=m" (u) : "m" (u));
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
250 /* If a1 + u.value is exact, the only rounding happens during
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
251 scaling down. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
252 if (j == 0)
111
kono
parents: 68
diff changeset
253 return v.value * 0x1p-228Q;
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
254 /* If result rounded to zero is not subnormal, no double
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
255 rounding will occur. */
111
kono
parents: 68
diff changeset
256 if (v.ieee.exponent > 228)
kono
parents: 68
diff changeset
257 return (a1 + u.value) * 0x1p-228Q;
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
258 /* If v.value * 0x1p-228L with round to zero is a subnormal above
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
259 or equal to FLT128_MIN / 2, then v.value * 0x1p-228L shifts mantissa
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
260 down just by 1 bit, which means v.ieee.mantissa3 |= j would
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
261 change the round bit, not sticky or guard bit.
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
262 v.value * 0x1p-228L never normalizes by shifting up,
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
263 so round bit plus sticky bit should be already enough
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
264 for proper rounding. */
111
kono
parents: 68
diff changeset
265 if (v.ieee.exponent == 228)
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
266 {
111
kono
parents: 68
diff changeset
267 /* If the exponent would be in the normal range when
kono
parents: 68
diff changeset
268 rounding to normal precision with unbounded exponent
kono
parents: 68
diff changeset
269 range, the exact result is known and spurious underflows
kono
parents: 68
diff changeset
270 must be avoided on systems detecting tininess after
kono
parents: 68
diff changeset
271 rounding. */
kono
parents: 68
diff changeset
272 if (TININESS_AFTER_ROUNDING)
kono
parents: 68
diff changeset
273 {
kono
parents: 68
diff changeset
274 w.value = a1 + u.value;
kono
parents: 68
diff changeset
275 if (w.ieee.exponent == 229)
kono
parents: 68
diff changeset
276 return w.value * 0x1p-228Q;
kono
parents: 68
diff changeset
277 }
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
278 /* v.ieee.mantissa3 & 2 is LSB bit of the result before rounding,
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
279 v.ieee.mantissa3 & 1 is the round bit and j is our sticky
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
280 bit. */
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
281 w.value = 0;
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
282 w.ieee.mantissa3 = ((v.ieee.mantissa3 & 3) << 1) | j;
111
kono
parents: 68
diff changeset
283 w.ieee.negative = v.ieee.negative;
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
284 v.ieee.mantissa3 &= ~3U;
111
kono
parents: 68
diff changeset
285 v.value *= 0x1p-228Q;
kono
parents: 68
diff changeset
286 w.value *= 0x1p-2Q;
kono
parents: 68
diff changeset
287 return v.value + w.value;
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
288 }
145
1830386684a0 gcc-9.2.0
anatofuz
parents: 111
diff changeset
289 v.ieee.mantissa3 |= j;
111
kono
parents: 68
diff changeset
290 return v.value * 0x1p-228Q;
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
291 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
292 }