V10/asm.libm/sqrt.s

Compare this file to the similar file:
Show the results in this format:
# double sqrt(arg):revised July 18,1980
# double arg
# if(arg<0.0) { _errno=EDOM; return(-sqrt(-arg)) }
# W. Kahan's magic sqrt
# coded by Heidi Stettner

	.set	EDOM,98
	.text
	.align	1
	.globl	_sqrt
	.globl	dsqrt_r5
	.globl	_errno
_sqrt:
	.word	0x003c          # save r5,r4,r3,r2
	movd	4(ap),r0
	bsbb	dsqrt_r5
	ret
dsqrt_r5:
	movd	r0,r4
	jleq	nonpos		#argument is not positive
	movzwl	r4,r2
	ashl	$-1,r2,r0
	addw2	$0x203c,r0	#r0 has magic initial appx

# Do two steps of Heron's rule

	divf3	r0,r4,r2	
	addf2	r2,r0
	subw2	$0x80,r0

	divf3	r0,r4,r2
	addf2	r2,r0
	subw2	$0x80,r0


# Scale argument and approximation to prevent over/underflow
# NOTE: The following four steps would not be necessary if underflow
#       were gentle.

	bicw3	$0xffff807f,r4,r1
	subw2	$0x4080,r1		# r1 contains scaling factor
	subw2	r1,r4
	movl	r0,r2
	subw2	r1,r2

# Cubic step

	clrl	r1
	clrl	r3
	muld2	r0,r2
	subd2	r2,r4
	addw2	$0x100,r2
	addd2	r4,r2
	muld2	r0,r4
	divd2	r2,r4
	addw2	$0x80,r4
	addd2	r4,r0
	rsb
nonpos:
	jneq	negarg
	clrd	r0		#argument is zero
	rsb
negarg:
	movl	$EDOM,_errno
	mnegd	r4,-(sp)
	calls	$2,_sqrt
	mnegd	r0,r0		# returns -sqrt(-arg)
	ret