1 | /* |
---|

2 | * ==================================================== |
---|

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

4 | * |
---|

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

6 | * Permission to use, copy, modify, and distribute this |
---|

7 | * software is freely granted, provided that this notice |
---|

8 | * is preserved. |
---|

9 | * ==================================================== |
---|

10 | */ |
---|

11 | |
---|

12 | /* |
---|

13 | * Modified for ALMOS-MKH OS at UPMC, France, August 2018. (Alain Greiner) |
---|

14 | */ |
---|

15 | |
---|

16 | /* |
---|

17 | * scalbln(double x, long n) |
---|

18 | * scalbln(x,n) returns x * 2**n computed by exponent |
---|

19 | * manipulation rather than by actually performing an |
---|

20 | * exponentiation or a multiplication. |
---|

21 | */ |
---|

22 | |
---|

23 | #include "math.h" |
---|

24 | #include "math_private.h" |
---|

25 | #include <limits.h> |
---|

26 | |
---|

27 | static const double |
---|

28 | two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ |
---|

29 | twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ |
---|

30 | huge = 1.0e+300, |
---|

31 | tiny = 1.0e-300; |
---|

32 | |
---|

33 | double scalbln(double x, long n) |
---|

34 | { |
---|

35 | int32_t k, hx, lx; |
---|

36 | |
---|

37 | EXTRACT_WORDS(hx, lx, x); |
---|

38 | k = (hx & 0x7ff00000) >> 20; /* extract exponent */ |
---|

39 | if (k == 0) { /* 0 or subnormal x */ |
---|

40 | if ((lx | (hx & 0x7fffffff)) == 0) |
---|

41 | return x; /* +-0 */ |
---|

42 | x *= two54; |
---|

43 | GET_HIGH_WORD(hx, x); |
---|

44 | k = ((hx & 0x7ff00000) >> 20) - 54; |
---|

45 | } |
---|

46 | if (k == 0x7ff) |
---|

47 | return x + x; /* NaN or Inf */ |
---|

48 | k = k + n; |
---|

49 | if (k > 0x7fe) |
---|

50 | return huge * copysign(huge, x); /* overflow */ |
---|

51 | if (n < -50000) |
---|

52 | return tiny * copysign(tiny, x); /* underflow */ |
---|

53 | if (k > 0) { /* normal result */ |
---|

54 | SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20)); |
---|

55 | return x; |
---|

56 | } |
---|

57 | if (k <= -54) { |
---|

58 | if (n > 50000) /* in case integer overflow in n+k */ |
---|

59 | return huge * copysign(huge, x); /* overflow */ |
---|

60 | return tiny * copysign(tiny, x); /* underflow */ |
---|

61 | } |
---|

62 | k += 54; /* subnormal result */ |
---|

63 | SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20)); |
---|

64 | return x * twom54; |
---|

65 | } |
---|

66 | |
---|

67 | double scalbn(double x, int n) |
---|

68 | { |
---|

69 | return scalbln(x, n); |
---|

70 | } |
---|

71 | |
---|