overby.me
this repo has no description
1// double scalbn( double x, int n );
2// double ldexp( double x, int n );
3// double scalb( double x, double n );
4// double scalbln( double x, long int n );
5//
6// Returns x * 2^n computed efficiently, without undue over/underflow.
7// Adapted from Ian Ollmann's scalbn implementation of 2005.
8//
9// -- Stephen Canon, January 2010
10
11.globl _scalbn
12.globl _ldexp
13.globl _scalb
14.globl _scalbln
15
16#if defined __x86_64__
17
18.literal8
19threek: .quad 0x40a7700000000000
20mthreek: .quad 0xc0a7700000000000
21
22 #define floatval -24(%rsp)
23 #define n %edi
24 #define ln %rdi
25 #define clamphi threek(%rip)
26 #define clamplo mthreek(%rip)
27
28.text
29.align 4
30_scalb:
31 movsd %xmm0, floatval
32 pcmpeqd %xmm0, %xmm0
33 fldl floatval // load x on x87
34 psllq $63, %xmm0 // conjure long dbl mantissa of 2^n
35 ucomisd %xmm1, %xmm1 // if isnan(n)
36 jp L_nIsNaN // branch to return n
37
38 minsd clamphi, %xmm1 // clamp n to [-3000, 3000]
39 maxsd clamplo, %xmm1
40 cvtsd2si %xmm1, n // convert clamped n to integer
41 jmp L_common // and jump to the shared scalbn path
42
43L_nIsNaN:
44 fstp %st(0)
45 movapd %xmm1, %xmm0 // since n is a NaN, we just return n
46 ret
47
48.align 4
49_scalbln:
50 movsd %xmm0, floatval
51 pcmpeqd %xmm0, %xmm0
52 fldl floatval // load x on x87
53 psllq $63, %xmm0 // conjure long dbl mantissa of 2^n
54 mov $3000, %rdx
55 cmp %rdx, ln // if n > 3000
56 cmovg %rdx, ln // n == 3000
57 neg %rdx
58 cmp %rdx, ln // if n < -3000
59 cmovl %rdx, ln // n == -3000
60 jmp L_common
61
62.align 4
63_scalbn:
64_ldexp:
65 movsd %xmm0, floatval
66 fldl floatval // load x on x87
67 pcmpeqd %xmm0, %xmm0
68 psllq $63, %xmm0 // conjure long dbl mantissa of 2^n
69 mov $3000, %edx
70 cmp %edx, n // if n > 3000
71 cmovg %edx, n // n == 3000
72 neg %edx
73 cmp %edx, n // if n < -3000
74 cmovl %edx, n // n == -3000
75
76L_common:
77 add $0x3fff, n // long dbl exponent field of 2^n
78 pinsrw $4, n, %xmm0 // insert exponent to make 2^n
79 movdqa %xmm0, floatval
80 fldt floatval
81 fmulp // x * 2^n
82 fstpl floatval // round to double
83 movsd floatval, %xmm0
84 ret
85
86#elif defined __i386__
87
88.literal8
89threek: .quad 0x40a7700000000000
90mthreek: .quad 0xc0a7700000000000
91
92 #define floatval 4(%esp)
93 #define exponent 12(%esp)
94 #define n %eax
95 #define nw %ax
96 #define clamphi (threek-0b)(%ecx)
97 #define clamplo (mthreek-0b)(%ecx)
98
99.text
100.align 4
101_scalb:
102 call 0f
1030: pop %ecx // conjure pic info
104 movsd 12(%esp), %xmm1 // load n
105 fldl floatval // load x on x87
106 ucomisd %xmm1, %xmm1 // if isnan(n)
107 jp L_nIsNaN // branch to return n
108
109 fld1
110 fstpt floatval // conjure 1.0L
111 minsd clamphi, %xmm1 // clamp n to [-3000, 3000]
112 maxsd clamplo, %xmm1
113 cvtsd2si %xmm1, n // convert clamped n to integer
114 jmp L_common // and jump to the shared scalbn path
115
116L_nIsNaN:
117 fstp %st(0)
118 fldl 12(%esp) // n is nan, return n
119 ret
120
121.align 4
122_scalbln:
123_scalbn:
124_ldexp:
125 mov 12(%esp), n // load n
126 fldl floatval // load x on x87
127 fld1
128 fstpt floatval // conjure 1.0L
129 mov $3000, %edx
130 cmp %edx, n // if n > 3000
131 cmovg %edx, n // n == 3000
132 mov $-3000, %edx
133 cmp %edx, n // if n < -3000
134 cmovl %edx, n // n == -3000
135
136L_common:
137 add nw, exponent // 2^n in memory
138 fldt floatval
139 fmulp // x * 2^n
140 fstpl floatval // round to double
141 fldl floatval
142 ret
143
144#else
145 #error "This implementation supports 32- and 64-bit Intel only"
146#endif
147