Re: [PATCH 2/2] MIPS: Hard-float rounding instructions support

[Xi Ruoyao <xry111@xry111.site>]

On Wed, 2024-01-31 at 11:28 +0800, Junxian Zhu wrote:
> 在 2023/12/26 16:29, Xi Ruoyao 写道:
> > On Tue, 2023-12-26 at 10:37 +0800, Junxian Zhu wrote:
> > > 在 2023/12/25 18:51, Xi Ruoyao 写道:
> > > > On Mon, 2023-12-25 at 18:35 +0800, Junxian Zhu wrote:
> > > > 
> > > > /* snip */
> > > > 
> > > > > +/*
> > > > > + * ceil(x)
> > > > > + * Return x rounded toward -inf to integral value
> > > > > + * Method:
> > > > > + *	Bit twiddling.
> > > > > + */
> > > > > +
> > > > > +#if ((__mips_fpr == 64) && (__mips_hard_float == 1) &&
> > > > > ((__mips == 32 && __mips_isa_rev > 1) || __mips == 64))
> > > > > +#include <sys/regdef.h>
> > > > > +#include <sysdep.h>
> > > > > +#include <libm-alias-double.h>
> > > > > +
> > > > > +ENTRY(__ceil)
> > > > > +	.set push
> > > > > +	.set noreorder
> > > > > +	.set noat
> > > > > +# $f0=ret, $f12=double, a0=int64/int32_h, a1=int32_l,
> > > > > a2=sign, a3=exp
> > > > > +#if __mips == 64
> > > > > +	dmfc1   a0, $f12 # assign int64
> > > > > +#else
> > > > > +	mfhc1   a0, $f12 # assign int64
> > > > > +#endif
> > > > > +	cfc1    t0, $f26
> > > > > +	ceil.l.d    $f0, $f12
> > > > No, C23 does not allow this function to raise an INEXACT
> > > > exception, but
> > > > ceil.l.d will do so.
> > > > 
> > > > Such optimizations should be performed in GCC which can be
> > > > controlled by
> > > > the programmer with -std=c23 and/or -f[no-]fp-int-builtin-
> > > > inexact, not
> > > > in Glibc where we cannot know if the programmer wants to deviate
> > > > from
> > > > C23.
> > > The cfc1 instruction will backup float point exception status
> > > before
> > > running ceil.l.d, and the following ctc1 will restore float point
> > > exception status to avoid INEXACT exception raised by ceil.l.d.
> > > It's the
> > > same way like what have been done in s_ceil.S for i386.
> > Still incorrect because when the Enable field of FCSR contains
> > INEXACT a
> > SIGFPE will be immediately delivered and there is no way to
> > recover.  A
> > demonstration:
> > 
> > #define _GNU_SOURCE
> > #include <stdio.h>
> > #include <fenv.h>
> > 
> > int main()
> > {
> >    printf("%d\n", feenableexcept(FE_INEXACT));
> > 
> >    double data = 114.514;
> >    long control;
> >    asm("cfc1\t%1,$f26\n\t"
> >        "ceil.l.d\t%0,%0\n\t"
> >        "cvt.d.l\t%0,%0\n\t"
> >        "ctc1\t%1,$f26": "+f"(data), "=r"(control));
> >    printf("%.15f\n", data);
> >    return 0;
> > }
> > 
> > On i386 the fnstenv instruction also masks out all the FP exceptions
> > so
> > this is not a problem.  See commit 26b0bf96000a.
> 
> I made some similar code to test and found that sqrt.fmt on MIPS will 
> raise INEXACT exception too.
> 
> So should we mask all the FP exceptions while compiling MIPS 
> SQRT{F}_BUILTIN in GCC ? (for C23)

Why?

    volatile double x = 2;
    sqrt(x); 

should definitely raise an INEXACT exception.  Unless explicitly stated
otherwise, all operations producing a result different than the
"infinite precision operation" should do it.  For the conversion
operations they are just "explicitly stated otherwise."

But if sqrt.fmt is raising it for

    volatile double x = 4;
    sqrt(x); 

it's obviously wrong, and I'd say such an exception is nonsense even
regardless of the standard so we should just not use this instruction.

Something not so obvious is: the standard explicitly says

    volatile double x = __builtin_inf();
    sqrt(x); 

should **not** raise an INEXACT exception.  So if you are raising
INEXACT for this you should only use sqrt.fmt when you can prove the
operand is finite (for example when -ffinite-math-only, or using some
information from the VRP pass).

-- 
Xi Ruoyao <xry111@xry111.site>
School of Aerospace Science and Technology, Xidian University