"Fossies" - the Fresh Open Source Software Archive

Member "pocl-1.8/lib/kernel/ocml/src/acosD.cl" (12 Oct 2021, 2216 Bytes) of package /linux/misc/pocl-1.8.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Lisp source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file.

    1 /*===--------------------------------------------------------------------------
    2  *                   ROCm Device Libraries
    3  *
    4  * This file is distributed under the University of Illinois Open Source
    5  * License. See LICENSE.TXT for details.
    6  *===------------------------------------------------------------------------*/
    7 
    8 #include "mathD.h"
    9 
   10 #define DOUBLE_SPECIALIZATION
   11 #include "ep.h"
   12 
   13 CONSTATTR double
   14 MATH_MANGLE(acos)(double x)
   15 {
   16     // Computes arccos(x).
   17     // The argument is first reduced by noting that arccos(x)
   18     // is invalid for abs(x) > 1. For denormal and small
   19     // arguments arccos(x) = pi/2 to machine accuracy.
   20     // Remaining argument ranges are handled as follows.
   21     // For abs(x) <= 0.5 use
   22     // arccos(x) = pi/2 - arcsin(x)
   23     // = pi/2 - (x + x^3*R(x^2))
   24     // where R(x^2) is a rational minimax approximation to
   25     // (arcsin(x) - x)/x^3.
   26     // For abs(x) > 0.5 exploit the identity:
   27     // arccos(x) = pi - 2*arcsin(sqrt(1-x)/2)
   28     // together with the above rational approximation, and
   29     // reconstruct the terms carefully.
   30 
   31     double y = BUILTIN_ABS_F64(x);
   32     bool transform = y >= 0.5;
   33 
   34     double rt = MATH_MAD(y, -0.5, 0.5);
   35     double y2 = y * y;
   36     double r = transform ? rt : y2;
   37 
   38     double u = r * MATH_MAD(r, MATH_MAD(r, MATH_MAD(r, MATH_MAD(r, 
   39                    MATH_MAD(r, MATH_MAD(r, MATH_MAD(r, MATH_MAD(r, 
   40                    MATH_MAD(r, MATH_MAD(r, MATH_MAD(r, 
   41                        0x1.059859fea6a70p-5, -0x1.0a5a378a05eafp-6), 0x1.4052137024d6ap-6), 0x1.ab3a098a70509p-8),
   42                        0x1.8ed60a300c8d2p-7), 0x1.c6fa84b77012bp-7), 0x1.1c6c111dccb70p-6), 0x1.6e89f0a0adacfp-6),
   43                        0x1.f1c72c668963fp-6), 0x1.6db6db41ce4bdp-5), 0x1.333333336fd5bp-4), 0x1.5555555555380p-3);
   44 
   45     double z = MATH_MAD(0x1.dd9ad336a0500p-1, 0x1.af154eeb562d6p+0, -MATH_MAD(x, u, x));
   46     if (transform) {
   47         double2 s = root2(r);
   48         double zm = MATH_MAD(0x1.dd9ad336a0500p+0, 0x1.af154eeb562d6p+0, -2.0*MATH_MAD(s.hi, u, s.hi));
   49         double zp = 2.0 * (s.hi + MATH_MAD(s.hi, u, s.lo));
   50         z = x < 0.0 ? zm : zp;
   51         z = x == -1.0 ? 0x1.921fb54442d18p+1 : z;
   52         z = x == 1.0 ? 0.0 : z;
   53     }
   54 
   55     return z;
   56 }
   57