1*20008Sdist# 2*20008Sdist# Copyright (c) 1980 Regents of the University of California. 3*20008Sdist# All rights reserved. The Berkeley software License Agreement 4*20008Sdist# specifies the terms and conditions for redistribution. 5*20008Sdist# 6*20008Sdist# @(#)sqrt.s 5.1 (Berkeley) 05/08/85 7*20008Sdist# 8*20008Sdist# 920000Sdist# double sqrt(arg):revised July 18,1980 1020000Sdist# double arg 1120000Sdist# if(arg<0.0) { _errno=EDOM; return(-sqrt(-arg)) } 1220000Sdist# W. Kahan's magic sqrt 1320000Sdist# coded by Heidi Stettner 1420000Sdist 1520000Sdist .set EDOM,98 1620000Sdist .text 1720000Sdist .align 1 1820000Sdist .globl _sqrt 1920000Sdist .globl dsqrt_r5 2020000Sdist .globl _errno 2120000Sdist_sqrt: 2220000Sdist .word 0x003c # save r5,r4,r3,r2 2320000Sdist bispsw $0xe0 2420000Sdist movd 4(ap),r0 2520000Sdist bsbb dsqrt_r5 2620000Sdist ret 2720000Sdistdsqrt_r5: 2820000Sdist movd r0,r4 2920000Sdist jleq nonpos #argument is not positive 3020000Sdist movzwl r4,r2 3120000Sdist ashl $-1,r2,r0 3220000Sdist addw2 $0x203c,r0 #r0 has magic initial appx 3320000Sdist 3420000Sdist# Do two steps of Heron's rule 3520000Sdist 3620000Sdist divf3 r0,r4,r2 3720000Sdist addf2 r2,r0 3820000Sdist subw2 $0x80,r0 3920000Sdist 4020000Sdist divf3 r0,r4,r2 4120000Sdist addf2 r2,r0 4220000Sdist subw2 $0x80,r0 4320000Sdist 4420000Sdist 4520000Sdist# Scale argument and approximation to prevent over/underflow 4620000Sdist# NOTE: The following four steps would not be necessary if underflow 4720000Sdist# were gentle. 4820000Sdist 4920000Sdist bicw3 $0xffff807f,r4,r1 5020000Sdist subw2 $0x4080,r1 # r1 contains scaling factor 5120000Sdist subw2 r1,r4 5220000Sdist movl r0,r2 5320000Sdist subw2 r1,r2 5420000Sdist 5520000Sdist# Cubic step 5620000Sdist 5720000Sdist clrl r1 5820000Sdist clrl r3 5920000Sdist muld2 r0,r2 6020000Sdist subd2 r2,r4 6120000Sdist addw2 $0x100,r2 6220000Sdist addd2 r4,r2 6320000Sdist muld2 r0,r4 6420000Sdist divd2 r2,r4 6520000Sdist addw2 $0x80,r4 6620000Sdist addd2 r4,r0 6720000Sdist rsb 6820000Sdistnonpos: 6920000Sdist jneq negarg 7020000Sdist clrd r0 #argument is zero 7120000Sdist rsb 7220000Sdistnegarg: 7320000Sdist movl $EDOM,_errno 7420000Sdist mnegd r4,-(sp) 7520000Sdist calls $2,_sqrt 7620000Sdist mnegd r0,r0 # returns -sqrt(-arg) 7720000Sdist ret 78