Merged
float类型接口性能优化 #6676
xuyiyangAdd on 26 days ago
float类型接口性能优化 #6676
Merged
Total: 1 files changed+117-114
| @@ -3680,46 +3680,61 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ static inline int32_t __asc_float_mantissa_expone | |||
| 3680 | // Return the exponent aligned to the mantissa integer returned above. | 3680 | // Return the exponent aligned to the mantissa integer returned above. |
| 3681 | // Subnormals are treated as exponent -149, normals as unbiased exponent minus 23. | 3681 | // Subnormals are treated as exponent -149, normals as unbiased exponent minus 23. |
| 3682 | const uint32_t exponent_bits = bits & ASCRT_EXP_BIT_FLOAT_U; | 3682 | const uint32_t exponent_bits = bits & ASCRT_EXP_BIT_FLOAT_U; |
| 3683 | - // -149: min subnormal exp, -150: bias + mantissa shift | 3683 | + // aligning the biased exponent field to the integer mantissa (value = mantissa * 2^exp) |
| 3684 | - return exponent_bits == 0U ? -149 : (static_cast<int32_t>(exponent_bits >> 23U) - 150); | 3684 | + return exponent_bits == 0U ? |
| 3685 | + -149 : | ||
| 3686 | + (static_cast<int32_t>(exponent_bits >> __internal_fp32_exponent_shift) - | ||
| 3687 | + 150); // -149: min subnormal exp, -150: bias + mantissa shift, 23: exponent field low bit | ||
| 3685 | } | 3688 | } |
| 3686 | 3689 | ||
| 3687 | __SIMT_DEVICE_FUNCTIONS_DECL__ static inline int32_t __asc_uint_floor_log2(uint32_t x) | 3690 | __SIMT_DEVICE_FUNCTIONS_DECL__ static inline int32_t __asc_uint_floor_log2(uint32_t x) |
| 3688 | { | 3691 | { |
| 3689 | - // Compute floor(log2(x)) for a nonzero integer mantissa. | 3692 | + // Branchless floor(log2(x)) for a nonzero integer mantissa below 2^24: the float32 |
| 3690 | - int32_t shift = 0; | 3693 | + // conversion of such integers is exact, so the biased exponent field of float(x) |
| 3691 | - while ((x >> 1U) != 0U) { | 3694 | + // equals floor(log2(x)) + 127. |
| 3692 | - x >>= 1U; | 3695 | + const uint32_t exp_field = (__float_as_uint(static_cast<float>(x)) >> __internal_fp32_exponent_shift) & |
| 3693 | - ++shift; | 3696 | + 0xFFU; // 0xFF: mask of the 8-bit biased exponent field |
| 3694 | - } | 3697 | + return static_cast<int32_t>(exp_field) - __internal_fp32_exponent_bias; |
| 3695 | - return shift; | ||
| 3696 | } | 3698 | } |
| 3697 | 3699 | ||
| 3698 | __SIMT_DEVICE_FUNCTIONS_DECL__ static inline float __asc_make_float_from_mantissa(uint32_t mantissa, int32_t exponent) | 3700 | __SIMT_DEVICE_FUNCTIONS_DECL__ static inline float __asc_make_float_from_mantissa(uint32_t mantissa, int32_t exponent) |
| 3699 | { | 3701 | { |
| 3700 | // Rebuild a float from an integer mantissa and an exponent aligned to that mantissa. | 3702 | // Rebuild a float from an integer mantissa and an exponent aligned to that mantissa. |
| 3701 | - float result = 0.0f; | 3703 | + const int32_t mantissa_log2 = __asc_uint_floor_log2(mantissa != 0U ? mantissa : 1U); |
| 3702 | - if (mantissa != 0U) { | 3704 | + const int32_t result_exponent = exponent + mantissa_log2; |
| 3703 | - const int32_t mantissa_log2 = __asc_uint_floor_log2(mantissa); | 3705 | + const bool is_normal = |
| 3704 | - const int32_t result_exponent = exponent + mantissa_log2; | 3706 | + (result_exponent >= -126) && (result_exponent <= 127); // -126: min normal exp, 127: max normal exp |
| 3705 | - // Clamp to fp32 range before reconstructing the final bit pattern. | 3707 | + const bool is_subnormal = |
| 3706 | - if (result_exponent > 127) { // 127: largest normal fp32 exponent, beyond it overflows to inf | 3708 | + (result_exponent < -126) && (result_exponent >= -149); // -126: min normal exp, -149: min subnormal exp |
| 3707 | - result = ASCRT_INF_F; | 3709 | + |
| 3708 | - } else if (result_exponent < -149) { // -149: exponent of the smallest fp32 subnormal, below it underflows to 0 | 3710 | + // Normal result: normalize mantissa and pack signless fp32 bits. |
| 3709 | - result = 0.0f; | 3711 | + const uint32_t normal_bits = |
| 3710 | - } else if (result_exponent >= -126) { // -126: smallest normal fp32 exponent | 3712 | + (static_cast<uint32_t>(result_exponent + __internal_fp32_exponent_bias) << __internal_fp32_exponent_shift) | |
| 3711 | - // Normal result: normalize mantissa and pack signless fp32 bits. | 3713 | + ((mantissa << static_cast<uint32_t>(23 - mantissa_log2)) & |
| 3712 | - const uint32_t normalized_mantissa = mantissa << static_cast<uint32_t>(23 - mantissa_log2); | 3714 | + ASCRT_MAN_BIT_FLOAT_U); // 23: mantissa bit width, implicit-one normalization |
| 3713 | - const uint32_t bits = | 3715 | + const int32_t raw_sub_shift = |
| 3714 | - (static_cast<uint32_t>(result_exponent + 127) << 23U) | (normalized_mantissa & ASCRT_MAN_BIT_FLOAT_U); | 3716 | + exponent + |
| 3715 | - result = __uint_as_float(bits); | 3717 | + 149; // 149: 126 (min normal exp offset) + 23 (mantissa bit width), lands bit 0 at the 2^-149 subnormal weight |
| 3716 | - } else { | 3718 | + const uint32_t left_shift = |
| 3717 | - // Subnormal result: shift mantissa directly into the denormal field. | 3719 | + raw_sub_shift < 0 ? |
| 3718 | - const uint32_t subnormal_shift = static_cast<uint32_t>(exponent + 149); | 3720 | + 0U : |
| 3719 | - result = __uint_as_float(mantissa << subnormal_shift); | 3721 | + (raw_sub_shift > 31 ? |
| 3720 | - } | 3722 | + 31U : |
| 3721 | - } | 3723 | + static_cast<uint32_t>(raw_sub_shift)); // 31: clamp, shift count must stay within the uint32 width |
| 3722 | - return result; | 3724 | + const int32_t raw_right = -raw_sub_shift; |
| 3725 | + const uint32_t right_shift = | ||
| 3726 | + raw_right < 0 ? | ||
| 3727 | + 0U : | ||
| 3728 | + (raw_right > 31 ? | ||
| 3729 | + 31U : | ||
| 3730 | + static_cast<uint32_t>(raw_right)); // 31: clamp, shift count must stay within the uint32 width | ||
| 3731 | + const uint32_t subnormal_bits = raw_sub_shift >= 0 ? (mantissa << left_shift) : (mantissa >> right_shift); | ||
| 3732 | + | ||
| 3733 | + uint32_t bits = is_normal ? normal_bits : subnormal_bits; | ||
| 3734 | + bits = is_subnormal ? subnormal_bits : bits; | ||
| 3735 | + bits = (result_exponent > 127) ? ASCRT_INF_U : bits; // beyond max normal exponent overflows to inf | ||
| 3736 | + bits = (mantissa == 0U || result_exponent < -149) ? 0U : bits; // zero mantissa or below min subnormal | ||
| 3737 | + return __uint_as_float(bits); | ||
| 3723 | } | 3738 | } |
| 3724 | 3739 | ||
| 3725 | __SIMT_DEVICE_FUNCTIONS_DECL__ static inline uint32_t __asc_double_abs_float_bits(uint32_t abs_bits) | 3740 | __SIMT_DEVICE_FUNCTIONS_DECL__ static inline uint32_t __asc_double_abs_float_bits(uint32_t abs_bits) |
| @@ -3752,39 +3767,42 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ static inline float __asc_remquo_core( | |||
| 3752 | } | 3767 | } |
| 3753 | 3768 | ||
| 3754 | // Split both operands into integer mantissas and aligned exponents. | 3769 | // Split both operands into integer mantissas and aligned exponents. |
| 3755 | - uint32_t quotient_low = 0U; | 3770 | + const uint32_t x_mantissa = __asc_float_mantissa_bits_from_abs_bits(abs_x_bits); |
| 3756 | - uint32_t remainder_mantissa = __asc_float_mantissa_bits_from_abs_bits(abs_x_bits); | 3771 | + uint32_t y_mantissa = __asc_float_mantissa_bits_from_abs_bits(abs_y_bits); |
| 3757 | - const uint32_t y_mantissa = __asc_float_mantissa_bits_from_abs_bits(abs_y_bits); | ||
| 3758 | const int32_t x_exponent = __asc_float_mantissa_exponent_from_abs_bits(abs_x_bits); | 3772 | const int32_t x_exponent = __asc_float_mantissa_exponent_from_abs_bits(abs_x_bits); |
| 3759 | - const int32_t y_exponent = __asc_float_mantissa_exponent_from_abs_bits(abs_y_bits); | 3773 | + int32_t y_exponent = __asc_float_mantissa_exponent_from_abs_bits(abs_y_bits); |
| 3760 | - const int32_t exponent_diff = x_exponent - y_exponent; | ||
| 3761 | 3774 | ||
| 3762 | - // Bit-by-bit long division over the mantissa domain. | 3775 | + // Rescale y's mantissa into [2^23, 2^24) (its value is preserved by lowering the |
| 3763 | - if (exponent_diff >= 0 && y_mantissa != 0U) { | 3776 | + // aligned exponent to compensate). With a normalized divisor the integer part of |
| 3764 | - uint32_t remainder_bits = 0U; | 3777 | + // the mantissa quotient is a single bit, so the whole mantissa pass collapses to |
| 3765 | - const int32_t x_mantissa_log2 = __asc_uint_floor_log2(remainder_mantissa); | 3778 | + // one compare-and-subtract step instead of a 24-step bit loop. The rescaled |
| 3766 | - for (int32_t bit = x_mantissa_log2; bit >= 0; --bit) { | 3779 | + // remainder stays exact: both modulus and remainder scale by the same power of two. |
| 3767 | - remainder_bits = (remainder_bits << 1U) | ((remainder_mantissa >> static_cast<uint32_t>(bit)) & 1U); | 3780 | + const uint32_t y_shift = static_cast<uint32_t>( |
| 3768 | - quotient_low = (quotient_low << 1U) & ASCRT_REMQUO_MASK_F; | 3781 | + 23 - __asc_uint_floor_log2(y_mantissa)); // 23: mantissa bit width, target of the [2^23, 2^24) rescale |
| 3769 | - if (remainder_bits >= y_mantissa) { | 3782 | + y_mantissa <<= y_shift; |
| 3770 | - remainder_bits -= y_mantissa; | 3783 | + y_exponent -= static_cast<int32_t>(y_shift); |
| 3771 | - quotient_low |= 1U; | 3784 | + int32_t exponent_diff = x_exponent - y_exponent; |
| 3772 | - } | 3785 | + exponent_diff = exponent_diff < 0 ? 0 : exponent_diff; // |x| >= |y| keeps this >= 0; guard only |
| 3773 | - } | 3786 | + |
| 3774 | - for (int32_t bit = 0; bit < exponent_diff; ++bit) { | 3787 | + // Pass 1: single step, floor(x_mantissa / y_mantissa) < 2 by construction. |
| 3775 | - remainder_bits <<= 1U; | 3788 | + const bool mantissa_ge = x_mantissa >= y_mantissa; |
| 3776 | - quotient_low = (quotient_low << 1U) & ASCRT_REMQUO_MASK_F; | 3789 | + uint32_t remainder_bits = mantissa_ge ? (x_mantissa - y_mantissa) : x_mantissa; |
| 3777 | - if (remainder_bits >= y_mantissa) { | 3790 | + uint32_t quotient_low = mantissa_ge ? 1U : 0U; |
| 3778 | - remainder_bits -= y_mantissa; | 3791 | + // Pass 2: fractional continuation, one step per exponent difference (inherently |
| 3779 | - quotient_low |= 1U; | 3792 | + // data dependent; select-based body). |
| 3780 | - } | 3793 | + for (int32_t bit = 0; bit < exponent_diff; ++bit) { |
| 3781 | - } | 3794 | + const uint32_t pulled = remainder_bits |
| 3782 | - remainder_mantissa = remainder_bits; | 3795 | + << 1U; // 1U: long-division shift, pulls the next remainder bit into the compare |
| 3796 | + const uint32_t sub = pulled >= y_mantissa ? y_mantissa : 0U; | ||
| 3797 | + const uint32_t quo_bit = pulled >= y_mantissa ? 1U : 0U; | ||
| 3798 | + remainder_bits = pulled - sub; | ||
| 3799 | + quotient_low = ((quotient_low << 1U) & ASCRT_REMQUO_MASK_F) | | ||
| 3800 | + quo_bit; // 1U: quotient bit shift, ASCRT_REMQUO_MASK_F keeps the low bits (remquo's low 3) | ||
| 3783 | } | 3801 | } |
| 3784 | 3802 | ||
| 3785 | // Reassemble the provisional remainder on y's exponent scale. | 3803 | // Reassemble the provisional remainder on y's exponent scale. |
| 3786 | *quo_abs_low = quotient_low; | 3804 | *quo_abs_low = quotient_low; |
| 3787 | - return __asc_make_float_from_mantissa(remainder_mantissa, y_exponent); | 3805 | + return __asc_make_float_from_mantissa(remainder_bits, y_exponent); |
| 3788 | } | 3806 | } |
| 3789 | 3807 | ||
| 3790 | __SIMT_DEVICE_FUNCTIONS_DECL__ static inline void __asc_remquo_tail( | 3808 | __SIMT_DEVICE_FUNCTIONS_DECL__ static inline void __asc_remquo_tail( |
| @@ -3996,8 +4014,8 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline void __internal_sincospif_core(float x, fl | |||
| 3996 | const float two_x = safe_x + safe_x; | 4014 | const float two_x = safe_x + safe_x; |
| 3997 | const int32_t quadrant = __float2int_rn(two_x); | 4015 | const int32_t quadrant = __float2int_rn(two_x); |
| 3998 | const float rounded_two_x = __int2float_rn(quadrant); | 4016 | const float rounded_two_x = __int2float_rn(quadrant); |
| 3999 | - const bool is_integer = !is_special && (truncf(x) == x); | ||
| 4000 | const float reduced = fmaf(-rounded_two_x, 0.5f, safe_x); | 4017 | const float reduced = fmaf(-rounded_two_x, 0.5f, safe_x); |
| 4018 | + const bool is_integer = !is_special && (reduced == 0.0f) && ((quadrant & 1) == 0); | ||
| 4001 | const float reduced2 = reduced * reduced; | 4019 | const float reduced2 = reduced * reduced; |
| 4002 | 4020 | ||
| 4003 | // sin(pi*r) polynomial on the reduced argument. | 4021 | // sin(pi*r) polynomial on the reduced argument. |
| @@ -4151,8 +4169,9 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float rhypotf(float x, float y) | |||
| 4151 | const float abs_y = fabsf(y); | 4169 | const float abs_y = fabsf(y); |
| 4152 | const uint32_t abs_x_bits = __float_as_uint(abs_x); | 4170 | const uint32_t abs_x_bits = __float_as_uint(abs_x); |
| 4153 | const uint32_t abs_y_bits = __float_as_uint(abs_y); | 4171 | const uint32_t abs_y_bits = __float_as_uint(abs_y); |
| 4154 | - const uint32_t min_bits = abs_x_bits < abs_y_bits ? abs_x_bits : abs_y_bits; | 4172 | + const bool x_is_greater = abs_x_bits > abs_y_bits; |
| 4155 | - const uint32_t max_bits = abs_x_bits > abs_y_bits ? abs_x_bits : abs_y_bits; | 4173 | + const uint32_t max_bits = x_is_greater ? abs_x_bits : abs_y_bits; |
| 4174 | + const uint32_t min_bits = x_is_greater ? abs_y_bits : abs_x_bits; | ||
| 4156 | const float min_abs = __uint_as_float(min_bits); | 4175 | const float min_abs = __uint_as_float(min_bits); |
| 4157 | const float max_abs = __uint_as_float(max_bits); | 4176 | const float max_abs = __uint_as_float(max_bits); |
| 4158 | 4177 | ||
| @@ -4163,7 +4182,7 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float rhypotf(float x, float y) | |||
| 4163 | const float scaled_max = max_abs * scale; | 4182 | const float scaled_max = max_abs * scale; |
| 4164 | // Evaluate the scaled square sum and take the reciprocal square root. | 4183 | // Evaluate the scaled square sum and take the reciprocal square root. |
| 4165 | const float square_sum = fmaf(scaled_max, scaled_max, scaled_min * scaled_min); | 4184 | const float square_sum = fmaf(scaled_max, scaled_max, scaled_min * scaled_min); |
| 4166 | - const float result = rsqrtf(square_sum) * scale; | 4185 | + const float result = 1.0f / __sqrtf(square_sum) * scale; |
| 4167 | 4186 | ||
| 4168 | // If the smaller magnitude is infinite, the reciprocal hypotenuse collapses to zero. | 4187 | // If the smaller magnitude is infinite, the reciprocal hypotenuse collapses to zero. |
| 4169 | return min_abs == ASCRT_INF_F ? 0.0f : result; | 4188 | return min_abs == ASCRT_INF_F ? 0.0f : result; |
| @@ -4193,19 +4212,10 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float rhypotf(float x, float y) | |||
| 4193 | */ | 4212 | */ |
| 4194 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float rnorm3df(float a, float b, float c) | 4213 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float rnorm3df(float a, float b, float c) |
| 4195 | { | 4214 | { |
| 4196 | - // Shared scaling constants used to keep the squared sum numerically safe. | ||
| 4197 | - constexpr float rsqrt_subnormal_scale = 4096.0f; // 2^12, compensates the temporary scaling after rsqrt refinement. | ||
| 4198 | - | ||
| 4199 | // Work with absolute values so the norm depends only on magnitudes. | 4215 | // Work with absolute values so the norm depends only on magnitudes. |
| 4200 | const float abs_a = fabsf(a); | 4216 | const float abs_a = fabsf(a); |
| 4201 | const float abs_b = fabsf(b); | 4217 | const float abs_b = fabsf(b); |
| 4202 | const float abs_c = fabsf(c); | 4218 | const float abs_c = fabsf(c); |
| 4203 | - if (isinf(abs_a) || isinf(abs_b) || isinf(abs_c)) { | ||
| 4204 | - return 0.0f; | ||
| 4205 | - } | ||
| 4206 | - if (isnan(abs_a) || isnan(abs_b) || isnan(abs_c)) { | ||
| 4207 | - return ASCRT_INF_F / ASCRT_INF_F; | ||
| 4208 | - } | ||
| 4209 | 4219 | ||
| 4210 | // Sort magnitudes so the largest term determines the scaling exponent. | 4220 | // Sort magnitudes so the largest term determines the scaling exponent. |
| 4211 | const float max_ab = fmaxf(abs_a, abs_b); | 4221 | const float max_ab = fmaxf(abs_a, abs_b); |
| @@ -4220,28 +4230,26 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float rnorm3df(float a, float b, float c) | |||
| 4220 | const float scaled_min = min_ab * scale; | 4230 | const float scaled_min = min_ab * scale; |
| 4221 | const float scaled_mid = mid_abs * scale; | 4231 | const float scaled_mid = mid_abs * scale; |
| 4222 | const float scaled_max = max_abs * scale; | 4232 | const float scaled_max = max_abs * scale; |
| 4223 | - float square_sum = scaled_mid * scaled_mid; | 4233 | + float square_sum = fmaf(scaled_min, scaled_min, fmaf(scaled_mid, scaled_mid, scaled_max * scaled_max)); |
| 4224 | - square_sum = fmaf(scaled_min, scaled_min, square_sum); | ||
| 4225 | - square_sum = fmaf(scaled_max, scaled_max, square_sum); | ||
| 4226 | - | ||
| 4227 | - // Subnormal squared sums are temporarily lifted before rsqrtf. | ||
| 4228 | - const bool is_subnormal_square_sum = fabsf(square_sum) < __internal_subnormal_boundary; | ||
| 4229 | - if (is_subnormal_square_sum) { | ||
| 4230 | - square_sum *= __internal_fp32_scale_2p24; | ||
| 4231 | - } | ||
| 4232 | - // Zero norm means the reciprocal norm is infinite. | ||
| 4233 | - if (square_sum == 0.0f) { | ||
| 4234 | - return ASCRT_INF_F; | ||
| 4235 | - } | ||
| 4236 | 4234 | ||
| 4237 | // Compute an approximate reciprocal norm and refine it with one Newton step. | 4235 | // Compute an approximate reciprocal norm and refine it with one Newton step. |
| 4238 | - float inv_norm = rsqrtf(square_sum); | 4236 | + float inv_norm = 1.0f / __sqrtf(square_sum); |
| 4239 | inv_norm = inv_norm * fmaf(-0.5f * square_sum, inv_norm * inv_norm, 1.5f); | 4237 | inv_norm = inv_norm * fmaf(-0.5f * square_sum, inv_norm * inv_norm, 1.5f); |
| 4240 | - if (is_subnormal_square_sum) { | ||
| 4241 | - inv_norm *= rsqrt_subnormal_scale; | ||
| 4242 | - } | ||
| 4243 | // Restore the outer scale factor. | 4238 | // Restore the outer scale factor. |
| 4244 | - return scale * inv_norm; | 4239 | + float result = scale * inv_norm; |
| 4240 | + // Zero norm means the reciprocal norm is infinite. | ||
| 4241 | + if (square_sum == 0.0f) { | ||
| 4242 | + result = ASCRT_INF_F; | ||
| 4243 | + } | ||
| 4244 | + | ||
| 4245 | + if (isnan(abs_a) || isnan(abs_b) || isnan(abs_c)) { | ||
| 4246 | + result = ASCRT_NAN_F; | ||
| 4247 | + } | ||
| 4248 | + | ||
| 4249 | + if (isinf(abs_a) || isinf(abs_b) || isinf(abs_c)) { | ||
| 4250 | + result = 0.0f; | ||
| 4251 | + } | ||
| 4252 | + return result; | ||
| 4245 | } | 4253 | } |
| 4246 | 4254 | ||
| 4247 | /** | 4255 | /** |
| @@ -4447,11 +4455,8 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float tanpif(float x) | |||
| 4447 | 4455 | ||
| 4448 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_euler_gamma_function(float x) | 4456 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_euler_gamma_function(float x) |
| 4449 | { | 4457 | { |
| 4450 | - // Reduce the argument around the nearest integer when x is nonzero. | 4458 | + // Reduce the argument around the nearest integer. |
| 4451 | - float reduced = x; | 4459 | + const float reduced = x - nearbyintf(x); |
| 4452 | - if (x != 0.0f) { | ||
| 4453 | - reduced = x - nearbyintf(x); | ||
| 4454 | - } | ||
| 4455 | 4460 | ||
| 4456 | // Polynomial approximation for the Euler gamma correction term. | 4461 | // Polynomial approximation for the Euler gamma correction term. |
| 4457 | float poly = fmaf(reduced, -0.0010728656779974699020f, 0.0071110534481704235077f); | 4462 | float poly = fmaf(reduced, -0.0010728656779974699020f, 0.0071110534481704235077f); |
| @@ -4469,13 +4474,12 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_euler_gamma_function(floa | |||
| 4469 | } | 4474 | } |
| 4470 | 4475 | ||
| 4471 | float denominator = fmaf(scale * poly, reduced, scale); | 4476 | float denominator = fmaf(scale * poly, reduced, scale); |
| 4472 | - // Fall back to a scaled reciprocal when the denominator becomes subnormal. | 4477 | + // The denominator stays below ~2 on this path, so scaling by 2^24 never |
| 4473 | - if (fabsf(denominator) < __internal_subnormal_boundary) { | 4478 | + // overflows and lifts a subnormal (or zero) denominator into the normal |
| 4474 | - int32_t e = 0; | 4479 | + // range. Dividing the lifted value and scaling back keeps the division |
| 4475 | - float m = frexpf(denominator, &e); | 4480 | + // correctly rounded (single rounding, exact power-of-two scaling) |
| 4476 | - return ldexpf(1.0f / m, 0 - e); | 4481 | + const float lifted_denominator = denominator * __internal_fp32_scale_2p24; |
| 4477 | - } | 4482 | + return (1.0f / lifted_denominator) * __internal_fp32_scale_2p24; |
| 4478 | - return 1.0f / denominator; | ||
| 4479 | } | 4483 | } |
| 4480 | 4484 | ||
| 4481 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_cal_abs_x(float x) | 4485 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_cal_abs_x(float x) |
| @@ -4613,17 +4617,16 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_stirling_and_euler_reflec | |||
| 4613 | const bool use_cos_poly = ((sinpi_quadrant & 1) != 0); | 4617 | const bool use_cos_poly = ((sinpi_quadrant & 1) != 0); |
| 4614 | float reduced_pi_arg = fmaf(rounded_two_abs_x, -0.5f, abs_x); | 4618 | float reduced_pi_arg = fmaf(rounded_two_abs_x, -0.5f, abs_x); |
| 4615 | const float reduced_pi_arg_square = reduced_pi_arg * reduced_pi_arg; | 4619 | const float reduced_pi_arg_square = reduced_pi_arg * reduced_pi_arg; |
| 4616 | - float sinpi_poly = 0.0f; | 4620 | + // Evaluate both quadrant polynomials unconditionally and select, so |
| 4617 | - if (use_cos_poly) { | 4621 | + // the quadrant split never diverges the warp. |
| 4618 | - sinpi_poly = fmaf(reduced_pi_arg_square, 0.22686031460762023926f, -1.334560394287109375f); | 4622 | + float cos_poly = fmaf(reduced_pi_arg_square, 0.22686031460762023926f, -1.334560394287109375f); |
| 4619 | - sinpi_poly = fmaf(reduced_pi_arg_square, sinpi_poly, 4.0586924552917480469f); | 4623 | + cos_poly = fmaf(reduced_pi_arg_square, cos_poly, 4.0586924552917480469f); |
| 4620 | - sinpi_poly = fmaf(reduced_pi_arg_square, sinpi_poly, -4.9348020553588867188f); | 4624 | + cos_poly = fmaf(reduced_pi_arg_square, cos_poly, -4.9348020553588867188f); |
| 4621 | - sinpi_poly = fmaf(reduced_pi_arg_square, sinpi_poly, 1.0f); | 4625 | + cos_poly = fmaf(reduced_pi_arg_square, cos_poly, 1.0f); |
| 4622 | - } else { | 4626 | + float sin_poly = fmaf(reduced_pi_arg_square, -0.59248024225234985352f, 2.550144195556640625f); |
| 4623 | - sinpi_poly = fmaf(reduced_pi_arg_square, -0.59248024225234985352f, 2.550144195556640625f); | 4627 | + sin_poly = fmaf(reduced_pi_arg_square, sin_poly, -5.1677198410034179688f); |
| 4624 | - sinpi_poly = fmaf(reduced_pi_arg_square, sinpi_poly, -5.1677198410034179688f); | 4628 | + sin_poly = fmaf(sin_poly, reduced_pi_arg * reduced_pi_arg_square, reduced_pi_arg * ASCRT_PI_F); |
| 4625 | - sinpi_poly = fmaf(sinpi_poly, reduced_pi_arg * reduced_pi_arg_square, reduced_pi_arg * ASCRT_PI_F); | 4629 | + const float sinpi_poly = use_cos_poly ? cos_poly : sin_poly; |
| 4626 | - } | ||
| 4627 | float sinpi = ((sinpi_quadrant & 2) != 0) ? -sinpi_poly : sinpi_poly; | 4630 | float sinpi = ((sinpi_quadrant & 2) != 0) ? -sinpi_poly : sinpi_poly; |
| 4628 | 4631 | ||
| 4629 | // Track the residual of the reflection denominator for the final correction. | 4632 | // Track the residual of the reflection denominator for the final correction. |
| @@ -4879,7 +4882,7 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float sinhf(float x) | |||
| 4879 | poly *= x2; | 4882 | poly *= x2; |
| 4880 | 4883 | ||
| 4881 | // Split the exponent for the large-argument exp2 reconstruction. | 4884 | // Split the exponent for the large-argument exp2 reconstruction. |
| 4882 | - float n = truncf(abs_x * 1.4426950216293334961f); | 4885 | + float n = floorf(abs_x * 1.4426950216293334961f); |
| 4883 | if (fabsf(n) > 126.0f) { | 4886 | if (fabsf(n) > 126.0f) { |
| 4884 | n = __internal_with_sign_bit(126.0f, n); | 4887 | n = __internal_with_sign_bit(126.0f, n); |
| 4885 | } | 4888 | } |
| @@ -4887,7 +4890,7 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float sinhf(float x) | |||
| 4887 | // Reconstruct the exponential tail with hi/lo corrections. | 4890 | // Reconstruct the exponential tail with hi/lo corrections. |
| 4888 | float r = fmaf(n, -0.69314718246459960938f, abs_x); | 4891 | float r = fmaf(n, -0.69314718246459960938f, abs_x); |
| 4889 | r = fmaf(n, 1.9046542121259335545e-09f, r); | 4892 | r = fmaf(n, 1.9046542121259335545e-09f, r); |
| 4890 | - const float exp2_residual = exp2f(r * 1.4426950216293334961f); | 4893 | + const float exp2_residual = __powf(2.0f, r * 1.4426950216293334961f); |
| 4891 | 4894 | ||
| 4892 | const float exponent_base = n + 12583037.0f; | 4895 | const float exponent_base = n + 12583037.0f; |
| 4893 | const uint32_t exponent_base_bits = __float_as_uint(exponent_base); | 4896 | const uint32_t exponent_base_bits = __float_as_uint(exponent_base); |