Merged
float类型性能优化 #7121
xuyiyangAdd on 8 days ago
float类型性能优化 #7121
Merged
Total: 1 files changed+49-53
| @@ -916,57 +916,50 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float atan2f(float y, float x) | |||
| 916 | // the first octant, so the magnitudes are reduced first and the octant is restored afterwards. | 916 | // the first octant, so the magnitudes are reduced first and the octant is restored afterwards. |
| 917 | float ay = __fabsf(y); | 917 | float ay = __fabsf(y); |
| 918 | float ax = __fabsf(x); | 918 | float ax = __fabsf(x); |
| 919 | - | ||
| 920 | - // Order the magnitudes so the ratio lo / hi never exceeds 1, which is the polynomial's valid range. | ||
| 921 | bool y_gt_x = ay > ax; | 919 | bool y_gt_x = ay > ax; |
| 922 | float hi = y_gt_x ? ay : ax; | 920 | float hi = y_gt_x ? ay : ax; |
| 923 | float lo = y_gt_x ? ax : ay; | 921 | float lo = y_gt_x ? ax : ay; |
| 922 | + bool nz = hi != 0.0f; | ||
| 924 | 923 | ||
| 925 | - float a = 0.0f; | 924 | + // Safe denominator: 1 instead of 0 keeps the shared path free of the 0/0 NaN; the garbage |
| 926 | - if (hi != 0.0f) { | 925 | + // result on that route is discarded by the tail select. |
| 927 | - // The ratio lies in [0, 1]. atan is odd, so with z = r^2 the tail is even in r: | 926 | + float safe_hi = nz ? hi : 1.0f; |
| 928 | - // atan(r) = r + r^3 * P(z), P(z) = -0.33333197 + 0.19993925 z - 0.14207722 z^2 + 0.10640416 z^3 | 927 | + float r = lo / safe_hi; |
| 929 | - // - 0.074792981 z^4 + 0.042200752 z^5 - 0.015681878 z^6 + 0.0027380611 z^7 | 928 | + float z = r * r; |
| 930 | - // P is minimax-fitted, not a Taylor truncation (leading term -0.33333197, not -1/3): the error is | ||
| 931 | - // flattened across all of [0, 1], so the coefficients cannot be extended by a series rule. | ||
| 932 | - float r = lo / hi; | ||
| 933 | - float z = r * r; | ||
| 934 | - float p = 0.0027380611281841993332f; | ||
| 935 | - p = __fma(z, p, -0.015681877732276916504f); | ||
| 936 | - p = __fma(z, p, 0.042200751602649688721f); | ||
| 937 | - p = __fma(z, p, -0.074792981147766113281f); | ||
| 938 | - p = __fma(z, p, 0.10640415549278259277f); | ||
| 939 | - p = __fma(z, p, -0.14207722246646881104f); | ||
| 940 | - p = __fma(z, p, 0.19993925094604492188f); | ||
| 941 | - p = __fma(z, p, -0.33333197236061096191f); | ||
| 942 | - // Assemble r + r^3 * P(z). r is the dominant term, so passing it as the fma addend keeps it exact | ||
| 943 | - // and confines rounding to the correction. | ||
| 944 | - float t = __fma(z * p, r, r); | ||
| 945 | 929 | ||
| 946 | - // Rebuild the full-circle angle from the octant. |y| == |x| is special-cased so the diagonals come | 930 | + // atan(r) = r + r^3 * P(r^2), P minimax-fitted on [0, 1]. |
| 947 | - // out as the exactly rounded pi/4 and 3pi/4, and because inf/inf above makes r NaN. | 931 | + float p = 0.0027380611281841993332f; |
| 948 | - if (ay == ax) { | 932 | + p = __fma(z, p, -0.015681877732276916504f); |
| 949 | - a = signbit(x) ? 2.35619449615478515625f : 0.78539818525314331055f; | 933 | + p = __fma(z, p, 0.042200751602649688721f); |
| 950 | - } else if (y_gt_x) { | 934 | + p = __fma(z, p, -0.074792981147766113281f); |
| 951 | - // Steeper than the diagonal: the ratio r = ax/ay is the reciprocal of y/x, so | 935 | + p = __fma(z, p, 0.10640415549278259277f); |
| 952 | - // atan(y/x) = pi/2 - atan(x/y) = pi/2 - t. A negative x lies in the second quadrant, so pi/2 + t. | 936 | + p = __fma(z, p, -0.14207722246646881104f); |
| 953 | - float pio2 = 1.57079637050628662109f; | 937 | + p = __fma(z, p, 0.19993925094604492188f); |
| 954 | - a = signbit(x) ? (pio2 + t) : (pio2 - t); | 938 | + p = __fma(z, p, -0.33333197236061096191f); |
| 955 | - } else { | 939 | + float t = __fma(z * p, r, r); |
| 956 | - // Shallow angles need no reflection; a negative x mirrors the first octant into the second. | 940 | + |
| 957 | - a = signbit(x) ? (3.14159274101257324219f - t) : t; | 941 | + // The only NaN source left is inf/inf; pi/4 makes the affine formula below produce every |
| 958 | - } | 942 | + // diagonal answer (0 + pi/4 = pi/4, pi - pi/4 = 3pi/4) with no special case. |
| 959 | - } else { | 943 | + t = (t != t) ? 0.78539818525314331055f : t; |
| 960 | - // Both arguments are zero. Sign of x picks pi or 0, matching atan2(+-0, -0) == +-pi. | 944 | + |
| 961 | - a = signbit(x) ? 3.14159274101257324219f : 0.0f; | 945 | + // Quadrant restore: every octant bucket is base + sgn*t, so the three-way select of the |
| 962 | - } | 946 | + // reference implementation is two bucket picks plus a single fma. |
| 963 | - float result = signbit(y) ? -a : a; | 947 | + bool xneg = __float_as_uint(x) >> 31; // 31: sign-bit offset of fp32 |
| 964 | - // Summing the magnitudes tests both arguments for NaN at once and propagates it as the result. | 948 | + float base = y_gt_x ? 1.57079637050628662109f : (xneg ? 3.14159274101257324219f : 0.0f); |
| 949 | + float sgn = (y_gt_x != xneg) ? -1.0f : 1.0f; | ||
| 950 | + float a_main = __fma(sgn, t, base); | ||
| 951 | + | ||
| 952 | + // hi == 0: sign of x picks pi or 0 (atan2(+-0, -0) == +-pi); nonzero lo is the y-axis (pi/2). | ||
| 953 | + float a0 = xneg ? 3.14159274101257324219f : 0.0f; | ||
| 954 | + a0 = ((hi == 0.0f) && (lo != 0.0f)) ? 1.57079637050628662109f : a0; | ||
| 955 | + float a = nz ? a_main : a0; | ||
| 956 | + | ||
| 957 | + // Odd in y: apply the sign (including -0.0) bit-wise, then let a NaN argument win. | ||
| 958 | + float result = __float_as_uint(y) >> 31 ? -a : a; // 31: sign-bit offset of fp32 | ||
| 965 | float sum = ax + ay; | 959 | float sum = ax + ay; |
| 966 | if (__isnan(sum)) { | 960 | if (__isnan(sum)) { |
| 967 | result = sum; | 961 | result = sum; |
| 968 | } | 962 | } |
| 969 | - // atan2 is odd in y, so the sign of y (including -0.0) is reapplied last. | ||
| 970 | return result; | 963 | return result; |
| 971 | } | 964 | } |
| 972 | 965 | ||
| @@ -4299,12 +4292,6 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float rnorm4df(float a, float b, float c, | |||
| 4299 | const float abs_b = fabsf(b); | 4292 | const float abs_b = fabsf(b); |
| 4300 | const float abs_c = fabsf(c); | 4293 | const float abs_c = fabsf(c); |
| 4301 | const float abs_d = fabsf(d); | 4294 | const float abs_d = fabsf(d); |
| 4302 | - if (isinf(abs_a) || isinf(abs_b) || isinf(abs_c) || isinf(abs_d)) { | ||
| 4303 | - return 0.0f; | ||
| 4304 | - } | ||
| 4305 | - if (isnan(abs_a) || isnan(abs_b) || isnan(abs_c) || isnan(abs_d)) { | ||
| 4306 | - return ASCRT_NAN_F; | ||
| 4307 | - } | ||
| 4308 | 4295 | ||
| 4309 | const float max_ab = (abs_a > abs_b) ? abs_a : abs_b; | 4296 | const float max_ab = (abs_a > abs_b) ? abs_a : abs_b; |
| 4310 | const float max_cd = (abs_c > abs_d) ? abs_c : abs_d; | 4297 | const float max_cd = (abs_c > abs_d) ? abs_c : abs_d; |
| @@ -4319,16 +4306,25 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float rnorm4df(float a, float b, float c, | |||
| 4319 | const float t_c = abs_c * scale; | 4306 | const float t_c = abs_c * scale; |
| 4320 | const float t_d = abs_d * scale; | 4307 | const float t_d = abs_d * scale; |
| 4321 | float square_sum = fmaf(t_a, t_a, fmaf(t_b, t_b, fmaf(t_c, t_c, t_d * t_d))); | 4308 | float square_sum = fmaf(t_a, t_a, fmaf(t_b, t_b, fmaf(t_c, t_c, t_d * t_d))); |
| 4322 | - // Zero norm means the reciprocal norm is infinite. | ||
| 4323 | - if (square_sum == 0.0f) { | ||
| 4324 | - return ASCRT_INF_F; | ||
| 4325 | - } | ||
| 4326 | 4309 | ||
| 4327 | // Compute an approximate reciprocal norm and refine it with one Newton step. | 4310 | // Compute an approximate reciprocal norm and refine it with one Newton step. |
| 4328 | float inv_norm = 1.0f / __sqrtf(square_sum); | 4311 | float inv_norm = 1.0f / __sqrtf(square_sum); |
| 4329 | inv_norm = inv_norm * fmaf(-0.5f * square_sum, inv_norm * inv_norm, 1.5f); | 4312 | inv_norm = inv_norm * fmaf(-0.5f * square_sum, inv_norm * inv_norm, 1.5f); |
| 4330 | // Restore the outer scale factor. | 4313 | // Restore the outer scale factor. |
| 4331 | - return scale * inv_norm; | 4314 | + float result = scale * inv_norm; |
| 4315 | + | ||
| 4316 | + // Zero norm means the reciprocal norm is infinite. | ||
| 4317 | + if (square_sum == 0.0f) { | ||
| 4318 | + result = ASCRT_INF_F; | ||
| 4319 | + } | ||
| 4320 | + if (isnan(abs_a) || isnan(abs_b) || isnan(abs_c) || isnan(abs_d)) { | ||
| 4321 | + result = ASCRT_NAN_F; | ||
| 4322 | + } | ||
| 4323 | + if (isinf(abs_a) || isinf(abs_b) || isinf(abs_c) || isinf(abs_d)) { | ||
| 4324 | + result = 0.0f; | ||
| 4325 | + } | ||
| 4326 | + | ||
| 4327 | + return result; | ||
| 4332 | } | 4328 | } |
| 4333 | 4329 | ||
| 4334 | /** | 4330 | /** |