"Fossies" - the Fresh Open Source Software Archive

Member "pocl-1.8/lib/kernel/ocml/src/acospiH.cl" (12 Oct 2021, 1533 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 "mathH.h"
    9 
   10 CONSTATTR UGEN(acospi)
   11 
   12 CONSTATTR half
   13 MATH_MANGLE(acospi)(half x)
   14 {
   15     // Computes arccos(x).
   16     // The argument is first reduced by noting that arccos(x)
   17     // is invalid for abs(x) > 1 and arccos(-x) = arccos(x).
   18     // For denormal and small arguments arccos(x) = pi/2 to machine
   19     // accuracy. Remaining argument ranges are handled as follows.
   20     // For abs(x) <= 0.5 use
   21     // arccos(x) = pi/2 - arcsin(x)
   22     // = pi/2 - (x + x^3*R(x^2))
   23     // where R(x^2) is a rational minimax approximation to
   24     // (arcsin(x) - x)/x^3.
   25     // For abs(x) > 0.5 exploit the identity:
   26     // arccos(x) = pi - 2*arcsin(sqrt(1-x)/2)
   27     // together with the above rational approximation, and
   28     // reconstruct the terms carefully.
   29 
   30     const half piinv = 0x1.46p-2h;
   31 
   32     half ax = BUILTIN_ABS_F16(x);
   33 
   34     half rt = MATH_MAD(-0.5h, ax, 0.5h);
   35     half x2 = ax * ax;
   36     half r = ax > 0.5h ? rt : x2;
   37 
   38     half u = r * MATH_MAD(r, 0x1.0b8p-5h, 0x1.a7cp-5h);
   39 
   40     half s = MATH_FAST_SQRT(r);
   41     half ztp = 2.0h * MATH_MAD(s, u, piinv*s);
   42     half ztn = 1.0h - ztp;
   43     half zt =  x < 0.0h ? ztn : ztp;
   44     half z = 0.5h - MATH_MAD(x, u, piinv*x);
   45     z = ax > 0.5h ? zt : z;
   46 
   47     return z;
   48 }
   49