From ae2d91698508701c83cab83714d42a1146dccf85 Mon Sep 17 00:00:00 2001 From: Keno Fischer Date: Wed, 19 Jan 2022 18:43:52 -0500 Subject: [PATCH] Correctly round double precision sqrt (#256) As discussed in https://github.com/JuliaLang/julia/pull/43786, openlibm's sqrt function is incorrectly rounded for i387. IEEE requires correct rounding for these functions and LLVM relies on it. Fix that by setting the precision in the FPU control word (see e.g. e_ceil.S for similar FPU modifications). --- i387/e_sqrt.S | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/i387/e_sqrt.S b/i387/e_sqrt.S index 7c36dfe5..5cd45b09 100644 --- a/i387/e_sqrt.S +++ b/i387/e_sqrt.S @@ -8,12 +8,29 @@ //__FBSDID("$FreeBSD: src/lib/msun/i387/e_sqrt.S,v 1.10 2011/01/07 16:13:12 kib Exp $") ENTRY(sqrt) - fldl 4(%esp) + pushl %ebp + movl %esp,%ebp + subl $8,%esp + + fstcw -4(%ebp) /* store fpu control word */ + movw -4(%ebp),%dx + andw $0xfeff,%dx /* Set precision field to 64 bits (53 bit mantissa). + We assume it's set to 0b11 (extended precision), + so zeroing out the low bit of the precision field, + will correctly set the precision */ + movw %dx,-8(%ebp) + fldcw -8(%ebp) /* load modfied control word */ + + fldl 8(%ebp) fsqrt + + fldcw -4(%ebp) /* restore original control word */ + + leave ret END(sqrt) - + /* Enable stack protection */ #if defined(__ELF__) .section .note.GNU-stack,"",%progbits