Remove sqrt from the volatile list

[Keno]

The LLVM IR spec now explicitly says that this intrinsic is required
to be rounded correctly. That means that without fasthmath flag
(which we do not set here), this intrinsic must have the bitwise
correctly rounded answer and should thus not differ between compile
and runtime. If there is still a case where it does differ, that is
likely some other underlying bug that we should fix instead.

Before:

julia> f() = sqrt(2)
f (generic function with 1 method)

julia> @code_typed f()
CodeInfo(
1 ─ %1 = Base.Math.sqrt_llvm(2.0)::Float64
└──      return %1
) => Float64

After:

julia> @code_typed f()
CodeInfo(
1 ─     return 1.4142135623730951
) => Float64

[oscardssmith]

Nice to get this fixed!

[Keno]

my guess is that the common way this could fail is if a libm uses fqsrt without accounting for double rounding. glibc has test vectors for this case at https://sourceware.org/bugzilla/show_bug.cgi?id=14032, so we could add those as a test to see if anybody reports issues

[oscardssmith]

If anyone fails this, it would probably be windows. Their libm isn't great in general.

[N5N3]

Random test passed locally on win64:

for I in 1:10000000
    r = reinterpret(Float64,rand(UInt))
    @test sqrt(abs(r)) === Base.sqrt_llvm(abs(r))
end

Edit: glibc test also passed.
So it should be OK for modern CPU? (with vsqrt support)

[Keno]

Base.sqrt is defined as sqrt_llvm, so unfortunately that test does not say anything. As I said, the interesting case are probably the test vectors in that glibc issue.

[N5N3]

IIUC, Base.sqrt_llvm always calls runtime sqrt like Base.fma_float. (#43530 use such feature to test julia_fma(f)'s precision).
Just for example: on offical 1.7.1 Win64

julia> fma(-1.9369631f13, 2.1513551f-7, -1.7354427f-24)
-4.1670958f6

julia> Base.fma_float(-1.9369631f13, 2.1513551f-7, -1.7354427f-24)
-4.1670955f6

[Keno]

Alright, I have added a test to detect if there's any platform with the fp80 double rounding problem. Otherwise, I think we should try this and see what happens.

[Keno]

Sigh, sqrt test failed on builtkite i686 Linux, so looks like whatever we're using there is buggy.

[Keno]

This might be our fault because it looks like openlibm has the double rounding bug on 32bit: https://github.com/JuliaMath/openlibm/blob/master/i387/e_sqrt.S#L12.

@oscardssmith wanna have a go at fixing this in openlibm?

That said, we assume a minimum of SSE2, so our minimum supported ISA has sqrtsd - we should probably just use that.

[Keno]

IIUC, Base.sqrt_llvm always calls runtime sqrt like Base.fma_float. (#43530 use such feature to test julia_fma(f)'s precision).

There's three cases that all need to be consistent:

1. The code generated by LLVM for the intrinsic
2. The code LLVM uses for constant folding the intrinsic
3. Our own runtime intrinsic used by the interpreter

Depending on how things are built, these can all pick up different version of the intrinsics.

[oscardssmith]

I'd honestly rather us move away from openlibm. can we just make llvm use the intrinsic?

[Keno]

Actually, looks like in this case LLVM doesn't actually constant fold sqrt, so we'll be ok once we fix runtime intrinsics (and since we require SSE2, it should be fine to do that using the SSE intrinsic). However, just for completeness, I do think we should also fix openlibm. Should be a 5 line change to set the precision flag in the control word appropriately.

[oscardssmith]

That makes sense. (I am, however unlikely to touch the openlibm version).

[Keno]

Alright, do you want to change the runtime intrinsics here and I'll make the patch to openlibm (but we won't block on an openlibm upgrade)?

[Keno]

Openlibm fix is JuliaMath/openlibm#256.

[Keno]

As discussed, for now we won't touch runtime_intrinsic and just pull in the openlibm update to enable this. #43869 tracks just getting rid of the intrinsic entirely to reduce the complexity.