// long double scalbnl( long double x, int n )// long double scalblnl( long double x, long int n )// Returns x * 2^n computed efficiently, without undue over/underflow.
//
// -- Stephen Canon, January 2010
.text
.globl _scalblnl
.globl _scalbnl
.globl _ldexpl
#if defined __x86_64__
#define floatval -24(%rsp)
#define exponent -16(%rsp)
#define n %edi
#define ln %rdi
#define nw %di
.align 4
_scalblnl:
fldt 8(%rsp) // load x
fld1
fstpt floatval // conjure 1.0L
cmp $16382, ln // if n is so large that 2^n is not
jg L_bigPositiveN // representable as a long double,
cmp $-16382, ln // branch to handle the scaling
jl L_bigNegativeN // carefully.
jmp L_scaleX
.align 4
_scalbnl:
_ldexpl:
fldt 8(%rsp) // load x
#elif defined __i386__
#define floatval 4(%esp)
#define exponent 12(%esp)
#define n %eax
#define nw %ax
.align 4
_scalblnl:
_scalbnl:
_ldexpl:
mov 20(%esp), n // load n
fldt floatval // load x
#else
#error "This implementation supports 32- and 64-bit Intel only"
#endif
fld1
fstpt floatval // conjure 1.0L
cmp $16382, n // if n is so large that 2^n is not
jg L_bigPositiveN // representable as a long double,
cmp $-16382, n // branch to handle the scaling
jl L_bigNegativeN // carefully.
L_scaleX:
add nw, exponent // 2^n
fldt floatval
fmulp // x * 2^n
ret
L_bigPositiveN:
addw $16382, exponent
fldt floatval // 2^16382
mov $(3*16382), %edx
cmp %edx, n
cmovg %edx, n // clamp n to 3 * 16382
0: fmul %st(0), %st(1) // do {
sub $16382, n // x *= 2^16382
cmp $16382, n // n -= 16382
jg 0b // } (while n > 16382)
fstp %st(0)
movw $0x3fff, exponent
jmp L_scaleX
L_bigNegativeN:
addw $-16382, exponent
fldt floatval // 2^-16382
mov $(-3*16382),%edx
cmp %edx, n
cmovl %edx, n // clamp n to -3 * 16382
0: fmul %st(0), %st(1) // do {
add $16382, n // x *= 2^-16382
cmp $-16382, n // n += 16382
jl 0b // } (while n < -16382)
fstp %st(0)
movw $0x3fff, exponent
jmp L_scaleX