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}