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