1 | */* e_sqrtf.c -- float version of e_sqrt.c.* |

2 | * */* |

3 | |

4 | */** |

5 | * * ====================================================* |

6 | * * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.* |

7 | * ** |

8 | * * Developed at SunPro, a Sun Microsystems, Inc. business.* |

9 | * * Permission to use, copy, modify, and distribute this* |

10 | * * software is freely granted, provided that this notice* |

11 | * * is preserved.* |

12 | * * ====================================================* |

13 | * */* |

14 | |

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

16 | __#include <math_private.h>__ |

17 | __#include <libm-alias-finite.h>__ |

18 | __#include <math-use-builtins.h>__ |

19 | |

20 | *float* |

21 | __ieee754_sqrtf(*float* x) |

22 | { |

23 | __#if USE_SQRTF_BUILTIN__ |

24 | **return** __builtin_sqrtf (x); |

25 | __#else__ |

26 | */* Use generic implementation. */* |

27 | *float* z; |

28 | int32_t sign = (*int*)`0x80000000`; |

29 | int32_t ix,s,q,m,t,i; |

30 | uint32_t r; |

31 | |

32 | GET_FLOAT_WORD(ix,x); |

33 | |

34 | */* take care of Inf and NaN */* |

35 | **if**((ix&`0x7f800000`)==`0x7f800000`) { |

36 | **return** x*x+x; */* sqrt(NaN)=NaN, sqrt(+inf)=+inf* |

37 | * sqrt(-inf)=sNaN */* |

38 | } |

39 | */* take care of zero */* |

40 | **if**(ix<=`0`) { |

41 | **if**((ix&(~sign))==`0`) **return** x;*/* sqrt(+-0) = +-0 */* |

42 | **else** **if**(ix<`0`) |

43 | **return** (x-x)/(x-x); */* sqrt(-ve) = sNaN */* |

44 | } |

45 | */* normalize x */* |

46 | m = (ix>>`23`); |

47 | **if**(m==`0`) { */* subnormal x */* |

48 | **for**(i=`0`;(ix&`0x00800000`)==`0`;i++) ix<<=`1`; |

49 | m -= i-`1`; |

50 | } |

51 | m -= `127`; */* unbias exponent */* |

52 | ix = (ix&`0x007fffff`)|`0x00800000`; |

53 | **if**(m&`1`) */* odd m, double x to make it even */* |

54 | ix += ix; |

55 | m >>= `1`; */* m = [m/2] */* |

56 | |

57 | */* generate sqrt(x) bit by bit */* |

58 | ix += ix; |

59 | q = s = `0`; */* q = sqrt(x) */* |

60 | r = `0x01000000`; */* r = moving bit from right to left */* |

61 | |

62 | **while**(r!=`0`) { |

63 | t = s+r; |

64 | **if**(t<=ix) { |

65 | s = t+r; |

66 | ix -= t; |

67 | q += r; |

68 | } |

69 | ix += ix; |

70 | r>>=`1`; |

71 | } |

72 | |

73 | */* use floating add to find out rounding direction */* |

74 | **if**(ix!=`0`) { |

75 | z = `0x1p0` - `0x1.4484cp-100`; */* trigger inexact flag. */* |

76 | **if** (z >= `0x1p0`) { |

77 | z = `0x1p0` + `0x1.4484cp-100`; |

78 | **if** (z > `0x1p0`) |

79 | q += `2`; |

80 | **else** |

81 | q += (q&`1`); |

82 | } |

83 | } |

84 | ix = (q>>`1`)+`0x3f000000`; |

85 | ix += (m <<`23`); |

86 | SET_FLOAT_WORD(z,ix); |

87 | **return** z; |

88 | __#endif /* ! USE_SQRTF_BUILTIN */__ |

89 | } |

90 | __#ifndef __ieee754_sqrtf__ |

91 | libm_alias_finite (__ieee754_sqrtf, __sqrtf) |

92 | __#endif__ |

93 | |