Changeset f250c5a in mainline
 Timestamp:
 20210618T16:05:56Z (2 years ago)
 Children:
 197d59d
 Parents:
 19520ba9
 gitauthor:
 Maurizio Lombardi <mlombard@…> (20181125 16:43:33)
 gitcommitter:
 Maurizio Lombardi <mlombard@…> (20210618 16:05:56)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

uspace/lib/math/generic/scalbn.c
r19520ba9 rf250c5a 18 18 #include <stdint.h> 19 19 20 #include "internal.h" 20 21 21 /* scalbn(x,n) returns x* 2**n computed by exponent 22 /* 23 * scalbn (double x, int n) 24 * scalbn(x,n) returns x* 2**n computed by exponent 22 25 * manipulation rather than by actually performing an 23 26 * exponentiation or a multiplication. 24 27 */ 25 28 26 double scalbn(double x, int n) 29 static const double 30 two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ 31 twom54 = 5.55111512312578270212e17, /* 0x3C900000, 0x00000000 */ 32 huge = 1.0e+300, 33 tiny = 1.0e300; 34 35 double 36 scalbn (double x, int n) 27 37 { 28 union {double f; uint64_t i;} u; 29 double_t y = x; 30 31 if (n > 1023) { 32 y *= 0x1p1023; 33 n = 1023; 34 if (n > 1023) { 35 y *= 0x1p1023; 36 n = 1023; 37 if (n > 1023) 38 n = 1023; 39 } 40 } else if (n < 1022) { 41 /* make sure final n < 53 to avoid double 42 rounding in the subnormal range */ 43 y *= 0x1p1022 * 0x1p53; 44 n += 1022  53; 45 if (n < 1022) { 46 y *= 0x1p1022 * 0x1p53; 47 n += 1022  53; 48 if (n < 1022) 49 n = 1022; 50 } 38 int32_t k,hx,lx; 39 EXTRACT_WORDS(hx,lx,x); 40 k = (hx&0x7ff00000)>>20; /* extract exponent */ 41 if (k==0) { /* 0 or subnormal x */ 42 if ((lx(hx&0x7fffffff))==0) return x; /* +0 */ 43 x *= two54; 44 GET_HIGH_WORD(hx,x); 45 k = ((hx&0x7ff00000)>>20)  54; 46 if (n< 50000) return tiny*x; /*underflow*/ 47 } 48 if (k==0x7ff) return x+x; /* NaN or Inf */ 49 k = k+n; 50 if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */ 51 if (k > 0) /* normal result */ 52 {SET_HIGH_WORD(x,(hx&0x800fffff)(k<<20)); return x;} 53 if (k <= 54) { 54 if (n > 50000) /* in case integer overflow in n+k */ 55 return huge*copysign(huge,x); /*overflow*/ 56 else 57 return tiny*copysign(tiny,x); /*underflow*/ 51 58 } 52 u.i = (uint64_t)(0x3ff+n)<<52; 53 x = y * u.f; 54 55 return x; 59 k += 54; /* subnormal result */ 60 SET_HIGH_WORD(x,(hx&0x800fffff)(k<<20)); 61 return x*twom54; 56 62 } 57 63
Note:
See TracChangeset
for help on using the changeset viewer.