![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to support.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 / national|
Return to support.s CVS log | Up to [CSRG BSD Unix] / 43BSDTahoe / usr.lib / libm / national |
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: ; @(#)support.s 5.3 (Berkeley) 6/30/88 ! 22: ; ! 23: ! 24: ; IEEE recommended functions ! 25: ; ! 26: ; double copysign(x,y) ! 27: ; double x,y; ! 28: ; IEEE 754 recommended function, return x*sign(y) ! 29: ; Coded by K.C. Ng in National 32k assembler, 11/9/85. ! 30: ; ! 31: .vers 2 ! 32: .text ! 33: .align 2 ! 34: .globl _copysign ! 35: _copysign: ! 36: movl 4(sp),f0 ! 37: movd 8(sp),r0 ! 38: movd 16(sp),r1 ! 39: xord r0,r1 ! 40: andd 0x80000000,r1 ! 41: cmpqd 0,r1 ! 42: beq end ! 43: negl f0,f0 ! 44: end: ret 0 ! 45: ! 46: ; ! 47: ; double logb(x) ! 48: ; double x; ! 49: ; IEEE p854 recommended function, return the exponent of x (return float(N) ! 50: ; such that 1 <= x*2**-N < 2, even for subnormal number. ! 51: ; Coded by K.C. Ng in National 32k assembler, 11/9/85. ! 52: ; Note: subnormal number (if implemented) will be taken care of. ! 53: ; ! 54: .vers 2 ! 55: .text ! 56: .align 2 ! 57: .globl _logb ! 58: _logb: ! 59: ; ! 60: ; extract the exponent of x ! 61: ; glossaries: r0 = high part of x ! 62: ; r1 = unbias exponent of x ! 63: ; r2 = 20 (first exponent bit position) ! 64: ; ! 65: movd 8(sp),r0 ! 66: movd 20,r2 ! 67: extd r2,r0,r1,11 ; extract the exponent of x ! 68: cmpqd 0,r1 ; if exponent bits = 0, goto L3 ! 69: beq L3 ! 70: cmpd 0x7ff,r1 ! 71: beq L2 ; if exponent bits = 0x7ff, goto L2 ! 72: L1: subd 1023,r1 ; unbias the exponent ! 73: movdl r1,f0 ; convert the exponent to floating value ! 74: ret 0 ! 75: ; ! 76: ; x is INF or NaN, simply return x ! 77: ; ! 78: L2: ! 79: movl 4(sp),f0 ; logb(+inf)=+inf, logb(NaN)=NaN ! 80: ret 0 ! 81: ; ! 82: ; x is 0 or subnormal ! 83: ; ! 84: L3: ! 85: movl 4(sp),f0 ! 86: cmpl 0f0,f0 ! 87: beq L5 ; x is 0 , goto L5 (return -inf) ! 88: ; ! 89: ; Now x is subnormal ! 90: ; ! 91: mull L64,f0 ; scale up f0 with 2**64 ! 92: movl f0,tos ! 93: movd tos,r0 ! 94: movd tos,r0 ; now r0 = new high part of x ! 95: extd r2,r0,r1,11 ; extract the exponent of x to r1 ! 96: subd 1087,r1 ; unbias the exponent with correction ! 97: movdl r1,f0 ; convert the exponent to floating value ! 98: ret 0 ! 99: ; ! 100: ; x is 0, return logb(0)= -INF ! 101: ; ! 102: L5: ! 103: movl 0f1.0e300,f0 ! 104: mull 0f-1.0e300,f0 ; multiply two big numbers to get -INF ! 105: ret 0 ! 106: ; ! 107: ; double rint(x) ! 108: ; double x; ! 109: ; ... delivers integer nearest x in direction of prevailing rounding ! 110: ; ... mode ! 111: ; Coded by K.C. Ng in National 32k assembler, 11/9/85. ! 112: ; Note: subnormal number (if implemented) will be taken care of. ! 113: ; ! 114: .vers 2 ! 115: .text ! 116: .align 2 ! 117: .globl _rint ! 118: _rint: ! 119: ; ! 120: movd 8(sp),r0 ! 121: movd 20,r2 ! 122: extd r2,r0,r1,11 ; extract the exponent of x ! 123: cmpd 0x433,r1 ! 124: ble itself ! 125: movl L52,f2 ; f2 = L = 2**52 ! 126: cmpqd 0,r0 ! 127: ble L1 ! 128: negl f2,f2 ; f2 = s = copysign(L,x) ! 129: L1: addl f2,f0 ; f0 = x + s ! 130: subl f2,f0 ; f0 = f0 - s ! 131: ret 0 ! 132: itself: movl 4(sp),f0 ! 133: ret 0 ! 134: L52: .double 0x0,0x43300000 ; L52=2**52 ! 135: ; ! 136: ; int finite(x) ! 137: ; double x; ! 138: ; IEEE 754 recommended function, return 0 if x is NaN or INF, else 0 ! 139: ; Coded by K.C. Ng in National 32k assembler, 11/9/85. ! 140: ; ! 141: .vers 2 ! 142: .text ! 143: .align 2 ! 144: .globl _finite ! 145: _finite: ! 146: movd 4(sp),r1 ! 147: andd 0x800fffff,r1 ! 148: cmpd 0x7ff00000,r1 ! 149: sned r0 ; r0=0 if exponent(x) = 0x7ff ! 150: ret 0 ! 151: ; ! 152: ; double scalb(x,N) ! 153: ; double x; int N; ! 154: ; IEEE 754 recommended function, return x*2**N by adjusting ! 155: ; exponent of x. ! 156: ; Coded by K.C. Ng in National 32k assembler, 11/9/85. ! 157: ; Note: subnormal number (if implemented) will be taken care of ! 158: ; ! 159: .vers 2 ! 160: .text ! 161: .align 2 ! 162: .globl _scalb ! 163: _scalb: ! 164: ; ! 165: ; if x=0 return 0 ! 166: ; ! 167: movl 4(sp),f0 ! 168: cmpl 0f0,f0 ! 169: beq end ; scalb(0,N) is x itself ! 170: ; ! 171: ; extract the exponent of x ! 172: ; glossaries: r0 = high part of x, ! 173: ; r1 = unbias exponent of x, ! 174: ; r2 = 20 (first exponent bit position). ! 175: ; ! 176: movd 8(sp),r0 ; r0 = high part of x ! 177: movd 20,r2 ; r2 = 20 ! 178: extd r2,r0,r1,11 ; extract the exponent of x in r1 ! 179: cmpd 0x7ff,r1 ! 180: ; ! 181: ; if exponent of x is 0x7ff, then x is NaN or INF; simply return x ! 182: ; ! 183: beq end ! 184: cmpqd 0,r1 ! 185: ; ! 186: ; if exponent of x is zero, then x is subnormal; goto L19 ! 187: ; ! 188: beq L19 ! 189: addd 12(sp),r1 ; r1 = (exponent of x) + N ! 190: bfs inof ; if integer overflows, goto inof ! 191: cmpqd 0,r1 ; if new exponent <= 0, goto underflow ! 192: bge underflow ! 193: cmpd 2047,r1 ; if new exponent >= 2047 goto overflow ! 194: ble overflow ! 195: insd r2,r1,r0,11 ; insert the new exponent ! 196: movd r0,tos ! 197: movd 8(sp),tos ! 198: movl tos,f0 ; return x*2**N ! 199: end: ret 0 ! 200: inof: bcs underflow ; negative int overflow if Carry bit is set ! 201: overflow: ! 202: andd 0x80000000,r0 ; keep the sign of x ! 203: ord 0x7fe00000,r0 ; set x to a huge number ! 204: movd r0,tos ! 205: movqd 0,tos ! 206: movl tos,f0 ! 207: mull 0f1.0e300,f0 ; multiply two huge number to get overflow ! 208: ret 0 ! 209: underflow: ! 210: addd 64,r1 ; add 64 to exonent to see if it is subnormal ! 211: cmpqd 0,r1 ! 212: bge zero ; underflow to zero ! 213: insd r2,r1,r0,11 ; insert the new exponent ! 214: movd r0,tos ! 215: movd 8(sp),tos ! 216: movl tos,f0 ! 217: mull L30,f0 ; readjust x by multiply it with 2**-64 ! 218: ret 0 ! 219: zero: andd 0x80000000,r0 ; keep the sign of x ! 220: ord 0x00100000,r0 ; set x to a tiny number ! 221: movd r0,tos ! 222: movqd 0,tos ! 223: movl tos,f0 ! 224: mull 0f1.0e-300,f0 ; underflow to 0 by multipling two tiny nos. ! 225: ret 0 ! 226: L19: ; subnormal number ! 227: mull L32,f0 ; scale up x by 2**64 ! 228: movl f0,tos ! 229: movd tos,r0 ! 230: movd tos,r0 ; get the high part of new x ! 231: extd r2,r0,r1,11 ; extract the exponent of x in r1 ! 232: addd 12(sp),r1 ; exponent of x + N ! 233: subd 64,r1 ; adjust it by subtracting 64 ! 234: cmpqd 0,r1 ! 235: bge underflow ! 236: cmpd 2047,r1 ! 237: ble overflow ! 238: insd r2,r1,r0,11 ; insert back the incremented exponent ! 239: movd r0,tos ! 240: movd 8(sp),tos ! 241: movl tos,f0 ! 242: end: ret 0 ! 243: L30: .double 0x0,0x3bf00000 ; floating point 2**-64 ! 244: L32: .double 0x0,0x43f00000 ; floating point 2**64 ! 245: ; ! 246: ; double drem(x,y) ! 247: ; double x,y; ! 248: ; IEEE double remainder function, return x-n*y, where n=x/y rounded to ! 249: ; nearest integer (half way case goes to even). Result exact. ! 250: ; Coded by K.C. Ng in National 32k assembly, 11/19/85. ! 251: ; ! 252: .vers 2 ! 253: .text ! 254: .align 2 ! 255: .globl _drem ! 256: _drem: ! 257: ; ! 258: ; glossaries: ! 259: ; r2 = high part of x ! 260: ; r3 = exponent of x ! 261: ; r4 = high part of y ! 262: ; r5 = exponent of y ! 263: ; r6 = sign of x ! 264: ; r7 = constant 0x7ff00000 ! 265: ; ! 266: ; 16(fp) : y ! 267: ; 8(fp) : x ! 268: ; -12(fp) : adjustment on y when y is subnormal ! 269: ; -16(fp) : fsr ! 270: ; -20(fp) : nx ! 271: ; -28(fp) : t ! 272: ; -36(fp) : t1 ! 273: ; -40(fp) : nf ! 274: ; ! 275: ; ! 276: enter [r3,r4,r5,r6,r7],40 ! 277: movl f6,tos ! 278: movl f4,tos ! 279: movl 0f0,-12(fp) ! 280: movd 0,-20(fp) ! 281: movd 0,-40(fp) ! 282: movd 0x7ff00000,r7 ; initialize r7=0x7ff00000 ! 283: movd 12(fp),r2 ; r2 = high(x) ! 284: movd r2,r3 ! 285: andd r7,r3 ; r3 = xexp ! 286: cmpd r7,r3 ! 287: ; if x is NaN or INF goto L1 ! 288: beq L1 ! 289: movd 20(fp),r4 ! 290: bicd [31],r4 ; r4 = high part of |y| ! 291: movd r4,20(fp) ; y = |y| ! 292: movd r4,r5 ! 293: andd r7,r5 ; r5 = yexp ! 294: cmpd r7,r5 ! 295: beq L2 ; if y is NaN or INF goto L2 ! 296: cmpd 0x04000000,r5 ; ! 297: bgt L3 ; if y is tiny goto L3 ! 298: ; ! 299: ; now y != 0 , x is finite ! 300: ; ! 301: L10: ! 302: movd r2,r6 ! 303: andd 0x80000000,r6 ; r6 = sign(x) ! 304: bicd [31],r2 ; x <- |x| ! 305: sfsr r1 ! 306: movd r1,-16(fp) ; save fsr in -16(fp) ! 307: bicd [5],r1 ! 308: lfsr r1 ; disable inexact interupt ! 309: movd 16(fp),r0 ; r0 = low part of y ! 310: movd r0,r1 ; r1 = r0 = low part of y ! 311: andd 0xf8000000,r1 ; mask off the lsb 27 bits of y ! 312: ! 313: movd r2,12(fp) ; update x to |x| ! 314: movd r0,-28(fp) ; ! 315: movd r4,-24(fp) ; t = y ! 316: movd r4,-32(fp) ; ! 317: movd r1,-36(fp) ; t1 = y with trialing 27 zeros ! 318: movd 0x01900000,r1 ; r1 = 25 in exponent field ! 319: LOOP: ! 320: movl 8(fp),f0 ; f0 = x ! 321: movl 16(fp),f2 ; f2 = y ! 322: cmpl f0,f2 ! 323: ble fnad ; goto fnad (final adjustment) if x <= y ! 324: movd r4,-32(fp) ! 325: movd r3,r0 ! 326: subd r5,r0 ; xexp - yexp ! 327: subd r1,r0 ; r0 = xexp - yexp - m25 ! 328: cmpqd 0,r0 ; r0 > 0 ? ! 329: bge 1f ! 330: addd r4,r0 ; scale up (high) y ! 331: movd r0,-24(fp) ; scale up t ! 332: movl -28(fp),f2 ; t ! 333: movd r0,-32(fp) ; scale up t1 ! 334: 1: ! 335: movl -36(fp),f4 ; t1 ! 336: movl f0,f6 ! 337: divl f2,f6 ; f6 = x/t ! 338: floorld f6,r0 ; r0 = [x/t] ! 339: movdl r0,f6 ; f6 = n ! 340: subl f4,f2 ; t = t - t1 (tail of t1) ! 341: mull f6,f4 ; f4 = n*t1 ...exact ! 342: subl f4,f0 ; x = x - n*t1 ! 343: mull f6,f2 ; n*(t-t1) ...exact ! 344: subl f2,f0 ; x = x - n*(t-t1) ! 345: ; update xexp ! 346: movl f0,8(fp) ! 347: movd 12(fp),r3 ! 348: andd r7,r3 ! 349: jump LOOP ! 350: fnad: ! 351: cmpqd 0,-20(fp) ; 0 = nx? ! 352: beq final ! 353: mull -12(fp),8(fp) ; scale up x the same amount as y ! 354: movd 0,-20(fp) ! 355: movd 12(fp),r2 ! 356: movd r2,r3 ! 357: andd r7,r3 ; update exponent of x ! 358: jump LOOP ! 359: ! 360: final: ! 361: movl 16(fp),f2 ; f2 = y (f0=x, r0=n) ! 362: subd 0x100000,r4 ; high y /2 ! 363: movd r4,-24(fp) ! 364: movl -28(fp),f4 ; f4 = y/2 ! 365: cmpl f0,f4 ; x > y/2 ? ! 366: bgt 1f ! 367: bne 2f ! 368: andd 1,r0 ; n is odd or even ! 369: cmpqd 0,r0 ! 370: beq 2f ! 371: 1: ! 372: subl f2,f0 ; x = x - y ! 373: 2: ! 374: cmpqd 0,-40(fp) ! 375: beq 3f ! 376: divl -12(fp),f0 ; scale down the answer ! 377: 3: ! 378: movl f0,tos ! 379: xord r6,tos ! 380: movl tos,f0 ! 381: movd -16(fp),r0 ! 382: lfsr r0 ; restore the fsr ! 383: ! 384: end: movl tos,f4 ! 385: movl tos,f6 ! 386: exit [r3,r4,r5,r6,r7] ! 387: ret 0 ! 388: ; ! 389: ; y is NaN or INF ! 390: ; ! 391: L2: ! 392: movd 16(fp),r0 ; r0 = low part of y ! 393: andd 0xfffff,r4 ; r4 = high part of y & 0x000fffff ! 394: ord r4,r0 ! 395: cmpqd 0,r0 ! 396: beq L4 ! 397: movl 16(fp),f0 ; y is NaN, return y ! 398: jump end ! 399: L4: movl 8(fp),f0 ; y is inf, return x ! 400: jump end ! 401: ; ! 402: ; exponent of y is less than 64, y may be zero or subnormal ! 403: ; ! 404: L3: ! 405: movl 16(fp),f0 ! 406: cmpl 0f0,f0 ! 407: bne L5 ! 408: divl f0,f0 ; y is 0, return NaN by doing 0/0 ! 409: jump end ! 410: ; ! 411: ; subnormal y or tiny y ! 412: ; ! 413: L5: ! 414: movd 0x04000000,-20(fp) ; nx = 64 in exponent field ! 415: movl L64,f2 ! 416: movl f2,-12(fp) ! 417: mull f2,f0 ! 418: cmpl f0,LTINY ! 419: bgt L6 ! 420: mull f2,f0 ! 421: addd 0x04000000,-20(fp) ; nx = nx + 64 in exponent field ! 422: mull f2,-12(fp) ! 423: L6: ! 424: movd -20(fp),-40(fp) ! 425: movl f0,16(fp) ! 426: movd 20(fp),r4 ! 427: movd r4,r5 ! 428: andd r7,r5 ; exponent of new y ! 429: jump L10 ! 430: ; ! 431: ; x is NaN or INF, return x-x ! 432: ; ! 433: L1: ! 434: movl 8(fp),f0 ! 435: subl f0,f0 ; if x is INF, then INF-INF is NaN ! 436: ret 0 ! 437: L64: .double 0x0,0x43f00000 ; L64 = 2**64 ! 438: LTINY: .double 0x0,0x04000000 ; LTINY = 2**-959
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.