Mercurial > hg > CbC > CbC_gcc

annotate libquadmath/math/fmaq.c @ 111:04ced10e8804

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gcc 7
author kono
date Fri, 27 Oct 2017 22:46:09 +0900
parents 561a7518be6b
children 1830386684a0
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.
111
kono
parents: 68
diff changeset
2 Copyright (C) 2010-2017 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 #include <math.h>
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
22 #include <float.h>
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
23 #ifdef HAVE_FENV_H
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
24 # include <fenv.h>
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
25 # if defined HAVE_FEHOLDEXCEPT && defined HAVE_FESETROUND \
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
26 && defined HAVE_FEUPDATEENV && defined HAVE_FETESTEXCEPT \
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
27 && defined FE_TOWARDZERO && defined FE_INEXACT
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
28 # define USE_FENV_H
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
29 # endif
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
30 #endif
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
31
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
32 /* 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
33 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
34 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
35
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
36 __float128
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
37 fmaq (__float128 x, __float128 y, __float128 z)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
38 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
39 ieee854_float128 u, v, w;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
40 int adjust = 0;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
41 u.value = x;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
42 v.value = y;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
43 w.value = z;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
44 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
45 >= 0x7fff + IEEE854_FLOAT128_BIAS
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
46 - FLT128_MANT_DIG, 0)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
47 || __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
48 || __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
49 || __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
50 || __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
51 <= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG, 0))
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
52 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
53 /* 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
54 z rather than NaN. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
55 if (w.ieee.exponent == 0x7fff
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
56 && 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 return (z + x) + y;
111
kono
parents: 68
diff changeset
59 /* If z is zero and x are y are nonzero, compute the result
kono
parents: 68
diff changeset
60 as x * y to avoid the wrong sign of a zero result if x * y
kono
parents: 68
diff changeset
61 underflows to 0. */
kono
parents: 68
diff changeset
62 if (z == 0 && x != 0 && y != 0)
kono
parents: 68
diff changeset
63 return x * y;
kono
parents: 68
diff changeset
64 /* If x or y or z is Inf/NaN, or if x * y is zero, compute as
kono
parents: 68
diff changeset
65 x * y + z. */
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
66 if (u.ieee.exponent == 0x7fff
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
67 || v.ieee.exponent == 0x7fff
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
68 || w.ieee.exponent == 0x7fff
111
kono
parents: 68
diff changeset
69 || x == 0
kono
parents: 68
diff changeset
70 || y == 0)
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
71 return x * y + z;
111
kono
parents: 68
diff changeset
72 /* If fma will certainly overflow, compute as x * y. */
kono
parents: 68
diff changeset
73 if (u.ieee.exponent + v.ieee.exponent
kono
parents: 68
diff changeset
74 > 0x7fff + IEEE854_FLOAT128_BIAS)
kono
parents: 68
diff changeset
75 return x * y;
kono
parents: 68
diff changeset
76 /* If x * y is less than 1/4 of FLT128_DENORM_MIN, neither the
kono
parents: 68
diff changeset
77 result nor whether there is underflow depends on its exact
kono
parents: 68
diff changeset
78 value, only on its sign. */
kono
parents: 68
diff changeset
79 if (u.ieee.exponent + v.ieee.exponent
kono
parents: 68
diff changeset
80 < IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG - 2)
kono
parents: 68
diff changeset
81 {
kono
parents: 68
diff changeset
82 int neg = u.ieee.negative ^ v.ieee.negative;
kono
parents: 68
diff changeset
83 __float128 tiny = neg ? -0x1p-16494Q : 0x1p-16494Q;
kono
parents: 68
diff changeset
84 if (w.ieee.exponent >= 3)
kono
parents: 68
diff changeset
85 return tiny + z;
kono
parents: 68
diff changeset
86 /* Scaling up, adding TINY and scaling down produces the
kono
parents: 68
diff changeset
87 correct result, because in round-to-nearest mode adding
kono
parents: 68
diff changeset
88 TINY has no effect and in other modes double rounding is
kono
parents: 68
diff changeset
89 harmless. But it may not produce required underflow
kono
parents: 68
diff changeset
90 exceptions. */
kono
parents: 68
diff changeset
91 v.value = z * 0x1p114Q + tiny;
kono
parents: 68
diff changeset
92 if (TININESS_AFTER_ROUNDING
kono
parents: 68
diff changeset
93 ? v.ieee.exponent < 115
kono
parents: 68
diff changeset
94 : (w.ieee.exponent == 0
kono
parents: 68
diff changeset
95 || (w.ieee.exponent == 1
kono
parents: 68
diff changeset
96 && w.ieee.negative != neg
kono
parents: 68
diff changeset
97 && w.ieee.mant_low == 0
kono
parents: 68
diff changeset
98 && w.ieee.mant_high == 0)))
kono
parents: 68
diff changeset
99 {
kono
parents: 68
diff changeset
100 __float128 force_underflow = x * y;
kono
parents: 68
diff changeset
101 math_force_eval (force_underflow);
kono
parents: 68
diff changeset
102 }
kono
parents: 68
diff changeset
103 return v.value * 0x1p-114Q;
kono
parents: 68
diff changeset
104 }
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
105 if (u.ieee.exponent + v.ieee.exponent
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
106 >= 0x7fff + IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
107 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
108 /* 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
109 at the end. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
110 if (u.ieee.exponent > v.ieee.exponent)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
111 u.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
112 else
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
113 v.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
114 /* 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
115 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
116 if (w.ieee.exponent > FLT128_MANT_DIG)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
117 w.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
118 adjust = 1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
119 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
120 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
121 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
122 /* Similarly.
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
123 If z exponent is very large and x and y exponents are
111
kono
parents: 68
diff changeset
124 very small, adjust them up to avoid spurious underflows,
kono
parents: 68
diff changeset
125 rather than down. */
kono
parents: 68
diff changeset
126 if (u.ieee.exponent + v.ieee.exponent
kono
parents: 68
diff changeset
127 <= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG)
kono
parents: 68
diff changeset
128 {
kono
parents: 68
diff changeset
129 if (u.ieee.exponent > v.ieee.exponent)
kono
parents: 68
diff changeset
130 u.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
kono
parents: 68
diff changeset
131 else
kono
parents: 68
diff changeset
132 v.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
kono
parents: 68
diff changeset
133 }
kono
parents: 68
diff changeset
134 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
135 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
136 if (u.ieee.exponent > FLT128_MANT_DIG)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
137 u.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
138 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
139 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
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 w.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
142 adjust = 1;
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 (u.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 u.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
147 if (v.ieee.exponent)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
148 v.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 v.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 (v.ieee.exponent >= 0x7fff - FLT128_MANT_DIG)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
153 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
154 v.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
155 if (u.ieee.exponent)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
156 u.ieee.exponent += FLT128_MANT_DIG;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
157 else
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
158 u.value *= 0x1p113Q;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
159 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
160 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
161 <= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG) */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
162 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
163 if (u.ieee.exponent > v.ieee.exponent)
111
kono
parents: 68
diff changeset
164 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
165 else
111
kono
parents: 68
diff changeset
166 v.ieee.exponent += 2 * FLT128_MANT_DIG + 2;
kono
parents: 68
diff changeset
167 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
168 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
169 if (w.ieee.exponent)
111
kono
parents: 68
diff changeset
170 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
171 else
111
kono
parents: 68
diff changeset
172 w.value *= 0x1p228Q;
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
173 adjust = -1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
174 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
175 /* Otherwise x * y should just affect inexact
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
176 and nothing else. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
177 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
178 x = u.value;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
179 y = v.value;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
180 z = w.value;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
181 }
111
kono
parents: 68
diff changeset
182
kono
parents: 68
diff changeset
183 /* Ensure correct sign of exact 0 + 0. */
kono
parents: 68
diff changeset
184 if (__builtin_expect ((x == 0 || y == 0) && z == 0, 0))
kono
parents: 68
diff changeset
185 {
kono
parents: 68
diff changeset
186 x = math_opt_barrier (x);
kono
parents: 68
diff changeset
187 return x * y + z;
kono
parents: 68
diff changeset
188 }
kono
parents: 68
diff changeset
189
kono
parents: 68
diff changeset
190 #ifdef USE_FENV_H
kono
parents: 68
diff changeset
191 fenv_t env;
kono
parents: 68
diff changeset
192 feholdexcept (&env);
kono
parents: 68
diff changeset
193 fesetround (FE_TONEAREST);
kono
parents: 68
diff changeset
194 #endif
kono
parents: 68
diff changeset
195
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
196 /* 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
197 #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
198 __float128 x1 = x * C;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
199 __float128 y1 = y * C;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
200 __float128 m1 = x * y;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
201 x1 = (x - x1) + x1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
202 y1 = (y - y1) + y1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
203 __float128 x2 = x - x1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
204 __float128 y2 = y - y1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
205 __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
206
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
207 /* 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
208 __float128 a1 = z + m1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
209 __float128 t1 = a1 - z;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
210 __float128 t2 = a1 - t1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
211 t1 = m1 - t1;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
212 t2 = z - t2;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
213 __float128 a2 = t1 + t2;
111
kono
parents: 68
diff changeset
214 /* Ensure the arithmetic is not scheduled after feclearexcept call. */
kono
parents: 68
diff changeset
215 math_force_eval (m2);
kono
parents: 68
diff changeset
216 math_force_eval (a2);
kono
parents: 68
diff changeset
217 #ifdef USE_FENV_H
kono
parents: 68
diff changeset
218 feclearexcept (FE_INEXACT);
kono
parents: 68
diff changeset
219 #endif
kono
parents: 68
diff changeset
220
kono
parents: 68
diff changeset
221 /* If the result is an exact zero, ensure it has the correct sign. */
kono
parents: 68
diff changeset
222 if (a1 == 0 && m2 == 0)
kono
parents: 68
diff changeset
223 {
kono
parents: 68
diff changeset
224 #ifdef USE_FENV_H
kono
parents: 68
diff changeset
225 feupdateenv (&env);
kono
parents: 68
diff changeset
226 #endif
kono
parents: 68
diff changeset
227 /* Ensure that round-to-nearest value of z + m1 is not reused. */
kono
parents: 68
diff changeset
228 z = math_opt_barrier (z);
kono
parents: 68
diff changeset
229 return z + m1;
kono
parents: 68
diff changeset
230 }
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
231
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
232 #ifdef USE_FENV_H
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
233 fesetround (FE_TOWARDZERO);
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
234 #endif
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
235 /* 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
236 u.value = a2 + m2;
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 if (__builtin_expect (adjust == 0, 1))
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
239 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
240 #ifdef USE_FENV_H
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
241 if ((u.ieee.mant_low & 1) == 0 && u.ieee.exponent != 0x7fff)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
242 u.ieee.mant_low |= fetestexcept (FE_INEXACT) != 0;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
243 feupdateenv (&env);
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
244 #endif
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
245 /* Result is a1 + u.value. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
246 return a1 + u.value;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
247 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
248 else if (__builtin_expect (adjust > 0, 1))
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
249 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
250 #ifdef USE_FENV_H
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
251 if ((u.ieee.mant_low & 1) == 0 && u.ieee.exponent != 0x7fff)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
252 u.ieee.mant_low |= fetestexcept (FE_INEXACT) != 0;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
253 feupdateenv (&env);
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
254 #endif
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
255 /* Result is a1 + u.value, scaled up. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
256 return (a1 + u.value) * 0x1p113Q;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
257 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
258 else
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
259 {
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
260 #ifdef USE_FENV_H
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
261 if ((u.ieee.mant_low & 1) == 0)
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
262 u.ieee.mant_low |= fetestexcept (FE_INEXACT) != 0;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
263 #endif
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
264 v.value = a1 + u.value;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
265 /* Ensure the addition is not scheduled after fetestexcept call. */
111
kono
parents: 68
diff changeset
266 asm volatile ("" : : "m" (v.value));
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
267 #ifdef USE_FENV_H
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
268 int j = fetestexcept (FE_INEXACT) != 0;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
269 feupdateenv (&env);
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
270 #else
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
271 int j = 0;
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
272 #endif
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
273 /* 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
274 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
275 asm volatile ("" : "=m" (u) : "m" (u));
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
276 /* 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
277 scaling down. */
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
278 if (j == 0)
111
kono
parents: 68
diff changeset
279 return v.value * 0x1p-228Q;
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
280 /* 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
281 rounding will occur. */
111
kono
parents: 68
diff changeset
282 if (v.ieee.exponent > 228)
kono
parents: 68
diff changeset
283 return (a1 + u.value) * 0x1p-228Q;
kono
parents: 68
diff changeset
284 /* If v.value * 0x1p-228Q with round to zero is a subnormal above
kono
parents: 68
diff changeset
285 or equal to FLT128_MIN / 2, then v.value * 0x1p-228Q shifts mantissa
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
286 down just by 1 bit, which means v.ieee.mant_low |= j would
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
287 change the round bit, not sticky or guard bit.
111
kono
parents: 68
diff changeset
288 v.value * 0x1p-228Q never normalizes by shifting up,
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
289 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
290 for proper rounding. */
111
kono
parents: 68
diff changeset
291 if (v.ieee.exponent == 228)
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
292 {
111
kono
parents: 68
diff changeset
293 /* If the exponent would be in the normal range when
kono
parents: 68
diff changeset
294 rounding to normal precision with unbounded exponent
kono
parents: 68
diff changeset
295 range, the exact result is known and spurious underflows
kono
parents: 68
diff changeset
296 must be avoided on systems detecting tininess after
kono
parents: 68
diff changeset
297 rounding. */
kono
parents: 68
diff changeset
298 if (TININESS_AFTER_ROUNDING)
kono
parents: 68
diff changeset
299 {
kono
parents: 68
diff changeset
300 w.value = a1 + u.value;
kono
parents: 68
diff changeset
301 if (w.ieee.exponent == 229)
kono
parents: 68
diff changeset
302 return w.value * 0x1p-228Q;
kono
parents: 68
diff changeset
303 }
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
304 /* v.ieee.mant_low & 2 is LSB bit of the result before rounding,
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
305 v.ieee.mant_low & 1 is the round bit and j is our sticky
111
kono
parents: 68
diff changeset
306 bit. */
kono
parents: 68
diff changeset
307 w.value = 0.0Q;
kono
parents: 68
diff changeset
308 w.ieee.mant_low = ((v.ieee.mant_low & 3) << 1) | j;
kono
parents: 68
diff changeset
309 w.ieee.negative = v.ieee.negative;
kono
parents: 68
diff changeset
310 v.ieee.mant_low &= ~3U;
kono
parents: 68
diff changeset
311 v.value *= 0x1p-228Q;
kono
parents: 68
diff changeset
312 w.value *= 0x1p-2Q;
kono
parents: 68
diff changeset
313 return v.value + w.value;
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
314 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
315 v.ieee.mant_low |= j;
111
kono
parents: 68
diff changeset
316 return v.value * 0x1p-228Q;
68
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
317 }
561a7518be6b update gcc-4.6
Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
parents:
diff changeset
318 }