github.com/afumu/libc@v0.0.6/musl/src/math/i386/sqrt.c (about)

     1  #include "libm.h"
     2  
     3  double sqrt(double x)
     4  {
     5  	union ldshape ux;
     6  	unsigned fpsr;
     7  	__asm__ ("fsqrt; fnstsw %%ax": "=t"(ux.f), "=a"(fpsr) : "0"(x));
     8  	if ((ux.i.m & 0x7ff) != 0x400)
     9  		return (double)ux.f;
    10  	/* Rounding to double would have encountered an exact halfway case.
    11  	   Adjust mantissa downwards if fsqrt rounded up, else upwards.
    12  	   (result of fsqrt could not have been exact) */
    13  	ux.i.m ^= (fpsr & 0x200) + 0x300;
    14  	return (double)ux.f;
    15  }