Merged
float类型数学接口性能优化 #6829
xuyiyangAdd on 19 days ago
float类型数学接口性能优化 #6829
Merged
Total: 1 files changed+130-127
| @@ -867,6 +867,27 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float cosf(float x) | |||
| 867 | return c; | 867 | return c; |
| 868 | } | 868 | } |
| 869 | 869 | ||
| 870 | +// Cody-Waite-only cosine: identical to cosf for |x| <= 71476.0625f (its | ||
| 871 | +// Payne-Hanek threshold). Callers that can prove the argument is inside the | ||
| 872 | +// range use this to keep the trig_radian_reduction threshold branch (and the | ||
| 873 | +// Payne-Hanek table-index branch) out of their instruction stream. | ||
| 874 | +__SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_cosf_cody_waite(float x) | ||
| 875 | +{ | ||
| 876 | + x = fmaf(x, 0.0f, x); // NaN for Inf | ||
| 877 | + int quadrant; | ||
| 878 | + float y = __internal_cody_waite_radian_reduction(x, &quadrant); | ||
| 879 | + float c = __internal_cos_poly(y); | ||
| 880 | + float s = __internal_sin_poly(y); | ||
| 881 | + if (quadrant & 2) { // 2: quadrants where cos and sin flip sign | ||
| 882 | + c = 0.0f - c; | ||
| 883 | + s = 0.0f - s; | ||
| 884 | + } | ||
| 885 | + if (quadrant & 1) { // 1: quadrants where cos(pi/2 + y) = -sin(y) | ||
| 886 | + c = 0.0f - s; | ||
| 887 | + } | ||
| 888 | + return c; | ||
| 889 | +} | ||
| 890 | + | ||
| 870 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float sinf(float x) | 891 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float sinf(float x) |
| 871 | { | 892 | { |
| 872 | int quadrant; | 893 | int quadrant; |
| @@ -2098,7 +2119,9 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_j0f_huge_range(float ax) | |||
| 2098 | } | 2119 | } |
| 2099 | 2120 | ||
| 2100 | float amp = beta * rsqrtf(ax) * 0.79788458347320556640625f; | 2121 | float amp = beta * rsqrtf(ax) * 0.79788458347320556640625f; |
| 2101 | - return amp * cosf(r - 0.785398185253143310546875f); | 2122 | + // PH leaves r in [0, 2*pi), so r - pi/4 is deep inside cosf's Cody-Waite |
| 2123 | + // range. | ||
| 2124 | + return amp * __internal_cosf_cody_waite(r - 0.785398185253143310546875f); | ||
| 2102 | } | 2125 | } |
| 2103 | 2126 | ||
| 2104 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_j1f_huge_range(float ax) | 2127 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_j1f_huge_range(float ax) |
| @@ -2125,7 +2148,9 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_j1f_huge_range(float ax) | |||
| 2125 | } | 2148 | } |
| 2126 | 2149 | ||
| 2127 | float amp = beta * rsqrtf(ax) * 0.79788458347320556640625f; | 2150 | float amp = beta * rsqrtf(ax) * 0.79788458347320556640625f; |
| 2128 | - return amp * cosf(r - 2.35619449615478515625f); | 2151 | + // cos(r - 3pi/4). PH leaves r in [0, 2*pi), so r - 3pi/4 is deep inside |
| 2152 | + // cosf's Cody-Waite range. | ||
| 2153 | + return amp * __internal_cosf_cody_waite(r - 2.35619449615478515625f); | ||
| 2129 | } | 2154 | } |
| 2130 | 2155 | ||
| 2131 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_j0f_less8(float x) | 2156 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_j0f_less8(float x) |
| @@ -2214,8 +2239,9 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_y0f_huge_range(float ax) | |||
| 2214 | } | 2239 | } |
| 2215 | 2240 | ||
| 2216 | float amp = beta * rsqrtf(ax) * 0.79788458347320556640625f; | 2241 | float amp = beta * rsqrtf(ax) * 0.79788458347320556640625f; |
| 2217 | - // cos(r - 3pi/4) = sin(r - pi/4) | 2242 | + // cos(r - 3pi/4) = sin(r - pi/4). PH leaves r in [0, 2*pi), so r - 3pi/4 is |
| 2218 | - return amp * cosf(r - 2.35619449615478515625f); | 2243 | + // deep inside cosf's Cody-Waite range. |
| 2244 | + return amp * __internal_cosf_cody_waite(r - 2.35619449615478515625f); | ||
| 2219 | } | 2245 | } |
| 2220 | 2246 | ||
| 2221 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_y1f_middle_range(float x) | 2247 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_y1f_middle_range(float x) |
| @@ -2253,7 +2279,7 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_y1f_huge_range(float ax) | |||
| 2253 | 2279 | ||
| 2254 | float amp = beta * rsqrtf(ax) * 0.79788458347320556640625f; | 2280 | float amp = beta * rsqrtf(ax) * 0.79788458347320556640625f; |
| 2255 | // cos(r - 5pi/4) = sin(r - 3pi/4) = -cos(r - pi/4) | 2281 | // cos(r - 5pi/4) = sin(r - 3pi/4) = -cos(r - pi/4) |
| 2256 | - return amp * cosf(r - 3.9269907474517822265625f); | 2282 | + return amp * __internal_cosf_cody_waite(r - 3.9269907474517822265625f); |
| 2257 | } | 2283 | } |
| 2258 | 2284 | ||
| 2259 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_jn_yn_asymptotic_bessel_amplitude(int n, float x, int index) | 2285 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float __internal_jn_yn_asymptotic_bessel_amplitude(int n, float x, int index) |
| @@ -2560,11 +2586,7 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float j1f(float x) | |||
| 2560 | } | 2586 | } |
| 2561 | float ax = __fabsf(x); | 2587 | float ax = __fabsf(x); |
| 2562 | if (__isinf(ax)) { | 2588 | if (__isinf(ax)) { |
| 2563 | - if (x < 0.0f) { | 2589 | + return __internal_with_sign_bit(0.0f, x); // +inf -> +0.0, -inf -> -0.0 |
| 2564 | - uint32_t neg_zero_bits = 0x80000000U; | ||
| 2565 | - return reinterpret_cast<float&>(neg_zero_bits); | ||
| 2566 | - } | ||
| 2567 | - return 0.0f; | ||
| 2568 | } | 2590 | } |
L | |||
| 2569 | float y; | 2591 | float y; |
| 2570 | if (ax <= 8.0f) { | 2592 | if (ax <= 8.0f) { |
| @@ -2589,7 +2611,7 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float jnf(int n, float x) | |||
| 2589 | return j1f(x); | 2611 | return j1f(x); |
| 2590 | } | 2612 | } |
| 2591 | if (n < 0 || isnan(x)) { | 2613 | if (n < 0 || isnan(x)) { |
| 2592 | - return ASCRT_INF_F / ASCRT_INF_F; | 2614 | + return ASCRT_NAN_F; |
| 2593 | } | 2615 | } |
| 2594 | if (x == 0) { | 2616 | if (x == 0) { |
| 2595 | return 0; | 2617 | return 0; |
| @@ -2620,7 +2642,7 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float y0f(float x) | |||
| 2620 | return x; | 2642 | return x; |
| 2621 | } | 2643 | } |
| 2622 | if (x < 0.0f) { | 2644 | if (x < 0.0f) { |
| 2623 | - return ASCRT_INF_F / ASCRT_INF_F; | 2645 | + return ASCRT_NAN_F; |
| 2624 | } | 2646 | } |
| 2625 | if (x == 0.0f) { | 2647 | if (x == 0.0f) { |
| 2626 | return -ASCRT_INF_F; | 2648 | return -ASCRT_INF_F; |
| @@ -2649,7 +2671,7 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float y1f(float x) | |||
| 2649 | return x; | 2671 | return x; |
| 2650 | } | 2672 | } |
| 2651 | if (x < 0.0f) { | 2673 | if (x < 0.0f) { |
| 2652 | - return ASCRT_INF_F / ASCRT_INF_F; | 2674 | + return ASCRT_NAN_F; |
| 2653 | } | 2675 | } |
| 2654 | if (x == 0.0f) { | 2676 | if (x == 0.0f) { |
| 2655 | return -ASCRT_INF_F; | 2677 | return -ASCRT_INF_F; |
| @@ -2678,15 +2700,6 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float y1f(float x) | |||
| 2678 | 2700 | ||
| 2679 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float ynf(int n, float x) | 2701 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float ynf(int n, float x) |
| 2680 | { | 2702 | { |
| 2681 | - if (n < 0 || x < 0 || isnan(x)) { | ||
| 2682 | - return ASCRT_INF_F / ASCRT_INF_F; | ||
| 2683 | - } | ||
| 2684 | - if (x == 0) { | ||
| 2685 | - return -ASCRT_INF_F; | ||
| 2686 | - } | ||
| 2687 | - if (isinf(x)) { | ||
| 2688 | - return 0; | ||
| 2689 | - } | ||
| 2690 | if (n == 0) { | 2703 | if (n == 0) { |
| 2691 | return y0f(x); | 2704 | return y0f(x); |
| 2692 | } | 2705 | } |
| @@ -2694,10 +2707,23 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float ynf(int n, float x) | |||
| 2694 | return y1f(x); | 2707 | return y1f(x); |
| 2695 | } | 2708 | } |
| 2696 | 2709 | ||
| 2710 | + float res; | ||
| 2697 | if (x > n * 10.0f) { | 2711 | if (x > n * 10.0f) { |
| 2698 | - return __internal_yn_asymptotic(n, x); // large x asymptotic | 2712 | + res = __internal_yn_asymptotic(n, x); // large x asymptotic |
| 2713 | + } else { | ||
| 2714 | + res = __internal_yn_recurrence(n, x); // small/medium x forward recurrence | ||
| 2699 | } | 2715 | } |
| 2700 | - return __internal_yn_recurrence(n, x); // small/medium x forward recurrence | 2716 | + |
| 2717 | + if (x == 0.0f) { | ||
| 2718 | + res = -ASCRT_INF_F; | ||
| 2719 | + } | ||
| 2720 | + if (isinf(x)) { | ||
| 2721 | + res = 0.0f; | ||
| 2722 | + } | ||
| 2723 | + if (n < 0 || x < 0.0f || isnan(x)) { | ||
| 2724 | + res = ASCRT_NAN_F; | ||
| 2725 | + } | ||
| 2726 | + return res; | ||
| 2701 | } | 2727 | } |
| 2702 | 2728 | ||
| 2703 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline long int labs(long int x) { return abs(x); } | 2729 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline long int labs(long int x) { return abs(x); } |
| @@ -3681,10 +3707,10 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ static inline int32_t __asc_float_mantissa_expone | |||
| 3681 | // Subnormals are treated as exponent -149, normals as unbiased exponent minus 23. | 3707 | // Subnormals are treated as exponent -149, normals as unbiased exponent minus 23. |
| 3682 | const uint32_t exponent_bits = bits & ASCRT_EXP_BIT_FLOAT_U; | 3708 | const uint32_t exponent_bits = bits & ASCRT_EXP_BIT_FLOAT_U; |
| 3683 | // aligning the biased exponent field to the integer mantissa (value = mantissa * 2^exp) | 3709 | // aligning the biased exponent field to the integer mantissa (value = mantissa * 2^exp) |
| 3684 | - return exponent_bits == 0U ? | 3710 | + //-149: min subnormal exp |
| 3685 | - -149 : | 3711 | + return exponent_bits == 0U ? -149 : |
| 3686 | - (static_cast<int32_t>(exponent_bits >> __internal_fp32_exponent_shift) - | 3712 | + (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 | 3713 | + 150); // -150: bias + mantissa shift |
| 3688 | } | 3714 | } |
| 3689 | 3715 | ||
| 3690 | __SIMT_DEVICE_FUNCTIONS_DECL__ static inline int32_t __asc_uint_floor_log2(uint32_t x) | 3716 | __SIMT_DEVICE_FUNCTIONS_DECL__ static inline int32_t __asc_uint_floor_log2(uint32_t x) |
| @@ -3732,8 +3758,8 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ static inline float __asc_make_float_from_mantiss | |||
| 3732 | 3758 | ||
| 3733 | uint32_t bits = is_normal ? normal_bits : subnormal_bits; | 3759 | uint32_t bits = is_normal ? normal_bits : subnormal_bits; |
| 3734 | bits = is_subnormal ? subnormal_bits : bits; | 3760 | bits = is_subnormal ? subnormal_bits : bits; |
| 3735 | - bits = (result_exponent > 127) ? ASCRT_INF_U : bits; // beyond max normal exponent overflows to inf | 3761 | + bits = (result_exponent > 127) ? ASCRT_INF_U : bits; // 127: max normal exp |
| 3736 | - bits = (mantissa == 0U || result_exponent < -149) ? 0U : bits; // zero mantissa or below min subnormal | 3762 | + bits = (mantissa == 0U || result_exponent < -149) ? 0U : bits; //-149: min subnormal exp |
| 3737 | return __uint_as_float(bits); | 3763 | return __uint_as_float(bits); |
| 3738 | } | 3764 | } |
| 3739 | 3765 | ||
| @@ -4268,58 +4294,41 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float rnorm3df(float a, float b, float c) | |||
| 4268 | */ | 4294 | */ |
| 4269 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float rnorm4df(float a, float b, float c, float d) | 4295 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float rnorm4df(float a, float b, float c, float d) |
| 4270 | { | 4296 | { |
| 4271 | - // Shared scaling constants used to keep the squared sum numerically safe. | ||
| 4272 | - constexpr float rsqrt_subnormal_scale = 4096.0f; // 2^12, compensates the temporary scaling after rsqrt refinement. | ||
| 4273 | - | ||
| 4274 | // Work with absolute values so the norm depends only on magnitudes. | 4297 | // Work with absolute values so the norm depends only on magnitudes. |
| 4275 | const float abs_a = fabsf(a); | 4298 | const float abs_a = fabsf(a); |
| 4276 | const float abs_b = fabsf(b); | 4299 | const float abs_b = fabsf(b); |
| 4277 | const float abs_c = fabsf(c); | 4300 | const float abs_c = fabsf(c); |
| 4278 | const float abs_d = fabsf(d); | 4301 | const float abs_d = fabsf(d); |
| 4302 | + | ||
| 4279 | if (isinf(abs_a) || isinf(abs_b) || isinf(abs_c) || isinf(abs_d)) { | 4303 | if (isinf(abs_a) || isinf(abs_b) || isinf(abs_c) || isinf(abs_d)) { |
| 4280 | return 0.0f; | 4304 | return 0.0f; |
| 4281 | } | 4305 | } |
| 4282 | if (isnan(abs_a) || isnan(abs_b) || isnan(abs_c) || isnan(abs_d)) { | 4306 | if (isnan(abs_a) || isnan(abs_b) || isnan(abs_c) || isnan(abs_d)) { |
| 4283 | - return ASCRT_INF_F / ASCRT_INF_F; | 4307 | + return ASCRT_NAN_F; |
| 4284 | } | 4308 | } |
| 4285 | 4309 | ||
| 4286 | - // Sort magnitudes so the largest term determines the scaling exponent. | 4310 | + const float max_ab = (abs_a > abs_b) ? abs_a : abs_b; |
| 4287 | - const float max_ab = fmaxf(abs_a, abs_b); | 4311 | + const float max_cd = (abs_c > abs_d) ? abs_c : abs_d; |
| 4288 | - const float min_ab = fminf(abs_a, abs_b); | 4312 | + const float max_abs = (max_ab > max_cd) ? max_ab : max_cd; |
| 4289 | - const float max_abc = fmaxf(abs_c, max_ab); | ||
| 4290 | - const float max_abs = fmaxf(abs_d, max_abc); | ||
| 4291 | - const float second_abs = fminf(abs_d, max_abc); | ||
| 4292 | - const float third_abs = fminf(abs_c, max_ab); | ||
| 4293 | const uint32_t max_bits = __float_as_uint(max_abs); | 4313 | const uint32_t max_bits = __float_as_uint(max_abs); |
| 4294 | const float scale = | 4314 | const float scale = |
| 4295 | __uint_as_float(__internal_norm_scale_anchor_bits - (max_bits & __internal_norm_exponent_bucket_mask)); | 4315 | __uint_as_float(__internal_norm_scale_anchor_bits - (max_bits & __internal_norm_exponent_bucket_mask)); |
| 4296 | 4316 | ||
| 4297 | - // Accumulate the scaled sum of squares in magnitude order. | 4317 | + // Accumulate the scaled sum of squares directly in input order. |
| 4298 | - const float scaled_min = min_ab * scale; | 4318 | + const float t_a = abs_a * scale; |
| 4299 | - const float scaled_second = second_abs * scale; | 4319 | + const float t_b = abs_b * scale; |
| 4300 | - const float scaled_third = third_abs * scale; | 4320 | + const float t_c = abs_c * scale; |
| 4301 | - const float scaled_max = max_abs * scale; | 4321 | + const float t_d = abs_d * scale; |
| 4302 | - float square_sum = scaled_second * scaled_second; | 4322 | + float square_sum = fmaf(t_a, t_a, fmaf(t_b, t_b, fmaf(t_c, t_c, t_d * t_d))); |
| 4303 | - square_sum = fmaf(scaled_third, scaled_third, square_sum); | ||
| 4304 | - square_sum = fmaf(scaled_min, scaled_min, square_sum); | ||
| 4305 | - square_sum = fmaf(scaled_max, scaled_max, square_sum); | ||
| 4306 | 4323 | ||
| 4307 | - // Subnormal squared sums are temporarily lifted before rsqrtf. | ||
| 4308 | - const bool is_subnormal_square_sum = fabsf(square_sum) < __internal_subnormal_boundary; | ||
| 4309 | - if (is_subnormal_square_sum) { | ||
| 4310 | - square_sum *= __internal_fp32_scale_2p24; | ||
| 4311 | - } | ||
| 4312 | // Zero norm means the reciprocal norm is infinite. | 4324 | // Zero norm means the reciprocal norm is infinite. |
| 4313 | if (square_sum == 0.0f) { | 4325 | if (square_sum == 0.0f) { |
| 4314 | return ASCRT_INF_F; | 4326 | return ASCRT_INF_F; |
| 4315 | } | 4327 | } |
| 4316 | 4328 | ||
| 4317 | // Compute an approximate reciprocal norm and refine it with one Newton step. | 4329 | // Compute an approximate reciprocal norm and refine it with one Newton step. |
| 4318 | - float inv_norm = rsqrtf(square_sum); | 4330 | + float inv_norm = 1.0f / __sqrtf(square_sum); |
| 4319 | inv_norm = inv_norm * fmaf(-0.5f * square_sum, inv_norm * inv_norm, 1.5f); | 4331 | inv_norm = inv_norm * fmaf(-0.5f * square_sum, inv_norm * inv_norm, 1.5f); |
| 4320 | - if (is_subnormal_square_sum) { | ||
| 4321 | - inv_norm *= rsqrt_subnormal_scale; | ||
| 4322 | - } | ||
| 4323 | // Restore the outer scale factor. | 4332 | // Restore the outer scale factor. |
| 4324 | return scale * inv_norm; | 4333 | return scale * inv_norm; |
| 4325 | } | 4334 | } |
| @@ -4669,7 +4678,7 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float tgammaf(float x) | |||
| 4669 | // Integer inputs on the negative side are poles of gamma. | 4678 | // Integer inputs on the negative side are poles of gamma. |
| 4670 | const float truncated_x = truncf(x); | 4679 | const float truncated_x = truncf(x); |
| 4671 | if (truncated_x == x) { | 4680 | if (truncated_x == x) { |
| 4672 | - result = x < 0.0f ? ASCRT_INF_F / ASCRT_INF_F : result; | 4681 | + result = x < 0.0f ? ASCRT_NAN_F : result; |
| 4673 | } | 4682 | } |
| 4674 | 4683 | ||
| 4675 | // For large negative non-integers, odd integer truncations force a sign correction to zero. | 4684 | // For large negative non-integers, odd integer truncations force a sign correction to zero. |
| @@ -4688,7 +4697,7 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float tgammaf(float x) | |||
| 4688 | result = ASCRT_INF_F; | 4697 | result = ASCRT_INF_F; |
| 4689 | } | 4698 | } |
| 4690 | if (x == -ASCRT_INF_F) { | 4699 | if (x == -ASCRT_INF_F) { |
| 4691 | - result = ASCRT_INF_F / ASCRT_INF_F; | 4700 | + result = ASCRT_NAN_F; |
| 4692 | } | 4701 | } |
| 4693 | if (x == 0.0f) { | 4702 | if (x == 0.0f) { |
| 4694 | result = 1.0f / x; | 4703 | result = 1.0f / x; |
| @@ -5191,91 +5200,85 @@ __SIMT_DEVICE_FUNCTIONS_DECL__ inline float fmodf(float x, float y) | |||
| 5191 | ux &= 0x7FFFFFFFU; | 5200 | ux &= 0x7FFFFFFFU; |
| 5192 | uy &= 0x7FFFFFFFU; | 5201 | uy &= 0x7FFFFFFFU; |
| 5193 | 5202 | ||
| 5194 | - // Invalid divisor, NaN, or infinity cases follow IEEE-style NaN propagation. | 5203 | + // IEEE special-value handling: a zero divisor, an infinite/NaN dividend, |
| 5204 | + // or a NaN divisor yields NaN. | ||
| 5195 | if (uy == 0U || ux >= 0x7F800000U || uy > 0x7F800000U) { | 5205 | if (uy == 0U || ux >= 0x7F800000U || uy > 0x7F800000U) { |
| 5196 | - return ASCRT_INF_F / ASCRT_INF_F; | 5206 | + return ASCRT_NAN_F; |
| 5197 | - } | ||
| 5198 | - // If |x| < |y|, the remainder is x itself. | ||
| 5199 | - if (ux < uy) { | ||
| 5200 | - return x; | ||
| 5201 | - } | ||
| 5202 | - // If |x| == |y|, the remainder is signed zero with the sign of x. | ||
| 5203 | - if (ux == uy) { | ||
| 5204 | - return reinterpret_cast<float&>(sx); | ||
| 5205 | } | 5207 | } |
| 5206 | 5208 | ||
| 5207 | - // Extract exponents and build mantissas in normalized form. | 5209 | + // From here on 1 <= uy <= ux < 0x7F800000, i.e. the plain finite domain: |
| 5208 | - int32_t ex = static_cast<int32_t>(ux >> 23); | 5210 | + // - |x| < |y| implies ex <= ey, so the alignment loop is skipped |
| 5209 | - int32_t ey = static_cast<int32_t>(uy >> 23); | 5211 | + // outright and the computation stays bounded. |
| 5210 | - uint32_t mx = ux; | 5212 | + // - |x| == |y| implies mx == my and ex == ey, so the post-loop |
| 5211 | - uint32_t my = uy; | 5213 | + // subtraction yields exactly zero and the exact-select produces the |
| 5214 | + // signed zero of x. | ||
| 5215 | + const bool x_normal = (ux >= 0x00800000U); | ||
| 5216 | + const int32_t x_log2 = __asc_uint_floor_log2(ux != 0U ? ux : 1U); | ||
| 5217 | + const uint32_t x_shift = x_normal ? 0U : static_cast<uint32_t>(__internal_fp32_exponent_shift - x_log2); | ||
| 5218 | + uint32_t mx = x_normal ? ((ux & 0x007FFFFFU) | 0x00800000U) : (ux << x_shift); | ||
| 5219 | + int32_t ex = x_normal ? static_cast<int32_t>(ux >> __internal_fp32_exponent_shift) : (x_log2 - 22); | ||
| 5212 | 5220 | ||
| 5213 | - // Normalize subnormal x so the hidden leading bit is restored. | 5221 | + // Same normalization for the divisor; uy >= 1 holds after the invalid |
| 5214 | - if (ex == 0) { | 5222 | + // check above, so the log2 input is never zero. |
| 5215 | - int32_t i = 0; | 5223 | + const bool y_normal = (uy >= 0x00800000U); |
| 5216 | - mx = ux; | 5224 | + const int32_t y_log2 = __asc_uint_floor_log2(uy); |
| 5217 | - while ((mx & 0x00800000U) == 0U) { | 5225 | + const uint32_t y_shift = y_normal ? 0U : static_cast<uint32_t>(__internal_fp32_exponent_shift - y_log2); |
| 5218 | - mx <<= 1; | 5226 | + uint32_t my = y_normal ? ((uy & 0x007FFFFFU) | 0x00800000U) : (uy << y_shift); |
| 5219 | - i++; | 5227 | + int32_t ey = y_normal ? static_cast<int32_t>(uy >> __internal_fp32_exponent_shift) : (y_log2 - 22); |
| 5220 | - } | ||
| 5221 | - ex = 1 - i; | ||
| 5222 | - } else { | ||
| 5223 | - mx = (ux & 0x007FFFFFU) | 0x00800000U; | ||
| 5224 | - } | ||
| 5225 | 5228 | ||
| 5226 | - // Normalize subnormal y in the same way. | 5229 | + // Align x to y: for each binary exponent step above ey, subtract the |
| 5227 | - if (ey == 0) { | 5230 | + // divisor mantissa when it fits (top bit of the difference stays clear) |
| 5228 | - int32_t i = 0; | 5231 | + // and shift the dividend left one place. A remainder that reaches zero |
| 5229 | - my = uy; | 5232 | + // keeps iterating with mx == 0 — every later subtraction borrows, so the |
| 5230 | - while ((my & 0x00800000U) == 0U) { | 5233 | + // value no longer changes — while the exact flag records the zero. |
| 5231 | - my <<= 1; | 5234 | + bool exact = false; |
| 5232 | - i++; | ||
| 5233 | - } | ||
| 5234 | - ey = 1 - i; | ||
| 5235 | - } else { | ||
| 5236 | - my = (uy & 0x007FFFFFU) | 0x00800000U; | ||
| 5237 | - } | ||
| 5238 | - | ||
| 5239 | - // Align x to y by repeatedly subtracting the divisor mantissa while | ||
| 5240 | - // shifting the dividend down one binary exponent step at a time. | ||
| 5241 | while (ex > ey) { | 5235 | while (ex > ey) { |
| 5242 | - uint32_t d = mx - my; | 5236 | + const uint32_t d = mx - my; |
| 5243 | - if ((d & 0x80000000U) == 0U) { | 5237 | + const bool fits = (d & 0x80000000U) == 0U; |
| 5244 | - if (d == 0U) { | 5238 | + mx = fits ? d : mx; |
| 5245 | - return reinterpret_cast<float&>(sx); | 5239 | + exact |= fits && (d == 0U); |
| 5246 | - } | ||
| 5247 | - mx = d; | ||
| 5248 | - } | ||
| 5249 | mx <<= 1; | 5240 | mx <<= 1; |
| 5250 | ex--; | 5241 | ex--; |
| 5251 | } | 5242 | } |
| 5252 | 5243 | ||
| 5253 | // Finish the aligned subtraction once exponents match. | 5244 | // Finish the aligned subtraction once exponents match. |
| 5254 | - uint32_t d = mx - my; | 5245 | + const uint32_t d = mx - my; |
| 5255 | - if ((d & 0x80000000U) == 0U) { | 5246 | + const bool fits = (d & 0x80000000U) == 0U; |
| 5256 | - if (d == 0U) { | 5247 | + exact |= fits && (d == 0U); |
| 5257 | - return reinterpret_cast<float&>(sx); | 5248 | + mx = fits ? d : mx; |
| 5258 | - } | ||
| 5259 | - mx = d; | ||
| 5260 | - } | ||
| 5261 | 5249 | ||
| 5262 | - // Renormalize the remainder mantissa so it can be packed back to float. | 5250 | + // Renormalize the remainder mantissa: floor(log2) gives the position of |
| 5263 | - while ((mx & 0x00800000U) == 0U) { | 5251 | + // the leading bit, so one shift brings it back to bit 23. A zero |
| 5264 | - mx <<= 1; | 5252 | + // significand (exact-zero remainder) is guarded to 1 so the shift count |
| 5265 | - ex--; | 5253 | + // stays in [0, 23]. |
| 5266 | - } | 5254 | + const int32_t r_log2 = __asc_uint_floor_log2(mx != 0U ? mx : 1U); |
| 5255 | + const int32_t r_shift = __internal_fp32_exponent_shift - r_log2; | ||
| 5256 | + mx <<= static_cast<uint32_t>(r_shift); | ||
| 5257 | + ex -= r_shift; | ||
| 5267 | 5258 | ||
| 5268 | - // Rebuild the float result, restore the sign of x, and return the remainder. | 5259 | + // Rebuild the float result: pack the renormalized mantissa with the |
| 5260 | + // adjusted exponent, or shift down when the exponent is subnormal. A | ||
| 5261 | + // shift of 32 or more yields zero. | ||
| 5269 | uint32_t out; | 5262 | uint32_t out; |
| 5270 | if (ex > 0) { | 5263 | if (ex > 0) { |
| 5271 | mx -= 0x00800000U; | 5264 | mx -= 0x00800000U; |
| 5272 | - out = mx | (static_cast<uint32_t>(ex) << 23); // 23: fp32 exponent field offset. | 5265 | + out = mx | (static_cast<uint32_t>(ex) << __internal_fp32_exponent_shift); |
| 5273 | } else { | 5266 | } else { |
| 5274 | - mx >>= static_cast<uint32_t>(1 - ex); | 5267 | + const int32_t sub_shift = 1 - ex; |
| 5275 | - out = mx; | 5268 | + out = (sub_shift < 32) ? (mx >> static_cast<uint32_t>(sub_shift)) : 0U; // 32: the width of uint32_t. |
| 5276 | } | 5269 | } |
| 5277 | out |= sx; | 5270 | out |= sx; |
| 5278 | - return reinterpret_cast<float&>(out); | 5271 | + float res = reinterpret_cast<float&>(out); |
| 5272 | + if (exact) { | ||
| 5273 | + // The remainder is exactly zero: return the signed zero of x. | ||
| 5274 | + res = reinterpret_cast<float&>(sx); | ||
| 5275 | + } | ||
| 5276 | + if (ux < uy) { | ||
| 5277 | + // |x| < |y|: the remainder is x itself. This check runs last so it | ||
| 5278 | + // wins over the exact-select; x == 0 also reaches it. | ||
| 5279 | + res = x; | ||
| 5280 | + } | ||
| 5281 | + return res; | ||
| 5279 | } | 5282 | } |
| 5280 | 5283 | ||
| 5281 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float acoshf(float x) | 5284 | __SIMT_DEVICE_FUNCTIONS_DECL__ inline float acoshf(float x) |
This code comment section is from+2581To+2590
这块放到最下面试试