; Copyright (c) 1985 Regents of the University of California.
; Redistribution and use in source and binary forms are permitted
; provided that the above copyright notice and this paragraph are
; duplicated in all such forms and that any documentation,
; advertising materials, and other materials related to such
; distribution and use acknowledge that the software was developed
; by the University of California, Berkeley. The name of the
; University may not be used to endorse or promote products derived
; from this software without specific prior written permission.
; THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
; IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
; WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
; All recipients should regard themselves as participants in an ongoing
; research project and hence should feel obligated to report their
; experiences (good or bad) with these elementary function codes, using
; the sendbug(8) program, to the authors.
; @(#)support.s 5.3 (Berkeley) %G%
; IEEE recommended functions
; IEEE 754 recommended function, return x*sign(y)
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
; IEEE p854 recommended function, return the exponent of x (return float(N)
; such that 1 <= x*2**-N < 2, even for subnormal number.
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
; Note: subnormal number (if implemented) will be taken care of.
; extract the exponent of x
; glossaries: r0 = high part of x
; r1 = unbias exponent of x
; r2 = 20 (first exponent bit position)
extd r2,r0,r1,11 ; extract the exponent of x
cmpqd 0,r1 ; if exponent bits = 0, goto L3
beq L2 ; if exponent bits = 0x7ff, goto L2
L1: subd 1023,r1 ; unbias the exponent
movdl r1,f0 ; convert the exponent to floating value
; x is INF or NaN, simply return x
movl 4(sp),f0 ; logb(+inf)=+inf, logb(NaN)=NaN
beq L5 ; x is 0 , goto L5 (return -inf)
mull L64,f0 ; scale up f0 with 2**64
movd tos,r0 ; now r0 = new high part of x
extd r2,r0,r1,11 ; extract the exponent of x to r1
subd 1087,r1 ; unbias the exponent with correction
movdl r1,f0 ; convert the exponent to floating value
; x is 0, return logb(0)= -INF
mull 0f-1.0e300,f0 ; multiply two big numbers to get -INF
; ... delivers integer nearest x in direction of prevailing rounding
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
; Note: subnormal number (if implemented) will be taken care of.
extd r2,r0,r1,11 ; extract the exponent of x
movl L52,f2 ; f2 = L = 2**52
negl f2,f2 ; f2 = s = copysign(L,x)
L1: addl f2,f0 ; f0 = x + s
L52: .double 0x0,0x43300000 ; L52=2**52
; IEEE 754 recommended function, return 0 if x is NaN or INF, else 0
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
sned r0 ; r0=0 if exponent(x) = 0x7ff
; IEEE 754 recommended function, return x*2**N by adjusting
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
; Note: subnormal number (if implemented) will be taken care of
beq end ; scalb(0,N) is x itself
; extract the exponent of x
; glossaries: r0 = high part of x,
; r1 = unbias exponent of x,
; r2 = 20 (first exponent bit position).
movd 8(sp),r0 ; r0 = high part of x
extd r2,r0,r1,11 ; extract the exponent of x in r1
; if exponent of x is 0x7ff, then x is NaN or INF; simply return x
; if exponent of x is zero, then x is subnormal; goto L19
addd 12(sp),r1 ; r1 = (exponent of x) + N
bfs inof ; if integer overflows, goto inof
cmpqd 0,r1 ; if new exponent <= 0, goto underflow
cmpd 2047,r1 ; if new exponent >= 2047 goto overflow
insd r2,r1,r0,11 ; insert the new exponent
movl tos,f0 ; return x*2**N
inof: bcs underflow ; negative int overflow if Carry bit is set
andd 0x80000000,r0 ; keep the sign of x
ord 0x7fe00000,r0 ; set x to a huge number
mull 0f1.0e300,f0 ; multiply two huge number to get overflow
addd 64,r1 ; add 64 to exonent to see if it is subnormal
bge zero ; underflow to zero
insd r2,r1,r0,11 ; insert the new exponent
mull L30,f0 ; readjust x by multiply it with 2**-64
zero: andd 0x80000000,r0 ; keep the sign of x
ord 0x00100000,r0 ; set x to a tiny number
mull 0f1.0e-300,f0 ; underflow to 0 by multipling two tiny nos.
mull L32,f0 ; scale up x by 2**64
movd tos,r0 ; get the high part of new x
extd r2,r0,r1,11 ; extract the exponent of x in r1
addd 12(sp),r1 ; exponent of x + N
subd 64,r1 ; adjust it by subtracting 64
insd r2,r1,r0,11 ; insert back the incremented exponent
L30: .double 0x0,0x3bf00000 ; floating point 2**-64
L32: .double 0x0,0x43f00000 ; floating point 2**64
; IEEE double remainder function, return x-n*y, where n=x/y rounded to
; nearest integer (half way case goes to even). Result exact.
; Coded by K.C. Ng in National 32k assembly, 11/19/85.
; r7 = constant 0x7ff00000
; -12(fp) : adjustment on y when y is subnormal
enter [r3,r4,r5,r6,r7],40
movd 0x7ff00000,r7 ; initialize r7=0x7ff00000
movd 12(fp),r2 ; r2 = high(x)
; if x is NaN or INF goto L1
bicd [31],r4 ; r4 = high part of |y|
beq L2 ; if y is NaN or INF goto L2
bgt L3 ; if y is tiny goto L3
; now y != 0 , x is finite
andd 0x80000000,r6 ; r6 = sign(x)
movd r1,-16(fp) ; save fsr in -16(fp)
lfsr r1 ; disable inexact interupt
movd 16(fp),r0 ; r0 = low part of y
movd r0,r1 ; r1 = r0 = low part of y
andd 0xf8000000,r1 ; mask off the lsb 27 bits of y
movd r2,12(fp) ; update x to |x|
movd r1,-36(fp) ; t1 = y with trialing 27 zeros
movd 0x01900000,r1 ; r1 = 25 in exponent field
ble fnad ; goto fnad (final adjustment) if x <= y
subd r1,r0 ; r0 = xexp - yexp - m25
addd r4,r0 ; scale up (high) y
movd r0,-24(fp) ; scale up t
movd r0,-32(fp) ; scale up t1
floorld f6,r0 ; r0 = [x/t]
subl f4,f2 ; t = t - t1 (tail of t1)
mull f6,f4 ; f4 = n*t1 ...exact
subl f4,f0 ; x = x - n*t1
mull f6,f2 ; n*(t-t1) ...exact
subl f2,f0 ; x = x - n*(t-t1)
cmpqd 0,-20(fp) ; 0 = nx?
mull -12(fp),8(fp) ; scale up x the same amount as y
andd r7,r3 ; update exponent of x
movl 16(fp),f2 ; f2 = y (f0=x, r0=n)
subd 0x100000,r4 ; high y /2
movl -28(fp),f4 ; f4 = y/2
andd 1,r0 ; n is odd or even
divl -12(fp),f0 ; scale down the answer
lfsr r0 ; restore the fsr
movd 16(fp),r0 ; r0 = low part of y
andd 0xfffff,r4 ; r4 = high part of y & 0x000fffff
movl 16(fp),f0 ; y is NaN, return y
L4: movl 8(fp),f0 ; y is inf, return x
; exponent of y is less than 64, y may be zero or subnormal
divl f0,f0 ; y is 0, return NaN by doing 0/0
movd 0x04000000,-20(fp) ; nx = 64 in exponent field
addd 0x04000000,-20(fp) ; nx = nx + 64 in exponent field
andd r7,r5 ; exponent of new y
; x is NaN or INF, return x-x
subl f0,f0 ; if x is INF, then INF-INF is NaN
L64: .double 0x0,0x43f00000 ; L64 = 2**64
LTINY: .double 0x0,0x04000000 ; LTINY = 2**-959