On 1/28/26 06:50, Richard Henderson wrote:
On 1/28/26 02:31, Ilya Leoshkevich wrote:
+ FloatParts128 a128, b128, *m128, m128_buf, n128, *q128,
q128_buf,
+ r128, *r128_precise;
+ int float_exception_flags = 0;
+ bool is_q128_smallish;
+ uint32_t r_flags;
+ int saved_flags;
+
+ /* Compute precise quotient */
+ parts_float_to_float_widen(&a128, a, status);
+ parts_float_to_float_widen(&b128, b, status);
+ q128_buf = a128;
+ q128 = parts_div(&q128_buf, &b128, status);
Why do you need FloatParts128?
You can see an inexact result from float64 with just FloatParts64.
C.f. soft_f64_div.
I thought I needed this extra precision, but turns out that in reality I
needed it to mask an issue with parts_round_to_int() - see below.
+
+ /* Final or partial case? */
+ is_q128_smallish = q128->exp < (fmt->frac_size + 1);
+
+ /*
+ * Final quotient is rounded using final-quotient-rounding
method, and
+ * partial quotient is rounded toward zero.
+ *
+ * Rounding of partial quotient may be inexact. This is the
whole point
+ * of distinguishing partial quotients, so ignore the
exception.
+ */
+ n128 = *q128;
+ saved_flags = status->float_exception_flags;
+ parts_round_to_int(&n128,
+ is_q128_smallish ?
final_quotient_rounding_mode :
+ float_round_to_zero,
+ 0, status, &float128_params);
float128_params is definitely wrong. The rounding is supposed to be
to the target format.
Hmm, yes, nothing in this code is about float128, that can't be right.
With some testcases I hit this condition in parts_round_to_int_normal():
if (a->exp >= frac_size) {
/* All integral */
return false;
}
which makes it a no-op.
I think the code assumes that FloatParts have just been unpacked and all
low fraction bits are zero, which is not the case for quotient here.
+ float_exception_flags = saved_flags;
+ parts_s390_precision_round_normal(&n128, fmt);
parts_round_to_int_normal already takes a rounding mode and frac_size.
It also returns whether or not the rounding was exact.
This appears to be trying to reinvent the wheel.
Apparently I use this to paper over the above issue, which is not great.
I guess improving parts_round_to_int() to work with non-zero low
fraction bits would be better.
What do you think?
There's still "... the two integers closest to this precise quotient
cannot be both be represented exactly in the precision of the quotent
..." to contend with. Is that your "smallish" test? If so, the
comment could use some improvement.
Ok.
+
+ /* Compute precise remainder */
+ m128_buf = b128;
+ m128 = parts_mul(&m128_buf, &n128, status);
+ m128->sign = !m128->sign;
+ status->float_exception_flags = 0;
+ r128_precise = parts_addsub(m128, &a128, status, false);
Surely parts_muladd_scalbn (with scale 0).
This works and is much shorter, thanks.
r~
