overby.me
this repo has no description
1// long double scalbnl( long double x, int n );
2// long double ldexpl( long double x, int n );
3// long double scalblnl( long double x, long int n );
4//
5// Returns x * 2^n computed efficiently, without undue over/underflow.
6//
7// -- Stephen Canon, January 2010
8
9
10.text
11.globl _scalblnl
12.globl _scalbnl
13.globl _ldexpl
14
15#if defined __x86_64__
16
17 #define floatval -24(%rsp)
18 #define exponent -16(%rsp)
19 #define n %edi
20 #define ln %rdi
21 #define nw %di
22
23.align 4
24_scalblnl:
25 fldt 8(%rsp) // load x
26 fld1
27 fstpt floatval // conjure 1.0L
28
29 cmp $16382, ln // if n is so large that 2^n is not
30 jg L_bigPositiveN // representable as a long double,
31 cmp $-16382, ln // branch to handle the scaling
32 jl L_bigNegativeN // carefully.
33 jmp L_scaleX
34
35.align 4
36_scalbnl:
37_ldexpl:
38 fldt 8(%rsp) // load x
39
40#elif defined __i386__
41
42 #define floatval 4(%esp)
43 #define exponent 12(%esp)
44 #define n %eax
45 #define nw %ax
46
47.align 4
48_scalblnl:
49_scalbnl:
50_ldexpl:
51 mov 20(%esp), n // load n
52 fldt floatval // load x
53
54#else
55 #error "This implementation supports 32- and 64-bit Intel only"
56#endif
57
58 fld1
59 fstpt floatval // conjure 1.0L
60
61 cmp $16382, n // if n is so large that 2^n is not
62 jg L_bigPositiveN // representable as a long double,
63 cmp $-16382, n // branch to handle the scaling
64 jl L_bigNegativeN // carefully.
65
66L_scaleX:
67 add nw, exponent // 2^n
68 fldt floatval
69 fmulp // x * 2^n
70 ret
71
72L_bigPositiveN:
73 addw $16382, exponent
74 fldt floatval // 2^16382
75 mov $(3*16382), %edx
76 cmp %edx, n
77 cmovg %edx, n // clamp n to 3 * 16382
780: fmul %st(0), %st(1) // do {
79 sub $16382, n // x *= 2^16382
80 cmp $16382, n // n -= 16382
81 jg 0b // } (while n > 16382)
82 fstp %st(0)
83 movw $0x3fff, exponent
84 jmp L_scaleX
85
86L_bigNegativeN:
87 addw $-16382, exponent
88 fldt floatval // 2^-16382
89 mov $(-3*16382),%edx
90 cmp %edx, n
91 cmovl %edx, n // clamp n to -3 * 16382
920: fmul %st(0), %st(1) // do {
93 add $16382, n // x *= 2^-16382
94 cmp $-16382, n // n += 16382
95 jl 0b // } (while n < -16382)
96 fstp %st(0)
97 movw $0x3fff, exponent
98 jmp L_scaleX