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