/*
 * Copyright (c) 2014 Advanced Micro Devices, Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

#include "math32.h"

/*
 * ====================================================
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunPro, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice
 * is preserved.
 * ====================================================
 */

#define pi   3.1415927410e+00f        /* 0x40490fdb */

#define a0   7.7215664089e-02f        /* 0x3d9e233f */
#define a1   3.2246702909e-01f        /* 0x3ea51a66 */
#define a2   6.7352302372e-02f        /* 0x3d89f001 */
#define a3   2.0580807701e-02f        /* 0x3ca89915 */
#define a4   7.3855509982e-03f        /* 0x3bf2027e */
#define a5   2.8905137442e-03f        /* 0x3b3d6ec6 */
#define a6   1.1927076848e-03f        /* 0x3a9c54a1 */
#define a7   5.1006977446e-04f        /* 0x3a05b634 */
#define a8   2.2086278477e-04f        /* 0x39679767 */
#define a9   1.0801156895e-04f        /* 0x38e28445 */
#define a10  2.5214456400e-05f        /* 0x37d383a2 */
#define a11  4.4864096708e-05f        /* 0x383c2c75 */

#define tc   1.4616321325e+00f        /* 0x3fbb16c3 */

#define tf  -1.2148628384e-01f        /* 0xbdf8cdcd */
/* tt -(tail of tf) */
#define tt   6.6971006518e-09f        /* 0x31e61c52 */

#define t0   4.8383611441e-01f        /* 0x3ef7b95e */
#define t1  -1.4758771658e-01f        /* 0xbe17213c */
#define t2   6.4624942839e-02f        /* 0x3d845a15 */
#define t3  -3.2788541168e-02f        /* 0xbd064d47 */
#define t4   1.7970675603e-02f        /* 0x3c93373d */
#define t5  -1.0314224288e-02f        /* 0xbc28fcfe */
#define t6   6.1005386524e-03f        /* 0x3bc7e707 */
#define t7  -3.6845202558e-03f        /* 0xbb7177fe */
#define t8   2.2596477065e-03f        /* 0x3b141699 */
#define t9  -1.4034647029e-03f        /* 0xbab7f476 */
#define t10  8.8108185446e-04f        /* 0x3a66f867 */
#define t11 -5.3859531181e-04f        /* 0xba0d3085 */
#define t12  3.1563205994e-04f        /* 0x39a57b6b */
#define t13 -3.1275415677e-04f        /* 0xb9a3f927 */
#define t14  3.3552918467e-04f        /* 0x39afe9f7 */

#define u0  -7.7215664089e-02f        /* 0xbd9e233f */
#define u1   6.3282704353e-01f        /* 0x3f2200f4 */
#define u2   1.4549225569e+00f        /* 0x3fba3ae7 */
#define u3   9.7771751881e-01f        /* 0x3f7a4bb2 */
#define u4   2.2896373272e-01f        /* 0x3e6a7578 */
#define u5   1.3381091878e-02f        /* 0x3c5b3c5e */

#define v1   2.4559779167e+00f        /* 0x401d2ebe */
#define v2   2.1284897327e+00f        /* 0x4008392d */
#define v3   7.6928514242e-01f        /* 0x3f44efdf */
#define v4   1.0422264785e-01f        /* 0x3dd572af */
#define v5   3.2170924824e-03f        /* 0x3b52d5db */

#define s0  -7.7215664089e-02f        /* 0xbd9e233f */
#define s1   2.1498242021e-01f        /* 0x3e5c245a */
#define s2   3.2577878237e-01f        /* 0x3ea6cc7a */
#define s3   1.4635047317e-01f        /* 0x3e15dce6 */
#define s4   2.6642270386e-02f        /* 0x3cda40e4 */
#define s5   1.8402845599e-03f        /* 0x3af135b4 */
#define s6   3.1947532989e-05f        /* 0x3805ff67 */

#define r1   1.3920053244e+00f        /* 0x3fb22d3b */
#define r2   7.2193557024e-01f        /* 0x3f38d0c5 */
#define r3   1.7193385959e-01f        /* 0x3e300f6e */
#define r4   1.8645919859e-02f        /* 0x3c98bf54 */
#define r5   7.7794247773e-04f        /* 0x3a4beed6 */
#define r6   7.3266842264e-06f        /* 0x36f5d7bd */

#define w0   4.1893854737e-01f        /* 0x3ed67f1d */
#define w1   8.3333335817e-02f        /* 0x3daaaaab */
#define w2  -2.7777778450e-03f        /* 0xbb360b61 */
#define w3   7.9365057172e-04f        /* 0x3a500cfd */
#define w4  -5.9518753551e-04f        /* 0xba1c065c */
#define w5   8.3633989561e-04f        /* 0x3a5b3dd2 */
#define w6  -1.6309292987e-03f        /* 0xbad5c4e8 */

__attribute__ ((overloadable, always_inline)) float
lgamma_r(float x, int *signp)
{
    int hx = as_int(x);
    int ix = hx & 0x7fffffff;
    float absx = as_float(ix);

    if (ix >= 0x7f800000) {
	*signp = 1;
        return x;
    }

    if (absx < 0x1.0p-70f) {
        *signp = hx < 0 ? -1 : 1;
	return -log(absx);
    }

    float r;

    if (absx == 1.0f | absx == 2.0f)
        r = 0.0f;

    else if (absx < 2.0f) {
	float y = 2.0f - absx;
	int i = 0;

	int c = absx < 0x1.bb4c30p+0f;
	float yt = absx - tc;
	y = c ? yt : y;
	i = c ? 1 : i;

	c = absx < 0x1.3b4c40p+0f;
	yt = absx - 1.0f;
	y = c ? yt : y;
	i = c ? 2 : i;

	r = -log(absx);
	yt = 1.0f - absx;
	c = absx <= 0x1.ccccccp-1f;
	r = c ? r : 0.0f;
	y = c ? yt : y;
	i = c ? 0 : i;

	c = absx < 0x1.769440p-1f;
	yt = absx - (tc - 1.0f);
	y = c ? yt : y;
	i = c ? 1 : i;

	c = absx < 0x1.da6610p-3f;
	y = c ? absx : y;
	i = c ? 2 : i;

	float z, w, p1, p2, p3, p;
        switch (i) {
        case 0:
            z = y * y;
	    p1 = mad(z, mad(z, mad(z, mad(z, mad(z, a10, a8), a6), a4), a2), a0);
	    p2 = z * mad(z, mad(z, mad(z, mad(z, mad(z, a11, a9), a7), a5), a3), a1);
	    p = mad(y, p1, p2);
	    r += mad(y, -0.5f, p);
            break;
        case 1:
            z = y * y;
            w = z * y;
	    p1 = mad(w, mad(w, mad(w, mad(w, t12, t9), t6), t3), t0);
	    p2 = mad(w, mad(w, mad(w, mad(w, t13, t10), t7), t4), t1);
            p3 = mad(w, mad(w, mad(w, mad(w, t14, t11), t8), t5), t2);
            p = mad(z, p1, -mad(w, -mad(y, p3, p2), tt));
            r += tf + p;
            break;
        case 2:
            p1 = y * mad(y, mad(y, mad(y, mad(y, mad(y, u5, u4), u3), u2), u1), u0);
            p2 = mad(y, mad(y, mad(y, mad(y, mad(y, v5, v4), v3), v2), v1), 1.0f);
            r += mad(y, -0.5f, MATH_DIVIDE(p1, p2));
	    break;
        }
    } else if (absx < 8.0f) {
        int i = (int)absx;
        float y = absx - (float) i;
        float p = y * mad(y, mad(y, mad(y, mad(y, mad(y, mad(y, s6, s5), s4), s3), s2), s1), s0);
        float q = mad(y, mad(y, mad(y, mad(y, mad(y, mad(y, r6, r5), r4), r3), r2), r1), 1.0f);
        r = mad(y, 0.5f, MATH_DIVIDE(p, q));

	float y6 = y + 6.0f;
	float y5 = y + 5.0f;
	float y4 = y + 4.0f;
	float y3 = y + 3.0f;
	float y2 = y + 2.0f;

        float z = 1.0f;
	z *= i > 6 ? y6 : 1.0f;
	z *= i > 5 ? y5 : 1.0f;
	z *= i > 4 ? y4 : 1.0f;
	z *= i > 3 ? y3 : 1.0f;
	z *= i > 2 ? y2 : 1.0f;

        r += log(z);
    } else if (absx < 0x1.0p+58f) {
        float z = 1.0f / absx;
        float y = z * z;
        float w = mad(z, mad(y, mad(y, mad(y, mad(y, mad(y, w6, w5), w4), w3), w2), w1), w0);
	r = mad(absx - 0.5f, log(absx) - 1.0f, w);
    } else
        // 2**58 <= x <= Inf
        r = absx * (log(absx) - 1.0f);

    int s = 1;

    if (x < 0.0f) {
        float t = sinpi(x);
        r = log(pi / fabs(t * x)) - r;
	r = t == 0.0f ? as_float(PINFBITPATT_SP32) : r;
	s = t < 0.0f ? -1 : s;
    }

    *signp = s;
    return r;
}

#if __OPENCL_C_VERSION__ < 200
__attribute__ ((overloadable, always_inline)) float
lgamma_r(float x, __local int *signp)
{
    int s;
    float l = lgamma_r(x, &s);
    *signp = s;
    return l;
}

__attribute__ ((overloadable, always_inline)) float
lgamma_r(float x, __global int *signp)
{
    int s;
    float l = lgamma_r(x, &s);
    *signp = s;
    return l;
}
#endif

__attribute__ ((overloadable, always_inline)) float
lgamma(float x)
{
    int s;
    float l = lgamma_r(x, &s);
    return l;
}