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