LXR uboot/post/lib_powerpc/fpu/darwin-ldouble.c

   1// SPDX-License-Identifier: GPL-2.0+
   2/*
   3 * Borrowed from GCC 4.2.2 (which still was GPL v2+)
   4 */
   5/* 128-bit long double support routines for Darwin.
   6   Copyright (C) 1993, 2003, 2004, 2005, 2006, 2007
   7   Free Software Foundation, Inc.
   8
   9This file is part of GCC.
  10 */
  11
  12/*
  13 * Implementations of floating-point long double basic arithmetic
  14 * functions called by the IBM C compiler when generating code for
  15 * PowerPC platforms.  In particular, the following functions are
  16 * implemented: __gcc_qadd, __gcc_qsub, __gcc_qmul, and __gcc_qdiv.
  17 * Double-double algorithms are based on the paper "Doubled-Precision
  18 * IEEE Standard 754 Floating-Point Arithmetic" by W. Kahan, February 26,
  19 * 1987.  An alternative published reference is "Software for
  20 * Doubled-Precision Floating-Point Computations", by Seppo Linnainmaa,
  21 * ACM TOMS vol 7 no 3, September 1981, pages 272-283.
  22 */
  23
  24/*
  25 * Each long double is made up of two IEEE doubles.  The value of the
  26 * long double is the sum of the values of the two parts.  The most
  27 * significant part is required to be the value of the long double
  28 * rounded to the nearest double, as specified by IEEE.  For Inf
  29 * values, the least significant part is required to be one of +0.0 or
  30 * -0.0.  No other requirements are made; so, for example, 1.0 may be
  31 * represented as (1.0, +0.0) or (1.0, -0.0), and the low part of a
  32 * NaN is don't-care.
  33 *
  34 * This code currently assumes big-endian.
  35 */
  36
  37#define fabs(x) __builtin_fabs(x)
  38#define isless(x, y) __builtin_isless(x, y)
  39#define inf() __builtin_inf()
  40#define unlikely(x) __builtin_expect((x), 0)
  41#define nonfinite(a) unlikely(!isless(fabs(a), inf()))
  42
  43typedef union {
  44        long double ldval;
  45        double dval[2];
  46} longDblUnion;
  47
  48/* Add two 'long double' values and return the result.  */
  49long double __gcc_qadd(double a, double aa, double c, double cc)
  50{
  51        longDblUnion x;
  52        double z, q, zz, xh;
  53
  54        z = a + c;
  55
  56        if (nonfinite(z)) {
  57                z = cc + aa + c + a;
  58                if (nonfinite(z))
  59                        return z;
  60                x.dval[0] = z;  /* Will always be DBL_MAX.  */
  61                zz = aa + cc;
  62                if (fabs(a) > fabs(c))
  63                        x.dval[1] = a - z + c + zz;
  64                else
  65                        x.dval[1] = c - z + a + zz;
  66        } else {
  67                q = a - z;
  68                zz = q + c + (a - (q + z)) + aa + cc;
  69
  70                /* Keep -0 result.  */
  71                if (zz == 0.0)
  72                        return z;
  73
  74                xh = z + zz;
  75                if (nonfinite(xh))
  76                        return xh;
  77
  78                x.dval[0] = xh;
  79                x.dval[1] = z - xh + zz;
  80        }
  81        return x.ldval;
  82}
  83
  84long double __gcc_qsub(double a, double b, double c, double d)
  85{
  86        return __gcc_qadd(a, b, -c, -d);
  87}
  88
  89long double __gcc_qmul(double a, double b, double c, double d)
  90{
  91        longDblUnion z;
  92        double t, tau, u, v, w;
  93
  94        t = a * c;              /* Highest order double term.  */
  95
  96        if (unlikely(t == 0)    /* Preserve -0.  */
  97            || nonfinite(t))
  98                return t;
  99
 100        /* Sum terms of two highest orders. */
 101
 102        /* Use fused multiply-add to get low part of a * c.  */
 103#ifndef __NO_FPRS__
 104        asm("fmsub %0,%1,%2,%3" : "=f"(tau) : "f"(a), "f"(c), "f"(t));
 105#else
 106        tau = fmsub(a, c, t);
 107#endif
 108        v = a * d;
 109        w = b * c;
 110        tau += v + w;           /* Add in other second-order terms.  */
 111        u = t + tau;
 112
 113        /* Construct long double result.  */
 114        if (nonfinite(u))
 115                return u;
 116        z.dval[0] = u;
 117        z.dval[1] = (t - u) + tau;
 118        return z.ldval;
 119}
 120