/* 128-bit long double support routines for Darwin.
|
/* 128-bit long double support routines for Darwin.
|
Copyright (C) 1993, 2003, 2004, 2005, 2006, 2007
|
Copyright (C) 1993, 2003, 2004, 2005, 2006, 2007
|
Free Software Foundation, Inc.
|
Free Software Foundation, Inc.
|
|
|
This file is part of GCC.
|
This file is part of GCC.
|
|
|
GCC is free software; you can redistribute it and/or modify it under
|
GCC is free software; you can redistribute it and/or modify it under
|
the terms of the GNU General Public License as published by the Free
|
the terms of the GNU General Public License as published by the Free
|
Software Foundation; either version 2, or (at your option) any later
|
Software Foundation; either version 2, or (at your option) any later
|
version.
|
version.
|
|
|
In addition to the permissions in the GNU General Public License, the
|
In addition to the permissions in the GNU General Public License, the
|
Free Software Foundation gives you unlimited permission to link the
|
Free Software Foundation gives you unlimited permission to link the
|
compiled version of this file into combinations with other programs,
|
compiled version of this file into combinations with other programs,
|
and to distribute those combinations without any restriction coming
|
and to distribute those combinations without any restriction coming
|
from the use of this file. (The General Public License restrictions
|
from the use of this file. (The General Public License restrictions
|
do apply in other respects; for example, they cover modification of
|
do apply in other respects; for example, they cover modification of
|
the file, and distribution when not linked into a combine
|
the file, and distribution when not linked into a combine
|
executable.)
|
executable.)
|
|
|
GCC is distributed in the hope that it will be useful, but WITHOUT ANY
|
GCC is distributed in the hope that it will be useful, but WITHOUT ANY
|
WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
for more details.
|
for more details.
|
|
|
You should have received a copy of the GNU General Public License
|
You should have received a copy of the GNU General Public License
|
along with GCC; see the file COPYING. If not, write to the Free
|
along with GCC; see the file COPYING. If not, write to the Free
|
Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
|
Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
|
02110-1301, USA. */
|
02110-1301, USA. */
|
|
|
/* Implementations of floating-point long double basic arithmetic
|
/* Implementations of floating-point long double basic arithmetic
|
functions called by the IBM C compiler when generating code for
|
functions called by the IBM C compiler when generating code for
|
PowerPC platforms. In particular, the following functions are
|
PowerPC platforms. In particular, the following functions are
|
implemented: __gcc_qadd, __gcc_qsub, __gcc_qmul, and __gcc_qdiv.
|
implemented: __gcc_qadd, __gcc_qsub, __gcc_qmul, and __gcc_qdiv.
|
Double-double algorithms are based on the paper "Doubled-Precision
|
Double-double algorithms are based on the paper "Doubled-Precision
|
IEEE Standard 754 Floating-Point Arithmetic" by W. Kahan, February 26,
|
IEEE Standard 754 Floating-Point Arithmetic" by W. Kahan, February 26,
|
1987. An alternative published reference is "Software for
|
1987. An alternative published reference is "Software for
|
Doubled-Precision Floating-Point Computations", by Seppo Linnainmaa,
|
Doubled-Precision Floating-Point Computations", by Seppo Linnainmaa,
|
ACM TOMS vol 7 no 3, September 1981, pages 272-283. */
|
ACM TOMS vol 7 no 3, September 1981, pages 272-283. */
|
|
|
/* Each long double is made up of two IEEE doubles. The value of the
|
/* Each long double is made up of two IEEE doubles. The value of the
|
long double is the sum of the values of the two parts. The most
|
long double is the sum of the values of the two parts. The most
|
significant part is required to be the value of the long double
|
significant part is required to be the value of the long double
|
rounded to the nearest double, as specified by IEEE. For Inf
|
rounded to the nearest double, as specified by IEEE. For Inf
|
values, the least significant part is required to be one of +0.0 or
|
values, the least significant part is required to be one of +0.0 or
|
-0.0. No other requirements are made; so, for example, 1.0 may be
|
-0.0. No other requirements are made; so, for example, 1.0 may be
|
represented as (1.0, +0.0) or (1.0, -0.0), and the low part of a
|
represented as (1.0, +0.0) or (1.0, -0.0), and the low part of a
|
NaN is don't-care.
|
NaN is don't-care.
|
|
|
This code currently assumes big-endian. */
|
This code currently assumes big-endian. */
|
|
|
#if ((!defined (__NO_FPRS__) || defined (_SOFT_FLOAT)) \
|
#if ((!defined (__NO_FPRS__) || defined (_SOFT_FLOAT)) \
|
&& !defined (__LITTLE_ENDIAN__) \
|
&& !defined (__LITTLE_ENDIAN__) \
|
&& (defined (__MACH__) || defined (__powerpc__) || defined (_AIX)))
|
&& (defined (__MACH__) || defined (__powerpc__) || defined (_AIX)))
|
|
|
#define fabs(x) __builtin_fabs(x)
|
#define fabs(x) __builtin_fabs(x)
|
#define isless(x, y) __builtin_isless (x, y)
|
#define isless(x, y) __builtin_isless (x, y)
|
#define inf() __builtin_inf()
|
#define inf() __builtin_inf()
|
|
|
#define unlikely(x) __builtin_expect ((x), 0)
|
#define unlikely(x) __builtin_expect ((x), 0)
|
|
|
#define nonfinite(a) unlikely (! isless (fabs (a), inf ()))
|
#define nonfinite(a) unlikely (! isless (fabs (a), inf ()))
|
|
|
/* Define ALIASNAME as a strong alias for NAME. */
|
/* Define ALIASNAME as a strong alias for NAME. */
|
# define strong_alias(name, aliasname) _strong_alias(name, aliasname)
|
# define strong_alias(name, aliasname) _strong_alias(name, aliasname)
|
# define _strong_alias(name, aliasname) \
|
# define _strong_alias(name, aliasname) \
|
extern __typeof (name) aliasname __attribute__ ((alias (#name)));
|
extern __typeof (name) aliasname __attribute__ ((alias (#name)));
|
|
|
/* All these routines actually take two long doubles as parameters,
|
/* All these routines actually take two long doubles as parameters,
|
but GCC currently generates poor code when a union is used to turn
|
but GCC currently generates poor code when a union is used to turn
|
a long double into a pair of doubles. */
|
a long double into a pair of doubles. */
|
|
|
long double __gcc_qadd (double, double, double, double);
|
long double __gcc_qadd (double, double, double, double);
|
long double __gcc_qsub (double, double, double, double);
|
long double __gcc_qsub (double, double, double, double);
|
long double __gcc_qmul (double, double, double, double);
|
long double __gcc_qmul (double, double, double, double);
|
long double __gcc_qdiv (double, double, double, double);
|
long double __gcc_qdiv (double, double, double, double);
|
|
|
#if defined __ELF__ && defined SHARED \
|
#if defined __ELF__ && defined SHARED \
|
&& (defined __powerpc64__ || !(defined __linux__ || defined __gnu_hurd__))
|
&& (defined __powerpc64__ || !(defined __linux__ || defined __gnu_hurd__))
|
/* Provide definitions of the old symbol names to satisfy apps and
|
/* Provide definitions of the old symbol names to satisfy apps and
|
shared libs built against an older libgcc. To access the _xlq
|
shared libs built against an older libgcc. To access the _xlq
|
symbols an explicit version reference is needed, so these won't
|
symbols an explicit version reference is needed, so these won't
|
satisfy an unadorned reference like _xlqadd. If dot symbols are
|
satisfy an unadorned reference like _xlqadd. If dot symbols are
|
not needed, the assembler will remove the aliases from the symbol
|
not needed, the assembler will remove the aliases from the symbol
|
table. */
|
table. */
|
__asm__ (".symver __gcc_qadd,_xlqadd@GCC_3.4\n\t"
|
__asm__ (".symver __gcc_qadd,_xlqadd@GCC_3.4\n\t"
|
".symver __gcc_qsub,_xlqsub@GCC_3.4\n\t"
|
".symver __gcc_qsub,_xlqsub@GCC_3.4\n\t"
|
".symver __gcc_qmul,_xlqmul@GCC_3.4\n\t"
|
".symver __gcc_qmul,_xlqmul@GCC_3.4\n\t"
|
".symver __gcc_qdiv,_xlqdiv@GCC_3.4\n\t"
|
".symver __gcc_qdiv,_xlqdiv@GCC_3.4\n\t"
|
".symver .__gcc_qadd,._xlqadd@GCC_3.4\n\t"
|
".symver .__gcc_qadd,._xlqadd@GCC_3.4\n\t"
|
".symver .__gcc_qsub,._xlqsub@GCC_3.4\n\t"
|
".symver .__gcc_qsub,._xlqsub@GCC_3.4\n\t"
|
".symver .__gcc_qmul,._xlqmul@GCC_3.4\n\t"
|
".symver .__gcc_qmul,._xlqmul@GCC_3.4\n\t"
|
".symver .__gcc_qdiv,._xlqdiv@GCC_3.4");
|
".symver .__gcc_qdiv,._xlqdiv@GCC_3.4");
|
#endif
|
#endif
|
|
|
typedef union
|
typedef union
|
{
|
{
|
long double ldval;
|
long double ldval;
|
double dval[2];
|
double dval[2];
|
} longDblUnion;
|
} longDblUnion;
|
|
|
/* Add two 'long double' values and return the result. */
|
/* Add two 'long double' values and return the result. */
|
long double
|
long double
|
__gcc_qadd (double a, double aa, double c, double cc)
|
__gcc_qadd (double a, double aa, double c, double cc)
|
{
|
{
|
longDblUnion x;
|
longDblUnion x;
|
double z, q, zz, xh;
|
double z, q, zz, xh;
|
|
|
z = a + c;
|
z = a + c;
|
|
|
if (nonfinite (z))
|
if (nonfinite (z))
|
{
|
{
|
z = cc + aa + c + a;
|
z = cc + aa + c + a;
|
if (nonfinite (z))
|
if (nonfinite (z))
|
return z;
|
return z;
|
x.dval[0] = z; /* Will always be DBL_MAX. */
|
x.dval[0] = z; /* Will always be DBL_MAX. */
|
zz = aa + cc;
|
zz = aa + cc;
|
if (fabs(a) > fabs(c))
|
if (fabs(a) > fabs(c))
|
x.dval[1] = a - z + c + zz;
|
x.dval[1] = a - z + c + zz;
|
else
|
else
|
x.dval[1] = c - z + a + zz;
|
x.dval[1] = c - z + a + zz;
|
}
|
}
|
else
|
else
|
{
|
{
|
q = a - z;
|
q = a - z;
|
zz = q + c + (a - (q + z)) + aa + cc;
|
zz = q + c + (a - (q + z)) + aa + cc;
|
|
|
/* Keep -0 result. */
|
/* Keep -0 result. */
|
if (zz == 0.0)
|
if (zz == 0.0)
|
return z;
|
return z;
|
|
|
xh = z + zz;
|
xh = z + zz;
|
if (nonfinite (xh))
|
if (nonfinite (xh))
|
return xh;
|
return xh;
|
|
|
x.dval[0] = xh;
|
x.dval[0] = xh;
|
x.dval[1] = z - xh + zz;
|
x.dval[1] = z - xh + zz;
|
}
|
}
|
return x.ldval;
|
return x.ldval;
|
}
|
}
|
|
|
long double
|
long double
|
__gcc_qsub (double a, double b, double c, double d)
|
__gcc_qsub (double a, double b, double c, double d)
|
{
|
{
|
return __gcc_qadd (a, b, -c, -d);
|
return __gcc_qadd (a, b, -c, -d);
|
}
|
}
|
|
|
#ifdef _SOFT_FLOAT
|
#ifdef _SOFT_FLOAT
|
static double fmsub (double, double, double);
|
static double fmsub (double, double, double);
|
#endif
|
#endif
|
|
|
long double
|
long double
|
__gcc_qmul (double a, double b, double c, double d)
|
__gcc_qmul (double a, double b, double c, double d)
|
{
|
{
|
longDblUnion z;
|
longDblUnion z;
|
double t, tau, u, v, w;
|
double t, tau, u, v, w;
|
|
|
t = a * c; /* Highest order double term. */
|
t = a * c; /* Highest order double term. */
|
|
|
if (unlikely (t == 0) /* Preserve -0. */
|
if (unlikely (t == 0) /* Preserve -0. */
|
|| nonfinite (t))
|
|| nonfinite (t))
|
return t;
|
return t;
|
|
|
/* Sum terms of two highest orders. */
|
/* Sum terms of two highest orders. */
|
|
|
/* Use fused multiply-add to get low part of a * c. */
|
/* Use fused multiply-add to get low part of a * c. */
|
#ifndef _SOFT_FLOAT
|
#ifndef _SOFT_FLOAT
|
asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
|
asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
|
#else
|
#else
|
tau = fmsub (a, c, t);
|
tau = fmsub (a, c, t);
|
#endif
|
#endif
|
v = a*d;
|
v = a*d;
|
w = b*c;
|
w = b*c;
|
tau += v + w; /* Add in other second-order terms. */
|
tau += v + w; /* Add in other second-order terms. */
|
u = t + tau;
|
u = t + tau;
|
|
|
/* Construct long double result. */
|
/* Construct long double result. */
|
if (nonfinite (u))
|
if (nonfinite (u))
|
return u;
|
return u;
|
z.dval[0] = u;
|
z.dval[0] = u;
|
z.dval[1] = (t - u) + tau;
|
z.dval[1] = (t - u) + tau;
|
return z.ldval;
|
return z.ldval;
|
}
|
}
|
|
|
long double
|
long double
|
__gcc_qdiv (double a, double b, double c, double d)
|
__gcc_qdiv (double a, double b, double c, double d)
|
{
|
{
|
longDblUnion z;
|
longDblUnion z;
|
double s, sigma, t, tau, u, v, w;
|
double s, sigma, t, tau, u, v, w;
|
|
|
t = a / c; /* highest order double term */
|
t = a / c; /* highest order double term */
|
|
|
if (unlikely (t == 0) /* Preserve -0. */
|
if (unlikely (t == 0) /* Preserve -0. */
|
|| nonfinite (t))
|
|| nonfinite (t))
|
return t;
|
return t;
|
|
|
/* Finite nonzero result requires corrections to the highest order term. */
|
/* Finite nonzero result requires corrections to the highest order term. */
|
|
|
s = c * t; /* (s,sigma) = c*t exactly. */
|
s = c * t; /* (s,sigma) = c*t exactly. */
|
w = -(-b + d * t); /* Written to get fnmsub for speed, but not
|
w = -(-b + d * t); /* Written to get fnmsub for speed, but not
|
numerically necessary. */
|
numerically necessary. */
|
|
|
/* Use fused multiply-add to get low part of c * t. */
|
/* Use fused multiply-add to get low part of c * t. */
|
#ifndef _SOFT_FLOAT
|
#ifndef _SOFT_FLOAT
|
asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s));
|
asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s));
|
#else
|
#else
|
sigma = fmsub (c, t, s);
|
sigma = fmsub (c, t, s);
|
#endif
|
#endif
|
v = a - s;
|
v = a - s;
|
|
|
tau = ((v-sigma)+w)/c; /* Correction to t. */
|
tau = ((v-sigma)+w)/c; /* Correction to t. */
|
u = t + tau;
|
u = t + tau;
|
|
|
/* Construct long double result. */
|
/* Construct long double result. */
|
if (nonfinite (u))
|
if (nonfinite (u))
|
return u;
|
return u;
|
z.dval[0] = u;
|
z.dval[0] = u;
|
z.dval[1] = (t - u) + tau;
|
z.dval[1] = (t - u) + tau;
|
return z.ldval;
|
return z.ldval;
|
}
|
}
|
|
|
#if defined (_SOFT_FLOAT) && defined (__LONG_DOUBLE_128__)
|
#if defined (_SOFT_FLOAT) && defined (__LONG_DOUBLE_128__)
|
|
|
long double __gcc_qneg (double, double);
|
long double __gcc_qneg (double, double);
|
int __gcc_qeq (double, double, double, double);
|
int __gcc_qeq (double, double, double, double);
|
int __gcc_qne (double, double, double, double);
|
int __gcc_qne (double, double, double, double);
|
int __gcc_qge (double, double, double, double);
|
int __gcc_qge (double, double, double, double);
|
int __gcc_qle (double, double, double, double);
|
int __gcc_qle (double, double, double, double);
|
int __gcc_qunord (double, double, double, double);
|
int __gcc_qunord (double, double, double, double);
|
long double __gcc_stoq (float);
|
long double __gcc_stoq (float);
|
long double __gcc_dtoq (double);
|
long double __gcc_dtoq (double);
|
float __gcc_qtos (double, double);
|
float __gcc_qtos (double, double);
|
double __gcc_qtod (double, double);
|
double __gcc_qtod (double, double);
|
int __gcc_qtoi (double, double);
|
int __gcc_qtoi (double, double);
|
unsigned int __gcc_qtou (double, double);
|
unsigned int __gcc_qtou (double, double);
|
long double __gcc_itoq (int);
|
long double __gcc_itoq (int);
|
long double __gcc_utoq (unsigned int);
|
long double __gcc_utoq (unsigned int);
|
|
|
extern int __eqdf2 (double, double);
|
extern int __eqdf2 (double, double);
|
extern int __ledf2 (double, double);
|
extern int __ledf2 (double, double);
|
extern int __gedf2 (double, double);
|
extern int __gedf2 (double, double);
|
extern int __unorddf2 (double, double);
|
extern int __unorddf2 (double, double);
|
|
|
/* Negate 'long double' value and return the result. */
|
/* Negate 'long double' value and return the result. */
|
long double
|
long double
|
__gcc_qneg (double a, double aa)
|
__gcc_qneg (double a, double aa)
|
{
|
{
|
longDblUnion x;
|
longDblUnion x;
|
|
|
x.dval[0] = -a;
|
x.dval[0] = -a;
|
x.dval[1] = -aa;
|
x.dval[1] = -aa;
|
return x.ldval;
|
return x.ldval;
|
}
|
}
|
|
|
/* Compare two 'long double' values for equality. */
|
/* Compare two 'long double' values for equality. */
|
int
|
int
|
__gcc_qeq (double a, double aa, double c, double cc)
|
__gcc_qeq (double a, double aa, double c, double cc)
|
{
|
{
|
if (__eqdf2 (a, c) == 0)
|
if (__eqdf2 (a, c) == 0)
|
return __eqdf2 (aa, cc);
|
return __eqdf2 (aa, cc);
|
return 1;
|
return 1;
|
}
|
}
|
|
|
strong_alias (__gcc_qeq, __gcc_qne);
|
strong_alias (__gcc_qeq, __gcc_qne);
|
|
|
/* Compare two 'long double' values for less than or equal. */
|
/* Compare two 'long double' values for less than or equal. */
|
int
|
int
|
__gcc_qle (double a, double aa, double c, double cc)
|
__gcc_qle (double a, double aa, double c, double cc)
|
{
|
{
|
if (__eqdf2 (a, c) == 0)
|
if (__eqdf2 (a, c) == 0)
|
return __ledf2 (aa, cc);
|
return __ledf2 (aa, cc);
|
return __ledf2 (a, c);
|
return __ledf2 (a, c);
|
}
|
}
|
|
|
strong_alias (__gcc_qle, __gcc_qlt);
|
strong_alias (__gcc_qle, __gcc_qlt);
|
|
|
/* Compare two 'long double' values for greater than or equal. */
|
/* Compare two 'long double' values for greater than or equal. */
|
int
|
int
|
__gcc_qge (double a, double aa, double c, double cc)
|
__gcc_qge (double a, double aa, double c, double cc)
|
{
|
{
|
if (__eqdf2 (a, c) == 0)
|
if (__eqdf2 (a, c) == 0)
|
return __gedf2 (aa, cc);
|
return __gedf2 (aa, cc);
|
return __gedf2 (a, c);
|
return __gedf2 (a, c);
|
}
|
}
|
|
|
strong_alias (__gcc_qge, __gcc_qgt);
|
strong_alias (__gcc_qge, __gcc_qgt);
|
|
|
/* Compare two 'long double' values for unordered. */
|
/* Compare two 'long double' values for unordered. */
|
int
|
int
|
__gcc_qunord (double a, double aa, double c, double cc)
|
__gcc_qunord (double a, double aa, double c, double cc)
|
{
|
{
|
if (__eqdf2 (a, c) == 0)
|
if (__eqdf2 (a, c) == 0)
|
return __unorddf2 (aa, cc);
|
return __unorddf2 (aa, cc);
|
return __unorddf2 (a, c);
|
return __unorddf2 (a, c);
|
}
|
}
|
|
|
/* Convert single to long double. */
|
/* Convert single to long double. */
|
long double
|
long double
|
__gcc_stoq (float a)
|
__gcc_stoq (float a)
|
{
|
{
|
longDblUnion x;
|
longDblUnion x;
|
|
|
x.dval[0] = (double) a;
|
x.dval[0] = (double) a;
|
x.dval[1] = 0.0;
|
x.dval[1] = 0.0;
|
|
|
return x.ldval;
|
return x.ldval;
|
}
|
}
|
|
|
/* Convert double to long double. */
|
/* Convert double to long double. */
|
long double
|
long double
|
__gcc_dtoq (double a)
|
__gcc_dtoq (double a)
|
{
|
{
|
longDblUnion x;
|
longDblUnion x;
|
|
|
x.dval[0] = a;
|
x.dval[0] = a;
|
x.dval[1] = 0.0;
|
x.dval[1] = 0.0;
|
|
|
return x.ldval;
|
return x.ldval;
|
}
|
}
|
|
|
/* Convert long double to single. */
|
/* Convert long double to single. */
|
float
|
float
|
__gcc_qtos (double a, double aa __attribute__ ((__unused__)))
|
__gcc_qtos (double a, double aa __attribute__ ((__unused__)))
|
{
|
{
|
return (float) a;
|
return (float) a;
|
}
|
}
|
|
|
/* Convert long double to double. */
|
/* Convert long double to double. */
|
double
|
double
|
__gcc_qtod (double a, double aa __attribute__ ((__unused__)))
|
__gcc_qtod (double a, double aa __attribute__ ((__unused__)))
|
{
|
{
|
return a;
|
return a;
|
}
|
}
|
|
|
/* Convert long double to int. */
|
/* Convert long double to int. */
|
int
|
int
|
__gcc_qtoi (double a, double aa)
|
__gcc_qtoi (double a, double aa)
|
{
|
{
|
double z = a + aa;
|
double z = a + aa;
|
return (int) z;
|
return (int) z;
|
}
|
}
|
|
|
/* Convert long double to unsigned int. */
|
/* Convert long double to unsigned int. */
|
unsigned int
|
unsigned int
|
__gcc_qtou (double a, double aa)
|
__gcc_qtou (double a, double aa)
|
{
|
{
|
double z = a + aa;
|
double z = a + aa;
|
return (unsigned int) z;
|
return (unsigned int) z;
|
}
|
}
|
|
|
/* Convert int to long double. */
|
/* Convert int to long double. */
|
long double
|
long double
|
__gcc_itoq (int a)
|
__gcc_itoq (int a)
|
{
|
{
|
return __gcc_dtoq ((double) a);
|
return __gcc_dtoq ((double) a);
|
}
|
}
|
|
|
/* Convert unsigned int to long double. */
|
/* Convert unsigned int to long double. */
|
long double
|
long double
|
__gcc_utoq (unsigned int a)
|
__gcc_utoq (unsigned int a)
|
{
|
{
|
return __gcc_dtoq ((double) a);
|
return __gcc_dtoq ((double) a);
|
}
|
}
|
|
|
#include "config/soft-fp/soft-fp.h"
|
#include "config/soft-fp/soft-fp.h"
|
#include "config/soft-fp/double.h"
|
#include "config/soft-fp/double.h"
|
#include "config/soft-fp/quad.h"
|
#include "config/soft-fp/quad.h"
|
|
|
/* Compute floating point multiply-subtract with higher (quad) precision. */
|
/* Compute floating point multiply-subtract with higher (quad) precision. */
|
static double
|
static double
|
fmsub (double a, double b, double c)
|
fmsub (double a, double b, double c)
|
{
|
{
|
FP_DECL_EX;
|
FP_DECL_EX;
|
FP_DECL_D(A);
|
FP_DECL_D(A);
|
FP_DECL_D(B);
|
FP_DECL_D(B);
|
FP_DECL_D(C);
|
FP_DECL_D(C);
|
FP_DECL_Q(X);
|
FP_DECL_Q(X);
|
FP_DECL_Q(Y);
|
FP_DECL_Q(Y);
|
FP_DECL_Q(Z);
|
FP_DECL_Q(Z);
|
FP_DECL_Q(U);
|
FP_DECL_Q(U);
|
FP_DECL_Q(V);
|
FP_DECL_Q(V);
|
FP_DECL_D(R);
|
FP_DECL_D(R);
|
double r;
|
double r;
|
long double u, v, x, y, z;
|
long double u, v, x, y, z;
|
|
|
FP_INIT_ROUNDMODE;
|
FP_INIT_ROUNDMODE;
|
FP_UNPACK_RAW_D (A, a);
|
FP_UNPACK_RAW_D (A, a);
|
FP_UNPACK_RAW_D (B, b);
|
FP_UNPACK_RAW_D (B, b);
|
FP_UNPACK_RAW_D (C, c);
|
FP_UNPACK_RAW_D (C, c);
|
|
|
/* Extend double to quad. */
|
/* Extend double to quad. */
|
#if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
|
#if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
|
FP_EXTEND(Q,D,4,2,X,A);
|
FP_EXTEND(Q,D,4,2,X,A);
|
FP_EXTEND(Q,D,4,2,Y,B);
|
FP_EXTEND(Q,D,4,2,Y,B);
|
FP_EXTEND(Q,D,4,2,Z,C);
|
FP_EXTEND(Q,D,4,2,Z,C);
|
#else
|
#else
|
FP_EXTEND(Q,D,2,1,X,A);
|
FP_EXTEND(Q,D,2,1,X,A);
|
FP_EXTEND(Q,D,2,1,Y,B);
|
FP_EXTEND(Q,D,2,1,Y,B);
|
FP_EXTEND(Q,D,2,1,Z,C);
|
FP_EXTEND(Q,D,2,1,Z,C);
|
#endif
|
#endif
|
FP_PACK_RAW_Q(x,X);
|
FP_PACK_RAW_Q(x,X);
|
FP_PACK_RAW_Q(y,Y);
|
FP_PACK_RAW_Q(y,Y);
|
FP_PACK_RAW_Q(z,Z);
|
FP_PACK_RAW_Q(z,Z);
|
FP_HANDLE_EXCEPTIONS;
|
FP_HANDLE_EXCEPTIONS;
|
|
|
/* Multiply. */
|
/* Multiply. */
|
FP_INIT_ROUNDMODE;
|
FP_INIT_ROUNDMODE;
|
FP_UNPACK_Q(X,x);
|
FP_UNPACK_Q(X,x);
|
FP_UNPACK_Q(Y,y);
|
FP_UNPACK_Q(Y,y);
|
FP_MUL_Q(U,X,Y);
|
FP_MUL_Q(U,X,Y);
|
FP_PACK_Q(u,U);
|
FP_PACK_Q(u,U);
|
FP_HANDLE_EXCEPTIONS;
|
FP_HANDLE_EXCEPTIONS;
|
|
|
/* Subtract. */
|
/* Subtract. */
|
FP_INIT_ROUNDMODE;
|
FP_INIT_ROUNDMODE;
|
FP_UNPACK_SEMIRAW_Q(U,u);
|
FP_UNPACK_SEMIRAW_Q(U,u);
|
FP_UNPACK_SEMIRAW_Q(Z,z);
|
FP_UNPACK_SEMIRAW_Q(Z,z);
|
FP_SUB_Q(V,U,Z);
|
FP_SUB_Q(V,U,Z);
|
FP_PACK_SEMIRAW_Q(v,V);
|
FP_PACK_SEMIRAW_Q(v,V);
|
FP_HANDLE_EXCEPTIONS;
|
FP_HANDLE_EXCEPTIONS;
|
|
|
/* Truncate quad to double. */
|
/* Truncate quad to double. */
|
FP_INIT_ROUNDMODE;
|
FP_INIT_ROUNDMODE;
|
FP_UNPACK_SEMIRAW_Q(V,v);
|
FP_UNPACK_SEMIRAW_Q(V,v);
|
#if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
|
#if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q
|
FP_TRUNC(D,Q,2,4,R,V);
|
FP_TRUNC(D,Q,2,4,R,V);
|
#else
|
#else
|
FP_TRUNC(D,Q,1,2,R,V);
|
FP_TRUNC(D,Q,1,2,R,V);
|
#endif
|
#endif
|
FP_PACK_SEMIRAW_D(r,R);
|
FP_PACK_SEMIRAW_D(r,R);
|
FP_HANDLE_EXCEPTIONS;
|
FP_HANDLE_EXCEPTIONS;
|
|
|
return r;
|
return r;
|
}
|
}
|
|
|
#endif
|
#endif
|
|
|
#endif
|
#endif
|
|
|
