1*24567Szliu# 2*24567Szliu# Copyright (c) 1985 Regents of the University of California. 3*24567Szliu# 4*24567Szliu# Use and reproduction of this software are granted in accordance with 5*24567Szliu# the terms and conditions specified in the Berkeley Software License 6*24567Szliu# Agreement (in particular, this entails acknowledgement of the programs' 7*24567Szliu# source, and inclusion of this notice) with the additional understanding 8*24567Szliu# that all recipients should regard themselves as participants in an 9*24567Szliu# ongoing research project and hence should feel obligated to report 10*24567Szliu# their experiences (good or bad) with these elementary function codes, 11*24567Szliu# using "sendbug 4bsd-bugs@BERKELEY", to the authors. 12*24567Szliu# 13*24567Szliu 14*24567Szliu# @(#)cabs.s 1.1 (ELEFUNT) 09/06/85 15*24567Szliu 16*24567Szliu# double precision complex absolute value 17*24567Szliu# CABS by W. Kahan, 9/7/80. 18*24567Szliu# Revised for reserved operands by E. LeBlanc, 8/18/82 19*24567Szliu# argument for complex absolute value by reference, *4(ap) 20*24567Szliu# argument for cabs and hypot (C fcns) by value, 4(ap) 21*24567Szliu# output is in r0:r1 (error less than 0.86 ulps) 22*24567Szliu 23*24567Szliu .text 24*24567Szliu .align 1 25*24567Szliu .globl _cabs 26*24567Szliu .globl _hypot 27*24567Szliu .globl _z_abs 28*24567Szliu .globl libm$cdabs_r6 29*24567Szliu .globl libm$dsqrt_r5 30*24567Szliu 31*24567Szliu# entry for c functions cabs and hypot 32*24567Szliu_cabs: 33*24567Szliu_hypot: 34*24567Szliu .word 0x807c # save r2-r6, enable floating overflow 35*24567Szliu movq 4(ap),r0 # r0:1 = x 36*24567Szliu movq 12(ap),r2 # r2:3 = y 37*24567Szliu jmp cabs2 38*24567Szliu# entry for Fortran use, call by: d = abs(z) 39*24567Szliu_z_abs: 40*24567Szliu .word 0x807c # save r2-r6, enable floating overflow 41*24567Szliu movl 4(ap),r2 # indirect addressing is necessary here 42*24567Szliu movq (r2)+,r0 # r0:1 = x 43*24567Szliu movq (r2),r2 # r2:3 = y 44*24567Szliu 45*24567Szliucabs2: 46*24567Szliu bicw3 $0x7f,r0,r4 # r4 has signed biased exp of x 47*24567Szliu cmpw $0x8000,r4 48*24567Szliu jeql return # x is a reserved operand, so return it 49*24567Szliu bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y 50*24567Szliu cmpw $0x8000,r5 51*24567Szliu jneq cont # y isn't a reserved operand 52*24567Szliu movq r2,r0 # return y if it's reserved 53*24567Szliu ret 54*24567Szliu 55*24567Szliucont: 56*24567Szliu bsbb regs_set # r0:1 = dsqrt(x^2+y^2)/2^r6 57*24567Szliu addw2 r6,r0 # unscaled cdabs in r0:1 58*24567Szliu jvc return # unless it overflows 59*24567Szliu subw2 $0x80,r0 # halve r0 to get meaningful overflow 60*24567Szliu addd2 r0,r0 # overflow; r0 is half of true abs value 61*24567Szliureturn: 62*24567Szliu ret 63*24567Szliu 64*24567Szliulibm$cdabs_r6: # ENTRY POINT for cdsqrt 65*24567Szliu # calculates a scaled (factor in r6) 66*24567Szliu # complex absolute value 67*24567Szliu 68*24567Szliu movq (r4)+,r0 # r0:r1 = x via indirect addressing 69*24567Szliu movq (r4),r2 # r2:r3 = y via indirect addressing 70*24567Szliu 71*24567Szliu bicw3 $0x7f,r0,r5 # r5 has signed biased exp of x 72*24567Szliu cmpw $0x8000,r5 73*24567Szliu jeql cdreserved # x is a reserved operand 74*24567Szliu bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y 75*24567Szliu cmpw $0x8000,r5 76*24567Szliu jneq regs_set # y isn't a reserved operand either? 77*24567Szliu 78*24567Szliucdreserved: 79*24567Szliu movl *4(ap),r4 # r4 -> (u,v), if x or y is reserved 80*24567Szliu movq r0,(r4)+ # copy u and v as is and return 81*24567Szliu movq r2,(r4) # (again addressing is indirect) 82*24567Szliu ret 83*24567Szliu 84*24567Szliuregs_set: 85*24567Szliu bicw2 $0x8000,r0 # r0:r1 = dabs(x) 86*24567Szliu bicw2 $0x8000,r2 # r2:r3 = dabs(y) 87*24567Szliu cmpw r0,r2 88*24567Szliu jgeq ordered 89*24567Szliu movq r0,r4 90*24567Szliu movq r2,r0 91*24567Szliu movq r4,r2 # force y's exp <= x's exp 92*24567Szliuordered: 93*24567Szliu bicw3 $0x7f,r0,r6 # r6 = exponent(x) + bias(129) 94*24567Szliu jeql retsb # if x = y = 0 then cdabs(x,y) = 0 95*24567Szliu subw2 $0x4780,r6 # r6 = exponent(x) - 14 96*24567Szliu subw2 r6,r0 # 2^14 <= scaled x < 2^15 97*24567Szliu bitw $0xff80,r2 98*24567Szliu jeql retsb # if y = 0 return dabs(x) 99*24567Szliu subw2 r6,r2 100*24567Szliu cmpw $0x3780,r2 # if scaled y < 2^-18 101*24567Szliu jgtr retsb # return dabs(x) 102*24567Szliu emodd r0,$0,r0,r4,r0 # r4 + r0:1 = scaled x^2 103*24567Szliu emodd r2,$0,r2,r5,r2 # r5 + r2:3 = scaled y^2 104*24567Szliu addd2 r2,r0 105*24567Szliu addl2 r5,r4 106*24567Szliu cvtld r4,r2 107*24567Szliu addd2 r2,r0 # r0:1 = scaled x^2 + y^2 108*24567Szliu jmp libm$dsqrt_r5 # r0:1 = dsqrt(x^2+y^2)/2^r6 109*24567Szliuretsb: 110*24567Szliu rsb # error < 0.86 ulp 111