1#include "libm.h" 2 3double 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} 16