124567Szliu# 224567Szliu# Copyright (c) 1985 Regents of the University of California. 324567Szliu# 424567Szliu# Use and reproduction of this software are granted in accordance with 524567Szliu# the terms and conditions specified in the Berkeley Software License 624567Szliu# Agreement (in particular, this entails acknowledgement of the programs' 724567Szliu# source, and inclusion of this notice) with the additional understanding 824567Szliu# that all recipients should regard themselves as participants in an 924567Szliu# ongoing research project and hence should feel obligated to report 1024567Szliu# their experiences (good or bad) with these elementary function codes, 1124567Szliu# using "sendbug 4bsd-bugs@BERKELEY", to the authors. 1224567Szliu# 13*24729Selefunt .data 14*24729Selefunt .align 2 15*24729Selefunt_sccsid: 16*24729Selefunt.asciz "@(#)cabs.s 1.2 (Berkeley) 8/21/85; 1.3 (ucb.elefunt) 09/12/85" 1724567Szliu 1824567Szliu# double precision complex absolute value 1924567Szliu# CABS by W. Kahan, 9/7/80. 2024567Szliu# Revised for reserved operands by E. LeBlanc, 8/18/82 2124567Szliu# argument for complex absolute value by reference, *4(ap) 2224567Szliu# argument for cabs and hypot (C fcns) by value, 4(ap) 2324567Szliu# output is in r0:r1 (error less than 0.86 ulps) 2424567Szliu 2524567Szliu .text 2624567Szliu .align 1 2724567Szliu .globl _cabs 2824567Szliu .globl _hypot 2924567Szliu .globl _z_abs 3024567Szliu .globl libm$cdabs_r6 3124567Szliu .globl libm$dsqrt_r5 3224567Szliu 3324567Szliu# entry for c functions cabs and hypot 3424567Szliu_cabs: 3524567Szliu_hypot: 3624567Szliu .word 0x807c # save r2-r6, enable floating overflow 3724567Szliu movq 4(ap),r0 # r0:1 = x 3824567Szliu movq 12(ap),r2 # r2:3 = y 3924567Szliu jmp cabs2 4024567Szliu# entry for Fortran use, call by: d = abs(z) 4124567Szliu_z_abs: 4224567Szliu .word 0x807c # save r2-r6, enable floating overflow 4324567Szliu movl 4(ap),r2 # indirect addressing is necessary here 4424567Szliu movq (r2)+,r0 # r0:1 = x 4524567Szliu movq (r2),r2 # r2:3 = y 4624567Szliu 4724567Szliucabs2: 4824567Szliu bicw3 $0x7f,r0,r4 # r4 has signed biased exp of x 4924567Szliu cmpw $0x8000,r4 5024567Szliu jeql return # x is a reserved operand, so return it 5124567Szliu bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y 5224567Szliu cmpw $0x8000,r5 5324567Szliu jneq cont # y isn't a reserved operand 5424567Szliu movq r2,r0 # return y if it's reserved 5524567Szliu ret 5624567Szliu 5724567Szliucont: 5824567Szliu bsbb regs_set # r0:1 = dsqrt(x^2+y^2)/2^r6 5924567Szliu addw2 r6,r0 # unscaled cdabs in r0:1 6024567Szliu jvc return # unless it overflows 6124567Szliu subw2 $0x80,r0 # halve r0 to get meaningful overflow 6224567Szliu addd2 r0,r0 # overflow; r0 is half of true abs value 6324567Szliureturn: 6424567Szliu ret 6524567Szliu 6624567Szliulibm$cdabs_r6: # ENTRY POINT for cdsqrt 6724567Szliu # calculates a scaled (factor in r6) 6824567Szliu # complex absolute value 6924567Szliu 7024567Szliu movq (r4)+,r0 # r0:r1 = x via indirect addressing 7124567Szliu movq (r4),r2 # r2:r3 = y via indirect addressing 7224567Szliu 7324567Szliu bicw3 $0x7f,r0,r5 # r5 has signed biased exp of x 7424567Szliu cmpw $0x8000,r5 7524567Szliu jeql cdreserved # x is a reserved operand 7624567Szliu bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y 7724567Szliu cmpw $0x8000,r5 7824567Szliu jneq regs_set # y isn't a reserved operand either? 7924567Szliu 8024567Szliucdreserved: 8124567Szliu movl *4(ap),r4 # r4 -> (u,v), if x or y is reserved 8224567Szliu movq r0,(r4)+ # copy u and v as is and return 8324567Szliu movq r2,(r4) # (again addressing is indirect) 8424567Szliu ret 8524567Szliu 8624567Szliuregs_set: 8724567Szliu bicw2 $0x8000,r0 # r0:r1 = dabs(x) 8824567Szliu bicw2 $0x8000,r2 # r2:r3 = dabs(y) 8924567Szliu cmpw r0,r2 9024567Szliu jgeq ordered 9124567Szliu movq r0,r4 9224567Szliu movq r2,r0 9324567Szliu movq r4,r2 # force y's exp <= x's exp 9424567Szliuordered: 9524567Szliu bicw3 $0x7f,r0,r6 # r6 = exponent(x) + bias(129) 9624567Szliu jeql retsb # if x = y = 0 then cdabs(x,y) = 0 9724567Szliu subw2 $0x4780,r6 # r6 = exponent(x) - 14 9824567Szliu subw2 r6,r0 # 2^14 <= scaled x < 2^15 9924567Szliu bitw $0xff80,r2 10024567Szliu jeql retsb # if y = 0 return dabs(x) 10124567Szliu subw2 r6,r2 10224567Szliu cmpw $0x3780,r2 # if scaled y < 2^-18 10324567Szliu jgtr retsb # return dabs(x) 10424567Szliu emodd r0,$0,r0,r4,r0 # r4 + r0:1 = scaled x^2 10524567Szliu emodd r2,$0,r2,r5,r2 # r5 + r2:3 = scaled y^2 10624567Szliu addd2 r2,r0 10724567Szliu addl2 r5,r4 10824567Szliu cvtld r4,r2 10924567Szliu addd2 r2,r0 # r0:1 = scaled x^2 + y^2 11024567Szliu jmp libm$dsqrt_r5 # r0:1 = dsqrt(x^2+y^2)/2^r6 11124567Szliuretsb: 11224567Szliu rsb # error < 0.86 ulp 113