source: trunk/libs/libmath/s_scalbn.c @ 623

Last change on this file since 623 was 469, checked in by alain, 6 years ago

1) Introduce the libsemaphore library.
2) Introduce a small libmath library, required by the "fft" application..
3) Introduce the multithreaded "fft" application.
4) Fix a bad synchronisation bug in the Copy-On-Write mechanism.

File size: 1.8 KB
Line 
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
27static const double
28two54  = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
29twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
30huge   = 1.0e+300,
31tiny   = 1.0e-300;
32
33double 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
67double scalbn(double x, int n)
68{
69        return scalbln(x, n);
70}
71
Note: See TracBrowser for help on using the repository browser.