| |
| # |
| # (C) 2008-2009 Advanced Micro Devices, Inc. All Rights Reserved. |
| # |
| # This file is part of libacml_mv. |
| # |
| # libacml_mv is free software; you can redistribute it and/or |
| # modify it under the terms of the GNU Lesser General Public |
| # License as published by the Free Software Foundation; either |
| # version 2.1 of the License, or (at your option) any later version. |
| # |
| # libacml_mv is distributed in the hope that it will be useful, |
| # but WITHOUT ANY WARRANTY; without even the implied warranty of |
| # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| # Lesser General Public License for more details. |
| # |
| # You should have received a copy of the GNU Lesser General Public |
| # License along with libacml_mv. If not, see |
| # <http://www.gnu.org/licenses/>. |
| # |
| # |
| |
| |
| # |
| # An implementation of the sinf function. |
| # |
| # Prototype: |
| # |
| # double sinf(double x); |
| # |
| # Computes sinf(x). |
| # It will provide proper C99 return values, |
| # but may not raise floating point status bits properly. |
| # Based on the NAG C implementation. |
| # |
| # |
| #ifdef __ELF__ |
| .section .note.GNU-stack,"",@progbits |
| #endif |
| |
| .data |
| .align 32 |
| .L__real_3ff0000000000000: .quad 0x03ff0000000000000 # 1.0 |
| .quad 0 # for alignment |
| .L__real_3fe0000000000000: .quad 0x03fe0000000000000 # 0.5 |
| .quad 0 |
| .L__real_3fc5555555555555: .quad 0x03fc5555555555555 # 0.166666666666 |
| .quad 0 |
| .L__real_3fe45f306dc9c883: .quad 0x03fe45f306dc9c883 # twobypi |
| .quad 0 |
| .L__real_3ff921fb54400000: .quad 0x03ff921fb54400000 # piby2_1 |
| .quad 0 |
| .L__real_3dd0b4611a626331: .quad 0x03dd0b4611a626331 # piby2_1tail |
| .quad 0 |
| .L__real_3dd0b4611a600000: .quad 0x03dd0b4611a600000 # piby2_2 |
| .quad 0 |
| .L__real_3ba3198a2e037073: .quad 0x03ba3198a2e037073 # piby2_2tail |
| .quad 0 |
| .L__real_411E848000000000: .quad 0x415312d000000000 # 5e6 0x0411E848000000000 # 5e5 |
| .quad 0 |
| |
| .align 32 |
| .Lcosfarray: |
| .quad 0x0bfe0000000000000 # -0.5 c0 |
| .quad 0 |
| .quad 0x03fa5555555555555 # 0.0416667 c1 |
| .quad 0 |
| .quad 0x0bf56c16c16c16c16 # -0.00138889 c2 |
| .quad 0 |
| .quad 0x03EFA01A01A01A019 # 2.48016e-005 c3 |
| .quad 0 |
| .quad 0x0be927e4fb7789f5c # -2.75573e-007 c4 |
| .quad 0 |
| |
| .align 32 |
| .Lsinfarray: |
| .quad 0x0bfc5555555555555 # -0.166667 s1 |
| .quad 0 |
| .quad 0x03f81111111111111 # 0.00833333 s2 |
| .quad 0 |
| .quad 0x0bf2a01a01a01a01a # -0.000198413 s3 |
| .quad 0 |
| .quad 0x03ec71de3a556c734 # 2.75573e-006 s4 |
| .quad 0 |
| |
| .text |
| .align 32 |
| .p2align 4,,15 |
| |
| #include "fn_macros.h" |
| #define fname FN_PROTOTYPE(sinf) |
| #define fname_special _sinf_special@PLT |
| |
| # define local variable storage offsets |
| .equ p_temp, 0x30 # temporary for get/put bits operation |
| .equ p_temp1, 0x40 # temporary for get/put bits operation |
| .equ r, 0x50 # pointer to r for amd_remainder_piby2 |
| .equ region, 0x60 # pointer to region for amd_remainder_piby2 |
| .equ stack_size, 0x88 |
| |
| .globl fname |
| .type fname,@function |
| |
| fname: |
| sub $stack_size, %rsp |
| xorpd %xmm2, %xmm2 # zeroed out for later use |
| |
| ## if NaN or inf |
| movd %xmm0, %edx |
| mov $0x07f800000, %eax |
| mov %eax, %r10d |
| and %edx, %r10d |
| cmp %eax, %r10d |
| jz .Lsinf_naninf |
| |
| # GET_BITS_DP64(x, ux); |
| # get the input value to an integer register. |
| cvtss2sd %xmm0, %xmm0 # convert input to double. |
| movsd %xmm0,p_temp(%rsp) # get the input value to an integer register. |
| |
| mov p_temp(%rsp), %rdx # rdx is ux |
| |
| # ax = (ux & ~SIGNBIT_DP64); |
| mov $0x07fffffffffffffff, %r10 |
| and %rdx, %r10 # r10 is ax |
| mov $1, %r8d # for determining region later on |
| |
| ## if (ax <= 0x3fe921fb54442d18) /* abs(x) <= pi/4 */ |
| mov $0x03fe921fb54442d18, %rax |
| cmp %rax, %r10 |
| jg .Lsinf_reduce |
| |
| ## if (ax < 0x3f80000000000000) /* abs(x) < 2.0^(-7) */ |
| mov $0x3f80000000000000, %rax |
| cmp %rax, %r10 |
| jge .Lsinf_small |
| |
| ## if (ax < 0x3f20000000000000) /* abs(x) < 2.0^(-13) */ |
| mov $0x3f20000000000000, %rax |
| cmp %rax, %r10 |
| jge .Lsinf_smaller |
| |
| # sinf = x; |
| jmp .Lsinf_cleanup # done |
| |
| ## else |
| |
| .Lsinf_smaller: |
| # sinf = x - x^3 * 0.1666666666666666666; |
| movsd %xmm0, %xmm2 |
| movsd .L__real_3fc5555555555555(%rip), %xmm4 # 0.1666666666666666666 |
| mulsd %xmm2, %xmm2 # x^2 |
| mulsd %xmm0, %xmm2 # x^3 |
| mulsd %xmm4, %xmm2 # x^3 * 0.1666666666666666666 |
| subsd %xmm2, %xmm0 # x - x^3 * 0.1666666666666666666 |
| jmp .Lsinf_cleanup |
| |
| .Lsinf_small: |
| movsd %xmm0, %xmm2 |
| mulsd %xmm0, %xmm2 # x2 |
| |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| # region 0 or 2 - do a sinf calculation |
| # zs = x + x3((s1 + x2 * s2) + x4(s3 + x2 * s4)); |
| movsd .Lsinfarray+0x30(%rip), %xmm1 # s4 |
| mulsd %xmm2, %xmm1 # s4x2 |
| movsd %xmm2, %xmm4 # move for x4 |
| movsd .Lsinfarray+0x10(%rip), %xmm5 # s2 |
| mulsd %xmm2, %xmm4 # x4 |
| movsd %xmm0, %xmm3 # move for x3 |
| mulsd %xmm2, %xmm5 # s2x2 |
| mulsd %xmm2, %xmm3 # x3 |
| addsd .Lsinfarray+0x20(%rip), %xmm1 # s3+s4x2 |
| mulsd %xmm4, %xmm1 # s3x4+s4x6 |
| addsd .Lsinfarray(%rip), %xmm5 # s1+s2x2 |
| addsd %xmm5, %xmm1 # s1+s2x2+s3x4+s4x6 |
| mulsd %xmm3, %xmm1 # x3(s1+s2x2+s3x4+s4x6) |
| addsd %xmm1, %xmm0 # x + x3(s1+s2x2+s3x4+s4x6) |
| jmp .Lsinf_cleanup |
| |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| .align 32 |
| .Lsinf_reduce: |
| |
| # xneg = (ax != ux); |
| cmp %r10, %rdx |
| mov $0, %r11d |
| |
| ## if (xneg) x = -x; |
| jz .L50e5 |
| mov $1, %r11d |
| subsd %xmm0, %xmm2 |
| movsd %xmm2, %xmm0 |
| |
| .L50e5: |
| ## if (x < 5.0e5) |
| cmp .L__real_411E848000000000(%rip), %r10 |
| jae .Lsinf_reduce_precise |
| |
| # reduce the argument to be in a range from -pi/4 to +pi/4 |
| # by subtracting multiples of pi/2 |
| movsd %xmm0, %xmm2 |
| movsd .L__real_3fe45f306dc9c883(%rip), %xmm3 # twobypi |
| movsd %xmm0, %xmm4 |
| movsd .L__real_3fe0000000000000(%rip), %xmm5 # .5 |
| mulsd %xmm3, %xmm2 |
| |
| #/* How many pi/2 is x a multiple of? */ |
| # xexp = ax >> EXPSHIFTBITS_DP64; |
| mov %r10, %r9 |
| shr $52, %r9 #>>EXPSHIFTBITS_DP64 |
| |
| # npi2 = (int)(x * twobypi + 0.5); |
| addsd %xmm5, %xmm2 # npi2 |
| |
| movsd .L__real_3ff921fb54400000(%rip), %xmm3 # piby2_1 |
| cvttpd2dq %xmm2, %xmm0 # convert to integer |
| movsd .L__real_3dd0b4611a626331(%rip), %xmm1 # piby2_1tail |
| cvtdq2pd %xmm0, %xmm2 # and back to double. |
| |
| # /* Subtract the multiple from x to get an extra-precision remainder */ |
| # rhead = x - npi2 * piby2_1; |
| mulsd %xmm2, %xmm3 |
| subsd %xmm3, %xmm4 # rhead |
| |
| # rtail = npi2 * piby2_1tail; |
| mulsd %xmm2, %xmm1 |
| movd %xmm0, %eax |
| |
| # GET_BITS_DP64(rhead-rtail, uy); ; originally only rhead |
| movsd %xmm4, %xmm0 |
| subsd %xmm1, %xmm0 |
| |
| movsd .L__real_3dd0b4611a600000(%rip), %xmm3 # piby2_2 |
| movsd %xmm0,p_temp(%rsp) |
| movsd .L__real_3ba3198a2e037073(%rip), %xmm5 # piby2_2tail |
| mov p_temp(%rsp), %rcx # rcx is rhead-rtail |
| |
| # xmm0=r, xmm4=rhead, xmm1=rtail, xmm2=npi2, xmm3=temp for calc, xmm5= temp for calc |
| # expdiff = xexp - ((uy & EXPBITS_DP64) >> EXPSHIFTBITS_DP64); |
| shl $1, %rcx # strip any sign bit |
| shr $53, %rcx #>> EXPSHIFTBITS_DP64 +1 |
| sub %rcx, %r9 #expdiff |
| |
| ## if (expdiff > 15) |
| cmp $15, %r9 |
| jle .Lexpdiff15 |
| |
| # /* The remainder is pretty small compared with x, which |
| # implies that x is a near multiple of pi/2 |
| # (x matches the multiple to at least 15 bits) */ |
| |
| # t = rhead; |
| movsd %xmm4, %xmm1 |
| |
| # rtail = npi2 * piby2_2; |
| mulsd %xmm2, %xmm3 |
| |
| # rhead = t - rtail; |
| mulsd %xmm2, %xmm5 # npi2 * piby2_2tail |
| subsd %xmm3, %xmm4 # rhead |
| |
| # rtail = npi2 * piby2_2tail - ((t - rhead) - rtail); |
| subsd %xmm4, %xmm1 # t - rhead |
| subsd %xmm3, %xmm1 # -rtail |
| subsd %xmm1, %xmm5 #rtail |
| |
| # r = rhead - rtail; |
| movsd %xmm4, %xmm0 |
| |
| #HARSHA |
| #xmm1=rtail |
| movsd %xmm5, %xmm1 |
| subsd %xmm5, %xmm0 |
| |
| # xmm0=r, xmm4=rhead, xmm1=rtail |
| .Lexpdiff15: |
| # region = npi2 & 3; |
| # No need rr for float case |
| |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| ## if the input was close to a pi/2 multiple |
| # The original NAG code missed this trick. If the input is very close to n*pi/2 after |
| # reduction, |
| # then the sinf is ~ 1.0 , to within 15 bits, when r is < 2^-13. We already |
| # have x at this point, so we can skip the sinf polynomials. |
| |
| cmp $0x03f2, %rcx ## if r small. |
| jge .Lsinf_piby4 # use taylor series if not |
| cmp $0x03de, %rcx ## if r really small. |
| jle .Lr_small # then sinf(r) = 0 |
| |
| movsd %xmm0, %xmm2 |
| mulsd %xmm2, %xmm2 #x^2 |
| |
| ## if region is 0 or 2 do a sinf calc. |
| and %eax, %r8d |
| jnz .Lcosfregion |
| |
| # region 0 or 2 do a sinf calculation |
| # use simply polynomial |
| # x - x*x*x*0.166666666666666666; |
| movsd .L__real_3fc5555555555555(%rip), %xmm3 # |
| mulsd %xmm0, %xmm3 # * x |
| mulsd %xmm2, %xmm3 # * x^2 |
| subsd %xmm3, %xmm0 # xs |
| jmp .Ladjust_region |
| |
| .align 32 |
| .Lcosfregion: |
| # region 1 or 3 do a cosf calculation |
| # use simply polynomial |
| # 1.0 - x*x*0.5; |
| movsd .L__real_3ff0000000000000(%rip), %xmm0 # 1.0 |
| mulsd .L__real_3fe0000000000000(%rip), %xmm2 # 0.5 *x^2 |
| subsd %xmm2, %xmm0 # xc |
| jmp .Ladjust_region |
| |
| .align 32 |
| .Lr_small: |
| ## if region is 1 or 3 do a cosf calc. |
| and %eax, %r8d |
| jz .Ladjust_region |
| |
| # odd |
| movsd .L__real_3ff0000000000000(%rip), %xmm0 # cosf(r) is a 1 |
| jmp .Ladjust_region |
| |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| .Lsinf_reduce_precise: |
| # // Reduce x into range [-pi/4,pi/4] |
| # __amd_remainder_piby2d2f(x, &r, ®ion); |
| |
| mov %r11,p_temp(%rsp) |
| lea region(%rsp), %rdx |
| lea r(%rsp), %rsi |
| movd %xmm0, %rdi |
| sub $0x20, %rsp |
| |
| call __amd_remainder_piby2d2f@PLT |
| |
| add $0x20, %rsp |
| mov p_temp(%rsp), %r11 |
| mov $1, %r8d # for determining region later on |
| movsd r(%rsp), %xmm1 #//x |
| mov region(%rsp), %eax #//region |
| |
| # xmm0 = x, xmm4 = xx, r8d = 1, eax= region |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| # perform taylor series to calc sinfx, cosfx |
| .Lsinf_piby4: |
| # x2 = r * r; |
| movsd %xmm0, %xmm2 |
| mulsd %xmm0, %xmm2 #x2 |
| |
| ## if region is 0 or 2 do a sinf calc. |
| and %eax, %r8d |
| jnz .Lcosfregion2 |
| |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| # region 0 or 2 do a sinf calculation |
| # zs = x + x3((s1 + x2 * s2) + x4(s3 + x2 * s4)); |
| movsd .Lsinfarray+0x30(%rip), %xmm1 # s4 |
| mulsd %xmm2, %xmm1 # s4x2 |
| movsd %xmm2, %xmm4 # move for x4 |
| mulsd %xmm2, %xmm4 # x4 |
| movsd .Lsinfarray+0x10(%rip), %xmm5 # s2 |
| mulsd %xmm2, %xmm5 # s2x2 |
| movsd %xmm0, %xmm3 # move for x3 |
| mulsd %xmm2, %xmm3 # x3 |
| addsd .Lsinfarray+0x20(%rip), %xmm1 # s3+s4x2 |
| mulsd %xmm4, %xmm1 # s3x4+s4x6 |
| addsd .Lsinfarray(%rip), %xmm5 # s1+s2x2 |
| addsd %xmm5, %xmm1 # s1+s2x2+s3x4+s4x6 |
| mulsd %xmm3, %xmm1 # x3(s1+s2x2+s3x4+s4x6) |
| addsd %xmm1, %xmm0 # x + x3(s1+s2x2+s3x4+s4x6) |
| |
| jmp .Ladjust_region |
| |
| .align 32 |
| .Lcosfregion2: |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| # region 1 or 3 - do a cosf calculation |
| # zc = 1-0.5*x2+ c1*x4 +c2*x6 +c3*x8 + c4*x10 for a higher precision |
| movsd .Lcosfarray+0x40(%rip), %xmm1 # c4 |
| movsd %xmm2, %xmm4 # move for x4 |
| mulsd %xmm2, %xmm1 # c4x2 |
| movsd .Lcosfarray+0x20(%rip), %xmm3 # c2 |
| mulsd %xmm2, %xmm4 # x4 |
| movsd .Lcosfarray(%rip), %xmm0 # c0 |
| mulsd %xmm2, %xmm3 # c2x2 |
| mulsd %xmm2, %xmm0 # c0x2 (=-0.5x2) |
| addsd .Lcosfarray+0x30(%rip), %xmm1 # c3+c4x2 |
| mulsd %xmm4, %xmm1 # c3x4 + c4x6 |
| addsd .Lcosfarray+0x10(%rip), %xmm3 # c1+c2x2 |
| addsd %xmm3, %xmm1 # c1 + c2x2 + c3x4 + c4x6 |
| mulsd %xmm4, %xmm1 # c1x4 + c2x6 + c3x8 + c4x10 |
| addsd .L__real_3ff0000000000000(%rip), %xmm0 # 1 - 0.5x2 |
| addsd %xmm1, %xmm0 # 1 - 0.5x2 + c1x4 + c2x6 + c3x8 + c4x10 |
| |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| |
| .align 32 |
| .Ladjust_region: # positive or negative |
| # switch (region) |
| shr $1, %eax |
| mov %eax, %ecx |
| and %r11d, %eax |
| |
| not %ecx |
| not %r11d |
| and %r11d, %ecx |
| |
| or %ecx, %eax |
| and $1, %eax |
| jnz .Lsinf_cleanup |
| |
| ## if the original region 0, 1 and arg is negative, then we negate the result. |
| ## if the original region 2, 3 and arg is positive, then we negate the result. |
| movsd %xmm0, %xmm2 |
| xorpd %xmm0, %xmm0 |
| subsd %xmm2, %xmm0 |
| |
| .align 32 |
| .Lsinf_cleanup: |
| cvtsd2ss %xmm0, %xmm0 |
| add $stack_size, %rsp |
| ret |
| |
| .align 32 |
| .Lsinf_naninf: |
| call fname_special |
| add $stack_size, %rsp |
| ret |
| |
| |
