Maneesh Gupta
1517 строки
43 KiB
C++
1517 строки
43 KiB
C++
/*
|
|
Copyright (c) 2015-2017 Advanced Micro Devices, Inc. All rights reserved.
|
|
|
|
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 <hc.hpp>
|
|
#include <grid_launch.h>
|
|
#include <hc_math.hpp>
|
|
#include "device_util.h"
|
|
#include "hip/hcc_detail/device_functions.h"
|
|
#include "hip/hip_runtime.h"
|
|
#include <atomic>
|
|
|
|
//=================================================================================================
|
|
/*
|
|
Implementation of malloc and free device functions.
|
|
|
|
This is the best place to put them because the device
|
|
global variables need to be initialized at the start.
|
|
*/
|
|
__device__ char gpuHeap[SIZE_OF_HEAP];
|
|
__device__ uint32_t gpuFlags[NUM_PAGES];
|
|
|
|
__device__ void* __hip_hc_malloc(size_t size) {
|
|
char* heap = (char*)gpuHeap;
|
|
if (size > SIZE_OF_HEAP) {
|
|
return (void*)nullptr;
|
|
}
|
|
uint32_t totalThreads =
|
|
blockDim.x * gridDim.x * blockDim.y * gridDim.y * blockDim.z * gridDim.z;
|
|
uint32_t currentWorkItem = threadIdx.x + blockDim.x * blockIdx.x;
|
|
|
|
uint32_t numHeapsPerWorkItem = NUM_PAGES / totalThreads;
|
|
uint32_t heapSizePerWorkItem = SIZE_OF_HEAP / totalThreads;
|
|
|
|
uint32_t stride = size / SIZE_OF_PAGE;
|
|
uint32_t start = numHeapsPerWorkItem * currentWorkItem;
|
|
|
|
uint32_t k = 0;
|
|
|
|
while (gpuFlags[k] > 0) {
|
|
k++;
|
|
}
|
|
|
|
for (uint32_t i = 0; i < stride - 1; i++) {
|
|
gpuFlags[i + start + k] = 1;
|
|
}
|
|
|
|
gpuFlags[start + stride - 1 + k] = 2;
|
|
|
|
void* ptr = (void*)(heap + heapSizePerWorkItem * currentWorkItem + k * SIZE_OF_PAGE);
|
|
|
|
return ptr;
|
|
}
|
|
|
|
__device__ void* __hip_hc_free(void* ptr) {
|
|
if (ptr == nullptr) {
|
|
return nullptr;
|
|
}
|
|
|
|
uint32_t offsetByte = (uint64_t)ptr - (uint64_t)gpuHeap;
|
|
uint32_t offsetPage = offsetByte / SIZE_OF_PAGE;
|
|
|
|
while (gpuFlags[offsetPage] != 0) {
|
|
if (gpuFlags[offsetPage] == 2) {
|
|
gpuFlags[offsetPage] = 0;
|
|
offsetPage++;
|
|
break;
|
|
} else {
|
|
gpuFlags[offsetPage] = 0;
|
|
offsetPage++;
|
|
}
|
|
}
|
|
|
|
return nullptr;
|
|
}
|
|
|
|
|
|
// loop unrolling
|
|
__device__ void* __hip_hc_memcpy(void* dst, const void* src, size_t size) {
|
|
auto dstPtr = static_cast<uint8_t*>(dst);
|
|
auto srcPtr = static_cast<const uint8_t*>(src);
|
|
|
|
while (size >= 4u) {
|
|
dstPtr[0] = srcPtr[0];
|
|
dstPtr[1] = srcPtr[1];
|
|
dstPtr[2] = srcPtr[2];
|
|
dstPtr[3] = srcPtr[3];
|
|
|
|
size -= 4u;
|
|
srcPtr += 4u;
|
|
dstPtr += 4u;
|
|
}
|
|
switch (size) {
|
|
case 3:
|
|
dstPtr[2] = srcPtr[2];
|
|
case 2:
|
|
dstPtr[1] = srcPtr[1];
|
|
case 1:
|
|
dstPtr[0] = srcPtr[0];
|
|
}
|
|
|
|
return dst;
|
|
}
|
|
|
|
__device__ void* __hip_hc_memset(void* dst, uint8_t val, size_t size) {
|
|
auto dstPtr = static_cast<uint8_t*>(dst);
|
|
|
|
while (size >= 4u) {
|
|
dstPtr[0] = val;
|
|
dstPtr[1] = val;
|
|
dstPtr[2] = val;
|
|
dstPtr[3] = val;
|
|
|
|
size -= 4u;
|
|
dstPtr += 4u;
|
|
}
|
|
switch (size) {
|
|
case 3:
|
|
dstPtr[2] = val;
|
|
case 2:
|
|
dstPtr[1] = val;
|
|
case 1:
|
|
dstPtr[0] = val;
|
|
}
|
|
|
|
return dst;
|
|
}
|
|
|
|
__device__ float __hip_erfinvf(float x) { return hc::precise_math::erfinvf(x); }
|
|
|
|
__device__ double __hip_erfinv(double x) { return hc::precise_math::erfinv(x); }
|
|
|
|
#define __hip_j0a1 57568490574.0
|
|
#define __hip_j0a2 -13362590354.0
|
|
#define __hip_j0a3 651619640.7
|
|
#define __hip_j0a4 -11214424.18
|
|
#define __hip_j0a5 77392.33017
|
|
#define __hip_j0a6 -184.9052456
|
|
|
|
#define __hip_j0b1 57568490411.0
|
|
#define __hip_j0b2 1029532985.0
|
|
#define __hip_j0b3 9494680.718
|
|
#define __hip_j0b4 59272.64853
|
|
#define __hip_j0b5 267.8532712
|
|
|
|
#define __hip_j0c 0.785398164
|
|
#define __hip_j0c1 -0.1098628627e-2
|
|
#define __hip_j0c2 0.2734510407e-4
|
|
#define __hip_j0c3 -0.2073370639e-5
|
|
#define __hip_j0c4 0.2093887211e-6
|
|
|
|
#define __hip_j0d1 -0.1562499995e-1
|
|
#define __hip_j0d2 0.1430488765e-3
|
|
#define __hip_j0d3 0.6911147651e-5
|
|
#define __hip_j0d4 0.7621095161e-6
|
|
#define __hip_j0d5 0.934935152e-7
|
|
|
|
#define __hip_j0e 0.636619772
|
|
|
|
__device__ double __hip_j0(double x) {
|
|
double ret, a = hc::precise_math::fabs(x);
|
|
if (a < 8.0) {
|
|
double y = x * x;
|
|
double y1 = __hip_j0a6 * y + __hip_j0a5;
|
|
double z1 = 1.0 * y + __hip_j0b5;
|
|
|
|
double y2 = y1 * y + __hip_j0a4;
|
|
double z2 = z1 * y + __hip_j0b4;
|
|
|
|
double y3 = y2 * y + __hip_j0a3;
|
|
double z3 = z2 * y + __hip_j0b3;
|
|
|
|
double y4 = y3 * y + __hip_j0a2;
|
|
double z4 = z3 * y + __hip_j0b2;
|
|
|
|
double y5 = y4 * y + __hip_j0a1;
|
|
double z5 = z4 * y + __hip_j0b1;
|
|
|
|
ret = y5 / z5;
|
|
|
|
} else {
|
|
double z = 8.0 / a;
|
|
double y = z * z;
|
|
double x1 = a - __hip_j0c;
|
|
|
|
double y1 = __hip_j0c4 * y + __hip_j0c3;
|
|
double z1 = __hip_j0d5 * y + __hip_j0d4;
|
|
|
|
double y2 = y1 * y + __hip_j0c2;
|
|
double z2 = z1 * z + __hip_j0d3;
|
|
|
|
double y3 = y2 * y + __hip_j0c1;
|
|
double z3 = z2 * y + __hip_j0d2;
|
|
|
|
double y4 = y3 * y + 1.0;
|
|
double z4 = z3 * y + __hip_j0d1;
|
|
|
|
ret = hc::precise_math::sqrt(__hip_j0e / a) *
|
|
(hc::precise_math::cos(x1) * y4 - z * hc::precise_math::sin(x1) * z4);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
__device__ float __hip_j0f(float x) {
|
|
float ret, a = hc::precise_math::fabsf(x);
|
|
if (a < 8.0) {
|
|
float y = x * x;
|
|
float y1 = __hip_j0a6 * y + __hip_j0a5;
|
|
float z1 = 1.0 * y + __hip_j0b5;
|
|
|
|
float y2 = y1 * y + __hip_j0a4;
|
|
float z2 = z1 * y + __hip_j0b4;
|
|
|
|
float y3 = y2 * y + __hip_j0a3;
|
|
float z3 = z2 * y + __hip_j0b3;
|
|
|
|
float y4 = y3 * y + __hip_j0a2;
|
|
float z4 = z3 * y + __hip_j0b2;
|
|
|
|
float y5 = y4 * y + __hip_j0a1;
|
|
float z5 = z4 * y + __hip_j0b1;
|
|
|
|
ret = y5 / z5;
|
|
|
|
} else {
|
|
float z = 8.0 / a;
|
|
float y = z * z;
|
|
float x1 = a - __hip_j0c;
|
|
|
|
float y1 = __hip_j0c4 * y + __hip_j0c3;
|
|
float z1 = __hip_j0d5 * y + __hip_j0d4;
|
|
|
|
float y2 = y1 * y + __hip_j0c2;
|
|
float z2 = z1 * z + __hip_j0d3;
|
|
|
|
float y3 = y2 * y + __hip_j0c1;
|
|
float z3 = z2 * y + __hip_j0d2;
|
|
|
|
float y4 = y3 * y + 1.0;
|
|
float z4 = z3 * y + __hip_j0d1;
|
|
|
|
ret = hc::precise_math::sqrtf(__hip_j0e / a) *
|
|
(hc::precise_math::cosf(x1) * y4 - z * hc::precise_math::sinf(x1) * z4);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
#define __hip_j1a1 -30.16036606
|
|
#define __hip_j1a2 15704.48260
|
|
#define __hip_j1a3 -2972611.439
|
|
#define __hip_j1a4 242396853.1
|
|
#define __hip_j1a5 -7895059235.0
|
|
#define __hip_j1a6 72362614232.0
|
|
|
|
#define __hip_j1b1 376.9991397
|
|
#define __hip_j1b2 99447.43394
|
|
#define __hip_j1b3 18583304.74
|
|
#define __hip_j1b4 2300535178.0
|
|
#define __hip_j1b5 144725228442.0
|
|
|
|
#define __hip_j1c 2.356194491
|
|
#define __hip_j1c1 -0.240337019e-6
|
|
#define __hip_j1c2 0.2457520174e-5
|
|
#define __hip_j1c3 -0.3516396496e-4
|
|
#define __hip_j1c4 0.183105e-2
|
|
|
|
#define __hip_j1d1 0.105787412e-6
|
|
#define __hip_j1d2 -0.88228987e-6
|
|
#define __hip_j1d3 0.8449199096e-5
|
|
#define __hip_j1d4 -0.2002690873e-3
|
|
#define __hip_j1d5 0.04687499995
|
|
|
|
#define __hip_j1e 0.636619772
|
|
|
|
__device__ double __hip_j1(double x) {
|
|
double ret, a = hc::precise_math::fabs(x);
|
|
if (a < 8.0) {
|
|
double y = x * x;
|
|
|
|
double y1 = __hip_j1a1 * y + __hip_j1a2;
|
|
double z1 = 1.0 * y + __hip_j1b1;
|
|
|
|
double y2 = y1 * y + __hip_j1a3;
|
|
double z2 = z1 * y + __hip_j1b2;
|
|
|
|
double y3 = y2 * y + __hip_j1a4;
|
|
double z3 = z2 * y + __hip_j1b3;
|
|
|
|
double y4 = y3 * y + __hip_j1a5;
|
|
double z4 = z3 * y + __hip_j1b4;
|
|
|
|
double y5 = y4 * y + __hip_j1a6;
|
|
double z5 = z4 * y + __hip_j1b5;
|
|
|
|
ret = x * y5 / z5;
|
|
|
|
} else {
|
|
double z = 8.0 / a;
|
|
double y = z * z;
|
|
double x1 = a - __hip_j1c;
|
|
|
|
double y1 = __hip_j1c1 * y + __hip_j1c2;
|
|
double y2 = y1 * y + __hip_j1c3;
|
|
double y3 = y2 * y + __hip_j1c4;
|
|
double y4 = y3 * y + 1.0;
|
|
|
|
double z1 = __hip_j1d1 * y + __hip_j1d2;
|
|
double z2 = z1 * y + __hip_j1d3;
|
|
double z3 = z2 * y + __hip_j1d4;
|
|
double z4 = z3 * y + __hip_j1d5;
|
|
|
|
ret = hc::precise_math::sqrt(__hip_j1e / a) *
|
|
(hc::precise_math::cos(x1) * y4 - z * hc::precise_math::sin(x1) * z4);
|
|
if (x < 0.0) ret = -ret;
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
__device__ float __hip_j1f(float x) {
|
|
double ret, a = hc::precise_math::fabsf(x);
|
|
if (a < 8.0) {
|
|
float y = x * x;
|
|
|
|
float y1 = __hip_j1a1 * y + __hip_j1a2;
|
|
float z1 = 1.0 * y + __hip_j1b1;
|
|
|
|
float y2 = y1 * y + __hip_j1a3;
|
|
float z2 = z1 * y + __hip_j1b2;
|
|
|
|
float y3 = y2 * y + __hip_j1a4;
|
|
float z3 = z2 * y + __hip_j1b3;
|
|
|
|
float y4 = y3 * y + __hip_j1a5;
|
|
float z4 = z3 * y + __hip_j1b4;
|
|
|
|
float y5 = y4 * y + __hip_j1a6;
|
|
float z5 = z4 * y + __hip_j1b5;
|
|
|
|
ret = x * y5 / z5;
|
|
|
|
} else {
|
|
float z = 8.0 / a;
|
|
float y = z * z;
|
|
float x1 = a - __hip_j1c;
|
|
|
|
float y1 = __hip_j1c1 * y + __hip_j1c2;
|
|
float y2 = y1 * y + __hip_j1c3;
|
|
float y3 = y2 * y + __hip_j1c4;
|
|
float y4 = y3 * y + 1.0;
|
|
|
|
float z1 = __hip_j1d1 * y + __hip_j1d2;
|
|
float z2 = z1 * y + __hip_j1d3;
|
|
float z3 = z2 * y + __hip_j1d4;
|
|
float z4 = z3 * y + __hip_j1d5;
|
|
|
|
ret = hc::precise_math::sqrtf(__hip_j1e / a) *
|
|
(hc::precise_math::cosf(x1) * y4 - z * hc::precise_math::sinf(x1) * z4);
|
|
if (x < 0.0) ret = -ret;
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
#define __hip_y0a1 228.4622733
|
|
#define __hip_y0a2 -86327.92757
|
|
#define __hip_y0a3 10879881.29
|
|
#define __hip_y0a4 -512359803.6
|
|
#define __hip_y0a5 7062834065.0
|
|
#define __hip_y0a6 -2957821389.0
|
|
|
|
#define __hip_y0b1 226.1030244
|
|
#define __hip_y0b2 47447.26470
|
|
#define __hip_y0b3 7189466.438
|
|
#define __hip_y0b4 745249964.8
|
|
#define __hip_y0b5 40076544269.0
|
|
|
|
#define __hip_y0c 0.636619772
|
|
|
|
#define __hip_y0d 0.785398164
|
|
|
|
#define __hip_y0e1 0.2093887211e-6
|
|
#define __hip_y0e2 -0.2073370639e-5
|
|
#define __hip_y0e3 0.2734510407e-4
|
|
#define __hip_y0e4 -0.1098628627e-2
|
|
|
|
|
|
#define __hip_y0f1 -0.934945152e-7
|
|
#define __hip_y0f2 0.7621095161e-6
|
|
#define __hip_y0f3 -0.6911147651e-5
|
|
#define __hip_y0f4 0.1430488765e-3
|
|
#define __hip_y0f5 -0.1562499995e-1
|
|
|
|
#define __hip_y1g 0.636619772
|
|
|
|
__device__ double __hip_y0(double x) {
|
|
double ret;
|
|
|
|
if (x < 8.0) {
|
|
double y = x * x;
|
|
double y1 = __hip_y0a1 * y + __hip_y0a2;
|
|
double y2 = y1 * y + __hip_y0a3;
|
|
double y3 = y2 * y + __hip_y0a4;
|
|
double y4 = y3 * y + __hip_y0a5;
|
|
double y5 = y4 * y + __hip_y0a6;
|
|
|
|
double z1 = 1.0 * y + __hip_y0b1;
|
|
double z2 = z1 * y + __hip_y0b2;
|
|
double z3 = z2 * y + __hip_y0b3;
|
|
double z4 = z3 * y + __hip_y0b4;
|
|
double z5 = z4 * y + __hip_y0b5;
|
|
|
|
|
|
ret = (y5 / z5) + __hip_y0c * __hip_j0(x) * hc::precise_math::log(x);
|
|
} else {
|
|
double z = 8.0 / x;
|
|
double y = z * z;
|
|
double x1 = x - __hip_y0d;
|
|
|
|
double y1 = __hip_y0e1 * y + __hip_y0e2;
|
|
double y2 = y1 * y + __hip_y0e3;
|
|
double y3 = y2 * y + __hip_y0e4;
|
|
double y4 = y3 * y + 1.0;
|
|
|
|
double z1 = __hip_y0f1 * y + __hip_y0f2;
|
|
double z2 = z1 * y + __hip_y0f3;
|
|
double z3 = z2 * y + __hip_y0f4;
|
|
double z4 = z3 * y + __hip_y0f5;
|
|
|
|
ret = hc::precise_math::sqrt(__hip_y1g / x) *
|
|
(hc::precise_math::sin(x1) * y4 + z * hc::precise_math::cos(x1) * z4);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
|
|
__device__ float __hip_y0f(float x) {
|
|
float ret;
|
|
|
|
if (x < 8.0) {
|
|
float y = x * x;
|
|
float y1 = __hip_y0a1 * y + __hip_y0a2;
|
|
float y2 = y1 * y + __hip_y0a3;
|
|
float y3 = y2 * y + __hip_y0a4;
|
|
float y4 = y3 * y + __hip_y0a5;
|
|
float y5 = y4 * y + __hip_y0a6;
|
|
|
|
float z1 = 1.0 * y + __hip_y0b1;
|
|
float z2 = z1 * y + __hip_y0b2;
|
|
float z3 = z2 * y + __hip_y0b3;
|
|
float z4 = z3 * y + __hip_y0b4;
|
|
float z5 = z4 * y + __hip_y0b5;
|
|
|
|
|
|
ret = (y5 / z5) + __hip_y0c * __hip_j0f(x) * hc::precise_math::logf(x);
|
|
} else {
|
|
float z = 8.0 / x;
|
|
float y = z * z;
|
|
float x1 = x - __hip_y0d;
|
|
|
|
float y1 = __hip_y0e1 * y + __hip_y0e2;
|
|
float y2 = y1 * y + __hip_y0e3;
|
|
float y3 = y2 * y + __hip_y0e4;
|
|
float y4 = y3 * y + 1.0;
|
|
|
|
float z1 = __hip_y0f1 * y + __hip_y0f2;
|
|
float z2 = z1 * y + __hip_y0f3;
|
|
float z3 = z2 * y + __hip_y0f4;
|
|
float z4 = z3 * y + __hip_y0f5;
|
|
|
|
ret = hc::precise_math::sqrtf(__hip_y1g / x) *
|
|
(hc::precise_math::sinf(x1) * y4 + z * hc::precise_math::cosf(x1) * z4);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
#define __hip_y1a1 0.8511937935e4
|
|
#define __hip_y1a2 -0.4237922726e7
|
|
#define __hip_y1a3 0.7349264551e9
|
|
#define __hip_y1a4 -0.5153438139e11
|
|
#define __hip_y1a5 0.1275274390e13
|
|
#define __hip_y1a6 -0.4900604943e13
|
|
|
|
#define __hip_y1b1 0.3549632885e3
|
|
#define __hip_y1b2 0.1020426050e6
|
|
#define __hip_y1b3 0.2245904002e8
|
|
#define __hip_y1b4 0.3733650367e10
|
|
#define __hip_y1b5 0.4244419664e12
|
|
#define __hip_y1b6 0.2499580570e14
|
|
|
|
#define __hip_y1c 0.636619772
|
|
|
|
#define __hip_y1d 2.356194491
|
|
|
|
#define __hip_y1e1 -0.240337019e-6
|
|
#define __hip_y1e2 0.2457520174e-5
|
|
#define __hip_y1e3 -0.3516396496e-4
|
|
#define __hip_y1e4 0.183105e-2
|
|
|
|
#define __hip_y1f1 0.105787412e-6
|
|
#define __hip_y1f2 -0.88228987e-6
|
|
#define __hip_y1f3 0.8449199096e-5
|
|
#define __hip_y1f4 -0.2002690873e-3
|
|
#define __hip_y1f5 0.04687499995
|
|
|
|
#define __hip_y1g 0.636619772
|
|
|
|
__device__ double __hip_y1(double x) {
|
|
double ret;
|
|
|
|
if (x < 8.0) {
|
|
double y = x * x;
|
|
|
|
double y1 = __hip_y1a1 * y + __hip_y1a2;
|
|
double y2 = y1 * y + __hip_y1a3;
|
|
double y3 = y2 * y + __hip_y1a4;
|
|
double y4 = y3 * y + __hip_y1a5;
|
|
double y5 = y4 * y + __hip_y1a6;
|
|
double y6 = x * y5;
|
|
|
|
double z1 = __hip_y1b1 + y;
|
|
double z2 = z1 * y + __hip_y1b2;
|
|
double z3 = z2 * y + __hip_y1b3;
|
|
double z4 = z3 * y + __hip_y1b4;
|
|
double z5 = z4 * y + __hip_y1b5;
|
|
double z6 = z5 * y + __hip_y1b6;
|
|
|
|
ret = (y6 / z6) + __hip_y1c * (__hip_j1(x) * hc::precise_math::log(x) - 1.0 / x);
|
|
} else {
|
|
double z = 8.0 / x;
|
|
double y = z * z;
|
|
double x1 = x - __hip_y1d;
|
|
|
|
double y1 = __hip_y1e1 * y + __hip_y1e2;
|
|
double y2 = y1 * y + __hip_y1e3;
|
|
double y3 = y2 * y + __hip_y1e4;
|
|
double y4 = y3 * y + 1.0;
|
|
|
|
double z1 = __hip_y1f1 * y + __hip_y1f2;
|
|
double z2 = z1 * y + __hip_y1f3;
|
|
double z3 = z2 * y + __hip_y1f4;
|
|
double z4 = z3 * y + __hip_y1f5;
|
|
|
|
ret = hc::precise_math::sqrt(__hip_y1g / x) *
|
|
(hc::precise_math::sin(x1) * y4 + z * hc::precise_math::cos(x1) * z4);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
__device__ float __hip_y1f(float x) {
|
|
float ret;
|
|
|
|
if (x < 8.0) {
|
|
float y = x * x;
|
|
|
|
float y1 = __hip_y1a1 * y + __hip_y1a2;
|
|
float y2 = y1 * y + __hip_y1a3;
|
|
float y3 = y2 * y + __hip_y1a4;
|
|
float y4 = y3 * y + __hip_y1a5;
|
|
float y5 = y4 * y + __hip_y1a6;
|
|
float y6 = x * y5;
|
|
|
|
float z1 = __hip_y1b1 + y;
|
|
float z2 = z1 * y + __hip_y1b2;
|
|
float z3 = z2 * y + __hip_y1b3;
|
|
float z4 = z3 * y + __hip_y1b4;
|
|
float z5 = z4 * y + __hip_y1b5;
|
|
float z6 = z5 * y + __hip_y1b6;
|
|
|
|
ret = (y6 / z6) + __hip_y1c * (__hip_j1f(x) * hc::precise_math::logf(x) - 1.0 / x);
|
|
} else {
|
|
float z = 8.0 / x;
|
|
float y = z * z;
|
|
float x1 = x - __hip_y1d;
|
|
|
|
float y1 = __hip_y1e1 * y + __hip_y1e2;
|
|
float y2 = y1 * y + __hip_y1e3;
|
|
float y3 = y2 * y + __hip_y1e4;
|
|
float y4 = y3 * y + 1.0;
|
|
|
|
float z1 = __hip_y1f1 * y + __hip_y1f2;
|
|
float z2 = z1 * y + __hip_y1f3;
|
|
float z3 = z2 * y + __hip_y1f4;
|
|
float z4 = z3 * y + __hip_y1f5;
|
|
|
|
ret = hc::precise_math::sqrtf(__hip_y1g / x) *
|
|
(hc::precise_math::sinf(x1) * y4 + z * hc::precise_math::cosf(x1) * z4);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
#define __hip_k1 40.0
|
|
#define __hip_k2 1.0e10
|
|
#define __hip_k3 1.0e-10
|
|
|
|
__device__ double __hip_jn(int n, double x) {
|
|
int sum = 0, m;
|
|
double a, b0, b1, b2, val, t, ret;
|
|
if (n < 0) {
|
|
return NAN;
|
|
}
|
|
a = hc::precise_math::fabs(x);
|
|
if (n == 0) {
|
|
return (__hip_j0(a));
|
|
}
|
|
if (n == 1) {
|
|
return (__hip_j1(a));
|
|
}
|
|
if (a == 0.0) {
|
|
return 0.0;
|
|
} else if (a > (double)n) {
|
|
t = 2.0 / a;
|
|
b1 = __hip_j0(a);
|
|
b0 = __hip_j1(a);
|
|
for (int i = 1; i < n; i++) {
|
|
b2 = i * t * b0 - b1;
|
|
b1 = b0;
|
|
b0 = b2;
|
|
}
|
|
ret = b0;
|
|
} else {
|
|
t = 2.0 / a;
|
|
m = 2 * ((n + (int)hc::precise_math::sqrt(__hip_k1 * n)) / 2);
|
|
b2 = ret = val = 0.0;
|
|
b0 = 1.0;
|
|
for (int i = m; i > 0; i--) {
|
|
b1 = i * t * b0 - b2;
|
|
b2 = b0;
|
|
b0 = b1;
|
|
if (hc::precise_math::fabs(b0) > __hip_k2) {
|
|
b0 *= __hip_k3;
|
|
b2 *= __hip_k3;
|
|
ret *= __hip_k3;
|
|
val *= __hip_k3;
|
|
}
|
|
if (sum) val += b0;
|
|
sum = !sum;
|
|
if (i == n) ret = b2;
|
|
}
|
|
val = 2.0 * val - b0;
|
|
ret /= val;
|
|
}
|
|
return x < 0.0 && n % 2 == 1 ? -ret : ret;
|
|
}
|
|
|
|
__device__ float __hip_jnf(int n, float x) {
|
|
int sum = 0, m;
|
|
float a, b0, b1, b2, val, t, ret;
|
|
if (n < 0) {
|
|
return NAN;
|
|
}
|
|
a = hc::precise_math::fabsf(x);
|
|
if (n == 0) {
|
|
return (__hip_j0f(a));
|
|
}
|
|
if (n == 1) {
|
|
return (__hip_j1f(a));
|
|
}
|
|
if (a == 0.0) {
|
|
return 0.0;
|
|
} else if (a > (float)n) {
|
|
t = 2.0 / a;
|
|
b1 = __hip_j0f(a);
|
|
b0 = __hip_j1f(a);
|
|
for (int i = 1; i < n; i++) {
|
|
b2 = i * t * b0 - b1;
|
|
b1 = b0;
|
|
b0 = b2;
|
|
}
|
|
ret = b0;
|
|
} else {
|
|
t = 2.0 / a;
|
|
m = 2 * ((n + (int)hc::precise_math::sqrtf(__hip_k1 * n)) / 2);
|
|
b2 = ret = val = 0.0;
|
|
b0 = 1.0;
|
|
for (int i = m; i > 0; i--) {
|
|
b1 = i * t * b0 - b2;
|
|
b2 = b0;
|
|
b0 = b1;
|
|
if (hc::precise_math::fabsf(b0) > __hip_k2) {
|
|
b0 *= __hip_k3;
|
|
b2 *= __hip_k3;
|
|
ret *= __hip_k3;
|
|
val *= __hip_k3;
|
|
}
|
|
if (sum) val += b0;
|
|
sum = !sum;
|
|
if (i == n) ret = b2;
|
|
}
|
|
val = 2.0 * val - b0;
|
|
ret /= val;
|
|
}
|
|
return x < 0.0 && n % 2 == 1 ? -ret : ret;
|
|
}
|
|
|
|
__device__ double __hip_yn(int n, double x) {
|
|
double b0, b1, b2, t;
|
|
|
|
if (n < 0 || x == 0.0) {
|
|
return NAN;
|
|
}
|
|
if (n == 0) {
|
|
return (__hip_y0(x));
|
|
}
|
|
if (n == 1) {
|
|
return (__hip_y1(x));
|
|
}
|
|
t = 2.0 / x;
|
|
b0 = __hip_y1(x);
|
|
b1 = __hip_y0(x);
|
|
for (int i = 1; i < n; i++) {
|
|
b2 = i * t * b0 - b1;
|
|
b1 = b0;
|
|
b0 = b2;
|
|
}
|
|
return b0;
|
|
}
|
|
|
|
__device__ float __hip_ynf(int n, float x) {
|
|
float b0, b1, b2, t;
|
|
|
|
if (n < 0 || x == 0.0) {
|
|
return NAN;
|
|
}
|
|
if (n == 0) {
|
|
return (__hip_y0f(x));
|
|
}
|
|
if (n == 1) {
|
|
return (__hip_y1f(x));
|
|
}
|
|
t = 2.0 / x;
|
|
b0 = __hip_y1f(x);
|
|
b1 = __hip_y0f(x);
|
|
for (int i = 1; i < n; i++) {
|
|
b2 = i * t * b0 - b1;
|
|
b1 = b0;
|
|
b0 = b2;
|
|
}
|
|
return b0;
|
|
}
|
|
|
|
__device__ long long int clock64() { return (long long int)hc::__cycle_u64(); };
|
|
__device__ clock_t clock() { return (clock_t)hc::__cycle_u64(); };
|
|
|
|
// abort
|
|
__device__ void abort() { return hc::abort(); }
|
|
|
|
// atomicAdd()
|
|
__device__ int atomicAdd(int* address, int val) { return hc::atomic_fetch_add(address, val); }
|
|
__device__ unsigned int atomicAdd(unsigned int* address, unsigned int val) {
|
|
return hc::atomic_fetch_add(address, val);
|
|
}
|
|
__device__ unsigned long long int atomicAdd(unsigned long long int* address,
|
|
unsigned long long int val) {
|
|
return (long long int)hc::atomic_fetch_add((uint64_t*)address, (uint64_t)val);
|
|
}
|
|
__device__ float atomicAdd(float* address, float val) { return hc::atomic_fetch_add(address, val); }
|
|
|
|
// atomicSub()
|
|
__device__ int atomicSub(int* address, int val) { return hc::atomic_fetch_sub(address, val); }
|
|
__device__ unsigned int atomicSub(unsigned int* address, unsigned int val) {
|
|
return hc::atomic_fetch_sub(address, val);
|
|
}
|
|
|
|
// atomicExch()
|
|
__device__ int atomicExch(int* address, int val) { return hc::atomic_exchange(address, val); }
|
|
__device__ unsigned int atomicExch(unsigned int* address, unsigned int val) {
|
|
return hc::atomic_exchange(address, val);
|
|
}
|
|
__device__ unsigned long long int atomicExch(unsigned long long int* address,
|
|
unsigned long long int val) {
|
|
return (long long int)hc::atomic_exchange((uint64_t*)address, (uint64_t)val);
|
|
}
|
|
__device__ float atomicExch(float* address, float val) { return hc::atomic_exchange(address, val); }
|
|
|
|
// atomicMin()
|
|
__device__ int atomicMin(int* address, int val) { return hc::atomic_fetch_min(address, val); }
|
|
__device__ unsigned int atomicMin(unsigned int* address, unsigned int val) {
|
|
return hc::atomic_fetch_min(address, val);
|
|
}
|
|
__device__ unsigned long long int atomicMin(unsigned long long int* address,
|
|
unsigned long long int val) {
|
|
return (long long int)hc::atomic_fetch_min((uint64_t*)address, (uint64_t)val);
|
|
}
|
|
|
|
// atomicMax()
|
|
__device__ int atomicMax(int* address, int val) { return hc::atomic_fetch_max(address, val); }
|
|
__device__ unsigned int atomicMax(unsigned int* address, unsigned int val) {
|
|
return hc::atomic_fetch_max(address, val);
|
|
}
|
|
__device__ unsigned long long int atomicMax(unsigned long long int* address,
|
|
unsigned long long int val) {
|
|
return (long long int)hc::atomic_fetch_max((uint64_t*)address, (uint64_t)val);
|
|
}
|
|
|
|
// atomicCAS()
|
|
template <typename T>
|
|
__device__ T atomicCAS_impl(T* address, T compare, T val) {
|
|
// the implementation assumes the atomic is lock-free and
|
|
// has the same size as the non-atmoic equivalent type
|
|
static_assert(sizeof(T) == sizeof(std::atomic<T>),
|
|
"size mismatch between atomic and non-atomic types");
|
|
|
|
union {
|
|
T* address;
|
|
std::atomic<T>* atomic_address;
|
|
} u;
|
|
u.address = address;
|
|
|
|
T expected = compare;
|
|
|
|
// hcc should generate a system scope atomic CAS
|
|
std::atomic_compare_exchange_weak_explicit(
|
|
u.atomic_address, &expected, val, std::memory_order_acq_rel, std::memory_order_relaxed);
|
|
return expected;
|
|
}
|
|
|
|
__device__ int atomicCAS(int* address, int compare, int val) {
|
|
return atomicCAS_impl(address, compare, val);
|
|
}
|
|
__device__ unsigned int atomicCAS(unsigned int* address, unsigned int compare, unsigned int val) {
|
|
return atomicCAS_impl(address, compare, val);
|
|
}
|
|
__device__ unsigned long long int atomicCAS(unsigned long long int* address,
|
|
unsigned long long int compare,
|
|
unsigned long long int val) {
|
|
return atomicCAS_impl(address, compare, val);
|
|
}
|
|
|
|
// atomicAnd()
|
|
__device__ int atomicAnd(int* address, int val) { return hc::atomic_fetch_and(address, val); }
|
|
__device__ unsigned int atomicAnd(unsigned int* address, unsigned int val) {
|
|
return hc::atomic_fetch_and(address, val);
|
|
}
|
|
__device__ unsigned long long int atomicAnd(unsigned long long int* address,
|
|
unsigned long long int val) {
|
|
return (long long int)hc::atomic_fetch_and((uint64_t*)address, (uint64_t)val);
|
|
}
|
|
|
|
// atomicOr()
|
|
__device__ int atomicOr(int* address, int val) { return hc::atomic_fetch_or(address, val); }
|
|
__device__ unsigned int atomicOr(unsigned int* address, unsigned int val) {
|
|
return hc::atomic_fetch_or(address, val);
|
|
}
|
|
__device__ unsigned long long int atomicOr(unsigned long long int* address,
|
|
unsigned long long int val) {
|
|
return (long long int)hc::atomic_fetch_or((uint64_t*)address, (uint64_t)val);
|
|
}
|
|
|
|
// atomicXor()
|
|
__device__ int atomicXor(int* address, int val) { return hc::atomic_fetch_xor(address, val); }
|
|
__device__ unsigned int atomicXor(unsigned int* address, unsigned int val) {
|
|
return hc::atomic_fetch_xor(address, val);
|
|
}
|
|
__device__ unsigned long long int atomicXor(unsigned long long int* address,
|
|
unsigned long long int val) {
|
|
return (long long int)hc::atomic_fetch_xor((uint64_t*)address, (uint64_t)val);
|
|
}
|
|
|
|
// atomicInc
|
|
__device__ unsigned int atomicInc(unsigned int* address, unsigned int val) {
|
|
return hc::__atomic_wrapinc(address, val);
|
|
}
|
|
|
|
// atomicDec
|
|
__device__ unsigned int atomicDec(unsigned int* address, unsigned int val) {
|
|
return hc::__atomic_wrapdec(address, val);
|
|
}
|
|
|
|
|
|
// warp vote function __all __any __ballot
|
|
__device__ int __all(int input) { return hc::__all(input); }
|
|
|
|
|
|
__device__ int __any(int input) {
|
|
#ifdef NVCC_COMPAT
|
|
if (hc::__any(input) != 0)
|
|
return 1;
|
|
else
|
|
return 0;
|
|
#else
|
|
return hc::__any(input);
|
|
#endif
|
|
}
|
|
|
|
__device__ unsigned long long int __ballot(int input) { return hc::__ballot(input); }
|
|
|
|
// warp shuffle functions
|
|
__device__ int __shfl(int input, int lane, int width) { return hc::__shfl(input, lane, width); }
|
|
|
|
__device__ int __shfl_up(int input, unsigned int lane_delta, int width) {
|
|
return hc::__shfl_up(input, lane_delta, width);
|
|
}
|
|
|
|
__device__ int __shfl_down(int input, unsigned int lane_delta, int width) {
|
|
return hc::__shfl_down(input, lane_delta, width);
|
|
}
|
|
|
|
__device__ int __shfl_xor(int input, int lane_mask, int width) {
|
|
return hc::__shfl_xor(input, lane_mask, width);
|
|
}
|
|
|
|
__device__ float __shfl(float input, int lane, int width) { return hc::__shfl(input, lane, width); }
|
|
|
|
__device__ float __shfl_up(float input, unsigned int lane_delta, int width) {
|
|
return hc::__shfl_up(input, lane_delta, width);
|
|
}
|
|
|
|
__device__ float __shfl_down(float input, unsigned int lane_delta, int width) {
|
|
return hc::__shfl_down(input, lane_delta, width);
|
|
}
|
|
|
|
__device__ float __shfl_xor(float input, int lane_mask, int width) {
|
|
return hc::__shfl_xor(input, lane_mask, width);
|
|
}
|
|
|
|
__host__ __device__ int min(int arg1, int arg2) {
|
|
return (int)(hc::precise_math::fmin((float)arg1, (float)arg2));
|
|
}
|
|
__host__ __device__ int max(int arg1, int arg2) {
|
|
return (int)(hc::precise_math::fmax((float)arg1, (float)arg2));
|
|
}
|
|
|
|
__device__ void* __get_dynamicgroupbaseptr() {
|
|
return hc::get_dynamic_group_segment_base_pointer();
|
|
}
|
|
|
|
__host__ void* __get_dynamicgroupbaseptr() { return nullptr; }
|
|
|
|
// Precise Math Functions
|
|
__device__ float __hip_precise_cosf(float x) { return hc::precise_math::cosf(x); }
|
|
|
|
__device__ float __hip_precise_exp10f(float x) { return hc::precise_math::exp10f(x); }
|
|
|
|
__device__ float __hip_precise_expf(float x) { return hc::precise_math::expf(x); }
|
|
|
|
__device__ float __hip_precise_frsqrt_rn(float x) { return hc::precise_math::rsqrt(x); }
|
|
|
|
__device__ float __hip_precise_fsqrt_rd(float x) { return hc::precise_math::sqrt(x); }
|
|
|
|
__device__ float __hip_precise_fsqrt_rn(float x) { return hc::precise_math::sqrt(x); }
|
|
|
|
__device__ float __hip_precise_fsqrt_ru(float x) { return hc::precise_math::sqrt(x); }
|
|
|
|
__device__ float __hip_precise_fsqrt_rz(float x) { return hc::precise_math::sqrt(x); }
|
|
|
|
__device__ float __hip_precise_log10f(float x) { return hc::precise_math::log10(x); }
|
|
|
|
__device__ float __hip_precise_log2f(float x) { return hc::precise_math::log2(x); }
|
|
|
|
__device__ float __hip_precise_logf(float x) { return hc::precise_math::logf(x); }
|
|
|
|
__device__ float __hip_precise_powf(float base, float exponent) {
|
|
return hc::precise_math::powf(base, exponent);
|
|
}
|
|
|
|
__device__ void __hip_precise_sincosf(float x, float* s, float* c) {
|
|
hc::precise_math::sincosf(x, s, c);
|
|
}
|
|
|
|
__device__ float __hip_precise_sinf(float x) { return hc::precise_math::sinf(x); }
|
|
|
|
__device__ float __hip_precise_tanf(float x) { return hc::precise_math::tanf(x); }
|
|
|
|
// Double Precision Math
|
|
__device__ double __hip_precise_dsqrt_rd(double x) { return hc::precise_math::sqrt(x); }
|
|
|
|
__device__ double __hip_precise_dsqrt_rn(double x) { return hc::precise_math::sqrt(x); }
|
|
|
|
__device__ double __hip_precise_dsqrt_ru(double x) { return hc::precise_math::sqrt(x); }
|
|
|
|
__device__ double __hip_precise_dsqrt_rz(double x) { return hc::precise_math::sqrt(x); }
|
|
|
|
#define LOG_BASE2_E 1.4426950408889634
|
|
#define LOG_BASE2_10 3.32192809488736
|
|
#define ONE_DIV_LOG_BASE2_E 0.69314718056
|
|
#define ONE_DIV_LOG_BASE2_10 0.30102999566
|
|
|
|
// Fast Math Intrinsics
|
|
__device__ float __hip_fast_exp10f(float x) { return __hip_fast_exp2f(x * LOG_BASE2_E); }
|
|
|
|
__device__ float __hip_fast_expf(float x) { return __hip_fast_exp2f(x * LOG_BASE2_10); }
|
|
|
|
__device__ float __hip_fast_frsqrt_rn(float x) {
|
|
return 1 / __hip_fast_fsqrt_rd(x);
|
|
;
|
|
}
|
|
|
|
__device__ float __hip_fast_fsqrt_rn(float x) { return __hip_fast_fsqrt_rd(x); }
|
|
|
|
__device__ float __hip_fast_fsqrt_ru(float x) { return __hip_fast_fsqrt_rd(x); }
|
|
|
|
__device__ float __hip_fast_fsqrt_rz(float x) { return __hip_fast_fsqrt_rd(x); }
|
|
|
|
__device__ float __hip_fast_log10f(float x) { return ONE_DIV_LOG_BASE2_E * __hip_fast_log2f(x); }
|
|
|
|
__device__ float __hip_fast_logf(float x) { return ONE_DIV_LOG_BASE2_10 * __hip_fast_log2f(x); }
|
|
|
|
__device__ float __hip_fast_powf(float base, float exponent) {
|
|
return hc::fast_math::powf(base, exponent);
|
|
}
|
|
|
|
__device__ void __hip_fast_sincosf(float x, float* s, float* c) {
|
|
*s = __hip_fast_sinf(x);
|
|
*c = __hip_fast_cosf(x);
|
|
}
|
|
|
|
__device__ float __hip_fast_tanf(float x) { return hc::fast_math::tanf(x); }
|
|
|
|
// Double Precision Math
|
|
// FIXME - HCC doesn't have a fast_math version double FP sqrt
|
|
// Another issue is that these intrinsics call for a specific rounding mode;
|
|
// however, their implementation all map to the same sqrt builtin
|
|
__device__ double __hip_fast_dsqrt_rd(double x) { return hc::precise_math::sqrt(x); }
|
|
|
|
__device__ double __hip_fast_dsqrt_rn(double x) { return hc::precise_math::sqrt(x); }
|
|
|
|
__device__ double __hip_fast_dsqrt_ru(double x) { return hc::precise_math::sqrt(x); }
|
|
|
|
__device__ double __hip_fast_dsqrt_rz(double x) { return hc::precise_math::sqrt(x); }
|
|
|
|
__device__ void __threadfence_system(void) { std::atomic_thread_fence(std::memory_order_seq_cst); }
|
|
|
|
double __hip_host_j0(double x) {
|
|
double ret, a = std::fabs(x);
|
|
if (a < 8.0) {
|
|
double y = x * x;
|
|
double y1 = __hip_j0a6 * y + __hip_j0a5;
|
|
double z1 = 1.0 * y + __hip_j0b5;
|
|
|
|
double y2 = y1 * y + __hip_j0a4;
|
|
double z2 = z1 * y + __hip_j0b4;
|
|
|
|
double y3 = y2 * y + __hip_j0a3;
|
|
double z3 = z2 * y + __hip_j0b3;
|
|
|
|
double y4 = y3 * y + __hip_j0a2;
|
|
double z4 = z3 * y + __hip_j0b2;
|
|
|
|
double y5 = y4 * y + __hip_j0a1;
|
|
double z5 = z4 * y + __hip_j0b1;
|
|
|
|
ret = y5 / z5;
|
|
|
|
} else {
|
|
double z = 8.0 / a;
|
|
double y = z * z;
|
|
double x1 = a - __hip_j0c;
|
|
|
|
double y1 = __hip_j0c4 * y + __hip_j0c3;
|
|
double z1 = __hip_j0d5 * y + __hip_j0d4;
|
|
|
|
double y2 = y1 * y + __hip_j0c2;
|
|
double z2 = z1 * z + __hip_j0d3;
|
|
|
|
double y3 = y2 * y + __hip_j0c1;
|
|
double z3 = z2 * y + __hip_j0d2;
|
|
|
|
double y4 = y3 * y + 1.0;
|
|
double z4 = z3 * y + __hip_j0d1;
|
|
|
|
ret = std::sqrt(__hip_j0e / a) * (std::cos(x1) * y4 - z * std::sin(x1) * z4);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
float __hip_host_j0f(float x) {
|
|
float ret, a = fabs(x);
|
|
if (a < 8.0) {
|
|
float y = x * x;
|
|
float y1 = __hip_j0a6 * y + __hip_j0a5;
|
|
float z1 = 1.0 * y + __hip_j0b5;
|
|
|
|
float y2 = y1 * y + __hip_j0a4;
|
|
float z2 = z1 * y + __hip_j0b4;
|
|
|
|
float y3 = y2 * y + __hip_j0a3;
|
|
float z3 = z2 * y + __hip_j0b3;
|
|
|
|
float y4 = y3 * y + __hip_j0a2;
|
|
float z4 = z3 * y + __hip_j0b2;
|
|
|
|
float y5 = y4 * y + __hip_j0a1;
|
|
float z5 = z4 * y + __hip_j0b1;
|
|
|
|
ret = y5 / z5;
|
|
|
|
} else {
|
|
float z = 8.0 / a;
|
|
float y = z * z;
|
|
float x1 = a - __hip_j0c;
|
|
|
|
float y1 = __hip_j0c4 * y + __hip_j0c3;
|
|
float z1 = __hip_j0d5 * y + __hip_j0d4;
|
|
|
|
float y2 = y1 * y + __hip_j0c2;
|
|
float z2 = z1 * z + __hip_j0d3;
|
|
|
|
float y3 = y2 * y + __hip_j0c1;
|
|
float z3 = z2 * y + __hip_j0d2;
|
|
|
|
float y4 = y3 * y + 1.0;
|
|
float z4 = z3 * y + __hip_j0d1;
|
|
|
|
ret = std::sqrt(__hip_j0e / a) * (std::cos(x1) * y4 - z * std::sin(x1) * z4);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
double __hip_host_j1(double x) {
|
|
double ret, a = std::fabs(x);
|
|
if (a < 8.0) {
|
|
double y = x * x;
|
|
|
|
double y1 = __hip_j1a1 * y + __hip_j1a2;
|
|
double z1 = 1.0 * y + __hip_j1b1;
|
|
|
|
double y2 = y1 * y + __hip_j1a3;
|
|
double z2 = z1 * y + __hip_j1b2;
|
|
|
|
double y3 = y2 * y + __hip_j1a4;
|
|
double z3 = z2 * y + __hip_j1b3;
|
|
|
|
double y4 = y3 * y + __hip_j1a5;
|
|
double z4 = z3 * y + __hip_j1b4;
|
|
|
|
double y5 = y4 * y + __hip_j1a6;
|
|
double z5 = z4 * y + __hip_j1b5;
|
|
|
|
ret = x * y5 / z5;
|
|
|
|
} else {
|
|
double z = 8.0 / a;
|
|
double y = z * z;
|
|
double x1 = a - __hip_j1c;
|
|
|
|
double y1 = __hip_j1c1 * y + __hip_j1c2;
|
|
double y2 = y1 * y + __hip_j1c3;
|
|
double y3 = y2 * y + __hip_j1c4;
|
|
double y4 = y3 * y + 1.0;
|
|
|
|
double z1 = __hip_j1d1 * y + __hip_j1d2;
|
|
double z2 = z1 * y + __hip_j1d3;
|
|
double z3 = z2 * y + __hip_j1d4;
|
|
double z4 = z3 * y + __hip_j1d5;
|
|
|
|
ret = std::sqrt(__hip_j1e / a) * (std::cos(x1) * y4 - z * std::sin(x1) * z4);
|
|
if (x < 0.0) ret = -ret;
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
float __hip_host_j1f(float x) {
|
|
double ret, a = fabs(x);
|
|
if (a < 8.0) {
|
|
float y = x * x;
|
|
|
|
float y1 = __hip_j1a1 * y + __hip_j1a2;
|
|
float z1 = 1.0 * y + __hip_j1b1;
|
|
|
|
float y2 = y1 * y + __hip_j1a3;
|
|
float z2 = z1 * y + __hip_j1b2;
|
|
|
|
float y3 = y2 * y + __hip_j1a4;
|
|
float z3 = z2 * y + __hip_j1b3;
|
|
|
|
float y4 = y3 * y + __hip_j1a5;
|
|
float z4 = z3 * y + __hip_j1b4;
|
|
|
|
float y5 = y4 * y + __hip_j1a6;
|
|
float z5 = z4 * y + __hip_j1b5;
|
|
|
|
ret = x * y5 / z5;
|
|
|
|
} else {
|
|
float z = 8.0 / a;
|
|
float y = z * z;
|
|
float x1 = a - __hip_j1c;
|
|
|
|
float y1 = __hip_j1c1 * y + __hip_j1c2;
|
|
float y2 = y1 * y + __hip_j1c3;
|
|
float y3 = y2 * y + __hip_j1c4;
|
|
float y4 = y3 * y + 1.0;
|
|
|
|
float z1 = __hip_j1d1 * y + __hip_j1d2;
|
|
float z2 = z1 * y + __hip_j1d3;
|
|
float z3 = z2 * y + __hip_j1d4;
|
|
float z4 = z3 * y + __hip_j1d5;
|
|
|
|
ret = std::sqrt(__hip_j1e / a) * (std::cos(x1) * y4 - z * std::sin(x1) * z4);
|
|
if (x < 0.0) ret = -ret;
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
double __hip_host_y0(double x) {
|
|
double ret;
|
|
|
|
if (x < 8.0) {
|
|
double y = x * x;
|
|
double y1 = __hip_y0a1 * y + __hip_y0a2;
|
|
double y2 = y1 * y + __hip_y0a3;
|
|
double y3 = y2 * y + __hip_y0a4;
|
|
double y4 = y3 * y + __hip_y0a5;
|
|
double y5 = y4 * y + __hip_y0a6;
|
|
|
|
double z1 = 1.0 * y + __hip_y0b1;
|
|
double z2 = z1 * y + __hip_y0b2;
|
|
double z3 = z2 * y + __hip_y0b3;
|
|
double z4 = z3 * y + __hip_y0b4;
|
|
double z5 = z4 * y + __hip_y0b5;
|
|
|
|
|
|
ret = (y5 / z5) + __hip_y0c * __hip_host_j0(x) * std::log(x);
|
|
} else {
|
|
double z = 8.0 / x;
|
|
double y = z * z;
|
|
double x1 = x - __hip_y0d;
|
|
|
|
double y1 = __hip_y0e1 * y + __hip_y0e2;
|
|
double y2 = y1 * y + __hip_y0e3;
|
|
double y3 = y2 * y + __hip_y0e4;
|
|
double y4 = y3 * y + 1.0;
|
|
|
|
double z1 = __hip_y0f1 * y + __hip_y0f2;
|
|
double z2 = z1 * y + __hip_y0f3;
|
|
double z3 = z2 * y + __hip_y0f4;
|
|
double z4 = z3 * y + __hip_y0f5;
|
|
|
|
ret = std::sqrt(__hip_y1g / x) * (std::sin(x1) * y4 + z * std::cos(x1) * z4);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
float __hip_host_y0f(float x) {
|
|
float ret;
|
|
|
|
if (x < 8.0) {
|
|
float y = x * x;
|
|
float y1 = __hip_y0a1 * y + __hip_y0a2;
|
|
float y2 = y1 * y + __hip_y0a3;
|
|
float y3 = y2 * y + __hip_y0a4;
|
|
float y4 = y3 * y + __hip_y0a5;
|
|
float y5 = y4 * y + __hip_y0a6;
|
|
|
|
float z1 = 1.0 * y + __hip_y0b1;
|
|
float z2 = z1 * y + __hip_y0b2;
|
|
float z3 = z2 * y + __hip_y0b3;
|
|
float z4 = z3 * y + __hip_y0b4;
|
|
float z5 = z4 * y + __hip_y0b5;
|
|
|
|
|
|
ret = (y5 / z5) + __hip_y0c * __hip_host_j0f(x) * log(x);
|
|
} else {
|
|
float z = 8.0 / x;
|
|
float y = z * z;
|
|
float x1 = x - __hip_y0d;
|
|
|
|
float y1 = __hip_y0e1 * y + __hip_y0e2;
|
|
float y2 = y1 * y + __hip_y0e3;
|
|
float y3 = y2 * y + __hip_y0e4;
|
|
float y4 = y3 * y + 1.0;
|
|
|
|
float z1 = __hip_y0f1 * y + __hip_y0f2;
|
|
float z2 = z1 * y + __hip_y0f3;
|
|
float z3 = z2 * y + __hip_y0f4;
|
|
float z4 = z3 * y + __hip_y0f5;
|
|
|
|
ret = std::sqrt(__hip_y1g / x) * (std::sin(x1) * y4 + z * std::cos(x1) * z4);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
double __hip_host_y1(double x) {
|
|
double ret;
|
|
|
|
if (x < 8.0) {
|
|
double y = x * x;
|
|
|
|
double y1 = __hip_y1a1 * y + __hip_y1a2;
|
|
double y2 = y1 * y + __hip_y1a3;
|
|
double y3 = y2 * y + __hip_y1a4;
|
|
double y4 = y3 * y + __hip_y1a5;
|
|
double y5 = y4 * y + __hip_y1a6;
|
|
double y6 = x * y5;
|
|
|
|
double z1 = __hip_y1b1 + y;
|
|
double z2 = z1 * y + __hip_y1b2;
|
|
double z3 = z2 * y + __hip_y1b3;
|
|
double z4 = z3 * y + __hip_y1b4;
|
|
double z5 = z4 * y + __hip_y1b5;
|
|
double z6 = z5 * y + __hip_y1b6;
|
|
|
|
ret = (y6 / z6) + __hip_y1c * (__hip_host_j1(x) * std::log(x) - 1.0 / x);
|
|
} else {
|
|
double z = 8.0 / x;
|
|
double y = z * z;
|
|
double x1 = x - __hip_y1d;
|
|
|
|
double y1 = __hip_y1e1 * y + __hip_y1e2;
|
|
double y2 = y1 * y + __hip_y1e3;
|
|
double y3 = y2 * y + __hip_y1e4;
|
|
double y4 = y3 * y + 1.0;
|
|
|
|
double z1 = __hip_y1f1 * y + __hip_y1f2;
|
|
double z2 = z1 * y + __hip_y1f3;
|
|
double z3 = z2 * y + __hip_y1f4;
|
|
double z4 = z3 * y + __hip_y1f5;
|
|
|
|
ret = std::sqrt(__hip_y1g / x) * (std::sin(x1) * y4 + z * std::cos(x1) * z4);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
float __hip_host_y1f(float x) {
|
|
float ret;
|
|
|
|
if (x < 8.0) {
|
|
float y = x * x;
|
|
|
|
float y1 = __hip_y1a1 * y + __hip_y1a2;
|
|
float y2 = y1 * y + __hip_y1a3;
|
|
float y3 = y2 * y + __hip_y1a4;
|
|
float y4 = y3 * y + __hip_y1a5;
|
|
float y5 = y4 * y + __hip_y1a6;
|
|
float y6 = x * y5;
|
|
|
|
float z1 = __hip_y1b1 + y;
|
|
float z2 = z1 * y + __hip_y1b2;
|
|
float z3 = z2 * y + __hip_y1b3;
|
|
float z4 = z3 * y + __hip_y1b4;
|
|
float z5 = z4 * y + __hip_y1b5;
|
|
float z6 = z5 * y + __hip_y1b6;
|
|
|
|
ret = (y6 / z6) + __hip_y1c * (__hip_host_j1f(x) * log(x) - 1.0 / x);
|
|
} else {
|
|
float z = 8.0 / x;
|
|
float y = z * z;
|
|
float x1 = x - __hip_y1d;
|
|
|
|
float y1 = __hip_y1e1 * y + __hip_y1e2;
|
|
float y2 = y1 * y + __hip_y1e3;
|
|
float y3 = y2 * y + __hip_y1e4;
|
|
float y4 = y3 * y + 1.0;
|
|
|
|
float z1 = __hip_y1f1 * y + __hip_y1f2;
|
|
float z2 = z1 * y + __hip_y1f3;
|
|
float z3 = z2 * y + __hip_y1f4;
|
|
float z4 = z3 * y + __hip_y1f5;
|
|
|
|
ret = std::sqrt(__hip_y1g / x) * (std::sin(x1) * y4 + z * std::cos(x1) * z4);
|
|
}
|
|
return ret;
|
|
}
|
|
|
|
double __hip_host_jn(int n, double x) {
|
|
int sum = 0, m;
|
|
double a, b0, b1, b2, val, t, ret;
|
|
if (n < 0) {
|
|
return NAN;
|
|
}
|
|
a = std::fabs(x);
|
|
if (n == 0) {
|
|
return (__hip_host_j0(a));
|
|
}
|
|
if (n == 1) {
|
|
return (__hip_host_j1(a));
|
|
}
|
|
if (a == 0.0) {
|
|
return 0.0;
|
|
} else if (a > (double)n) {
|
|
t = 2.0 / a;
|
|
b1 = __hip_host_j0(a);
|
|
b0 = __hip_host_j1(a);
|
|
for (int i = 1; i < n; i++) {
|
|
b2 = i * t * b0 - b1;
|
|
b1 = b0;
|
|
b0 = b2;
|
|
}
|
|
ret = b0;
|
|
} else {
|
|
t = 2.0 / a;
|
|
m = 2 * ((n + (int)std::sqrt(__hip_k1 * n)) / 2);
|
|
b2 = ret = val = 0.0;
|
|
b0 = 1.0;
|
|
for (int i = m; i > 0; i--) {
|
|
b1 = i * t * b0 - b2;
|
|
b2 = b0;
|
|
b0 = b1;
|
|
if (std::fabs(b0) > __hip_k2) {
|
|
b0 *= __hip_k3;
|
|
b2 *= __hip_k3;
|
|
ret *= __hip_k3;
|
|
val *= __hip_k3;
|
|
}
|
|
if (sum) val += b0;
|
|
sum = !sum;
|
|
if (i == n) ret = b2;
|
|
}
|
|
val = 2.0 * val - b0;
|
|
ret /= val;
|
|
}
|
|
return x < 0.0 && n % 2 == 1 ? -ret : ret;
|
|
}
|
|
|
|
float __hip_host_jnf(int n, float x) {
|
|
int sum = 0, m;
|
|
float a, b0, b1, b2, val, t, ret;
|
|
if (n < 0) {
|
|
return NAN;
|
|
}
|
|
a = fabs(x);
|
|
if (n == 0) {
|
|
return (__hip_host_j0f(a));
|
|
}
|
|
if (n == 1) {
|
|
return (__hip_host_j1f(a));
|
|
}
|
|
if (a == 0.0) {
|
|
return 0.0;
|
|
} else if (a > (float)n) {
|
|
t = 2.0 / a;
|
|
b1 = __hip_host_j0f(a);
|
|
b0 = __hip_host_j1f(a);
|
|
for (int i = 1; i < n; i++) {
|
|
b2 = i * t * b0 - b1;
|
|
b1 = b0;
|
|
b0 = b2;
|
|
}
|
|
ret = b0;
|
|
} else {
|
|
t = 2.0 / a;
|
|
m = 2 * ((n + (int)std::sqrt(__hip_k1 * n)) / 2);
|
|
b2 = ret = val = 0.0;
|
|
b0 = 1.0;
|
|
for (int i = m; i > 0; i--) {
|
|
b1 = i * t * b0 - b2;
|
|
b2 = b0;
|
|
b0 = b1;
|
|
if (std::fabs(b0) > __hip_k2) {
|
|
b0 *= __hip_k3;
|
|
b2 *= __hip_k3;
|
|
ret *= __hip_k3;
|
|
val *= __hip_k3;
|
|
}
|
|
if (sum) val += b0;
|
|
sum = !sum;
|
|
if (i == n) ret = b2;
|
|
}
|
|
val = 2.0 * val - b0;
|
|
ret /= val;
|
|
}
|
|
return x < 0.0 && n % 2 == 1 ? -ret : ret;
|
|
}
|
|
|
|
double __hip_host_yn(int n, double x) {
|
|
double b0, b1, b2, t;
|
|
|
|
if (n < 0 || x == 0.0) {
|
|
return NAN;
|
|
}
|
|
if (n == 0) {
|
|
return (__hip_host_y0(x));
|
|
}
|
|
if (n == 1) {
|
|
return (__hip_host_y1(x));
|
|
}
|
|
t = 2.0 / x;
|
|
b0 = __hip_host_y1(x);
|
|
b1 = __hip_host_y0(x);
|
|
for (int i = 1; i < n; i++) {
|
|
b2 = i * t * b0 - b1;
|
|
b1 = b0;
|
|
b0 = b2;
|
|
}
|
|
return b0;
|
|
}
|
|
|
|
float __hip_host_ynf(int n, float x) {
|
|
float b0, b1, b2, t;
|
|
|
|
if (n < 0 || x == 0.0) {
|
|
return NAN;
|
|
}
|
|
if (n == 0) {
|
|
return (__hip_host_y0f(x));
|
|
}
|
|
if (n == 1) {
|
|
return (__hip_host_y1f(x));
|
|
}
|
|
t = 2.0 / x;
|
|
b0 = __hip_host_y1f(x);
|
|
b1 = __hip_host_y0f(x);
|
|
for (int i = 1; i < n; i++) {
|
|
b2 = i * t * b0 - b1;
|
|
b1 = b0;
|
|
b0 = b2;
|
|
}
|
|
return b0;
|
|
}
|