projects / unix-history / blob
? search:
summary | tags | clone url | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
add Berkeley specific headers
[unix-history] / usr / src / lib / libm / national / support.s
;
; Copyright (c) 1985 Regents of the University of California.
; All rights reserved.
;
; Redistribution and use in source and binary forms are permitted
; provided that this notice is preserved and that due credit is given
; to the University of California at 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'' without express or implied warranty.
;
; 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.2 (Berkeley) %G%
;
; IEEE recommended functions
;
; double copysign(x,y)
; double x,y;
; IEEE 754 recommended function, return x*sign(y)
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
;
.vers 2
.text
.align 2
.globl _copysign
_copysign:
movl 4(sp),f0
movd 8(sp),r0
movd 16(sp),r1
xord r0,r1
andd 0x80000000,r1
cmpqd 0,r1
beq end
negl f0,f0
end: ret 0
;
; double logb(x)
; double x;
; 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.
;
.vers 2
.text
.align 2
.globl _logb
_logb:
;
; 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
movd 20,r2
extd r2,r0,r1,11 ; extract the exponent of x
cmpqd 0,r1 ; if exponent bits = 0, goto L3
beq L3
cmpd 0x7ff,r1
beq L2 ; if exponent bits = 0x7ff, goto L2
L1: subd 1023,r1 ; unbias the exponent
movdl r1,f0 ; convert the exponent to floating value
ret 0
;
; x is INF or NaN, simply return x
;
L2:
movl 4(sp),f0 ; logb(+inf)=+inf, logb(NaN)=NaN
ret 0
;
; x is 0 or subnormal
;
L3:
movl 4(sp),f0
cmpl 0f0,f0
beq L5 ; x is 0 , goto L5 (return -inf)
;
; Now x is subnormal
;
mull L64,f0 ; scale up f0 with 2**64
movl f0,tos
movd tos,r0
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
ret 0
;
; x is 0, return logb(0)= -INF
;
L5:
movl 0f1.0e300,f0
mull 0f-1.0e300,f0 ; multiply two big numbers to get -INF
ret 0
;
; double rint(x)
; double x;
; ... delivers integer nearest x in direction of prevailing rounding
; ... mode
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
; Note: subnormal number (if implemented) will be taken care of.
;
.vers 2
.text
.align 2
.globl _rint
_rint:
;
movd 8(sp),r0
movd 20,r2
extd r2,r0,r1,11 ; extract the exponent of x
cmpd 0x433,r1
ble itself
movl L52,f2 ; f2 = L = 2**52
cmpqd 0,r0
ble L1
negl f2,f2 ; f2 = s = copysign(L,x)
L1: addl f2,f0 ; f0 = x + s
subl f2,f0 ; f0 = f0 - s
ret 0
itself: movl 4(sp),f0
ret 0
L52: .double 0x0,0x43300000 ; L52=2**52
;
; int finite(x)
; double x;
; 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.
;
.vers 2
.text
.align 2
.globl _finite
_finite:
movd 4(sp),r1
andd 0x800fffff,r1
cmpd 0x7ff00000,r1
sned r0 ; r0=0 if exponent(x) = 0x7ff
ret 0
;
; double scalb(x,N)
; double x; int N;
; IEEE 754 recommended function, return x*2**N by adjusting
; exponent of x.
; Coded by K.C. Ng in National 32k assembler, 11/9/85.
; Note: subnormal number (if implemented) will be taken care of
;
.vers 2
.text
.align 2
.globl _scalb
_scalb:
;
; if x=0 return 0
;
movl 4(sp),f0
cmpl 0f0,f0
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
movd 20,r2 ; r2 = 20
extd r2,r0,r1,11 ; extract the exponent of x in r1
cmpd 0x7ff,r1
;
; if exponent of x is 0x7ff, then x is NaN or INF; simply return x
;
beq end
cmpqd 0,r1
;
; if exponent of x is zero, then x is subnormal; goto L19
;
beq 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
bge underflow
cmpd 2047,r1 ; if new exponent >= 2047 goto overflow
ble overflow
insd r2,r1,r0,11 ; insert the new exponent
movd r0,tos
movd 8(sp),tos
movl tos,f0 ; return x*2**N
end: ret 0
inof: bcs underflow ; negative int overflow if Carry bit is set
overflow:
andd 0x80000000,r0 ; keep the sign of x
ord 0x7fe00000,r0 ; set x to a huge number
movd r0,tos
movqd 0,tos
movl tos,f0
mull 0f1.0e300,f0 ; multiply two huge number to get overflow
ret 0
underflow:
addd 64,r1 ; add 64 to exonent to see if it is subnormal
cmpqd 0,r1
bge zero ; underflow to zero
insd r2,r1,r0,11 ; insert the new exponent
movd r0,tos
movd 8(sp),tos
movl tos,f0
mull L30,f0 ; readjust x by multiply it with 2**-64
ret 0
zero: andd 0x80000000,r0 ; keep the sign of x
ord 0x00100000,r0 ; set x to a tiny number
movd r0,tos
movqd 0,tos
movl tos,f0
mull 0f1.0e-300,f0 ; underflow to 0 by multipling two tiny nos.
ret 0
L19: ; subnormal number
mull L32,f0 ; scale up x by 2**64
movl f0,tos
movd tos,r0
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
cmpqd 0,r1
bge underflow
cmpd 2047,r1
ble overflow
insd r2,r1,r0,11 ; insert back the incremented exponent
movd r0,tos
movd 8(sp),tos
movl tos,f0
end: ret 0
L30: .double 0x0,0x3bf00000 ; floating point 2**-64
L32: .double 0x0,0x43f00000 ; floating point 2**64
;
; double drem(x,y)
; double x,y;
; 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.
;
.vers 2
.text
.align 2
.globl _drem
_drem:
;
; glossaries:
; r2 = high part of x
; r3 = exponent of x
; r4 = high part of y
; r5 = exponent of y
; r6 = sign of x
; r7 = constant 0x7ff00000
;
; 16(fp) : y
; 8(fp) : x
; -12(fp) : adjustment on y when y is subnormal
; -16(fp) : fsr
; -20(fp) : nx
; -28(fp) : t
; -36(fp) : t1
; -40(fp) : nf
;
;
enter [r3,r4,r5,r6,r7],40
movl f6,tos
movl f4,tos
movl 0f0,-12(fp)
movd 0,-20(fp)
movd 0,-40(fp)
movd 0x7ff00000,r7 ; initialize r7=0x7ff00000
movd 12(fp),r2 ; r2 = high(x)
movd r2,r3
andd r7,r3 ; r3 = xexp
cmpd r7,r3
; if x is NaN or INF goto L1
beq L1
movd 20(fp),r4
bicd [31],r4 ; r4 = high part of |y|
movd r4,20(fp) ; y = |y|
movd r4,r5
andd r7,r5 ; r5 = yexp
cmpd r7,r5
beq L2 ; if y is NaN or INF goto L2
cmpd 0x04000000,r5 ;
bgt L3 ; if y is tiny goto L3
;
; now y != 0 , x is finite
;
L10:
movd r2,r6
andd 0x80000000,r6 ; r6 = sign(x)
bicd [31],r2 ; x <- |x|
sfsr r1
movd r1,-16(fp) ; save fsr in -16(fp)
bicd [5],r1
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 r0,-28(fp) ;
movd r4,-24(fp) ; t = y
movd r4,-32(fp) ;
movd r1,-36(fp) ; t1 = y with trialing 27 zeros
movd 0x01900000,r1 ; r1 = 25 in exponent field
LOOP:
movl 8(fp),f0 ; f0 = x
movl 16(fp),f2 ; f2 = y
cmpl f0,f2
ble fnad ; goto fnad (final adjustment) if x <= y
movd r4,-32(fp)
movd r3,r0
subd r5,r0 ; xexp - yexp
subd r1,r0 ; r0 = xexp - yexp - m25
cmpqd 0,r0 ; r0 > 0 ?
bge 1f
addd r4,r0 ; scale up (high) y
movd r0,-24(fp) ; scale up t
movl -28(fp),f2 ; t
movd r0,-32(fp) ; scale up t1
1:
movl -36(fp),f4 ; t1
movl f0,f6
divl f2,f6 ; f6 = x/t
floorld f6,r0 ; r0 = [x/t]
movdl r0,f6 ; f6 = n
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)
; update xexp
movl f0,8(fp)
movd 12(fp),r3
andd r7,r3
jump LOOP
fnad:
cmpqd 0,-20(fp) ; 0 = nx?
beq final
mull -12(fp),8(fp) ; scale up x the same amount as y
movd 0,-20(fp)
movd 12(fp),r2
movd r2,r3
andd r7,r3 ; update exponent of x
jump LOOP
final:
movl 16(fp),f2 ; f2 = y (f0=x, r0=n)
subd 0x100000,r4 ; high y /2
movd r4,-24(fp)
movl -28(fp),f4 ; f4 = y/2
cmpl f0,f4 ; x > y/2 ?
bgt 1f
bne 2f
andd 1,r0 ; n is odd or even
cmpqd 0,r0
beq 2f
1:
subl f2,f0 ; x = x - y
2:
cmpqd 0,-40(fp)
beq 3f
divl -12(fp),f0 ; scale down the answer
3:
movl f0,tos
xord r6,tos
movl tos,f0
movd -16(fp),r0
lfsr r0 ; restore the fsr
end: movl tos,f4
movl tos,f6
exit [r3,r4,r5,r6,r7]
ret 0
;
; y is NaN or INF
;
L2:
movd 16(fp),r0 ; r0 = low part of y
andd 0xfffff,r4 ; r4 = high part of y & 0x000fffff
ord r4,r0
cmpqd 0,r0
beq L4
movl 16(fp),f0 ; y is NaN, return y
jump end
L4: movl 8(fp),f0 ; y is inf, return x
jump end
;
; exponent of y is less than 64, y may be zero or subnormal
;
L3:
movl 16(fp),f0
cmpl 0f0,f0
bne L5
divl f0,f0 ; y is 0, return NaN by doing 0/0
jump end
;
; subnormal y or tiny y
;
L5:
movd 0x04000000,-20(fp) ; nx = 64 in exponent field
movl L64,f2
movl f2,-12(fp)
mull f2,f0
cmpl f0,LTINY
bgt L6
mull f2,f0
addd 0x04000000,-20(fp) ; nx = nx + 64 in exponent field
mull f2,-12(fp)
L6:
movd -20(fp),-40(fp)
movl f0,16(fp)
movd 20(fp),r4
movd r4,r5
andd r7,r5 ; exponent of new y
jump L10
;
; x is NaN or INF, return x-x
;
L1:
movl 8(fp),f0
subl f0,f0 ; if x is INF, then INF-INF is NaN
ret 0
L64: .double 0x0,0x43f00000 ; L64 = 2**64
LTINY: .double 0x0,0x04000000 ; LTINY = 2**-959
Continuous source code history from original PDP-7 UNIX to FreeBSD 2, the first fully unencumbered FreeBSD release.