Merged
float类型数学接口性能优化 #6500
xuyiyangAdd on Sep 29
float类型数学接口性能优化 #6500
Merged
Total: 1 files changed+55-52
| @@ -3708,7 +3708,7 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ static inline float __asc_make_float_from_mantiss | |||
| 3708 | } else if (result_exponent < -149) { // -149: exponent of the smallest fp32 subnormal, below it underflows to 0 | 3708 | } else if (result_exponent < -149) { // -149: exponent of the smallest fp32 subnormal, below it underflows to 0 |
| 3709 | result = 0.0f; | 3709 | result = 0.0f; |
| 3710 | } else if (result_exponent >= -126) { // -126: smallest normal fp32 exponent | 3710 | } else if (result_exponent >= -126) { // -126: smallest normal fp32 exponent |
| 3711 | - // Normal result: normalize mantissa and pack signless fp32 bits. | 3711 | + // Normal result: normalize mantissa and pack signless fp32 bits. |
| 3712 | const uint32_t normalized_mantissa = mantissa << static_cast<uint32_t>(23 - mantissa_log2); | 3712 | const uint32_t normalized_mantissa = mantissa << static_cast<uint32_t>(23 - mantissa_log2); |
| 3713 | const uint32_t bits = | 3713 | const uint32_t bits = |
| 3714 | (static_cast<uint32_t>(result_exponent + 127) << 23U) | (normalized_mantissa & ASCRT_MAN_BIT_FLOAT_U); | 3714 | (static_cast<uint32_t>(result_exponent + 127) << 23U) | (normalized_mantissa & ASCRT_MAN_BIT_FLOAT_U); |
| @@ -4958,18 +4958,18 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float log1pf(float x) | |||
| 4958 | // Add the separate exponent contribution back as ln(2) * exponent. | 4958 | // Add the separate exponent contribution back as ln(2) * exponent. |
| 4959 | output = fmaf(exponent, 0.69314718246459960938f, output); | 4959 | output = fmaf(exponent, 0.69314718246459960938f, output); |
| 4960 | 4960 | ||
| 4961 | - // Positive infinity/NaN and all negative inputs enter the tail correction. | 4961 | + // Compute the candidate unconditionally; it is only meaningful for special |
| 4962 | - if (x_bits >= fp32_positive_inf_bits) { | 4962 | + // inputs and discarded by the select on the normal path. |
| 4963 | - // x is +inf/NaN or x < -1: generate the result via x * +inf + +inf. | 4963 | + const float inf_path = fmaf(x, ASCRT_INF_F, ASCRT_INF_F); |
| 4964 | - // x == -1 keeps the -inf produced by the reduction and polynomial path above. | 4964 | + |
| 4965 | - if (!(x_bits >= fp32_sign_bit && x_bits < log1p_lower_bound_bits)) { | 4965 | + // Merge the outer and inner conditions: the inf path is taken when both hold. |
| 4966 | - output = fmaf(x, ASCRT_INF_F, ASCRT_INF_F); | 4966 | + const bool need_inf = |
| 4967 | - } | 4967 | + (x_bits >= fp32_positive_inf_bits) && !(x_bits >= fp32_sign_bit && x_bits < log1p_lower_bound_bits); |
| 4968 | - // Inputs comparing equal to 0 select -0 so the sign of zero is preserved. | 4968 | + // Only x == -0 matches here (the outer condition already excludes +0). |
| 4969 | - if (x == 0.0f) { | 4969 | + const bool need_neg_zero = (x_bits >= fp32_positive_inf_bits) && (x == 0.0f); |
| 4970 | - output = ASCRT_NEG_ZERO_F; | 4970 | + |
| 4971 | - } | 4971 | + output = need_inf ? inf_path : output; |
| 4972 | - } | 4972 | + output = need_neg_zero ? ASCRT_NEG_ZERO_F : output; |
| 4973 | return output; | 4973 | return output; |
| 4974 | } | 4974 | } |
| 4975 | 4975 | ||
| @@ -5532,15 +5532,6 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float erfcf(float x) | |||
| 5532 | 5532 | ||
| 5533 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float exp10f(float x) | 5533 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float exp10f(float x) |
| 5534 | { | 5534 | { |
| 5535 | - // 10^x = 2^(x * log2(10)). Split x*log2(10) into an integer k and a small remainder r, then compute | ||
| 5536 | - // 2^k * 2^r with 2^k constructed exactly from the bit pattern. | ||
| 5537 | - if (__isnan(x)) { | ||
| 5538 | - return x; | ||
| 5539 | - } | ||
| 5540 | - if (__isinf(x)) { | ||
| 5541 | - return x > 0.0f ? x : 0.0f; | ||
| 5542 | - } | ||
| 5543 | - | ||
| 5544 | // t = x * (log2(10) / 252) + 0.5, a scaled-and-biased form of x * log2(10) used for floor reduction. | 5535 | // t = x * (log2(10) / 252) + 0.5, a scaled-and-biased form of x * log2(10) used for floor reduction. |
| 5545 | // log2(10) / 252 is precomputed so that floor(t * 252) directly yields the integer index k in the | 5536 | // log2(10) / 252 is precomputed so that floor(t * 252) directly yields the integer index k in the |
| 5546 | // range [0, 252] (covering the full float32 overflow/underflow domain of 10^x). | 5537 | // range [0, 252] (covering the full float32 overflow/underflow domain of 10^x). |
| @@ -5563,25 +5554,22 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float exp10f(float x) | |||
| 5563 | // drops from ~1e-8 to ~1e-15. r is confined to [-0.5, 0.5], where __powf(2, r) is most accurate. | 5554 | // drops from ~1e-8 to ~1e-15. r is confined to [-0.5, 0.5], where __powf(2, r) is most accurate. |
| 5564 | float r = __fma(x, 3.3219280242919921875f, -k); | 5555 | float r = __fma(x, 3.3219280242919921875f, -k); |
| 5565 | r = __fma(x, 7.0595369550119357882e-08f, r); | 5556 | r = __fma(x, 7.0595369550119357882e-08f, r); |
| 5566 | - return scale * __powf(2.0f, r); | 5557 | + float result = scale * __powf(2.0f, r); |
| 5558 | + | ||
| 5559 | + // 10^x = 2^(x * log2(10)). Split x*log2(10) into an integer k and a small remainder r, then compute | ||
| 5560 | + // 2^k * 2^r with 2^k constructed exactly from the bit pattern. | ||
| 5561 | + if (__isinf(x)) { | ||
| 5562 | + result = x > 0.0f ? x : 0.0f; | ||
| 5563 | + } | ||
| 5564 | + | ||
| 5565 | + if (__isnan(x)) { | ||
| 5566 | + result = x; | ||
| 5567 | + } | ||
| 5568 | + return result; | ||
| 5567 | } | 5569 | } |
| 5568 | 5570 | ||
| 5569 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float expm1f(float x) | 5571 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float expm1f(float x) |
| 5570 | { | 5572 | { |
| 5571 | - // expm1f(x) = e^x - 1, computed accurately for small |x| where e^x-1 would otherwise lose | ||
| 5572 | - // precision to catastrophic cancellation. The approach follows the same range-reduction + | ||
| 5573 | - // 2^k * 2^r decomposition as expf: split x * log2(e) into integer k and remainder r, build | ||
| 5574 | - // 2^k exactly from its bit pattern, then reconstruct (2^k * 2^r) - 1 in a cancellation-safe form. | ||
| 5575 | - if (x == 0.0f) { | ||
| 5576 | - return x; | ||
| 5577 | - } | ||
| 5578 | - if (__isnan(x)) { | ||
| 5579 | - return x; | ||
| 5580 | - } | ||
| 5581 | - if (__isinf(x)) { | ||
| 5582 | - return x > 0.0f ? x : -1.0f; | ||
| 5583 | - } | ||
| 5584 | - | ||
| 5585 | // Clamp |x| to the float32 overflow threshold (log2(e) * 128 ~= 88.72) so the reduction below | 5573 | // Clamp |x| to the float32 overflow threshold (log2(e) * 128 ~= 88.72) so the reduction below |
| 5586 | // stays within the representable range; the clamped value still selects the correct overflow path. | 5574 | // stays within the representable range; the clamped value still selects the correct overflow path. |
| 5587 | float ax = __fabsf(x); | 5575 | float ax = __fabsf(x); |
| @@ -5629,6 +5617,18 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float expm1f(float x) | |||
| 5629 | if (large_k) { | 5617 | if (large_k) { |
| 5630 | y = y + y; | 5618 | y = y + y; |
| 5631 | } | 5619 | } |
| 5620 | + | ||
| 5621 | + // expm1f(x) = e^x - 1, computed accurately for small |x| where e^x-1 would otherwise lose | ||
| 5622 | + // precision to catastrophic cancellation. The approach follows the same range-reduction + | ||
| 5623 | + // 2^k * 2^r decomposition as expf: split x * log2(e) into integer k and remainder r, build | ||
| 5624 | + // 2^k exactly from its bit pattern, then reconstruct (2^k * 2^r) - 1 in a cancellation-safe form. | ||
| 5625 | + if (__isinf(x)) { | ||
| 5626 | + y = x > 0.0f ? x : -1.0f; | ||
| 5627 | + } | ||
| 5628 | + | ||
| 5629 | + if (__isnan(x) || x == 0.0f) { | ||
| 5630 | + y = x; | ||
| 5631 | + } | ||
| 5632 | return y; | 5632 | return y; |
| 5633 | } | 5633 | } |
| 5634 | 5634 | ||
| @@ -5657,35 +5657,38 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline int32_t __internal_ilogbf_finite_abs(float | |||
| 5657 | // logbf(x) = (float) floor(log2(|x|)) = the unbiased exponent of |x| as a float. | 5657 | // logbf(x) = (float) floor(log2(|x|)) = the unbiased exponent of |x| as a float. |
| 5658 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float logbf(float x) | 5658 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float logbf(float x) |
| 5659 | { | 5659 | { |
| 5660 | - if (__isnan(x)) { | ||
| 5661 | - return x; // NaN propagates. | ||
| 5662 | - } | ||
| 5663 | - | ||
| 5664 | float ax = __fabsf(x); | 5660 | float ax = __fabsf(x); |
| 5665 | - if (ax == 0.0f) { | 5661 | + float result = static_cast<float>(__internal_ilogbf_finite_abs(ax)); |
| 5666 | - return -ASCRT_INF_F; // logb(0) = -inf, matching the C standard. | 5662 | + |
| 5667 | - } | ||
| 5668 | if (ax == ASCRT_INF_F) { | 5663 | if (ax == ASCRT_INF_F) { |
| 5669 | - return ASCRT_INF_F; // logb(inf) = +inf. | 5664 | + result = ASCRT_INF_F; // logb(inf) = +inf. |
| 5670 | } | 5665 | } |
| 5671 | 5666 | ||
| 5672 | - return static_cast<float>(__internal_ilogbf_finite_abs(ax)); | 5667 | + if (ax == 0.0f) { |
| 5668 | + result = -ASCRT_INF_F; // logb(0) = -inf, matching the C standard. | ||
| 5669 | + } | ||
| 5670 | + | ||
| 5671 | + if (__isnan(x)) { | ||
| 5672 | + result = x; // NaN propagates. | ||
| 5673 | + } | ||
| 5674 | + return result; | ||
| 5673 | } | 5675 | } |
| 5674 | 5676 | ||
| 5675 | // ilogbf(x): integer variant of logbf. Returns the unbiased exponent as int32_t, with the IEEE-754 | 5677 | // ilogbf(x): integer variant of logbf. Returns the unbiased exponent as int32_t, with the IEEE-754 |
| 5676 | // special-value encoding: NaN or 0 -> INT_MIN (0x80000000), +inf / -inf -> INT_MAX (0x7FFFFFFF). | 5678 | // special-value encoding: NaN or 0 -> INT_MIN (0x80000000), +inf / -inf -> INT_MAX (0x7FFFFFFF). |
| 5677 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline int32_t ilogbf(float x) | 5679 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline int32_t ilogbf(float x) |
| 5678 | { | 5680 | { |
| 5679 | - if (__isnan(x) || x == 0.0f) { | ||
| 5680 | - return static_cast<int32_t>(0x80000000U); // INT_MIN: NaN and 0 map to this sentinel. | ||
| 5681 | - } | ||
| 5682 | - | ||
| 5683 | float ax = __fabsf(x); | 5681 | float ax = __fabsf(x); |
| 5682 | + int32_t result = __internal_ilogbf_finite_abs(ax); | ||
| 5683 | + | ||
| 5684 | if (ax == ASCRT_INF_F) { | 5684 | if (ax == ASCRT_INF_F) { |
| 5685 | - return static_cast<int32_t>(0x7FFFFFFFU); // INT_MAX: infinity maps to this sentinel. | 5685 | + result = ASCRT_MAX_VAL_S; // INT_MAX: infinity maps to this sentinel. |
| 5686 | + } | ||
| 5687 | + if (__isnan(x) || x == 0.0f) { | ||
| 5688 | + result = ASCRT_MIN_VAL_S; // INT_MIN: NaN and 0 map to this sentinel. | ||
| 5686 | } | 5689 | } |
| 5687 | 5690 | ||
| 5688 | - return __internal_ilogbf_finite_abs(ax); | 5691 | + return result; |
| 5689 | } | 5692 | } |
| 5690 | 5693 | ||
| 5691 | /* | 5694 | /* |