1# Copyright (c) 1985 Regents of the University of California. 2# All rights reserved. 3# 4# Redistribution and use in source and binary forms are permitted 5# provided that this notice is preserved and that due credit is given 6# to the University of California at Berkeley. The name of the University 7# may not be used to endorse or promote products derived from this 8# software without specific prior written permission. This software 9# is provided ``as is'' without express or implied warranty. 10# 11# All recipients should regard themselves as participants in an ongoing 12# research project and hence should feel obligated to report their 13# experiences (good or bad) with these elementary function codes, using 14# the sendbug(8) program, to the authors. 15# 16# @(#)cabs.s 5.2 (Berkeley) 04/29/88 17# 18 .data 19 .align 2 20_sccsid: 21.asciz "@(#)cabs.s 1.2 (Berkeley) 8/21/85; 5.2 (ucb.elefunt) 04/29/88" 22 23# double precision complex absolute value 24# CABS by W. Kahan, 9/7/80. 25# Revised for reserved operands by E. LeBlanc, 8/18/82 26# argument for complex absolute value by reference, *4(ap) 27# argument for cabs and hypot (C fcns) by value, 4(ap) 28# output is in r0:r1 (error less than 0.86 ulps) 29 30 .text 31 .align 1 32 .globl _cabs 33 .globl _hypot 34 .globl _z_abs 35 .globl libm$cdabs_r6 36 .globl libm$dsqrt_r5 37 38# entry for c functions cabs and hypot 39_cabs: 40_hypot: 41 .word 0x807c # save r2-r6, enable floating overflow 42 movq 4(ap),r0 # r0:1 = x 43 movq 12(ap),r2 # r2:3 = y 44 jmp cabs2 45# entry for Fortran use, call by: d = abs(z) 46_z_abs: 47 .word 0x807c # save r2-r6, enable floating overflow 48 movl 4(ap),r2 # indirect addressing is necessary here 49 movq (r2)+,r0 # r0:1 = x 50 movq (r2),r2 # r2:3 = y 51 52cabs2: 53 bicw3 $0x7f,r0,r4 # r4 has signed biased exp of x 54 cmpw $0x8000,r4 55 jeql return # x is a reserved operand, so return it 56 bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y 57 cmpw $0x8000,r5 58 jneq cont # y isn't a reserved operand 59 movq r2,r0 # return y if it's reserved 60 ret 61 62cont: 63 bsbb regs_set # r0:1 = dsqrt(x^2+y^2)/2^r6 64 addw2 r6,r0 # unscaled cdabs in r0:1 65 jvc return # unless it overflows 66 subw2 $0x80,r0 # halve r0 to get meaningful overflow 67 addd2 r0,r0 # overflow; r0 is half of true abs value 68return: 69 ret 70 71libm$cdabs_r6: # ENTRY POINT for cdsqrt 72 # calculates a scaled (factor in r6) 73 # complex absolute value 74 75 movq (r4)+,r0 # r0:r1 = x via indirect addressing 76 movq (r4),r2 # r2:r3 = y via indirect addressing 77 78 bicw3 $0x7f,r0,r5 # r5 has signed biased exp of x 79 cmpw $0x8000,r5 80 jeql cdreserved # x is a reserved operand 81 bicw3 $0x7f,r2,r5 # r5 has signed biased exp of y 82 cmpw $0x8000,r5 83 jneq regs_set # y isn't a reserved operand either? 84 85cdreserved: 86 movl *4(ap),r4 # r4 -> (u,v), if x or y is reserved 87 movq r0,(r4)+ # copy u and v as is and return 88 movq r2,(r4) # (again addressing is indirect) 89 ret 90 91regs_set: 92 bicw2 $0x8000,r0 # r0:r1 = dabs(x) 93 bicw2 $0x8000,r2 # r2:r3 = dabs(y) 94 cmpw r0,r2 95 jgeq ordered 96 movq r0,r4 97 movq r2,r0 98 movq r4,r2 # force y's exp <= x's exp 99ordered: 100 bicw3 $0x7f,r0,r6 # r6 = exponent(x) + bias(129) 101 jeql retsb # if x = y = 0 then cdabs(x,y) = 0 102 subw2 $0x4780,r6 # r6 = exponent(x) - 14 103 subw2 r6,r0 # 2^14 <= scaled x < 2^15 104 bitw $0xff80,r2 105 jeql retsb # if y = 0 return dabs(x) 106 subw2 r6,r2 107 cmpw $0x3780,r2 # if scaled y < 2^-18 108 jgtr retsb # return dabs(x) 109 emodd r0,$0,r0,r4,r0 # r4 + r0:1 = scaled x^2 110 emodd r2,$0,r2,r5,r2 # r5 + r2:3 = scaled y^2 111 addd2 r2,r0 112 addl2 r5,r4 113 cvtld r4,r2 114 addd2 r2,r0 # r0:1 = scaled x^2 + y^2 115 jmp libm$dsqrt_r5 # r0:1 = dsqrt(x^2+y^2)/2^r6 116retsb: 117 rsb # error < 0.86 ulp 118