fp.s

早期freebsd实现

	.word	cmp_d		# 55
	.word	ill		# 55
	.word	ill		# 55
	.word	ill		# 55
	.word	ill		# 55
	.word	ill		# 55
	.word	ill		# 55
	.word	cmp_s		# 56
	.word	cmp_d		# 56
	.word	ill		# 56
	.word	ill		# 56
	.word	ill		# 56
	.word	ill		# 56
	.word	ill		# 56
	.word	ill		# 56
	.word	cmp_s		# 57
	.word	cmp_d		# 57
	.word	ill		# 57
	.word	ill		# 57
	.word	ill		# 57
	.word	ill		# 57
	.word	ill		# 57
	.word	ill		# 57
	.word	cmp_s		# 58
	.word	cmp_d		# 58
	.word	ill		# 58
	.word	ill		# 58
	.word	ill		# 58
	.word	ill		# 58
	.word	ill		# 58
	.word	ill		# 58
	.word	cmp_s		# 59
	.word	cmp_d		# 59
	.word	ill		# 59
	.word	ill		# 59
	.word	ill		# 59
	.word	ill		# 59
	.word	ill		# 59
	.word	ill		# 59
	.word	cmp_s		# 60
	.word	cmp_d		# 60
	.word	ill		# 60
	.word	ill		# 60
	.word	ill		# 60
	.word	ill		# 60
	.word	ill		# 60
	.word	ill		# 60
	.word	cmp_s		# 61
	.word	cmp_d		# 61
	.word	ill		# 61
	.word	ill		# 61
	.word	ill		# 61
	.word	ill		# 61
	.word	ill		# 61
	.word	ill		# 61
	.word	cmp_s		# 62
	.word	cmp_d		# 62
	.word	ill		# 62
	.word	ill		# 62
	.word	ill		# 62
	.word	ill		# 62
	.word	ill		# 62
	.word	ill		# 62
	.word	cmp_s		# 63
	.word	cmp_d		# 63
	.word	ill		# 63
	.word	ill		# 63
	.word	ill		# 63
	.word	ill		# 63
	.word	ill		# 63
	.word	ill		# 63
	.text

/*
 * Single precision subtract.
 */
sub_s:
	jal	get_ft_fs_s
	xor	t4, t4, 1			# negate FT sign bit
	b	add_sub_s
/*
 * Single precision add.
 */
add_s:
	jal	get_ft_fs_s
add_sub_s:
	bne	t1, SEXP_INF, 1f		# is FS an infinity?
	bne	t5, SEXP_INF, result_fs_s	# if FT is not inf, result=FS
	bne	t2, zero, result_fs_s		# if FS is NAN, result is FS
	bne	t6, zero, result_ft_s		# if FT is NAN, result is FT
	bne	t0, t4, invalid_s		# both infinities same sign?
	b	result_fs_s			# result is in FS
1:
	beq	t5, SEXP_INF, result_ft_s	# if FT is inf, result=FT
	bne	t1, zero, 4f			# is FS a denormalized num?
	beq	t2, zero, 3f			# is FS zero?
	bne	t5, zero, 2f			# is FT a denormalized num?
	beq	t6, zero, result_fs_s		# FT is zero, result=FS
	jal	renorm_fs_s
	jal	renorm_ft_s
	b	5f
2:
	jal	renorm_fs_s
	subu	t5, t5, SEXP_BIAS		# unbias FT exponent
	or	t6, t6, SIMPL_ONE		# set implied one bit
	b	5f
3:
	bne	t5, zero, result_ft_s		# if FT != 0, result=FT
	bne	t6, zero, result_ft_s
	and	v0, a1, MACH_FPC_ROUNDING_BITS	# get rounding mode
	bne	v0, MACH_FPC_ROUND_RM, 1f	# round to -infinity?
	or	t0, t0, t4			# compute result sign
	b	result_fs_s
1:
	and	t0, t0, t4			# compute result sign
	b	result_fs_s
4:
	bne	t5, zero, 2f			# is FT a denormalized num?
	beq	t6, zero, result_fs_s		# FT is zero, result=FS
	subu	t1, t1, SEXP_BIAS		# unbias FS exponent
	or	t2, t2, SIMPL_ONE		# set implied one bit
	jal	renorm_ft_s
	b	5f
2:
	subu	t1, t1, SEXP_BIAS		# unbias FS exponent
	or	t2, t2, SIMPL_ONE		# set implied one bit
	subu	t5, t5, SEXP_BIAS		# unbias FT exponent
	or	t6, t6, SIMPL_ONE		# set implied one bit
/*
 * Perform the addition.
 */
5:
	move	t8, zero			# no shifted bits (sticky reg)
	beq	t1, t5, 4f			# no shift needed
	subu	v0, t1, t5			# v0 = difference of exponents
	move	v1, v0				# v1 = abs(difference)
	bge	v0, zero, 1f
	negu	v1
1:
	ble	v1, SFRAC_BITS+2, 2f		# is difference too great?
	li	t8, STICKYBIT			# set the sticky bit
	bge	v0, zero, 1f			# check which exp is larger
	move	t1, t5				# result exp is FTs
	move	t2, zero			# FSs fraction shifted is zero
	b	4f
1:
	move	t6, zero			# FTs fraction shifted is zero
	b	4f
2:
	li	t9, 32				# compute 32 - abs(exp diff)
	subu	t9, t9, v1
	bgt	v0, zero, 3f			# if FS > FT, shift FTs frac
	move	t1, t5				# FT > FS, result exp is FTs
	sll	t8, t2, t9			# save bits shifted out
	srl	t2, t2, v1			# shift FSs fraction
	b	4f
3:
	sll	t8, t6, t9			# save bits shifted out
	srl	t6, t6, v1			# shift FTs fraction
4:
	bne	t0, t4, 1f			# if signs differ, subtract
	addu	t2, t2, t6			# add fractions
	b	norm_s
1:
	blt	t2, t6, 3f			# subtract larger from smaller
	bne	t2, t6, 2f			# if same, result=0
	move	t1, zero			# result=0
	move	t2, zero
	and	v0, a1, MACH_FPC_ROUNDING_BITS	# get rounding mode
	bne	v0, MACH_FPC_ROUND_RM, 1f	# round to -infinity?
	or	t0, t0, t4			# compute result sign
	b	result_fs_s
1:
	and	t0, t0, t4			# compute result sign
	b	result_fs_s
2:
	sltu	t9, zero, t8			# compute t2:zero - t6:t8
	subu	t8, zero, t8
	subu	t2, t2, t6			# subtract fractions
	subu	t2, t2, t9			# subtract barrow
	b	norm_s
3:
	move	t0, t4				# sign of result = FTs
	sltu	t9, zero, t8			# compute t6:zero - t2:t8
	subu	t8, zero, t8
	subu	t2, t6, t2			# subtract fractions
	subu	t2, t2, t9			# subtract barrow
	b	norm_s

/*
 * Double precision subtract.
 */
sub_d:
	jal	get_ft_fs_d
	xor	t4, t4, 1			# negate sign bit
	b	add_sub_d
/*
 * Double precision add.
 */
add_d:
	jal	get_ft_fs_d
add_sub_d:
	bne	t1, DEXP_INF, 1f		# is FS an infinity?
	bne	t5, DEXP_INF, result_fs_d	# if FT is not inf, result=FS
	bne	t2, zero, result_fs_d		# if FS is NAN, result is FS
	bne	t3, zero, result_fs_d
	bne	t6, zero, result_ft_d		# if FT is NAN, result is FT
	bne	t7, zero, result_ft_d
	bne	t0, t4, invalid_d		# both infinities same sign?
	b	result_fs_d			# result is in FS
1:
	beq	t5, DEXP_INF, result_ft_d	# if FT is inf, result=FT
	bne	t1, zero, 4f			# is FS a denormalized num?
	bne	t2, zero, 1f			# is FS zero?
	beq	t3, zero, 3f
1:
	bne	t5, zero, 2f			# is FT a denormalized num?
	bne	t6, zero, 1f
	beq	t7, zero, result_fs_d		# FT is zero, result=FS
1:
	jal	renorm_fs_d
	jal	renorm_ft_d
	b	5f
2:
	jal	renorm_fs_d
	subu	t5, t5, DEXP_BIAS		# unbias FT exponent
	or	t6, t6, DIMPL_ONE		# set implied one bit
	b	5f
3:
	bne	t5, zero, result_ft_d		# if FT != 0, result=FT
	bne	t6, zero, result_ft_d
	bne	t7, zero, result_ft_d
	and	v0, a1, MACH_FPC_ROUNDING_BITS	# get rounding mode
	bne	v0, MACH_FPC_ROUND_RM, 1f	# round to -infinity?
	or	t0, t0, t4			# compute result sign
	b	result_fs_d
1:
	and	t0, t0, t4			# compute result sign
	b	result_fs_d
4:
	bne	t5, zero, 2f			# is FT a denormalized num?
	bne	t6, zero, 1f
	beq	t7, zero, result_fs_d		# FT is zero, result=FS
1:
	subu	t1, t1, DEXP_BIAS		# unbias FS exponent
	or	t2, t2, DIMPL_ONE		# set implied one bit
	jal	renorm_ft_d
	b	5f
2:
	subu	t1, t1, DEXP_BIAS		# unbias FS exponent
	or	t2, t2, DIMPL_ONE		# set implied one bit
	subu	t5, t5, DEXP_BIAS		# unbias FT exponent
	or	t6, t6, DIMPL_ONE		# set implied one bit
/*
 * Perform the addition.
 */
5:
	move	t8, zero			# no shifted bits (sticky reg)
	beq	t1, t5, 4f			# no shift needed
	subu	v0, t1, t5			# v0 = difference of exponents
	move	v1, v0				# v1 = abs(difference)
	bge	v0, zero, 1f
	negu	v1
1:
	ble	v1, DFRAC_BITS+2, 2f		# is difference too great?
	li	t8, STICKYBIT			# set the sticky bit
	bge	v0, zero, 1f			# check which exp is larger
	move	t1, t5				# result exp is FTs
	move	t2, zero			# FSs fraction shifted is zero
	move	t3, zero
	b	4f
1:
	move	t6, zero			# FTs fraction shifted is zero
	move	t7, zero
	b	4f
2:
	li	t9, 32
	bge	v0, zero, 3f			# if FS > FT, shift FTs frac
	move	t1, t5				# FT > FS, result exp is FTs
	blt	v1, t9, 1f			# shift right by < 32?
	subu	v1, v1, t9
	subu	t9, t9, v1
	sll	t8, t2, t9			# save bits shifted out
	sltu	t9, zero, t3			# dont lose any one bits
	or	t8, t8, t9			# save sticky bit
	srl	t3, t2, v1			# shift FSs fraction
	move	t2, zero
	b	4f
1:
	subu	t9, t9, v1
	sll	t8, t3, t9			# save bits shifted out
	srl	t3, t3, v1			# shift FSs fraction
	sll	t9, t2, t9			# save bits shifted out of t2
	or	t3, t3, t9			# and put into t3
	srl	t2, t2, v1
	b	4f
3:
	blt	v1, t9, 1f			# shift right by < 32?
	subu	v1, v1, t9
	subu	t9, t9, v1
	sll	t8, t6, t9			# save bits shifted out
	srl	t7, t6, v1			# shift FTs fraction
	move	t6, zero
	b	4f
1:
	subu	t9, t9, v1
	sll	t8, t7, t9			# save bits shifted out
	srl	t7, t7, v1			# shift FTs fraction
	sll	t9, t6, t9			# save bits shifted out of t2
	or	t7, t7, t9			# and put into t3
	srl	t6, t6, v1
4:
	bne	t0, t4, 1f			# if signs differ, subtract
	addu	t3, t3, t7			# add fractions
	sltu	t9, t3, t7			# compute carry
	addu	t2, t2, t6			# add fractions
	addu	t2, t2, t9			# add carry
	b	norm_d
1:
	blt	t2, t6, 3f			# subtract larger from smaller
	bne	t2, t6, 2f
	bltu	t3, t7, 3f
	bne	t3, t7, 2f			# if same, result=0
	move	t1, zero			# result=0
	move	t2, zero
	move	t3, zero
	and	v0, a1, MACH_FPC_ROUNDING_BITS	# get rounding mode
	bne	v0, MACH_FPC_ROUND_RM, 1f	# round to -infinity?
	or	t0, t0, t4			# compute result sign
	b	result_fs_d
1:
	and	t0, t0, t4			# compute result sign
	b	result_fs_d
2:
	beq	t8, zero, 1f			# compute t2:t3:zero - t6:t7:t8
	subu	t8, zero, t8
	sltu	v0, t3, 1			# compute barrow out
	subu	t3, t3, 1			# subtract barrow
	subu	t2, t2, v0
1:
	sltu	v0, t3, t7
	subu	t3, t3, t7			# subtract fractions
	subu	t2, t2, t6			# subtract fractions
	subu	t2, t2, v0			# subtract barrow
	b	norm_d
3:
	move	t0, t4				# sign of result = FTs
	beq	t8, zero, 1f			# compute t6:t7:zero - t2:t3:t8
	subu	t8, zero, t8
	sltu	v0, t7, 1			# compute barrow out
	subu	t7, t7, 1			# subtract barrow
	subu	t6, t6, v0
1:
	sltu	v0, t7, t3
	subu	t3, t7, t3			# subtract fractions
	subu	t2, t6, t2			# subtract fractions
	subu	t2, t2, v0			# subtract barrow
	b	norm_d

/*
 * Single precision multiply.
 */
mul_s:
	jal	get_ft_fs_s
	xor	t0, t0, t4			# compute sign of result
	move	t4, t0
	bne	t1, SEXP_INF, 2f		# is FS an infinity?
	bne	t2, zero, result_fs_s		# if FS is a NAN, result=FS
	bne	t5, SEXP_INF, 1f		# FS is inf, is FT an infinity?
	bne	t6, zero, result_ft_s		# if FT is a NAN, result=FT
	b	result_fs_s			# result is infinity
1:
	bne	t5, zero, result_fs_s		# inf * zero? if no, result=FS
	bne	t6, zero, result_fs_s
	b	invalid_s			# infinity * zero is invalid
2:
	bne	t5, SEXP_INF, 1f		# FS != inf, is FT an infinity?
	bne	t1, zero, result_ft_s		# zero * inf? if no, result=FT
	bne	t2, zero, result_ft_s
	bne	t6, zero, result_ft_s		# if FT is a NAN, result=FT
	b	invalid_s			# zero * infinity is invalid
1:
	bne	t1, zero, 1f			# is FS zero?
	beq	t2, zero, result_fs_s		# result is zero
	jal	renorm_fs_s
	b	2f
1:
	subu	t1, t1, SEXP_BIAS		# unbias FS exponent
	or	t2, t2, SIMPL_ONE		# set implied one bit
2:
	bne	t5, zero, 1f			# is FT zero?
	beq	t6, zero, result_ft_s		# result is zero
	jal	renorm_ft_s
	b	2f
1:
	subu	t5, t5, SEXP_BIAS		# unbias FT exponent
	or	t6, t6, SIMPL_ONE		# set implied one bit
2:
	addu	t1, t1, t5			# compute result exponent
	addu	t1, t1, 9			# account for binary point
	multu	t2, t6				# multiply fractions
	mflo	t8
	mfhi	t2
	b	norm_s

/*
 * Double precision multiply.
 */
mul_d:
	jal	get_ft_fs_d
	xor	t0, t0, t4			# compute sign of result
	move	t4, t0
	bne	t1, DEXP_INF, 2f		# is FS an infinity?
	bne	t2, zero, result_fs_d		# if FS is a NAN, result=FS
	bne	t3, zero, result_fs_d
	bne	t5, DEXP_INF, 1f		# FS is inf, is FT an infinity?
	bne	t6, zero, result_ft_d		# if FT is a NAN, result=FT
	bne	t7, zero, result_ft_d
	b	result_fs_d			# result is infinity
1:
	bne	t5, zero, result_fs_d		# inf * zero? if no, result=FS
	bne	t6, zero, result_fs_d
	bne	t7, zero, result_fs_d
	b	invalid_d			# infinity * zero is invalid
2:
	bne	t5, DEXP_INF, 1f		# FS != inf, is FT an infinity?
	bne	t1, zero, result_ft_d		# zero * inf? if no, result=FT
	bne	t2, zero, result_ft_d		# if FS is a NAN, result=FS
	bne	t3, zero, result_ft_d
	bne	t6, zero, result_ft_d		# if FT is a NAN, result=FT
	bne	t7, zero, result_ft_d
	b	invalid_d			# zero * infinity is invalid
1:
	bne	t1, zero, 2f			# is FS zero?
	bne	t2, zero, 1f
	beq	t3, zero, result_fs_d		# result is zero
1:
	jal	renorm_fs_d
	b	3f
2:
	subu	t1, t1, DEXP_BIAS		# unbias FS exponent
	or	t2, t2, DIMPL_ONE		# set implied one bit
3:
	bne	t5, zero, 2f			# is FT zero?
	bne	t6, zero, 1f
	beq	t7, zero, result_ft_d		# result is zero
1:
	jal	renorm_ft_d
	b	3f
2:
	subu	t5, t5, DEXP_BIAS		# unbias FT exponent
	or	t6, t6, DIMPL_ONE		# set implied one bit
3:
	addu	t1, t1, t5			# compute result exponent
	addu	t1, t1, 12			# ???
	multu	t3, t7				# multiply fractions (low * low)
	move	t4, t2				# free up t2,t3 for result
	move	t5, t3
	mflo	a3				# save low order bits
	mfhi	t8
	not	v0, t8
	multu	t4, t7				# multiply FS(high) * FT(low)
	mflo	v1
	mfhi	t3				# init low result
	sltu	v0, v0, v1			# compute carry
	addu	t8, v1
	multu	t5, t6				# multiply FS(low) * FT(high)
	addu	t3, t3, v0			# add carry
	not	v0, t8
	mflo	v1
	mfhi	t2
	sltu	v0, v0, v1
	addu	t8, v1
	multu	t4, t6				# multiply FS(high) * FT(high)
	addu	t3, v0
	not	v1, t3
	sltu	v1, v1, t2
	addu	t3, t2
	not	v0, t3
	mfhi	t2
	addu	t2, v1
	mflo	v1
	sltu	v0, v0, v1
	addu	t2, v0
	addu	t3, v1
	sltu	a3, zero, a3			# reduce t8,a3 to just t8
	or	t8, a3
	b	norm_d

/*
 * Single precision divide.
 */
div_s:
	jal	get_ft_fs_s
	xor	t0, t0, t4			# compute sign of result
	move	t4, t0
	bne	t1, SEXP_INF, 1f		# is FS an infinity?
	bne	t2, zero, result_fs_s		# if FS is NAN, result is FS
	bne	t5, SEXP_INF, result_fs_s	# is FT an infinity?
	bne	t6, zero, result_ft_s		# if FT is NAN, result is FT
	b	invalid_s			# infinity/infinity is invalid
1:
	bne	t5, SEXP_INF, 1f		# is FT an infinity?
	bne	t6, zero, result_ft_s		# if FT is NAN, result is FT
	move	t1, zero			# x / infinity is zero
	move	t2, zero
	b	result_fs_s
1:
	bne	t1, zero, 2f			# is FS zero?
	bne	t2, zero, 1f
	bne	t5, zero, result_fs_s		# FS=zero, is FT zero?
	beq	t6, zero, invalid_s		# 0 / 0
	b	result_fs_s			# result = zero
1:
	jal	renorm_fs_s
	b	3f
2:
	subu	t1, t1, SEXP_BIAS		# unbias FS exponent
	or	t2, t2, SIMPL_ONE		# set implied one bit
3:
	bne	t5, zero, 2f			# is FT zero?
	bne	t6, zero, 1f
	or	a1, a1, MACH_FPC_EXCEPTION_DIV0 | MACH_FPC_STICKY_DIV0
	and	v0, a1, MACH_FPC_ENABLE_DIV0	# trap enabled?
	bne	v0, zero, fpe_trap
	ctc1	a1, MACH_FPC_CSR		# save exceptions
	li	t1, SEXP_INF			# result is infinity
	move	t2, zero
	b	result_fs_s
1:
	jal	renorm_ft_s
	b	3f
2:
	subu	t5, t5, SEXP_BIAS		# unbias FT exponent
	or	t6, t6, SIMPL_ONE		# set implied one bit
3:
	subu	t1, t1, t5			# compute exponent
	subu	t1, t1, 3			# compensate for result position
	li	v0, SFRAC_BITS+3		# number of bits to divide
	move	t8, t2				# init dividend
	move	t2, zero			# init result
1:
	bltu	t8, t6, 3f			# is dividend >= divisor?
2:
	subu	t8, t8, t6			# subtract divisor from dividend
	or	t2, t2, 1			# remember that we did
	bne	t8, zero, 3f			# if not done, continue
	sll	t2, t2, v0			# shift result to final position
	b	norm_s
3:
	sll	t8, t8, 1			# shift dividend
	sll	t2, t2, 1			# shift result
	subu	v0, v0, 1			# are we done?
	bne	v0, zero, 1b			# no, continue
	b	norm_s

/*
 * Double precision divide.
 */