1# 2# Copyright (c) 1980 Regents of the University of California. 3# All rights reserved. The Berkeley software License Agreement 4# specifies the terms and conditions for redistribution. 5# 6# @(#)sqrt.s 5.1 (Berkeley) 05/08/85 7# 8# 9# double sqrt(arg):revised July 18,1980 10# double arg 11# if(arg<0.0) { _errno=EDOM; return(-sqrt(-arg)) } 12# W. Kahan's magic sqrt 13# coded by Heidi Stettner 14 15 .set EDOM,98 16 .text 17 .align 1 18 .globl _sqrt 19 .globl dsqrt_r5 20 .globl _errno 21_sqrt: 22 .word 0x003c # save r5,r4,r3,r2 23 bispsw $0xe0 24 movd 4(ap),r0 25 bsbb dsqrt_r5 26 ret 27dsqrt_r5: 28 movd r0,r4 29 jleq nonpos #argument is not positive 30 movzwl r4,r2 31 ashl $-1,r2,r0 32 addw2 $0x203c,r0 #r0 has magic initial appx 33 34# Do two steps of Heron's rule 35 36 divf3 r0,r4,r2 37 addf2 r2,r0 38 subw2 $0x80,r0 39 40 divf3 r0,r4,r2 41 addf2 r2,r0 42 subw2 $0x80,r0 43 44 45# Scale argument and approximation to prevent over/underflow 46# NOTE: The following four steps would not be necessary if underflow 47# were gentle. 48 49 bicw3 $0xffff807f,r4,r1 50 subw2 $0x4080,r1 # r1 contains scaling factor 51 subw2 r1,r4 52 movl r0,r2 53 subw2 r1,r2 54 55# Cubic step 56 57 clrl r1 58 clrl r3 59 muld2 r0,r2 60 subd2 r2,r4 61 addw2 $0x100,r2 62 addd2 r4,r2 63 muld2 r0,r4 64 divd2 r2,r4 65 addw2 $0x80,r4 66 addd2 r4,r0 67 rsb 68nonpos: 69 jneq negarg 70 clrd r0 #argument is zero 71 rsb 72negarg: 73 movl $EDOM,_errno 74 mnegd r4,-(sp) 75 calls $2,_sqrt 76 mnegd r0,r0 # returns -sqrt(-arg) 77 ret 78