124567Szliu# Copyright (c) 1985 Regents of the University of California. 2*34125Sbostic# All rights reserved. 324567Szliu# 4*34125Sbostic# Redistribution and use in source and binary forms are permitted 5*34125Sbostic# provided that this notice is preserved and that due credit is given 6*34125Sbostic# to the University of California at Berkeley. The name of the University 7*34125Sbostic# may not be used to endorse or promote products derived from this 8*34125Sbostic# software without specific prior written permission. This software 9*34125Sbostic# is provided ``as is'' without express or implied warranty. 10*34125Sbostic# 11*34125Sbostic# All recipients should regard themselves as participants in an ongoing 12*34125Sbostic# research project and hence should feel obligated to report their 13*34125Sbostic# experiences (good or bad) with these elementary function codes, using 14*34125Sbostic# the sendbug(8) program, to the authors. 15*34125Sbostic# 16*34125Sbostic# @(#)cabs.s 5.2 (Berkeley) 04/29/88 17*34125Sbostic# 1824729Selefunt .data 1924729Selefunt .align 2 2024729Selefunt_sccsid: 21*34125Sbostic.asciz "@(#)cabs.s 1.2 (Berkeley) 8/21/85; 5.2 (ucb.elefunt) 04/29/88" 2224567Szliu 2324567Szliu# double precision complex absolute value 2424567Szliu# CABS by W. Kahan, 9/7/80. 2524567Szliu# Revised for reserved operands by E. LeBlanc, 8/18/82 2624567Szliu# argument for complex absolute value by reference, *4(ap) 2724567Szliu# argument for cabs and hypot (C fcns) by value, 4(ap) 2824567Szliu# output is in r0:r1 (error less than 0.86 ulps) 2924567Szliu 3024567Szliu .text 3124567Szliu .align 1 3224567Szliu .globl _cabs 3324567Szliu .globl _hypot 3424567Szliu .globl _z_abs 3524567Szliu .globl libm$cdabs_r6 3624567Szliu .globl libm$dsqrt_r5 3724567Szliu 3824567Szliu# entry for c functions cabs and hypot 3924567Szliu_cabs: 4024567Szliu_hypot: 4124567Szliu .word 0x807c # save r2-r6, enable floating overflow 4224567Szliu movq 4(ap),r0 # r0:1 = x 4324567Szliu movq 12(ap),r2 # r2:3 = y 4424567Szliu jmp cabs2 4524567Szliu# entry for Fortran use, call by: d = abs(z) 4624567Szliu_z_abs: 4724567Szliu .word 0x807c # save r2-r6, enable floating overflow 4824567Szliu movl 4(ap),r2 # indirect addressing is necessary here 4924567Szliu movq (r2)+,r0 # r0:1 = x 5024567Szliu movq (r2),r2 # r2:3 = y 5124567Szliu 5224567Szliucabs2: 5324567Szliu bicw3 $0x7f,r0,r4 # r4 has signed biased exp of x 5424567Szliu cmpw $0x8000,r4 5524567Szliu jeql return # x is a reserved operand, so return it 5624567Szliu bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y 5724567Szliu cmpw $0x8000,r5 5824567Szliu jneq cont # y isn't a reserved operand 5924567Szliu movq r2,r0 # return y if it's reserved 6024567Szliu ret 6124567Szliu 6224567Szliucont: 6324567Szliu bsbb regs_set # r0:1 = dsqrt(x^2+y^2)/2^r6 6424567Szliu addw2 r6,r0 # unscaled cdabs in r0:1 6524567Szliu jvc return # unless it overflows 6624567Szliu subw2 $0x80,r0 # halve r0 to get meaningful overflow 6724567Szliu addd2 r0,r0 # overflow; r0 is half of true abs value 6824567Szliureturn: 6924567Szliu ret 7024567Szliu 7124567Szliulibm$cdabs_r6: # ENTRY POINT for cdsqrt 7224567Szliu # calculates a scaled (factor in r6) 7324567Szliu # complex absolute value 7424567Szliu 7524567Szliu movq (r4)+,r0 # r0:r1 = x via indirect addressing 7624567Szliu movq (r4),r2 # r2:r3 = y via indirect addressing 7724567Szliu 7824567Szliu bicw3 $0x7f,r0,r5 # r5 has signed biased exp of x 7924567Szliu cmpw $0x8000,r5 8024567Szliu jeql cdreserved # x is a reserved operand 8124567Szliu bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y 8224567Szliu cmpw $0x8000,r5 8324567Szliu jneq regs_set # y isn't a reserved operand either? 8424567Szliu 8524567Szliucdreserved: 8624567Szliu movl *4(ap),r4 # r4 -> (u,v), if x or y is reserved 8724567Szliu movq r0,(r4)+ # copy u and v as is and return 8824567Szliu movq r2,(r4) # (again addressing is indirect) 8924567Szliu ret 9024567Szliu 9124567Szliuregs_set: 9224567Szliu bicw2 $0x8000,r0 # r0:r1 = dabs(x) 9324567Szliu bicw2 $0x8000,r2 # r2:r3 = dabs(y) 9424567Szliu cmpw r0,r2 9524567Szliu jgeq ordered 9624567Szliu movq r0,r4 9724567Szliu movq r2,r0 9824567Szliu movq r4,r2 # force y's exp <= x's exp 9924567Szliuordered: 10024567Szliu bicw3 $0x7f,r0,r6 # r6 = exponent(x) + bias(129) 10124567Szliu jeql retsb # if x = y = 0 then cdabs(x,y) = 0 10224567Szliu subw2 $0x4780,r6 # r6 = exponent(x) - 14 10324567Szliu subw2 r6,r0 # 2^14 <= scaled x < 2^15 10424567Szliu bitw $0xff80,r2 10524567Szliu jeql retsb # if y = 0 return dabs(x) 10624567Szliu subw2 r6,r2 10724567Szliu cmpw $0x3780,r2 # if scaled y < 2^-18 10824567Szliu jgtr retsb # return dabs(x) 10924567Szliu emodd r0,$0,r0,r4,r0 # r4 + r0:1 = scaled x^2 11024567Szliu emodd r2,$0,r2,r5,r2 # r5 + r2:3 = scaled y^2 11124567Szliu addd2 r2,r0 11224567Szliu addl2 r5,r4 11324567Szliu cvtld r4,r2 11424567Szliu addd2 r2,r0 # r0:1 = scaled x^2 + y^2 11524567Szliu jmp libm$dsqrt_r5 # r0:1 = dsqrt(x^2+y^2)/2^r6 11624567Szliuretsb: 11724567Szliu rsb # error < 0.86 ulp 118