| |
| # |
| # (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 sin function. |
| # |
| # Prototype: |
| # |
| # double sin(double x); |
| # |
| # Computes sin(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_3fe0000000000000: .quad 0x03fe0000000000000 # 0.5 |
| .quad 0 # for alignment |
| .L__real_3ff0000000000000: .quad 0x03ff0000000000000 # 1.0 |
| .quad 0 |
| .L__real_3fc5555555555555: .quad 0x03fc5555555555555 # 0.166666666666 |
| .quad 0 |
| .L__real_3fe45f306dc9c883: .quad 0x03fe45f306dc9c883 # twobypi |
| .quad 0 |
| .L__real_411E848000000000: .quad 0x415312d000000000 # 5e6 0x0411E848000000000 # 5e5 |
| .quad 0 |
| .L__real_7fffffffffffffff: .quad 0x07fffffffffffffff # Sign bit zero |
| .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 |
| |
| .align 32 |
| .Lcosarray: |
| .quad 0x03fa5555555555555 # 0.0416667 c1 |
| .quad 0 |
| .quad 0x0bf56c16c16c16967 # -0.00138889 c2 |
| .quad 0 |
| .quad 0x03EFA01A019F4EC91 # 2.48016e-005 c3 |
| .quad 0 |
| .quad 0x0bE927E4FA17F667B # -2.75573e-007 c4 |
| .quad 0 |
| .quad 0x03E21EEB690382EEC # 2.08761e-009 c5 |
| .quad 0 |
| .quad 0x0bDA907DB47258AA7 # -1.13826e-011 c6 |
| .quad 0 |
| |
| .align 32 |
| .Lsinarray: |
| .quad 0x0bfc5555555555555 # -0.166667 s1 |
| .quad 0 |
| .quad 0x03f81111111110bb3 # 0.00833333 s2 |
| .quad 0 |
| .quad 0x0bf2a01a019e83e5c # -0.000198413 s3 |
| .quad 0 |
| .quad 0x03ec71de3796cde01 # 2.75573e-006 s4 |
| .quad 0 |
| .quad 0x0be5ae600b42fdfa7 # -2.50511e-008 s5 |
| .quad 0 |
| .quad 0x03de5e0b2f9a43bb8 # 1.59181e-010 s6 |
| .quad 0 |
| |
| .text |
| .align 32 |
| .p2align 4,,15 |
| |
| #include "fn_macros.h" |
| #define fname FN_PROTOTYPE(sin) |
| #define fname_special _sin_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 rr, 0x60 # pointer to rr for amd_remainder_piby2 |
| .equ region, 0x70 # pointer to region for amd_remainder_piby2 |
| .equ stack_size, 0x98 |
| |
| .globl fname |
| .type fname,@function |
| |
| fname: |
| sub $stack_size, %rsp |
| xorpd %xmm2, %xmm2 # zeroed out for later use |
| |
| # GET_BITS_DP64(x, ux); |
| # get the input value to an integer register. |
| movsd %xmm0, p_temp(%rsp) |
| mov p_temp(%rsp), %rdx # rdx is ux |
| |
| ## if NaN or inf |
| mov $0x07ff0000000000000, %rax |
| mov %rax, %r10 |
| and %rdx, %r10 |
| cmp %rax, %r10 |
| jz .Lsin_naninf |
| |
| # 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 .Lsin_reduce |
| |
| ## if (ax < 0x3f20000000000000) /* abs(x) < 2.0^(-13) */ |
| mov $0x03f20000000000000, %rax |
| cmp %rax, %r10 |
| jge .Lsin_small |
| |
| ## if (ax < 0x3e40000000000000) /* abs(x) < 2.0^(-27) */ |
| mov $0x03e40000000000000, %rax |
| cmp %rax, %r10 |
| jge .Lsin_smaller |
| |
| # sin = 1.0; |
| jmp .Lsin_cleanup |
| |
| .align 32 |
| .Lsin_smaller: |
| # sin = 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 .Lsin_cleanup |
| |
| .align 32 |
| .Lsin_small: |
| # sin = sin_piby4(x, 0.0); |
| movsd .L__real_3fe0000000000000(%rip), %xmm5 # .5 |
| |
| .Lsin_piby4_noreduce: |
| movsd %xmm0, %xmm2 |
| mulsd %xmm0, %xmm2 # x2 |
| |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| # region 0 or 2 - do a sin calculation |
| # zs = (s2 + x2 * (s3 + x2 * (s4 + x2 * (s5 + x2 * s6)))); |
| movsd .Lsinarray+0x50(%rip), %xmm3 # s6 |
| mulsd %xmm2, %xmm3 # x2s6 |
| movsd .Lsinarray+0x20(%rip), %xmm5 # s3 |
| movsd %xmm2, %xmm1 # move for x4 |
| mulsd %xmm2, %xmm1 # x4 |
| mulsd %xmm2, %xmm5 # x2s3 |
| movsd %xmm0, %xmm4 # move for x3 |
| addsd .Lsinarray+0x40(%rip), %xmm3 # s5+x2s6 |
| mulsd %xmm2, %xmm1 # x6 |
| mulsd %xmm2, %xmm3 # x2(s5+x2s6) |
| mulsd %xmm2, %xmm4 # x3 |
| addsd .Lsinarray+0x10(%rip), %xmm5 # s2+x2s3 |
| mulsd %xmm2, %xmm5 # x2(s2+x2s3) |
| addsd .Lsinarray+0x30(%rip), %xmm3 # s4 + x2(s5+x2s6) |
| mulsd %xmm1, %xmm3 # x6(s4 + x2(s5+x2s6)) |
| addsd .Lsinarray(%rip), %xmm5 # s1+x2(s2+x2s3) |
| addsd %xmm5, %xmm3 # zs |
| mulsd %xmm3, %xmm4 # *x3 |
| addsd %xmm4, %xmm0 # +x |
| jmp .Lsin_cleanup |
| |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| .align 16 |
| .Lsin_reduce: |
| # xneg = (ax != ux); |
| cmp %r10, %rdx |
| mov $0, %r11d |
| |
| ## if (xneg) x = -x; |
| jz .Lpositive |
| mov $1, %r11d |
| subsd %xmm0, %xmm2 |
| movsd %xmm2, %xmm0 |
| |
| .align 16 |
| .Lpositive: |
| ## if (x < 5.0e5) |
| cmp .L__real_411E848000000000(%rip), %r10 |
| jae .Lsin_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 float. |
| |
| # /* 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); |
| 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 .Lexplediff15 |
| |
| # /* 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 |
| .Lexplediff15: |
| # region = npi2 & 3; |
| |
| subsd %xmm0, %xmm4 # rhead-r |
| subsd %xmm1, %xmm4 # rr = (rhead-r) - rtail |
| |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| ## 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 sin is ~ 1.0 , to within 53 bits, when r is < 2^-27. We already |
| # have x at this point, so we can skip the sin polynomials. |
| |
| cmp $0x03f2, %rcx # if r small. |
| jge .Lsin_piby4 # use taylor series if not |
| cmp $0x03de, %rcx # if r really small. |
| jle .Lr_small # then sin(r) = 0 |
| |
| movsd %xmm0, %xmm2 |
| mulsd %xmm2, %xmm2 # x^2 |
| |
| ## if region is 0 or 2 do a sin calc. |
| and %eax, %r8d |
| jnz .Lcossmall |
| |
| # region 0 or 2 do a sin 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 16 |
| .Lcossmall: |
| # region 1 or 3 do a cos 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 16 |
| .Lr_small: |
| ## if region is 1 or 3 do a cos calc. |
| and %eax, %r8d |
| jz .Ladjust_region |
| |
| # odd |
| movsd .L__real_3ff0000000000000(%rip), %xmm0 # cos(r) is a 1 |
| jmp .Ladjust_region |
| |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| .align 32 |
| .Lsin_reduce_precise: |
| # // Reduce x into range [-pi/4,pi/4] |
| # __amd_remainder_piby2(x, &r, &rr, ®ion); |
| |
| mov %r11,p_temp(%rsp) |
| lea region(%rsp), %rdx |
| lea rr(%rsp), %rsi |
| lea r(%rsp), %rdi |
| |
| call __amd_remainder_piby2@PLT |
| |
| mov p_temp(%rsp), %r11 |
| mov $1, %r8d # for determining region later on |
| movsd r(%rsp), %xmm0 # x |
| movsd rr(%rsp), %xmm4 # xx |
| mov region(%rsp), %eax # region |
| |
| # xmm0 = x, xmm4 = xx, r8d = 1, eax= region |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| .align 16 |
| # perform taylor series to calc sinx, sinx |
| .Lsin_piby4: |
| # x2 = r * r; |
| |
| #xmm4 = a part of rr for the sin path, xmm4 is overwritten in the sin path |
| #instead use xmm3 because that was freed up in the sin path, xmm3 is overwritten in sin path |
| movsd %xmm0, %xmm3 |
| movsd %xmm0, %xmm2 |
| mulsd %xmm0, %xmm2 # x2 |
| |
| ## if region is 0 or 2 do a sin calc. |
| and %eax, %r8d |
| jnz .Lcosregion |
| |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| # region 0 or 2 do a sin calculation |
| movsd .Lsinarray+0x50(%rip), %xmm3 # s6 |
| mulsd %xmm2, %xmm3 # x2s6 |
| movsd .Lsinarray+0x20(%rip), %xmm5 # s3 |
| movsd %xmm4,p_temp(%rsp) # store xx |
| movsd %xmm2, %xmm1 # move for x4 |
| mulsd %xmm2, %xmm1 # x4 |
| movsd %xmm0,p_temp1(%rsp) # store x |
| mulsd %xmm2, %xmm5 # x2s3 |
| movsd %xmm0, %xmm4 # move for x3 |
| addsd .Lsinarray+0x40(%rip), %xmm3 # s5+x2s6 |
| mulsd %xmm2, %xmm1 # x6 |
| mulsd %xmm2, %xmm3 # x2(s5+x2s6) |
| mulsd %xmm2, %xmm4 # x3 |
| addsd .Lsinarray+0x10(%rip), %xmm5 # s2+x2s3 |
| mulsd %xmm2, %xmm5 # x2(s2+x2s3) |
| addsd .Lsinarray+0x30(%rip), %xmm3 # s4 + x2(s5+x2s6) |
| mulsd .L__real_3fe0000000000000(%rip), %xmm2 # 0.5 *x2 |
| movsd p_temp(%rsp), %xmm0 # load xx |
| mulsd %xmm1, %xmm3 # x6(s4 + x2(s5+x2s6)) |
| addsd .Lsinarray(%rip), %xmm5 # s1+x2(s2+x2s3) |
| mulsd %xmm0, %xmm2 # 0.5 * x2 *xx |
| addsd %xmm5, %xmm3 # zs |
| mulsd %xmm3, %xmm4 # *x3 |
| subsd %xmm2, %xmm4 # x3*zs - 0.5 * x2 *xx |
| addsd %xmm4, %xmm0 # +xx |
| addsd p_temp1(%rsp), %xmm0 # +x |
| jmp .Ladjust_region |
| |
| .align 16 |
| .Lcosregion: |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| # region 1 or 3 - do a cos calculation |
| # zc = (c2 + x2 * (c3 + x2 * (c4 + x2 * (c5 + x2 * c6)))); |
| mulsd %xmm0, %xmm4 # x*xx |
| movsd .L__real_3fe0000000000000(%rip), %xmm5 |
| movsd .Lcosarray+0x50(%rip), %xmm1 # c6 |
| movsd .Lcosarray+0x20(%rip), %xmm0 # c3 |
| mulsd %xmm2, %xmm5 # r = 0.5 *x2 |
| movsd %xmm2, %xmm3 # copy of x2 |
| movsd %xmm4,p_temp(%rsp) # store x*xx |
| mulsd %xmm2, %xmm1 # c6*x2 |
| mulsd %xmm2, %xmm0 # c3*x2 |
| subsd .L__real_3ff0000000000000(%rip), %xmm5 # -t=r-1.0 ;trash r |
| mulsd %xmm2, %xmm3 # x4 |
| addsd .Lcosarray+0x40(%rip), %xmm1 # c5+x2c6 |
| addsd .Lcosarray+0x10(%rip), %xmm0 # c2+x2C3 |
| addsd .L__real_3ff0000000000000(%rip), %xmm5 # 1 + (-t) ;trash t |
| mulsd %xmm2, %xmm3 # x6 |
| mulsd %xmm2, %xmm1 # x2(c5+x2c6) |
| mulsd %xmm2, %xmm0 # x2(c2+x2C3) |
| movsd %xmm2, %xmm4 # copy of x2 |
| mulsd .L__real_3fe0000000000000(%rip), %xmm4 # r recalculate |
| addsd .Lcosarray+0x30(%rip), %xmm1 # c4 + x2(c5+x2c6) |
| addsd .Lcosarray(%rip), %xmm0 # c1+x2(c2+x2C3) |
| mulsd %xmm2, %xmm2 # x4 recalculate |
| subsd %xmm4, %xmm5 # (1 + (-t)) - r |
| mulsd %xmm3, %xmm1 # x6(c4 + x2(c5+x2c6)) |
| addsd %xmm1, %xmm0 # zc |
| subsd .L__real_3ff0000000000000(%rip), %xmm4 # t relaculate |
| subsd p_temp(%rsp), %xmm5 # ((1 + (-t)) - r) - x*xx |
| mulsd %xmm2, %xmm0 # x4 * zc |
| addsd %xmm5, %xmm0 # x4 * zc + ((1 + (-t)) - r -x*xx) |
| subsd %xmm4, %xmm0 # result - (-t) |
| |
| #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; |
| |
| .align 16 |
| .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 .Lsin_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 16 |
| .Lsin_cleanup: |
| add $stack_size, %rsp |
| ret |
| |
| .align 16 |
| .Lsin_naninf: |
| call fname_special |
| add $stack_size, %rsp |
| ret |
| |
| |
| |
