| 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) %G% |
| 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 |