
#
#  (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/>.
#
#


#
# __vrs4_expf.s
#
# A vector implementation of the expf libm function.
#  This routine implemented in single precision.  It is slightly
#  less accurate than the double precision version, but it will
#  be better for vectorizing.
#
# Prototype:
#
#     __m128 __vrs4_expf(__m128 x);
#
#   Computes e raised to the x power for 4 floats at a time.
# Does not perform error handling, but does return C99 values for error
# inputs.   Denormal results are truncated to 0.
#
#

#ifdef __ELF__
.section .note.GNU-stack,"",@progbits
#endif

        .text
        .align 16
        .p2align 4,,15

# define local variable storage offsets
.equ	p_temp,0	# temporary for get/put bits operation
.equ	p_ux,0x10		# local storage for ux array
.equ	p_m,0x20		# local storage for m array
.equ	p_j,0x30		# local storage for m array
.equ    save_rbx,0x040          #qword
.equ	stack_size,0x48


.globl __vrs4_expf
       .type   __vrs4_expf,@function
__vrs4_expf:
	sub	$stack_size,%rsp
        mov             %rbx,save_rbx(%rsp)


	movaps	 %xmm0,p_ux(%rsp)
        maxps   .L__real_m8192(%rip),%xmm0   # protect against small input values


#        /* Find m, z1 and z2 such that exp(x) = 2**m * (z1 + z2) */
#      Step 1. Reduce the argument.
	#    r = x * thirtytwo_by_logbaseof2;
	movaps	.L__real_thirtytwo_by_log2(%rip),%xmm2	#
	mulps	%xmm0,%xmm2
	xor		%rax,%rax
        minps   .L__real_8192(%rip),%xmm2   # protect against large input values

#    /* Set n = nearest integer to r */
	cvtps2dq	%xmm2,%xmm3
	lea		.L__two_to_jby32_table(%rip),%rdi
	cvtdq2ps	%xmm3,%xmm1


#    r1 = x - n * logbaseof2_by_32_lead;
	movaps	.L__real_log2_by_32_head(%rip),%xmm2
	mulps	%xmm1,%xmm2
	subps	%xmm2,%xmm0				# r1 in xmm0,

#    r2 = - n * logbaseof2_by_32_lead;
	mulps	.L__real_log2_by_32_tail(%rip),%xmm1

#    j = n & 0x0000001f;
	movdqa	%xmm3,%xmm4
	movdqa	.L__int_mask_1f(%rip),%xmm2
	pand	%xmm4,%xmm2
	movdqa	%xmm2,p_j(%rsp)
#    f1 = two_to_jby32_lead_table[j];

#    *m = (n - j) / 32;
	psubd	%xmm2,%xmm4
	psrad	$5,%xmm4
	movdqa	%xmm4,p_m(%rsp)

	movaps	%xmm0,%xmm3
	addps	%xmm1,%xmm3

	mov		p_j(%rsp),%eax 			# get an individual index
	mov		(%rdi,%rax,4),%edx		# get the f1 value
	mov		%edx,p_j(%rsp) 			# save the f1 value

#      Step 2. Compute the polynomial.
#     q = r1 +
#              r*r*( 5.00000000000000008883e-01 +
#                      r*( 1.66666666665260878863e-01 +
#                      r*( 4.16666666662260795726e-02 +
#                      r*( 8.33336798434219616221e-03 +
#                      r*( 1.38889490863777199667e-03 )))));
#    q = r + r^2/2 + r^3/6 + r^4/24 + r^5/120 + r^6/720
#    q = r + r^2/2 + r^3/6 + r^4/24 good enough for single precision
	movaps	%xmm3,%xmm4
	movaps	%xmm3,%xmm2			# x*x
	mulps	%xmm2,%xmm2
	mulps	.L__real_1_24(%rip),%xmm4	# /24

	mov		p_j+4(%rsp),%eax 			# get an individual index
	mov		(%rdi,%rax,4),%edx		# get the f1 value
	mov		%edx,p_j+4(%rsp) 			# save the f1 value

	addps 	.L__real_1_6(%rip),%xmm4	# +1/6

	mulps	%xmm2,%xmm3			# x^3
	mov		p_j+8(%rsp),%eax 			# get an individual index
	mov		(%rdi,%rax,4),%edx		# get the f1 value
	mov		%edx,p_j+8(%rsp) 			# save the f1 value
	mulps	.L__real_half(%rip),%xmm2	# x^2/2
	mulps	%xmm3,%xmm4			# *x^3

	mov		p_j+12(%rsp),%eax 			# get an individual index
	mov		(%rdi,%rax,4),%edx		# get the f1 value
	mov		%edx,p_j+12(%rsp) 			# save the f1 value
	addps	%xmm4,%xmm1			# +r2

	addps	%xmm2,%xmm1			# + x^2/2
	addps	%xmm1,%xmm0			# +r1

# deal with infinite or denormal results
        movdqa  p_m(%rsp),%xmm1
        movdqa  p_m(%rsp),%xmm2
        pcmpgtd .L__int_127(%rip),%xmm2
        pminsw  .L__int_128(%rip),%xmm1 # ceil at 128
        movmskps        %xmm2,%eax
        test            $0x0f,%eax

        paddd   .L__int_127(%rip),%xmm1 # add bias

#    *z2 = f2 + ((f1 + f2) * q);
        mulps   p_j(%rsp),%xmm0         # * f1
        addps   p_j(%rsp),%xmm0         # + f1
        jnz             .L__exp_largef
.L__check1:

	pxor	%xmm2,%xmm2			# floor at 0
	pmaxsw	%xmm2,%xmm1

	pslld	$23,%xmm1			# build 2^n

	movaps	%xmm1,%xmm2


# check for infinity or nan
	movaps	p_ux(%rsp),%xmm1
	andps	.L__real_infinity(%rip),%xmm1
	cmpps	$0,.L__real_infinity(%rip),%xmm1
	movmskps	%xmm1,%eax
	test		$0xf,%eax


# end of splitexp
#        /* Scale (z1 + z2) by 2.0**m */
#      Step 3. Reconstitute.

	mulps	%xmm2,%xmm0		# result *= 2^n

# we'd like to avoid a branch, and can use cmp's and and's to
# eliminate them.  But it adds cycles for normal cases
# to handle events that are supposed to be exceptions.
#  Using this branch with the
# check above results in faster code for the normal cases.
# And branch mispredict penalties should only come into
# play for nans and infinities.
	jnz		.L__exp_naninf

#
.L__final_check:
        mov             save_rbx(%rsp),%rbx             # restore rbx
	add		$stack_size,%rsp
	ret

#  deal with nans and infinities

.L__exp_naninf:
	movaps	%xmm0,p_temp(%rsp)		# save the computed values
	mov		%eax,%ecx				# save the error mask
	test	$1,%ecx					# first value?
	jz		.__Lni2
	mov		p_ux(%rsp),%edx			# get the input
	call	.L__naninf
	mov		%edx,p_temp(%rsp)		# save the new result
.__Lni2:
	test	$2,%ecx					# first value?
	jz		.__Lni3
	mov		p_ux+4(%rsp),%edx		# get the input
	call	.L__naninf
	mov		%edx,p_temp+4(%rsp)		# save the new result
.__Lni3:
	test	$4,%ecx					# first value?
	jz		.__Lni4
	mov		p_ux+8(%rsp),%edx		# get the input
	call	.L__naninf
	mov		%edx,p_temp+8(%rsp)		# save the new result
.__Lni4:
	test	$8,%ecx					# first value?
	jz		.__Lnie
	mov		p_ux+12(%rsp),%edx		# get the input
	call	.L__naninf
	mov		%edx,p_temp+12(%rsp)		# save the new result
.__Lnie:
	movaps	p_temp(%rsp),%xmm0		# get the answers
	jmp		.L__final_check


# a simple subroutine to check a scalar input value for infinity
# or NaN and return the correct result
# expects input in edx, and returns value in edx.  Destroys eax.
.L__naninf:
	mov	$0x0007FFFFF,%eax
	test	%eax,%edx
	jnz		.L__enan		# jump if mantissa not zero, so it's a NaN
# inf
	mov		%edx,%eax
	rcl		$1,%eax
	jnc		.L__r				# exp(+inf) = inf
	xor		%edx,%edx				# exp(-inf) = 0
	jmp		.L__r

#NaN
.L__enan:
	mov		$0x000400000,%eax	# convert to quiet
	or		%eax,%edx
.L__r:
	ret

        .align 16
#  deal with m > 127.  In some instances, rounding during calculations
#  can result in infinity when it shouldn't.  For these cases, we scale
#  m down, and scale the mantissa up.

.L__exp_largef:
        movdqa    %xmm0,p_temp(%rsp)    # save the mantissa portion
        movdqa    %xmm1,p_m(%rsp)       # save the exponent portion
        mov       %eax,%ecx              # save the error mask
        test    $1,%ecx                  # first value?
        jz       .L__Lf2
        mov      p_m(%rsp),%edx # get the exponent
        sub      $1,%edx                # scale it down
        mov      %edx,p_m(%rsp)       # save the exponent
        movss   p_temp(%rsp),%xmm3     # get the mantissa
        mulss   .L__real_two(%rip),%xmm3        # scale it up
        movss    %xmm3,p_temp(%rsp)   # save the mantissa
.L__Lf2:
        test    $2,%ecx                 # second value?
        jz       .L__Lf3
        mov      p_m+4(%rsp),%edx # get the exponent
        sub      $1,%edx                # scale it down
        mov      %edx,p_m+4(%rsp)       # save the exponent
        movss   p_temp+4(%rsp),%xmm3     # get the mantissa
        mulss   .L__real_two(%rip),%xmm3        # scale it up
        movss     %xmm3,p_temp+4(%rsp)   # save the mantissa
.L__Lf3:
        test    $4,%ecx                 # third value?
        jz       .L__Lf4
        mov      p_m+8(%rsp),%edx # get the exponent
        sub      $1,%edx                # scale it down
        mov      %edx,p_m+8(%rsp)       # save the exponent
        movss   p_temp+8(%rsp),%xmm3     # get the mantissa
        mulss   .L__real_two(%rip),%xmm3        # scale it up
        movss    %xmm3,p_temp+8(%rsp)   # save the mantissa
.L__Lf4:
        test    $8,%ecx                                 # fourth value?
        jz       .L__Lfe
        mov      p_m+12(%rsp),%edx        # get the exponent
        sub      $1,%edx                 # scale it down
        mov      %edx,p_m+12(%rsp)      # save the exponent
        movss   p_temp+12(%rsp),%xmm3    # get the mantissa
        mulss   .L__real_two(%rip),%xmm3        # scale it up
        movss     %xmm3,p_temp+12(%rsp)  # save the mantissa
.L__Lfe:
        movaps  p_temp(%rsp),%xmm0      # restore the mantissa portion back
        movdqa  p_m(%rsp),%xmm1         # restore the exponent portion
        jmp             .L__check1

	.data
        .align 64

.L__real_half:			.long 0x3f000000	# 1/2
				.long 0x3f000000
				.long 0x3f000000
				.long 0x3f000000

.L__real_two:			.long 0x40000000	# 2
				.long 0x40000000
				.long 0x40000000
				.long 0x40000000

.L__real_8192:			.long 0x46000000	# 8192, to protect against really large numbers
				.long 0x46000000
				.long 0x46000000
				.long 0x46000000
.L__real_m8192:			.long 0xC6000000	# -8192, to protect against really small numbers
				.long 0xC6000000
				.long 0xC6000000
				.long 0xC6000000

.L__real_thirtytwo_by_log2: 	.long 0x4238AA3B	# thirtytwo_by_log2
				.long 0x4238AA3B
				.long 0x4238AA3B
				.long 0x4238AA3B

.L__real_log2_by_32:		.long 0x3CB17218	# log2_by_32
				.long 0x3CB17218
				.long 0x3CB17218
				.long 0x3CB17218

.L__real_log2_by_32_head:	.long 0x3CB17000	# log2_by_32
				.long 0x3CB17000
				.long 0x3CB17000
				.long 0x3CB17000

.L__real_log2_by_32_tail:	.long 0xB585FDF4	# log2_by_32
				.long 0xB585FDF4
				.long 0xB585FDF4
				.long 0xB585FDF4

.L__real_1_6:			.long 0x3E2AAAAB	# 0.16666666666 used in polynomial
				.long 0x3E2AAAAB
				.long 0x3E2AAAAB
				.long 0x3E2AAAAB

.L__real_1_24:			.long 0x3D2AAAAB	# 0.041666668 used in polynomial
				.long 0x3D2AAAAB
				.long 0x3D2AAAAB
				.long 0x3D2AAAAB

.L__real_infinity:		.long 0x7f800000	# infinity
				.long 0x7f800000
				.long 0x7f800000
				.long 0x7f800000
.L__int_mask_1f:		.long 0x00000001f
				.long 0x00000001f
				.long 0x00000001f
				.long 0x00000001f
.L__int_128:			.long 0x000000080
				.long 0x000000080
				.long 0x000000080
				.long 0x000000080
.L__int_127:			.long 0x00000007f
				.long 0x00000007f
				.long 0x00000007f
				.long 0x00000007f


.L__two_to_jby32_table:
   .long  0x3F800000  #          1
   .long  0x3F82CD87  #  1.0218972
   .long  0x3F85AAC3  #  1.0442737
   .long  0x3F88980F  #  1.0671405
   .long  0x3F8B95C2  #  1.0905077
   .long  0x3F8EA43A  #  1.1143868
   .long  0x3F91C3D3  #  1.1387886
   .long  0x3F94F4F0  #  1.1637249
   .long  0x3F9837F0  #  1.1892071
   .long  0x3F9B8D3A  #  1.2152474
   .long  0x3F9EF532  #  1.2418578
   .long  0x3FA27043  #   1.269051
   .long  0x3FA5FED7  #  1.2968396
   .long  0x3FA9A15B  #  1.3252367
   .long  0x3FAD583F  #  1.3542556
   .long  0x3FB123F6  #  1.3839099
   .long  0x3FB504F3  #  1.4142135
   .long  0x3FB8FBAF  #  1.4451808
   .long  0x3FBD08A4  #  1.4768262
   .long  0x3FC12C4D  #  1.5091645
   .long  0x3FC5672A  #  1.5422108
   .long  0x3FC9B9BE  #  1.5759809
   .long  0x3FCE248C  #  1.6104903
   .long  0x3FD2A81E  #  1.6457555
   .long  0x3FD744FD  #  1.6817929
   .long  0x3FDBFBB8  #  1.7186193
   .long  0x3FE0CCDF  #  1.7562522
   .long  0x3FE5B907  #  1.7947091
   .long  0x3FEAC0C7  #  1.8340081
   .long  0x3FEFE4BA  #  1.8741677
   .long  0x3FF5257D  #  1.9152066
   .long  0x3FFA83B3  #  1.9571441
   .long 0					# for alignment