1 | __#include "quadmath-imp.h"__ |

2 | __#include <math.h>__ |

3 | __#include <float.h>__ |

4 | |

5 | **__float128** |

6 | sqrtq (*const* **__float128** x) |

7 | { |

8 | **__float128** y; |

9 | *int* exp; |

10 | |

11 | **if** (isnanq (x) || (isinfq (x) && x > `0`)) |

12 | **return** x; |

13 | |

14 | **if** (x == `0`) |

15 | **return** x; |

16 | |

17 | **if** (x < `0`) |

18 | { |

19 | */* Return NaN with invalid signal. */* |

20 | **return** (x - x) / (x - x); |

21 | } |

22 | |

23 | **if** (x <= DBL_MAX && x >= DBL_MIN) |

24 | { |

25 | */* Use double result as starting point. */* |

26 | y = sqrt ((*double*) x); |

27 | |

28 | */* Two Newton iterations. */* |

29 | y -= `0.5q` * (y - x / y); |

30 | y -= `0.5q` * (y - x / y); |

31 | **return** y; |

32 | } |

33 | |

34 | __#ifdef HAVE_SQRTL__ |

35 | { |

36 | *long* *double* xl = (*long* *double*) x; |

37 | **if** (xl <= LDBL_MAX && xl >= LDBL_MIN) |

38 | { |

39 | */* Use long double result as starting point. */* |

40 | y = (**__float128**) sqrtl (xl); |

41 | |

42 | */* One Newton iteration. */* |

43 | y -= `0.5q` * (y - x / y); |

44 | **return** y; |

45 | } |

46 | } |

47 | __#endif__ |

48 | |

49 | */* If we're outside of the range of C types, we have to compute* |

50 | * the initial guess the hard way. */* |

51 | y = frexpq (x, &exp); |

52 | **if** (exp % `2`) |

53 | y *= `2`, exp--; |

54 | |

55 | y = sqrt (y); |

56 | y = scalbnq (y, exp / `2`); |

57 | |

58 | */* Two Newton iterations. */* |

59 | y -= `0.5q` * (y - x / y); |

60 | y -= `0.5q` * (y - x / y); |

61 | **return** y; |

62 | } |

63 | |

64 | |