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