blob: 40813600117f275a7d9c07fa94e63917cbd7e3bb (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
|
/***********************************************************************
** File: scalb.c
**
** Contains: C source code for implementations of floating-point
** scalb functions defined in header <fp.h>. In
** particular, this file contains implementations of
** functions scalb and scalbl for double and long double
** formats on PowerPC platforms.
**
** Written by: Jon Okada, SANEitation Engineer, ext. 4-4838
**
** Copyright: © 1992 by Apple Computer, Inc., all rights reserved
**
** Change History ( most recent first ):
**
** 28 May 97 ali made an speed improvement for large n,
** removed scalbl.
** 12 Dec 92 JPO First created.
**
***********************************************************************/
#include <math.h>
#include <endian.h>
typedef union
{
struct {
#if (__BYTE_ORDER == __BIG_ENDIAN)
unsigned long int hi;
unsigned long int lo;
#else
unsigned long int lo;
unsigned long int hi;
#endif
} words;
double dbl;
} DblInHex;
static const double twoTo1023 = 8.988465674311579539e307; // 0x1p1023
static const double twoToM1022 = 2.225073858507201383e-308; // 0x1p-1022
/***********************************************************************
double scalb( double x, long int n ) returns its argument x scaled
by the factor 2^m. NaNs, signed zeros, and infinities are propagated
by this function regardless of the value of n.
Exceptions: OVERFLOW/INEXACT or UNDERFLOW inexact may occur;
INVALID for signaling NaN inputs ( quiet NaN returned ).
Calls: none.
***********************************************************************/
libm_hidden_proto(scalb)
#ifdef _SCALB_INT
double scalb ( double x, int n )
#else
double scalb ( double x, double n )
#endif
{
DblInHex xInHex;
xInHex.words.lo = 0UL; // init. low half of xInHex
if ( n > 1023 )
{ // large positive scaling
if ( n > 2097 ) // huge scaling
return ( ( x * twoTo1023 ) * twoTo1023 ) * twoTo1023;
while ( n > 1023 )
{ // scale reduction loop
x *= twoTo1023; // scale x by 2^1023
n -= 1023; // reduce n by 1023
}
}
else if ( n < -1022 )
{ // large negative scaling
if ( n < -2098 ) // huge negative scaling
return ( ( x * twoToM1022 ) * twoToM1022 ) * twoToM1022;
while ( n < -1022 )
{ // scale reduction loop
x *= twoToM1022; // scale x by 2^( -1022 )
n += 1022; // incr n by 1022
}
}
/*******************************************************************************
* -1022 <= n <= 1023; convert n to double scale factor. *
*******************************************************************************/
xInHex.words.hi = ( ( unsigned long ) ( n + 1023 ) ) << 20;
return ( x * xInHex.dbl );
}
libm_hidden_def(scalb)
|
