Bell 32V development

Work on file usr/src/libnm/cbrt.s
Work on file usr/src/libnm/sqrt.s

Co-Authored-By: John Reiser <jfr@research.uucp>
Synthesized-from: 32v

diff --git a/usr/src/libnm/cbrt.s b/usr/src/libnm/cbrt.s
new file mode 100644 (file)
index 0000000..007f133
--- /dev/null
+++ b/usr/src/libnm/cbrt.s
@@ -0,0 +1,53 @@
+# double cbrt(arg)
+# double arg
+# no error exits
+#method: range reduction to [1/8,1], poly appox, newtons method
+# J F Jarvis, August 10,1978
+.globl _cbrt
+.text
+.align 1
+_cbrt:
+       .word   0x00c0
+       clrl    r3
+       movd    4(ap),r4
+       jgtr    range
+       jeql    retz
+       mnegd   r4,r4   # |arg| in r0,r1
+       movl    $0x100,r3       # sign bit of result
+range:
+       extzv   $7,$8,r4,r6
+       insv    $128,$7,$8,r4   # 0.5<= frac: r4,r5 <1.0
+       clrl    r7
+       ediv    $3,r6,r6,r7     # r6= expnt/3; r7= expnt%3
+       addb2   $86,r6
+       bisl2   r3,r6   # sign,exponent of result
+       polyf   r4,$3,pcoef     # initial estimate is Hart&Cheney CBRT 0642
+                                               # D=4.1
+       muld3   r0,r0,r2        # Newtons method, iteration 1, H&C 6.1.10
+       divd3   r2,r4,r2
+       subd3   r2,r0,r2
+       muld2   third,r2
+       subd2   r2,r0   # D=8.2
+       muld3   r0,r0,r2        # iteration 2
+       divd3   r2,r4,r2
+       subd3   r2,r0,r2
+       muld2   third,r2
+       subd2   r2,r0
+       muld2   hc[r7],r0       # set range
+       insv    r6,$7,$9,r0     # set sign,exponent
+       ret
+retz:
+       clrd    r0
+       ret
+.data
+.align 3
+third: .double 0d0.33333333333333333333e+0
+hc:
+       .double 0d1.25992104989487316476e+0
+       .double 0d1.58740105196819947475e+0
+       .double 0d1.0e+0
+pcoef:
+       .float 0f0.1467073818e+0
+       .float 0f-0.5173964673e+0
+       .float 0f0.9319858515e+0
+       .float 0f0.4387762363e+0

diff --git a/usr/src/libnm/sqrt.s b/usr/src/libnm/sqrt.s
new file mode 100644 (file)
index 0000000..22b9f3e
--- /dev/null
+++ b/usr/src/libnm/sqrt.s
@@ -0,0 +1,46 @@
+# double sqrt(arg)
+# double arg
+# if(arg<0.0) { _errno=EDOM; return(0.0) }
+# J. F. Jarvis August 2, 1978
+.set   EDOM,98
+.text
+.align 1
+.globl _sqrt
+.globl _errno
+_sqrt:
+       .word   0x0c0
+       movd    4(ap),r4
+       jgtr    range
+       jeql    retz
+       movl    $EDOM,_errno
+retz:
+       clrd    r0
+       ret
+#
+range:
+       extzv   $7,$8,r4,r6     # r6 is exponent of arg
+       insv    $128,$7,$8,r4   # r2: 0.5<=fraction(arg)<1.0
+       incl    r6
+       clrl    r7
+       ediv    $2,r6,r6,r7     # r6=(exp+1)/2; r7=(exp+1)%2
+       addb2   $64,r6  # r6 is correct exponent for result
+       polyf   r4,$4,pcoef     # init estimate of sqrt(frac)
+                                               # Hart&Cheney SQRT 0132 D=5.1
+       divd3   r0,r4,r2        # Newtons method, 2 iterations
+       addd2   r2,r0
+       muld2   $0d0.5e+0,r0
+       divd3   r0,r4,r2        # Hart&Cheney 6.1.7
+       addd2   r2,r0   #d=21 at exit
+       muld2   hc[r7],r0       # *sqrt(2) requ for even org exp.
+       insv    r6,$7,$8,r0     # insert correct exp.
+       ret
+.data
+.align 3
+pcoef:
+       .float 0f-0.1214683825e+0
+       .float 0f0.5010420763e+0
+       .float 0f-0.9093210498e+0
+       .float 0f0.1300669049e+1
+       .float 0f0.2290699453e+0
+hc:            .double 0d0.35355339059327376220e+0     # sqrt(2)/4
+               .double 0d0.5e+0