1*61318Sbostic# Copyright (c) 1985, 1993 2*61318Sbostic# The Regents of the University of California. All rights reserved. 324567Szliu# 445308Sbostic# %sccs.include.redist.sh% 534125Sbostic# 6*61318Sbostic# @(#)cabs.s 8.1 (Berkeley) 06/04/93 734125Sbostic# 824729Selefunt .data 924729Selefunt .align 2 1024729Selefunt_sccsid: 11*61318Sbostic.asciz "@(#)cabs.s 1.2 (Berkeley) 8/21/85; 8.1 (ucb.elefunt) 06/04/93" 1224567Szliu 1324567Szliu# double precision complex absolute value 1424567Szliu# CABS by W. Kahan, 9/7/80. 1524567Szliu# Revised for reserved operands by E. LeBlanc, 8/18/82 1624567Szliu# argument for complex absolute value by reference, *4(ap) 1724567Szliu# argument for cabs and hypot (C fcns) by value, 4(ap) 1824567Szliu# output is in r0:r1 (error less than 0.86 ulps) 1924567Szliu 2024567Szliu .text 2124567Szliu .align 1 2224567Szliu .globl _cabs 2324567Szliu .globl _hypot 2424567Szliu .globl _z_abs 2524567Szliu .globl libm$cdabs_r6 2624567Szliu .globl libm$dsqrt_r5 2724567Szliu 2824567Szliu# entry for c functions cabs and hypot 2924567Szliu_cabs: 3024567Szliu_hypot: 3124567Szliu .word 0x807c # save r2-r6, enable floating overflow 3224567Szliu movq 4(ap),r0 # r0:1 = x 3324567Szliu movq 12(ap),r2 # r2:3 = y 3424567Szliu jmp cabs2 3524567Szliu# entry for Fortran use, call by: d = abs(z) 3624567Szliu_z_abs: 3724567Szliu .word 0x807c # save r2-r6, enable floating overflow 3824567Szliu movl 4(ap),r2 # indirect addressing is necessary here 3924567Szliu movq (r2)+,r0 # r0:1 = x 4024567Szliu movq (r2),r2 # r2:3 = y 4124567Szliu 4224567Szliucabs2: 4324567Szliu bicw3 $0x7f,r0,r4 # r4 has signed biased exp of x 4424567Szliu cmpw $0x8000,r4 4524567Szliu jeql return # x is a reserved operand, so return it 4624567Szliu bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y 4724567Szliu cmpw $0x8000,r5 4824567Szliu jneq cont # y isn't a reserved operand 4924567Szliu movq r2,r0 # return y if it's reserved 5024567Szliu ret 5124567Szliu 5224567Szliucont: 5324567Szliu bsbb regs_set # r0:1 = dsqrt(x^2+y^2)/2^r6 5424567Szliu addw2 r6,r0 # unscaled cdabs in r0:1 5524567Szliu jvc return # unless it overflows 5624567Szliu subw2 $0x80,r0 # halve r0 to get meaningful overflow 5724567Szliu addd2 r0,r0 # overflow; r0 is half of true abs value 5824567Szliureturn: 5924567Szliu ret 6024567Szliu 6124567Szliulibm$cdabs_r6: # ENTRY POINT for cdsqrt 6224567Szliu # calculates a scaled (factor in r6) 6324567Szliu # complex absolute value 6424567Szliu 6524567Szliu movq (r4)+,r0 # r0:r1 = x via indirect addressing 6624567Szliu movq (r4),r2 # r2:r3 = y via indirect addressing 6724567Szliu 6824567Szliu bicw3 $0x7f,r0,r5 # r5 has signed biased exp of x 6924567Szliu cmpw $0x8000,r5 7024567Szliu jeql cdreserved # x is a reserved operand 7124567Szliu bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y 7224567Szliu cmpw $0x8000,r5 7324567Szliu jneq regs_set # y isn't a reserved operand either? 7424567Szliu 7524567Szliucdreserved: 7624567Szliu movl *4(ap),r4 # r4 -> (u,v), if x or y is reserved 7724567Szliu movq r0,(r4)+ # copy u and v as is and return 7824567Szliu movq r2,(r4) # (again addressing is indirect) 7924567Szliu ret 8024567Szliu 8124567Szliuregs_set: 8224567Szliu bicw2 $0x8000,r0 # r0:r1 = dabs(x) 8324567Szliu bicw2 $0x8000,r2 # r2:r3 = dabs(y) 8424567Szliu cmpw r0,r2 8524567Szliu jgeq ordered 8624567Szliu movq r0,r4 8724567Szliu movq r2,r0 8824567Szliu movq r4,r2 # force y's exp <= x's exp 8924567Szliuordered: 9024567Szliu bicw3 $0x7f,r0,r6 # r6 = exponent(x) + bias(129) 9124567Szliu jeql retsb # if x = y = 0 then cdabs(x,y) = 0 9224567Szliu subw2 $0x4780,r6 # r6 = exponent(x) - 14 9324567Szliu subw2 r6,r0 # 2^14 <= scaled x < 2^15 9424567Szliu bitw $0xff80,r2 9524567Szliu jeql retsb # if y = 0 return dabs(x) 9624567Szliu subw2 r6,r2 9724567Szliu cmpw $0x3780,r2 # if scaled y < 2^-18 9824567Szliu jgtr retsb # return dabs(x) 9924567Szliu emodd r0,$0,r0,r4,r0 # r4 + r0:1 = scaled x^2 10024567Szliu emodd r2,$0,r2,r5,r2 # r5 + r2:3 = scaled y^2 10124567Szliu addd2 r2,r0 10224567Szliu addl2 r5,r4 10324567Szliu cvtld r4,r2 10424567Szliu addd2 r2,r0 # r0:1 = scaled x^2 + y^2 10524567Szliu jmp libm$dsqrt_r5 # r0:1 = dsqrt(x^2+y^2)/2^r6 10624567Szliuretsb: 10724567Szliu rsb # error < 0.86 ulp 108