1# Copyright (c) 1985, 1993 2# The Regents of the University of California. All rights reserved. 3# 4# %sccs.include.redist.sh% 5# 6# @(#)cabs.s 8.1 (Berkeley) 06/04/93 7# 8 .data 9 .align 2 10_sccsid: 11.asciz "@(#)cabs.s 1.2 (Berkeley) 8/21/85; 8.1 (ucb.elefunt) 06/04/93" 12 13# double precision complex absolute value 14# CABS by W. Kahan, 9/7/80. 15# Revised for reserved operands by E. LeBlanc, 8/18/82 16# argument for complex absolute value by reference, *4(ap) 17# argument for cabs and hypot (C fcns) by value, 4(ap) 18# output is in r0:r1 (error less than 0.86 ulps) 19 20 .text 21 .align 1 22 .globl _cabs 23 .globl _hypot 24 .globl _z_abs 25 .globl libm$cdabs_r6 26 .globl libm$dsqrt_r5 27 28# entry for c functions cabs and hypot 29_cabs: 30_hypot: 31 .word 0x807c # save r2-r6, enable floating overflow 32 movq 4(ap),r0 # r0:1 = x 33 movq 12(ap),r2 # r2:3 = y 34 jmp cabs2 35# entry for Fortran use, call by: d = abs(z) 36_z_abs: 37 .word 0x807c # save r2-r6, enable floating overflow 38 movl 4(ap),r2 # indirect addressing is necessary here 39 movq (r2)+,r0 # r0:1 = x 40 movq (r2),r2 # r2:3 = y 41 42cabs2: 43 bicw3 $0x7f,r0,r4 # r4 has signed biased exp of x 44 cmpw $0x8000,r4 45 jeql return # x is a reserved operand, so return it 46 bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y 47 cmpw $0x8000,r5 48 jneq cont # y isn't a reserved operand 49 movq r2,r0 # return y if it's reserved 50 ret 51 52cont: 53 bsbb regs_set # r0:1 = dsqrt(x^2+y^2)/2^r6 54 addw2 r6,r0 # unscaled cdabs in r0:1 55 jvc return # unless it overflows 56 subw2 $0x80,r0 # halve r0 to get meaningful overflow 57 addd2 r0,r0 # overflow; r0 is half of true abs value 58return: 59 ret 60 61libm$cdabs_r6: # ENTRY POINT for cdsqrt 62 # calculates a scaled (factor in r6) 63 # complex absolute value 64 65 movq (r4)+,r0 # r0:r1 = x via indirect addressing 66 movq (r4),r2 # r2:r3 = y via indirect addressing 67 68 bicw3 $0x7f,r0,r5 # r5 has signed biased exp of x 69 cmpw $0x8000,r5 70 jeql cdreserved # x is a reserved operand 71 bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y 72 cmpw $0x8000,r5 73 jneq regs_set # y isn't a reserved operand either? 74 75cdreserved: 76 movl *4(ap),r4 # r4 -> (u,v), if x or y is reserved 77 movq r0,(r4)+ # copy u and v as is and return 78 movq r2,(r4) # (again addressing is indirect) 79 ret 80 81regs_set: 82 bicw2 $0x8000,r0 # r0:r1 = dabs(x) 83 bicw2 $0x8000,r2 # r2:r3 = dabs(y) 84 cmpw r0,r2 85 jgeq ordered 86 movq r0,r4 87 movq r2,r0 88 movq r4,r2 # force y's exp <= x's exp 89ordered: 90 bicw3 $0x7f,r0,r6 # r6 = exponent(x) + bias(129) 91 jeql retsb # if x = y = 0 then cdabs(x,y) = 0 92 subw2 $0x4780,r6 # r6 = exponent(x) - 14 93 subw2 r6,r0 # 2^14 <= scaled x < 2^15 94 bitw $0xff80,r2 95 jeql retsb # if y = 0 return dabs(x) 96 subw2 r6,r2 97 cmpw $0x3780,r2 # if scaled y < 2^-18 98 jgtr retsb # return dabs(x) 99 emodd r0,$0,r0,r4,r0 # r4 + r0:1 = scaled x^2 100 emodd r2,$0,r2,r5,r2 # r5 + r2:3 = scaled y^2 101 addd2 r2,r0 102 addl2 r5,r4 103 cvtld r4,r2 104 addd2 r2,r0 # r0:1 = scaled x^2 + y^2 105 jmp libm$dsqrt_r5 # r0:1 = dsqrt(x^2+y^2)/2^r6 106retsb: 107 rsb # error < 0.86 ulp 108