Fix x86 sqrt rounding (bug 14032)

[Rich Felker]

On Thu, Nov 28, 2013 at 10:34:13AM +1300, Richard Henderson wrote:
> On 11/28/2013 07:06 AM, Joseph S. Myers wrote:
> >  ENTRY(__ieee754_sqrt)
> >  	fldl	4(%esp)
> >  	fsqrt
> > +	fstsw	%ax
> > +	subl	$12, %esp
> > +	cfi_adjust_cfa_offset (12)
> > +	fld	%st
> > +	fstpt	(%esp)
> > +	movl	(%esp), %edx
> > +	andl	$0x7ff, %edx
> > +	cmpl	$0x400, %edx
> > +	jne	1f
> > +	andl	$0x200, %eax
> > +	subl	$0x100, %eax
> > +	subl	%eax, (%esp)
> > +	fstp	%st
> > +	fldt	(%esp)
> > +1:	fstpl	(%esp)
> > +	fldl	(%esp)
> > +	addl	$12,%esp
> > +	cfi_adjust_cfa_offset (-12)
> >  	ret
> 
> Surely this is much more complex than simply setting the rounding precision to
> double before using fsqrt.

Setting the rounding precision is likely to be much slower, and has
issues with signals (signal handlers could be invoked with wrong
rounding precision, resulting in completely wrong results). Also, in
general, setting the rounding precision is not a valid solution to
excess-precision problems, since it doesn't handle denormals
correctly (the exponent still has full ld80 range), but for sqrt it
would work since the result of sqrt can never be denormal.

Rich