CbC/CbC_gcc: libquadmath/math/fmaq.c annotate

author	Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
date	Sun, 21 Aug 2011 07:07:55 +0900
parents
children	04ced10e8804

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.
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	2 Copyright (C) 2010 Free Software Foundation, Inc.
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
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	17 License along with the GNU C Library; if not, write to the Free
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	18 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	19 02111-1307 USA. */
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	20
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	21 #include "quadmath-imp.h"
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	22 #include <math.h>
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	23 #include <float.h>
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	24 #ifdef HAVE_FENV_H
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	25 # include <fenv.h>
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	26 # if defined HAVE_FEHOLDEXCEPT && defined HAVE_FESETROUND \
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	27 && defined HAVE_FEUPDATEENV && defined HAVE_FETESTEXCEPT \
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	28 && defined FE_TOWARDZERO && defined FE_INEXACT
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	29 # define USE_FENV_H
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 #endif
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	32
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	33 /* 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	34 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	35 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	36
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	37 __float128
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	38 fmaq (__float128 x, __float128 y, __float128 z)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	39 {
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	40 ieee854_float128 u, v, w;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	41 int adjust = 0;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	42 u.value = x;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	43 v.value = y;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	44 w.value = z;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	45 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	46 >= 0x7fff + IEEE854_FLOAT128_BIAS
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	47 - FLT128_MANT_DIG, 0)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	48 \|\| __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	49 \|\| __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	50 \|\| __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	51 \|\| __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	52 <= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG, 0))
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	53 {
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	54 /* 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	55 z rather than NaN. */
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	56 if (w.ieee.exponent == 0x7fff
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	57 && u.ieee.exponent != 0x7fff
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	58 && v.ieee.exponent != 0x7fff)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	59 return (z + x) + y;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	60 /* If x or y or z is Inf/NaN, or if fma will certainly overflow,
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	61 or if x * y is less than half of FLT128_DENORM_MIN,
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	62 compute as x * y + z. */
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	63 if (u.ieee.exponent == 0x7fff
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	64 \|\| v.ieee.exponent == 0x7fff
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	65 \|\| w.ieee.exponent == 0x7fff
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	66 \|\| u.ieee.exponent + v.ieee.exponent
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	67 > 0x7fff + IEEE854_FLOAT128_BIAS
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	68 \|\| u.ieee.exponent + v.ieee.exponent
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	69 < IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG - 2)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	70 return x * y + z;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	71 if (u.ieee.exponent + v.ieee.exponent
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	72 >= 0x7fff + IEEE854_FLOAT128_BIAS - FLT128_MANT_DIG)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	73 {
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	74 /* 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	75 at the end. */
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	76 if (u.ieee.exponent > v.ieee.exponent)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	77 u.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	78 else
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	79 v.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	80 /* 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	81 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	82 if (w.ieee.exponent > FLT128_MANT_DIG)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	83 w.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	84 adjust = 1;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	85 }
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	86 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	87 {
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	88 /* Similarly.
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	89 If z exponent is very large and x and y exponents are
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	90 very small, 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	91 if (u.ieee.exponent > v.ieee.exponent)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	92 {
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	93 if (u.ieee.exponent > FLT128_MANT_DIG)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	94 u.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	95 }
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	96 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	97 v.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	98 w.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	99 adjust = 1;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	100 }
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	101 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	102 {
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 if (v.ieee.exponent)
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 else
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	107 v.value *= 0x1p113Q;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	108 }
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	109 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	110 {
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	111 v.ieee.exponent -= FLT128_MANT_DIG;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	112 if (u.ieee.exponent)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	113 u.ieee.exponent += FLT128_MANT_DIG;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	114 else
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	115 u.value *= 0x1p113Q;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	116 }
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	117 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	118 <= IEEE854_FLOAT128_BIAS + FLT128_MANT_DIG) */
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 if (u.ieee.exponent > v.ieee.exponent)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	121 u.ieee.exponent += 2 * FLT128_MANT_DIG;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	122 else
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	123 v.ieee.exponent += 2 * FLT128_MANT_DIG;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	124 if (w.ieee.exponent <= 4 * FLT128_MANT_DIG + 4)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	125 {
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	126 if (w.ieee.exponent)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	127 w.ieee.exponent += 2 * FLT128_MANT_DIG;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	128 else
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	129 w.value *= 0x1p226Q;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	130 adjust = -1;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	131 }
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	132 /* Otherwise x * y should just affect inexact
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	133 and nothing else. */
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	134 }
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	135 x = u.value;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	136 y = v.value;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	137 z = w.value;
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 /* 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	140 #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	141 __float128 x1 = x * C;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	142 __float128 y1 = y * C;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	143 __float128 m1 = x * y;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	144 x1 = (x - x1) + x1;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	145 y1 = (y - y1) + y1;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	146 __float128 x2 = x - x1;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	147 __float128 y2 = y - y1;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	148 __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	149
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	150 /* 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	151 __float128 a1 = z + m1;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	152 __float128 t1 = a1 - z;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	153 __float128 t2 = a1 - t1;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	154 t1 = m1 - t1;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	155 t2 = z - t2;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	156 __float128 a2 = t1 + t2;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	157
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	158 #ifdef USE_FENV_H
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	159 fenv_t env;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	160 feholdexcept (&env);
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	161 fesetround (FE_TOWARDZERO);
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	162 #endif
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	163 /* 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	164 u.value = a2 + m2;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	165
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	166 if (__builtin_expect (adjust == 0, 1))
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	167 {
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	168 #ifdef USE_FENV_H
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	169 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	170 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	171 feupdateenv (&env);
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	172 #endif
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	173 /* Result is a1 + u.value. */
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	174 return a1 + u.value;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	175 }
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	176 else if (__builtin_expect (adjust > 0, 1))
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 #ifdef USE_FENV_H
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	179 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	180 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	181 feupdateenv (&env);
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	182 #endif
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	183 /* Result is a1 + u.value, scaled up. */
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	184 return (a1 + u.value) * 0x1p113Q;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	185 }
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	186 else
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	187 {
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	188 #ifdef USE_FENV_H
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	189 if ((u.ieee.mant_low & 1) == 0)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	190 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	191 #endif
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	192 v.value = a1 + u.value;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	193 /* Ensure the addition is not scheduled after fetestexcept call. */
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	194 asm volatile ("" : : "m" (v));
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	195 #ifdef USE_FENV_H
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	196 int j = fetestexcept (FE_INEXACT) != 0;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	197 feupdateenv (&env);
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	198 #else
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	199 int j = 0;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	200 #endif
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	201 /* 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	202 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	203 asm volatile ("" : "=m" (u) : "m" (u));
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	204 /* 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	205 scaling down. */
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	206 if (j == 0)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	207 return v.value * 0x1p-226Q;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	208 /* 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	209 rounding will occur. */
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	210 if (v.ieee.exponent > 226)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	211 return (a1 + u.value) * 0x1p-226Q;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	212 /* If v.value * 0x1p-226Q with round to zero is a subnormal above
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	213 or equal to FLT128_MIN / 2, then v.value * 0x1p-226Q shifts mantissa
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	214 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	215 change the round bit, not sticky or guard bit.
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	216 v.value * 0x1p-226Q never normalizes by shifting up,
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	217 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	218 for proper rounding. */
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	219 if (v.ieee.exponent == 226)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	220 {
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	221 /* 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	222 v.ieee.mant_low & 1 is the round bit and j is our sticky
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	223 bit. In round-to-nearest 001 rounds down like 00,
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	224 011 rounds up, even though 01 rounds down (thus we need
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	225 to adjust), 101 rounds down like 10 and 111 rounds up
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	226 like 11. */
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	227 if ((v.ieee.mant_low & 3) == 1)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	228 {
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	229 v.value *= 0x1p-226Q;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	230 if (v.ieee.negative)
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	231 return v.value - 0x1p-16494Q /* __FLT128_DENORM_MIN__ */;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	232 else
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	233 return v.value + 0x1p-16494Q /* __FLT128_DENORM_MIN__ */;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	234 }
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	235 else
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	236 return v.value * 0x1p-226Q;
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 v.ieee.mant_low \|= j;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	239 return v.value * 0x1p-226Q;
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	240 }
561a7518be6b update gcc-4.6 Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp> parents: diff changeset	241 }

68

561a7518be6b update gcc-4.6