fpu: Reorganize partsN(muladd)
Check the likely case of normal product and normal or zero addend first; shift NaN and infinity detection down; end with zero product + addend. Reviewed-by: Peter Maydell <peter.maydell@linaro.org> Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
Richard Henderson committed
Apr 28, 2026 at 10:13 UTC
8b44b423ff039ce9098a9243e21ec32f594d19f9
1 file changed
+68
-82
fpu/softfloat-parts.c.inc
+68
-82
@@ -681,11 +681,47 @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b,
681
FloatPartsN *c,
682
int flags, float_status *s)
683
{
684
- int ab_mask, abc_mask;
685
- FloatPartsW p_widen, c_widen;
684
+ int ab_mask = float_cmask(a->cls) | float_cmask(b->cls);
685
+ int c_mask = float_cmask(c->cls);
686
+ int abc_mask = ab_mask | c_mask;
687
+ bool c_sign = c->sign ^ !!(flags & float_muladd_negate_c);
688
+ bool p_sign = a->sign ^ b->sign ^ !!(flags & float_muladd_negate_product);
689
+
690
+ /*
691
+ * The "likely" case is A and B normal, so that the product is normal,
692
+ * and C normal or zero so that the result is normal.
693
+ */
694
+ int likely_mask = ab_mask | (c_mask & ~float_cmask_zero);
695
+ if (likely(cmask_is_only_normals(likely_mask))) {
696
+ record_denormals_used(abc_mask, s);
697
+
698
+ /* Perform the multiplication step. */
699
+ FloatPartsW p_widen = { .sign = p_sign, .exp = a->exp + b->exp + 1 };
700
+ fracN(mulw)(&p_widen, a, b);
701
+ if (!(p_widen.frac_hi & DECOMPOSED_IMPLICIT_BIT)) {
702
+ fracW(add)(&p_widen, &p_widen, &p_widen);
703
+ p_widen.exp -= 1;
704
+ }
705
+
706
+ /* Perform the addition step. */
707
+ if (!(c_mask & float_cmask_zero)) {
708
+ /* Zero-extend C to less significant bits. */
709
+ FloatPartsW c_widen = { .sign = c_sign, .exp = c->exp };
710
+ fracN(widen)(&c_widen, c);
711
687
- ab_mask = float_cmask(a->cls) | float_cmask(b->cls);
688
- abc_mask = float_cmask(c->cls) | ab_mask;
712
+ if (p_sign == c_sign) {
713
+ partsW(add_normal)(&p_widen, &c_widen);
714
+ } else if (!partsW(sub_normal)(&p_widen, &c_widen)) {
715
+ goto return_sub_zero;
716
+ }
717
+ }
718
+
719
+ /* Narrow with sticky bit, for proper rounding later. */
720
+ fracN(truncjam)(a, &p_widen);
721
+ a->sign = p_widen.sign;
722
+ a->exp = p_widen.exp;
723
+ return a;
724
+ }
725
726
/*
727
* It is implementation-defined whether the cases of (0,inf,qnan)
@@ -698,97 +734,47 @@ static FloatPartsN *partsN(muladd)(FloatPartsN *a, FloatPartsN *b,
734
return a;
735
}
736
701
- if (flags & float_muladd_negate_c) {
702
- c->sign ^= 1;
703
- }
704
-
705
- /* Compute the sign of the product into A. */
706
- a->sign ^= b->sign;
707
- if (flags & float_muladd_negate_product) {
708
- a->sign ^= 1;
737
+ if (unlikely(ab_mask == float_cmask_infzero)) {
738
+ /* Inf * Zero == NaN */
739
+ float_raise(float_flag_invalid | float_flag_invalid_imz, s);
740
+ goto d_nan;
741
}
742
711
- if (unlikely(!cmask_is_only_normals(ab_mask))) {
712
- if (unlikely(ab_mask == float_cmask_infzero)) {
713
- float_raise(float_flag_invalid | float_flag_invalid_imz, s);
743
+ if (unlikely(ab_mask & float_cmask_inf)) {
744
+ if ((c_mask & float_cmask_inf) && p_sign != c_sign) {
745
+ /* Inf - Inf == NaN */
746
+ float_raise(float_flag_invalid | float_flag_invalid_isi, s);
747
goto d_nan;
748
}
716
-
717
- if (ab_mask & float_cmask_inf) {
718
- if (c->cls == float_class_inf && a->sign != c->sign) {
719
- float_raise(float_flag_invalid | float_flag_invalid_isi, s);
720
- goto d_nan;
721
- }
722
- goto return_inf;
723
- }
724
-
725
- g_assert(ab_mask & float_cmask_zero);
726
- if (is_anynorm(c->cls)) {
727
- *a = *c;
728
- goto finish_sign;
729
- }
730
- if (c->cls == float_class_zero) {
731
- if (flags & float_muladd_suppress_add_product_zero) {
732
- a->sign = c->sign;
733
- } else if (a->sign != c->sign) {
734
- goto return_sub_zero;
735
- }
736
- goto return_zero;
737
- }
738
- g_assert(c->cls == float_class_inf);
739
- }
740
-
741
- if (unlikely(c->cls == float_class_inf)) {
742
- a->sign = c->sign;
743
- goto return_inf;
744
- }
745
-
746
- /* Perform the multiplication step. */
747
- p_widen.sign = a->sign;
748
- p_widen.exp = a->exp + b->exp + 1;
749
- fracN(mulw)(&p_widen, a, b);
750
- if (!(p_widen.frac_hi & DECOMPOSED_IMPLICIT_BIT)) {
751
- fracW(add)(&p_widen, &p_widen, &p_widen);
752
- p_widen.exp -= 1;
753
- }
754
-
755
- /* Perform the addition step. */
756
- if (c->cls != float_class_zero) {
757
- /* Zero-extend C to less significant bits. */
758
- fracN(widen)(&c_widen, c);
759
- c_widen.exp = c->exp;
760
-
761
- if (a->sign == c->sign) {
762
- partsW(add_normal)(&p_widen, &c_widen);
763
- } else if (!partsW(sub_normal)(&p_widen, &c_widen)) {
764
- goto return_sub_zero;
765
- }
749
+ /* Inf + C == Inf */
750
+ record_denormals_used(abc_mask, s);
751
+ a->sign = p_sign;
752
+ a->cls = float_class_inf;
753
+ return a;
754
}
755
+ record_denormals_used(abc_mask, s);
756
768
- /* Narrow with sticky bit, for proper rounding later. */
769
- fracN(truncjam)(a, &p_widen);
770
- a->sign = p_widen.sign;
771
- a->exp = p_widen.exp;
757
+ /* Only remaining cases are zero product or inf addend. */
758
+ assert((ab_mask & float_cmask_zero) | (c_mask & float_cmask_inf));
759
773
- finish_sign:
760
/*
775
- * All result types except for "return the default NaN
776
- * because this is an Invalid Operation" go through here;
777
- * this matches the set of cases where we consumed a
778
- * denormal input.
761
+ * P + Inf == Inf, or
762
+ * 0 + C == C,
763
+ * except for 0 - 0, which needs special rounding,
764
+ * except for when we want to suppress this addition step.
765
*/
780
- record_denormals_used(abc_mask, s);
781
- return a;
766
+ if (!(c_mask & float_cmask_zero)
767
+ || p_sign == c_sign
768
+ || (flags & float_muladd_suppress_add_product_zero)) {
769
+ c->sign = c_sign;
770
+ return c;
771
+ }
772
773
return_sub_zero:
774
+ /* 0 - 0 == -0 for round_down, +0 otherwise. */
775
a->sign = s->float_rounding_mode == float_round_down;
785
- return_zero:
776
a->cls = float_class_zero;
787
- goto finish_sign;
788
-
789
- return_inf:
790
- a->cls = float_class_inf;
791
- goto finish_sign;
777
+ return a;
778
779
d_nan:
780
*a = partsN(default_nan)(s);