Subversion Repositories openrisc_me

/* 128-bit long double support routines for Darwin.

   Copyright (C) 1993, 2003, 2004, 2005, 2006, 2007, 2008, 2009

   Free Software Foundation, Inc.

This file is part of GCC.

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

Software Foundation; either version 3, or (at your option) any later

version.

GCC is distributed in the hope that it will be useful, but WITHOUT ANY

WARRANTY; without even the implied warranty of MERCHANTABILITY or

FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License

for more details.

Under Section 7 of GPL version 3, you are granted additional

permissions described in the GCC Runtime Library Exception, version

3.1, as published by the Free Software Foundation.

You should have received a copy of the GNU General Public License and

a copy of the GCC Runtime Library Exception along with this program;

see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see

<http://www.gnu.org/licenses/>.  */

/* Implementations of floating-point long double basic arithmetic

   functions called by the IBM C compiler when generating code for

   PowerPC platforms.  In particular, the following functions are

   implemented: __gcc_qadd, __gcc_qsub, __gcc_qmul, and __gcc_qdiv.

   Double-double algorithms are based on the paper "Doubled-Precision

   IEEE Standard 754 Floating-Point Arithmetic" by W. Kahan, February 26,

   1987.  An alternative published reference is "Software for

   Doubled-Precision Floating-Point Computations", by Seppo Linnainmaa,

   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

   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

   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

   -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

   NaN is don't-care.

   This code currently assumes big-endian.  */

#if (!defined (__LITTLE_ENDIAN__) \

     && (defined (__MACH__) || defined (__powerpc__) || defined (_AIX)))

#define fabs(x) __builtin_fabs(x)

#define isless(x, y) __builtin_isless (x, y)

#define inf() __builtin_inf()

#define unlikely(x) __builtin_expect ((x), 0)

#define nonfinite(a) unlikely (! isless (fabs (a), inf ()))

/* Define ALIASNAME as a strong alias for NAME.  */

# define strong_alias(name, aliasname) _strong_alias(name, aliasname)

# define _strong_alias(name, aliasname) \

  extern __typeof (name) aliasname __attribute__ ((alias (#name)));

/* All these routines actually take two long doubles as parameters,

   but GCC currently generates poor code when a union is used to turn

   a long double into a pair of doubles.  */

long double __gcc_qadd (double, double, double, double);

long double __gcc_qsub (double, double, double, double);

long double __gcc_qmul (double, double, double, double);

long double __gcc_qdiv (double, double, double, double);

#if defined __ELF__ && defined SHARED \

    && (defined __powerpc64__ || !(defined __linux__ || defined __gnu_hurd__))

/* Provide definitions of the old symbol names to satisfy apps and

   shared libs built against an older libgcc.  To access the _xlq

   symbols an explicit version reference is needed, so these won't

   satisfy an unadorned reference like _xlqadd.  If dot symbols are

   not needed, the assembler will remove the aliases from the symbol

   table.  */

__asm__ (".symver __gcc_qadd,_xlqadd@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_qdiv,_xlqdiv@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_qmul,._xlqmul@GCC_3.4\n\t"

         ".symver .__gcc_qdiv,._xlqdiv@GCC_3.4");

#endif

typedef union

{

  long double ldval;

  double dval[2];

} longDblUnion;

/* Add two 'long double' values and return the result.  */

long double

__gcc_qadd (double a, double aa, double c, double cc)

{

  longDblUnion x;

  double z, q, zz, xh;

  z = a + c;

  if (nonfinite (z))

    {

      z = cc + aa + c + a;

      if (nonfinite (z))

        return z;

      x.dval[0] = z;  /* Will always be DBL_MAX.  */

      zz = aa + cc;

      if (fabs(a) > fabs(c))

        x.dval[1] = a - z + c + zz;

      else

        x.dval[1] = c - z + a + zz;

    }

  else

    {

      q = a - z;

      zz = q + c + (a - (q + z)) + aa + cc;

      /* Keep -0 result.  */

      if (zz == 0.0)

        return z;

      xh = z + zz;

      if (nonfinite (xh))

        return xh;

      x.dval[0] = xh;

      x.dval[1] = z - xh + zz;

    }

  return x.ldval;

}

long double

__gcc_qsub (double a, double b, double c, double d)

{

  return __gcc_qadd (a, b, -c, -d);

}

#ifdef __NO_FPRS__

static double fmsub (double, double, double);

#endif

long double

__gcc_qmul (double a, double b, double c, double d)

{

  longDblUnion z;

  double t, tau, u, v, w;

  t = a * c;                    /* Highest order double term.  */

  if (unlikely (t == 0)          /* Preserve -0.  */

      || nonfinite (t))

    return t;

  /* Sum terms of two highest orders. */

  /* Use fused multiply-add to get low part of a * c.  */

#ifndef __NO_FPRS__

  asm ("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));

#else

  tau = fmsub (a, c, t);

#endif

  v = a*d;

  w = b*c;

  tau += v + w;     /* Add in other second-order terms.  */

  u = t + tau;

  /* Construct long double result.  */

  if (nonfinite (u))

    return u;

  z.dval[0] = u;

  z.dval[1] = (t - u) + tau;

  return z.ldval;

}

long double

__gcc_qdiv (double a, double b, double c, double d)

{

  longDblUnion z;

  double s, sigma, t, tau, u, v, w;

  t = a / c;                    /* highest order double term */

  if (unlikely (t == 0)          /* Preserve -0.  */

      || nonfinite (t))

    return t;

  /* Finite nonzero result requires corrections to the highest order term.  */

  s = c * t;                    /* (s,sigma) = c*t exactly.  */

  w = -(-b + d * t);    /* Written to get fnmsub for speed, but not

                           numerically necessary.  */

  /* Use fused multiply-add to get low part of c * t.    */

#ifndef __NO_FPRS__

  asm ("fmsub %0,%1,%2,%3" : "=f"(sigma) : "f"(c), "f"(t), "f"(s));

#else

  sigma = fmsub (c, t, s);

#endif

  v = a - s;

  tau = ((v-sigma)+w)/c;   /* Correction to t.  */

  u = t + tau;

  /* Construct long double result.  */

  if (nonfinite (u))

    return u;

  z.dval[0] = u;

  z.dval[1] = (t - u) + tau;

  return z.ldval;

}

#if defined (_SOFT_DOUBLE) && defined (__LONG_DOUBLE_128__)

long double __gcc_qneg (double, double);

int __gcc_qeq (double, double, double, double);

int __gcc_qne (double, double, double, double);

int __gcc_qge (double, double, double, double);

int __gcc_qle (double, double, double, double);

long double __gcc_stoq (float);

long double __gcc_dtoq (double);

float __gcc_qtos (double, double);

double __gcc_qtod (double, double);

int __gcc_qtoi (double, double);

unsigned int __gcc_qtou (double, double);

long double __gcc_itoq (int);

long double __gcc_utoq (unsigned int);

extern int __eqdf2 (double, double);

extern int __ledf2 (double, double);

extern int __gedf2 (double, double);

/* Negate 'long double' value and return the result.    */

long double

__gcc_qneg (double a, double aa)

{

  longDblUnion x;

  x.dval[0] = -a;

  x.dval[1] = -aa;

  return x.ldval;

}

/* Compare two 'long double' values for equality.  */

int

__gcc_qeq (double a, double aa, double c, double cc)

{

  if (__eqdf2 (a, c) == 0)

    return __eqdf2 (aa, cc);

  return 1;

}

strong_alias (__gcc_qeq, __gcc_qne);

/* Compare two 'long double' values for less than or equal.  */

int

__gcc_qle (double a, double aa, double c, double cc)

{

  if (__eqdf2 (a, c) == 0)

    return __ledf2 (aa, cc);

  return __ledf2 (a, c);

}

strong_alias (__gcc_qle, __gcc_qlt);

/* Compare two 'long double' values for greater than or equal.  */

int

__gcc_qge (double a, double aa, double c, double cc)

{

  if (__eqdf2 (a, c) == 0)

    return __gedf2 (aa, cc);

  return __gedf2 (a, c);

}

strong_alias (__gcc_qge, __gcc_qgt);

/* Convert single to long double.  */

long double

__gcc_stoq (float a)

{

  longDblUnion x;

  x.dval[0] = (double) a;

  x.dval[1] = 0.0;

  return x.ldval;

}

/* Convert double to long double.  */

long double

__gcc_dtoq (double a)

{

  longDblUnion x;

  x.dval[0] = a;

  x.dval[1] = 0.0;

  return x.ldval;

}

/* Convert long double to single.  */

float

__gcc_qtos (double a, double aa __attribute__ ((__unused__)))

{

  return (float) a;

}

/* Convert long double to double.  */

double

__gcc_qtod (double a, double aa __attribute__ ((__unused__)))

{

  return a;

}

/* Convert long double to int.  */

int

__gcc_qtoi (double a, double aa)

{

  double z = a + aa;

  return (int) z;

}

/* Convert long double to unsigned int.  */

unsigned int

__gcc_qtou (double a, double aa)

{

  double z = a + aa;

  return (unsigned int) z;

}

/* Convert int to long double.  */

long double

__gcc_itoq (int a)

{

  return __gcc_dtoq ((double) a);

}

/* Convert unsigned int to long double.  */

long double

__gcc_utoq (unsigned int a)

{

  return __gcc_dtoq ((double) a);

}

#endif

#ifdef __NO_FPRS__

int __gcc_qunord (double, double, double, double);

extern int __eqdf2 (double, double);

extern int __unorddf2 (double, double);

/* Compare two 'long double' values for unordered.  */

int

__gcc_qunord (double a, double aa, double c, double cc)

{

  if (__eqdf2 (a, c) == 0)

    return __unorddf2 (aa, cc);

  return __unorddf2 (a, c);

}

#include "config/soft-fp/soft-fp.h"

#include "config/soft-fp/double.h"

#include "config/soft-fp/quad.h"

/* Compute floating point multiply-subtract with higher (quad) precision.  */

static double

fmsub (double a, double b, double c)

{

    FP_DECL_EX;

    FP_DECL_D(A);

    FP_DECL_D(B);

    FP_DECL_D(C);

    FP_DECL_Q(X);

    FP_DECL_Q(Y);

    FP_DECL_Q(Z);

    FP_DECL_Q(U);

    FP_DECL_Q(V);

    FP_DECL_D(R);

    double r;

    long double u, x, y, z;

    FP_INIT_ROUNDMODE;

    FP_UNPACK_RAW_D (A, a);

    FP_UNPACK_RAW_D (B, b);

    FP_UNPACK_RAW_D (C, c);

    /* Extend double to quad.  */

#if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q

    FP_EXTEND(Q,D,4,2,X,A);

    FP_EXTEND(Q,D,4,2,Y,B);

    FP_EXTEND(Q,D,4,2,Z,C);

#else

    FP_EXTEND(Q,D,2,1,X,A);

    FP_EXTEND(Q,D,2,1,Y,B);

    FP_EXTEND(Q,D,2,1,Z,C);

#endif

    FP_PACK_RAW_Q(x,X);

    FP_PACK_RAW_Q(y,Y);

    FP_PACK_RAW_Q(z,Z);

    FP_HANDLE_EXCEPTIONS;

    /* Multiply.  */

    FP_INIT_ROUNDMODE;

    FP_UNPACK_Q(X,x);

    FP_UNPACK_Q(Y,y);

    FP_MUL_Q(U,X,Y);

    FP_PACK_Q(u,U);

    FP_HANDLE_EXCEPTIONS;

    /* Subtract.  */

    FP_INIT_ROUNDMODE;

    FP_UNPACK_SEMIRAW_Q(U,u);

    FP_UNPACK_SEMIRAW_Q(Z,z);

    FP_SUB_Q(V,U,Z);

    /* Truncate quad to double.  */

#if (2 * _FP_W_TYPE_SIZE) < _FP_FRACBITS_Q

    V_f[3] &= 0x0007ffff;

    FP_TRUNC(D,Q,2,4,R,V);

#else

    V_f1 &= 0x0007ffffffffffffL;

    FP_TRUNC(D,Q,1,2,R,V);

#endif

    FP_PACK_SEMIRAW_D(r,R);

    FP_HANDLE_EXCEPTIONS;

    return r;

}

#endif

#endif