![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to atan2.s CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 43BSDTahoe / usr.lib / libm / vax|
Return to atan2.s CVS log | Up to [CSRG BSD Unix] / 43BSDTahoe / usr.lib / libm / vax |
1.1 ! root 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: # @(#)atan2.s 5.3 (Berkeley) 6/30/88 ! 22: # ! 23: .data ! 24: .align 2 ! 25: _sccsid: ! 26: .asciz "@(#)atan2.s 1.2 (Berkeley) 8/21/85; 5.3 (ucb.elefunt) 6/30/88" ! 27: ! 28: # ATAN2(Y,X) ! 29: # RETURN ARG (X+iY) ! 30: # VAX D FORMAT (56 BITS PRECISION) ! 31: # CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/16/85; ! 32: # ! 33: # ! 34: # Method : ! 35: # 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x). ! 36: # 2. Reduce x to positive by (if x and y are unexceptional): ! 37: # ARG (x+iy) = arctan(y/x) ... if x > 0, ! 38: # ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0, ! 39: # 3. According to the integer k=4t+0.25 truncated , t=y/x, the argument ! 40: # is further reduced to one of the following intervals and the ! 41: # arctangent of y/x is evaluated by the corresponding formula: ! 42: # ! 43: # [0,7/16] atan(y/x) = t - t^3*(a1+t^2*(a2+...(a10+t^2*a11)...) ! 44: # [7/16,11/16] atan(y/x) = atan(1/2) + atan( (y-x/2)/(x+y/2) ) ! 45: # [11/16.19/16] atan(y/x) = atan( 1 ) + atan( (y-x)/(x+y) ) ! 46: # [19/16,39/16] atan(y/x) = atan(3/2) + atan( (y-1.5x)/(x+1.5y) ) ! 47: # [39/16,INF] atan(y/x) = atan(INF) + atan( -x/y ) ! 48: # ! 49: # Special cases: ! 50: # Notations: atan2(y,x) == ARG (x+iy) == ARG(x,y). ! 51: # ! 52: # ARG( NAN , (anything) ) is NaN; ! 53: # ARG( (anything), NaN ) is NaN; ! 54: # ARG(+(anything but NaN), +-0) is +-0 ; ! 55: # ARG(-(anything but NaN), +-0) is +-PI ; ! 56: # ARG( 0, +-(anything but 0 and NaN) ) is +-PI/2; ! 57: # ARG( +INF,+-(anything but INF and NaN) ) is +-0 ; ! 58: # ARG( -INF,+-(anything but INF and NaN) ) is +-PI; ! 59: # ARG( +INF,+-INF ) is +-PI/4 ; ! 60: # ARG( -INF,+-INF ) is +-3PI/4; ! 61: # ARG( (anything but,0,NaN, and INF),+-INF ) is +-PI/2; ! 62: # ! 63: # Accuracy: ! 64: # atan2(y,x) returns the exact ARG(x+iy) nearly rounded. ! 65: # ! 66: .text ! 67: .align 1 ! 68: .globl _atan2 ! 69: _atan2 : ! 70: .word 0x0ff4 ! 71: movq 4(ap),r2 # r2 = y ! 72: movq 12(ap),r4 # r4 = x ! 73: bicw3 $0x7f,r2,r0 ! 74: bicw3 $0x7f,r4,r1 ! 75: cmpw r0,$0x8000 # y is the reserved operand ! 76: jeql resop ! 77: cmpw r1,$0x8000 # x is the reserved operand ! 78: jeql resop ! 79: subl2 $8,sp ! 80: bicw3 $0x7fff,r2,-4(fp) # copy y sign bit to -4(fp) ! 81: bicw3 $0x7fff,r4,-8(fp) # copy x sign bit to -8(fp) ! 82: cmpd r4,$0x4080 # x = 1.0 ? ! 83: bneq xnot1 ! 84: movq r2,r0 ! 85: bicw2 $0x8000,r0 # t = |y| ! 86: movq r0,r2 # y = |y| ! 87: brb begin ! 88: xnot1: ! 89: bicw3 $0x807f,r2,r11 # yexp ! 90: jeql yeq0 # if y=0 goto yeq0 ! 91: bicw3 $0x807f,r4,r10 # xexp ! 92: jeql pio2 # if x=0 goto pio2 ! 93: subw2 r10,r11 # k = yexp - xexp ! 94: cmpw r11,$0x2000 # k >= 64 (exp) ? ! 95: jgeq pio2 # atan2 = +-pi/2 ! 96: divd3 r4,r2,r0 # t = y/x never overflow ! 97: bicw2 $0x8000,r0 # t > 0 ! 98: bicw2 $0xff80,r2 # clear the exponent of y ! 99: bicw2 $0xff80,r4 # clear the exponent of x ! 100: bisw2 $0x4080,r2 # normalize y to [1,2) ! 101: bisw2 $0x4080,r4 # normalize x to [1,2) ! 102: subw2 r11,r4 # scale x so that yexp-xexp=k ! 103: begin: ! 104: cmpw r0,$0x411c # t : 39/16 ! 105: jgeq L50 ! 106: addl3 $0x180,r0,r10 # 8*t ! 107: cvtrfl r10,r10 # [8*t] rounded to int ! 108: ashl $-1,r10,r10 # [8*t]/2 ! 109: casel r10,$0,$4 ! 110: L1: ! 111: .word L20-L1 ! 112: .word L20-L1 ! 113: .word L30-L1 ! 114: .word L40-L1 ! 115: .word L40-L1 ! 116: L10: ! 117: movq $0xb4d9940f985e407b,r6 # Hi=.98279372324732906796d0 ! 118: movq $0x21b1879a3bc2a2fc,r8 # Lo=-.17092002525602665777d-17 ! 119: subd3 r4,r2,r0 # y-x ! 120: addw2 $0x80,r0 # 2(y-x) ! 121: subd2 r4,r0 # 2(y-x)-x ! 122: addw2 $0x80,r4 # 2x ! 123: movq r2,r10 ! 124: addw2 $0x80,r10 # 2y ! 125: addd2 r10,r2 # 3y ! 126: addd2 r4,r2 # 3y+2x ! 127: divd2 r2,r0 # (2y-3x)/(2x+3y) ! 128: brw L60 ! 129: L20: ! 130: cmpw r0,$0x3280 # t : 2**(-28) ! 131: jlss L80 ! 132: clrq r6 # Hi=r6=0, Lo=r8=0 ! 133: clrq r8 ! 134: brw L60 ! 135: L30: ! 136: movq $0xda7b2b0d63383fed,r6 # Hi=.46364760900080611433d0 ! 137: movq $0xf0ea17b2bf912295,r8 # Lo=.10147340032515978826d-17 ! 138: movq r2,r0 ! 139: addw2 $0x80,r0 # 2y ! 140: subd2 r4,r0 # 2y-x ! 141: addw2 $0x80,r4 # 2x ! 142: addd2 r2,r4 # 2x+y ! 143: divd2 r4,r0 # (2y-x)/(2x+y) ! 144: brb L60 ! 145: L50: ! 146: movq $0x68c2a2210fda40c9,r6 # Hi=1.5707963267948966135d1 ! 147: movq $0x06e0145c26332326,r8 # Lo=.22517417741562176079d-17 ! 148: cmpw r0,$0x5100 # y : 2**57 ! 149: bgeq L90 ! 150: divd3 r2,r4,r0 ! 151: bisw2 $0x8000,r0 # -x/y ! 152: brb L60 ! 153: L40: ! 154: movq $0x68c2a2210fda4049,r6 # Hi=.78539816339744830676d0 ! 155: movq $0x06e0145c263322a6,r8 # Lo=.11258708870781088040d-17 ! 156: subd3 r4,r2,r0 # y-x ! 157: addd2 r4,r2 # y+x ! 158: divd2 r2,r0 # (y-x)/(y+x) ! 159: L60: ! 160: movq r0,r10 ! 161: muld2 r0,r0 ! 162: polyd r0,$12,ptable ! 163: muld2 r10,r0 ! 164: subd2 r0,r8 ! 165: addd3 r8,r10,r0 ! 166: addd2 r6,r0 ! 167: L80: ! 168: movw -8(fp),r2 ! 169: bneq pim ! 170: bisw2 -4(fp),r0 # return sign(y)*r0 ! 171: ret ! 172: L90: # x >= 2**25 ! 173: movq r6,r0 ! 174: brb L80 ! 175: pim: ! 176: subd3 r0,$0x68c2a2210fda4149,r0 # pi-t ! 177: bisw2 -4(fp),r0 ! 178: ret ! 179: yeq0: ! 180: movw -8(fp),r2 ! 181: beql zero # if sign(x)=1 return pi ! 182: movq $0x68c2a2210fda4149,r0 # pi=3.1415926535897932270d1 ! 183: ret ! 184: zero: ! 185: clrq r0 # return 0 ! 186: ret ! 187: pio2: ! 188: movq $0x68c2a2210fda40c9,r0 # pi/2=1.5707963267948966135d1 ! 189: bisw2 -4(fp),r0 # return sign(y)*pi/2 ! 190: ret ! 191: resop: ! 192: movq $0x8000,r0 # propagate the reserved operand ! 193: ret ! 194: .align 2 ! 195: ptable: ! 196: .quad 0xb50f5ce96e7abd60 ! 197: .quad 0x51e44a42c1073e02 ! 198: .quad 0x3487e3289643be35 ! 199: .quad 0xdb62066dffba3e54 ! 200: .quad 0xcf8e2d5199abbe70 ! 201: .quad 0x26f39cb884883e88 ! 202: .quad 0x135117d18998be9d ! 203: .quad 0x602ce9742e883eba ! 204: .quad 0xa35ad0be8e38bee3 ! 205: .quad 0xffac922249243f12 ! 206: .quad 0x7f14ccccccccbf4c ! 207: .quad 0xaa8faaaaaaaa3faa ! 208: .quad 0x0000000000000000
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.