Clean up & prepare for ELF. Don't define isnan since it's in libc. Add a

powf (really a wrapper for pow).
This commit is contained in:
matt 2000-07-14 04:50:58 +00:00
parent 77c2009ea0
commit 584a2f53ee
10 changed files with 133 additions and 173 deletions

View File

@ -1,4 +1,4 @@
/* $NetBSD: n_argred.S,v 1.4 1999/07/02 15:37:35 simonb Exp $ */
/* $NetBSD: n_argred.S,v 1.5 2000/07/14 04:50:58 matt Exp $ */
/*
* Copyright (c) 1985, 1993
* The Regents of the University of California. All rights reserved.
@ -34,6 +34,8 @@
* @(#)argred.s 8.1 (Berkeley) 6/4/93
*/
#include <machine/asm.h>
/*
* libm$argred implements Bob Corbett's argument reduction and
* libm$sincos implements Peter Tang's double precision sin/cos.
@ -44,14 +46,8 @@
* method: true range reduction to [-pi/4,pi/4], P. Tang & B. Corbett
* S. McDonald, April 4, 1985
*/
.text
.globl libm$argred
.type libm$argred,@label
.globl libm$sincos
.type libm$sincos,@label
.align 1
libm$argred:
ENTRY(__libm_argred, 0)
/*
* Compare the argument with the largest possible that can
* be reduced by table lookup. r3 := |x| will be used in table_lookup .
@ -74,7 +70,8 @@ small_arg:
* r3 contains a F-format extension to the reduced argument;
* r4 contains a 0 or 1 corresponding to a sin or cos entry.
*/
libm$sincos:
ENTRY(__libm_sincos, 0)
/*
* Compensate for a cosine entry by adding one to the quadrant number.
*/
@ -109,7 +106,7 @@ sine:
/* cvtfd r4,r4 ... r5 = 0 after a polyd. */
addd2 r4,r0 # S(X) = beta + S(X)
addd2 r6,r0 # S(X) = X + S(X)
brb done
jbr done
cosine:
muld2 r6,r6 # Xsq = X * X
beql zero_arg
@ -126,8 +123,7 @@ done:
even:
rsb
.data
.align 2
_ALIGN_TEXT
sin_coef:
.double 0d-7.53080332264191085773e-13 # s7 = 2^-29 -1.a7f2504ffc49f8..
@ -264,8 +260,9 @@ leading:
twoOverPi:
.double 0d+6.36619772367581343076e-01
.text
.align 1
_ALIGN_TEXT
table_lookup:
muld3 r3,twoOverPi,r0
@ -298,7 +295,7 @@ abs2:
* p.0
*/
.text
.align 2
_ALIGN_TEXT
/*
* Only 256 (actually 225) bits of 2/pi are needed for VAX double
* precision; this was determined by enumerating all the nearest
@ -450,7 +447,7 @@ more1:
addl2 r1,r5
emul r0,r6,r5,r5
addl2 r0,r6
brb addbits1
jbr addbits1
signoff4:
emul r1,r6,$0,r4
@ -484,7 +481,7 @@ more2:
emul r0,r6,r4,r5
addl2 r0,r6
brb addbits2
jbr addbits2
signoff5:
emul r0,r6,r4,r5
@ -557,7 +554,7 @@ small:
*/
ashq r6,r9,r10
ashq r6,r8,r9
brb mult
jbr mult
/* p.13
*
* Test to see if the sign bit of r9 is on.
@ -571,7 +568,7 @@ tiny:
movl $32,r6
movq r8,r10
tstl r10
brb mult
jbr mult
/* p.14
*
* Test whether r9 is zero. It is probably
@ -597,7 +594,7 @@ tinier:
*/
ashq r1,r8,r10
ashq r1,r7,r9
brb mult
jbr mult
/* p.15
*
* The following code sets the reduced
@ -607,7 +604,7 @@ zero:
clrl r1
clrl r2
clrl r3
brw return
jbr return
/* p.16
*
* At this point, r0 contains the octant number,
@ -735,7 +732,7 @@ chop:
/* subw2 $217,r6 -S.McD */
subw2 $64,r6
insv r6,$7,$8,r3
brb return
jbr return
/* p.21
*
* The following code generates the appropriate

View File

@ -1,4 +1,4 @@
/* $NetBSD: n_atan2.S,v 1.3 1999/07/02 15:37:35 simonb Exp $ */
/* $NetBSD: n_atan2.S,v 1.4 2000/07/14 04:50:58 matt Exp $ */
/*
* Copyright (c) 1985, 1993
* The Regents of the University of California. All rights reserved.
@ -34,6 +34,8 @@
* @(#)atan2.s 8.1 (Berkeley) 6/4/93
*/
#include <machine/asm.h>
/*
* ATAN2(Y,X)
* RETURN ARG (X+iY)
@ -74,12 +76,7 @@
* atan2(y,x) returns the exact ARG(x+iy) nearly rounded.
*/
.text
.align 1
.globl _atan2
.type _atan2,@function
_atan2 :
.word 0x0ff4
ENTRY(atan2, 0x0fc0)
movq 4(ap),r2 # r2 = y
movq 12(ap),r4 # r4 = x
bicw3 $0x7f,r2,r0
@ -96,7 +93,7 @@ _atan2 :
movq r2,r0
bicw2 $0x8000,r0 # t = |y|
movq r0,r2 # y = |y|
brb begin
jbr begin
xnot1:
bicw3 $0x807f,r2,r11 # yexp
jeql yeq0 # if y=0 goto yeq0
@ -137,13 +134,13 @@ L10:
addd2 r10,r2 # 3y
addd2 r4,r2 # 3y+2x
divd2 r2,r0 # (2y-3x)/(2x+3y)
brw L60
jbr L60
L20:
cmpw r0,$0x3280 # t : 2**(-28)
jlss L80
clrq r6 # Hi=r6=0, Lo=r8=0
clrq r8
brw L60
jbr L60
L30:
movq $0xda7b2b0d63383fed,r6 # Hi=.46364760900080611433d0
movq $0xf0ea17b2bf912295,r8 # Lo=.10147340032515978826d-17
@ -153,7 +150,7 @@ L30:
addw2 $0x80,r4 # 2x
addd2 r2,r4 # 2x+y
divd2 r4,r0 # (2y-x)/(2x+y)
brb L60
jbr L60
L50:
movq $0x68c2a2210fda40c9,r6 # Hi=1.5707963267948966135d1
movq $0x06e0145c26332326,r8 # Lo=.22517417741562176079d-17
@ -161,7 +158,7 @@ L50:
bgeq L90
divd3 r2,r4,r0
bisw2 $0x8000,r0 # -x/y
brb L60
jbr L60
L40:
movq $0x68c2a2210fda4049,r6 # Hi=.78539816339744830676d0
movq $0x06e0145c263322a6,r8 # Lo=.11258708870781088040d-17
@ -183,7 +180,7 @@ L80:
ret
L90: # x >= 2**25
movq r6,r0
brb L80
jbr L80
pim:
subd3 r0,$0x68c2a2210fda4149,r0 # pi-t
bisw2 -4(fp),r0
@ -203,7 +200,8 @@ pio2:
resop:
movq $0x8000,r0 # propagate the reserved operand
ret
.align 2
_ALIGN_TEXT
ptable:
.quad 0xb50f5ce96e7abd60
.quad 0x51e44a42c1073e02

View File

@ -1,4 +1,4 @@
/* $NetBSD: n_cabs.S,v 1.2 1998/10/31 02:06:02 matt Exp $ */
/* $NetBSD: n_cabs.S,v 1.3 2000/07/14 04:50:58 matt Exp $ */
/*
* Copyright (c) 1985, 1993
* The Regents of the University of California. All rights reserved.
@ -34,6 +34,9 @@
* @(#)cabs.s 8.1 (Berkeley) 6/4/93
*/
#include <machine/asm.h>
.globl _C_LABEL(__libm_dsqrt_r5)
/*
* double precision complex absolute value
* CABS by W. Kahan, 9/7/80.
@ -43,29 +46,15 @@
* output is in r0:r1 (error less than 0.86 ulps)
*/
.text
.align 1
.globl _cabs
.type _cabs,@function
.globl _hypot
.type _hypot,@function
.globl _z_abs
.type _z_abs,@function
.globl libm$cdabs_r6
.type libm$cdabs_r6,@label
.globl libm$dsqrt_r5
.type libm$dsqrt_r5,@label
/* entry for c functions cabs and hypot */
_cabs:
_hypot:
.word 0x807c # save r2-r6, enable floating overflow
ALTENTRY(cabs)
ENTRY(hypot, 0x8040) # save r6, enable floating overflow
movq 4(ap),r0 # r0:1 = x
movq 12(ap),r2 # r2:3 = y
jmp cabs2
jbr cabs2
/* entry for Fortran use, call by: d = abs(z) */
_z_abs:
.word 0x807c # save r2-r6, enable floating overflow
ENTRY(z_abs, 0x8040) # save r6, enable floating overflow
movl 4(ap),r2 # indirect addressing is necessary here
movq (r2)+,r0 # r0:1 = x
movq (r2),r2 # r2:3 = y
@ -89,7 +78,7 @@ cont:
return:
ret
libm$cdabs_r6: # ENTRY POINT for cdsqrt
ENTRY(__libm_cdabs_r6,0) # ENTRY POINT for cdsqrt
# calculates a scaled (factor in r6)
# complex absolute value
@ -133,6 +122,7 @@ ordered:
addl2 r5,r4
cvtld r4,r2
addd2 r2,r0 # r0:1 = scaled x^2 + y^2
jmp libm$dsqrt_r5 # r0:1 = dsqrt(x^2+y^2)/2^r6
jmp _C_LABEL(__libm_dsqrt_r5)+2
# r0:1 = dsqrt(x^2+y^2)/2^r6
retsb:
rsb # error < 0.86 ulp

View File

@ -1,4 +1,4 @@
/* $NetBSD: n_cbrt.S,v 1.3 1999/07/02 15:37:35 simonb Exp $ */
/* $NetBSD: n_cbrt.S,v 1.4 2000/07/14 04:50:58 matt Exp $ */
/*
* Copyright (c) 1985, 1993
* The Regents of the University of California. All rights reserved.
@ -34,6 +34,8 @@
* @(#)cbrt.s 8.1 (Berkeley) 6/4/93
*/
#include <machine/asm.h>
/*
* double cbrt(double arg)
* W. Kahan, 10/13/80. revised 1/13/84 for keeping sign symmetry
@ -42,22 +44,12 @@
* Max error less than 0.667 ulps (unit in the last places)
*/
.globl _cbrt
.type _cbrt,@function
.globl _d_cbrt
.type _d_cbrt,@function
.globl _dcbrt_
.type _dcbrt_,@function
.text
.align 1
_cbrt:
_d_cbrt:
.word 0x00fc # save r2 to r7
ALTENTRY(cbrt)
ENTRY(d_cbrt, 0x00c0) # save r6 & r7
movq 4(ap),r0 # r0 = argument x
jmp dcbrt2
_dcbrt_:
.word 0x00fc # save r2 to r7
jbr dcbrt2
ENTRY(dcbrt_, 0x00c0) # save r6 & r7
movq *4(ap),r0 # r0 = argument x
dcbrt2: bicw3 $0x807f,r0,r2 # biased exponent of x
@ -93,12 +85,10 @@ dcbrt2: bicw3 $0x807f,r0,r2 # biased exponent of x
return:
ret # error less than 0.667 ulps
.data
.align 2
_ALIGN_TEXT
B : .long 721142941 # (86-0.03306235651)*(2^23)
C : .float 0f0.5428571429 # 19/35
D : .float 0f-0.7053061224 # -864/1225
E : .float 0f1.414285714 # 99/70
F : .float 0f1.607142857 # 45/28
G : .float 0f0.3571428571 # 5/14

View File

@ -1,4 +1,4 @@
/* $NetBSD: n_infnan.S,v 1.3 1999/07/02 15:37:35 simonb Exp $ */
/* $NetBSD: n_infnan.S,v 1.4 2000/07/14 04:50:58 matt Exp $ */
/*
* Copyright (c) 1985, 1993
* The Regents of the University of California. All rights reserved.
@ -33,10 +33,11 @@
*
* @(#)infnan.s 8.1 (Berkeley) 6/4/93
*/
.data
.align 2
#include <machine/asm.h>
.text
_sccsid:
.asciz "@(#)infnan.s 1.1 (Berkeley) 8/21/85; 8.1 (ucb.elefunt) 6/4/93"
.asciz "@(#)infnan.s\t1.1 (Berkeley) 8/21/85; 8.1 (ucb.elefunt) 6/4/93"
/*
* infnan(arg) int arg;
@ -48,16 +49,12 @@ _sccsid:
*/
.set EDOM,33
.set ERANGE,34
.text
.align 1
.globl _infnan
.type _infnan,@function
_infnan:
.word 0x0
ENTRY(infnan, 0)
cmpl 4(ap),$ERANGE
bneq 1f
movl $ERANGE,_errno
movl $ERANGE,_C_LABEL(errno)
brb 2f
1: movl $EDOM,_errno
1: movl $EDOM,_C_LABEL(errno)
2: emodd $0,$0,$0x8000,r0,r0 # generates the reserved operand fault
ret

View File

@ -1,4 +1,4 @@
/* $NetBSD: n_sincos.S,v 1.3 1999/07/02 15:37:35 simonb Exp $ */
/* $NetBSD: n_sincos.S,v 1.4 2000/07/14 04:50:58 matt Exp $ */
/*
* Copyright (c) 1985, 1993
* The Regents of the University of California. All rights reserved.
@ -49,12 +49,8 @@
* S. McDonald, April 4, 1985
*/
#include <machine/asm.h>
.text
.align 1
.globl _sin
.type _sin@function
_sin:
.word 0xffc # save r2-r11
ENTRY(sin, 0xfc0)
movq 4(ap),r0
bicw3 $0x807f,r0,r2
beql 1f # if x is zero or reserved operand then return x
@ -70,9 +66,9 @@ _sin:
/*
* Entered by sine ; save 0 in r4 .
*/
jsb libm$argred
jsb _C_LABEL(__libm_argred)+2
movl $0,r4
jsb libm$sincos
jsb _C_LABEL(__libm_sincos)+2
bispsw (sp)+
1: ret
@ -82,12 +78,8 @@ _sin:
* method: true range reduction to [-pi/4,pi/4], P. Tang & B. Corbett
* S. McDonald, April 4, 1985
*/
.text
.align 1
.globl _cos
.type _cos,@function
_cos:
.word 0xffc # save r2-r11
ENTRY(cos, 0x0fc0)
movq 4(ap),r0
bicw3 $0x7f,r0,r2
cmpw $0x8000,r2
@ -104,8 +96,8 @@ _cos:
/*
* Entered by cosine ; save 1 in r4 .
*/
jsb libm$argred
jsb _C_LABEL(__libm_argred)+2
movl $1,r4
jsb libm$sincos
jsb _C_LABEL(__libm_sincos)+2
bispsw (sp)+
1: ret

View File

@ -1,4 +1,4 @@
/* $NetBSD: n_sqrt.S,v 1.2 1998/10/31 02:06:02 matt Exp $ */
/* $NetBSD: n_sqrt.S,v 1.3 2000/07/14 04:50:58 matt Exp $ */
/*
* Copyright (c) 1985, 1993
* The Regents of the University of California. All rights reserved.
@ -34,6 +34,8 @@
* @(#)sqrt.s 8.1 (Berkeley) 6/4/93
*/
#include <machine/asm.h>
/*
* double sqrt(arg) revised August 15,1982
* double arg;
@ -45,31 +47,27 @@
* entry points:_d_sqrt address of double arg is on the stack
* _sqrt double arg is on the stack
*/
.text
.align 1
.globl _sqrt
.type _sqrt,@function
.globl _d_sqrt
.type _d_sqrt,@function
.globl libm$dsqrt_r5
.type libm$dsqrt_r5,@label
.set EDOM,33
_d_sqrt:
.word 0x003c # save r5,r4,r3,r2
ENTRY(d_sqrt, 0x003c) # save r5,r4,r3,r2
movq *4(ap),r0
jmp dsqrt2
_sqrt:
.word 0x003c # save r5,r4,r3,r2
jbr dsqrt2
ENTRY(sqrt, 0x003c) # save r5,r4,r3,r2
movq 4(ap),r0
dsqrt2: bicw3 $0x807f,r0,r2 # check exponent of input
jeql noexp # biased exponent is zero -> 0.0 or reserved
bsbb libm$dsqrt_r5
bsbb _C_LABEL(__libm_dsqrt_r5)+2
noexp: ret
/* **************************** internal procedure */
libm$dsqrt_r5: /* ENTRY POINT FOR cdabs and cdsqrt */
ALTENTRY(__libm_dsqrt_r5)
nop
nop
/* ENTRY POINT FOR cdabs and cdsqrt */
/* returns double square root scaled by */
/* 2^r6 */
@ -117,8 +115,8 @@ libm$dsqrt_r5: /* ENTRY POINT FOR cdabs and cdsqrt */
rsb # DONE!
nonpos:
jneq negarg
ret # argument and root are zero
ret # argument and root are zero
negarg:
pushl $EDOM
calls $1,_infnan # generate the reserved op fault
calls $1,_C_LABEL(infnan) # generate the reserved op fault
ret

View File

@ -1,4 +1,4 @@
/* $NetBSD: n_support.S,v 1.2 1999/07/02 15:37:35 simonb Exp $ */
/* $NetBSD: n_support.S,v 1.3 2000/07/14 04:51:00 matt Exp $ */
/*
* Copyright (c) 1985, 1993
* The Regents of the University of California. All rights reserved.
@ -33,10 +33,11 @@
*
* @(#)support.s 8.1 (Berkeley) 6/4/93
*/
.data
.align 2
#include <machine/asm.h>
.text
_sccsid:
.asciz "@(#)support.s 1.3 (Berkeley) 8/21/85; 8.1 (ucb.elefunt) 6/4/93"
.asciz "@(#)support.s\t1.3 (Berkeley) 8/21/85; 8.1 (ucb.elefunt) 6/4/93"
/*
* copysign(x,y),
@ -49,14 +50,10 @@ _sccsid:
*/
/*
* double copysign(x,y)
* double x,y;
* double copysign(double x,double y)
*/
.globl _copysign
.text
.align 1
_copysign:
.word 0x4
ENTRY(copysign, 0)
movq 4(ap),r0 # load x into r0
bicw3 $0x807f,r0,r2 # mask off the exponent of x
beql Lz # if zero or reserved op then return x
@ -66,14 +63,9 @@ _copysign:
Lz: ret
/*
* double logb(x)
* double x;
* double logb(double x);
*/
.globl _logb
.text
.align 1
_logb:
.word 0x0
ENTRY(logb, 0)
bicl3 $0xffff807f,4(ap),r0 # mask off the exponent of x
beql Ln
ashl $-7,r0,r0 # get the bias exponent
@ -86,14 +78,9 @@ Ln: movq 4(ap),r0 # r0:1 = x (zero or reserved op)
1: ret
/*
* long finite(x)
* double x;
* long finite(double x);
*/
.globl _finite
.text
.align 1
_finite:
.word 0x0000
ENTRY(finite, 0)
bicw3 $0x7f,4(ap),r0 # mask off the mantissa
cmpw r0,$0x8000 # to see if x is the reserved op
beql 1f # if so, return FALSE (0)
@ -102,16 +89,27 @@ _finite:
1: clrl r0
ret
/* int isnan(double x);
*/
#if 0
ENTRY(isnan, 0)
clrl r0
ret
#endif
/* int isnanf(float x);
*/
ENTRY(isnanf, 0)
clrl r0
ret
/*
* double scalb(x,N)
* double x; double N;
*/
.globl _scalb
.set ERANGE,34
.text
.align 1
_scalb:
.word 0x3c
ENTRY(scalb, 0)
movq 4(ap),r0
bicl3 $0xffff807f,r0,r3
beql ret1 # 0 or reserved operand
@ -129,7 +127,7 @@ _scalb:
addl2 r2,r0
ret
ovfl: pushl $ERANGE
calls $1,_infnan # if it returns
calls $1,_C_LABEL(infnan) # if it returns
bicw3 $0x7fff,4(ap),r2 # get the sign of input arg
bisw2 r2,r0 # re-attach the sign to r0/1
ret
@ -142,12 +140,9 @@ ret1: ret
* DOUBLE PRECISION (VAX D format 56 bits)
* CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/8/85.
*/
.globl _drem
.set EDOM,33
.text
.align 1
_drem:
.word 0xffc
ENTRY(drem, 0x0fc0)
subl2 $12,sp
movq 4(ap),r0 #r0=x
movq 12(ap),r2 #r2=y
@ -192,13 +187,13 @@ C2:
muld2 r8,r6 #n*t1
subd2 r6,r0 #x-n*t1
subd2 r4,r0 #(x-n*t1)-n*(t-t1)
brb loop
jbr loop
E1:
movw -4(fp),r6 #r6=nx
beql C3 #if nx=0 goto C3
addw2 r6,r0 #x:=x*2**57 scale up x by nx
movw $0,-4(fp) #clear nx
brb loop
jbr loop
C3:
movq r2,r4 #r4 = y
subw2 $0x80,r4 #r4 = y/2
@ -223,7 +218,7 @@ C5:
ret
Rop: #Reserved operand
pushl $EDOM
calls $1,_infnan #generate reserved op fault
calls $1,_C_LABEL(infnan) #generate reserved op fault
ret
Ret:
movq $0x8000,r0 #propagate reserved op

View File

@ -1,4 +1,4 @@
/* $NetBSD: n_tan.S,v 1.3 1999/07/02 15:37:35 simonb Exp $ */
/* $NetBSD: n_tan.S,v 1.4 2000/07/14 04:51:00 matt Exp $ */
/*
* Copyright (c) 1985, 1993
* The Regents of the University of California. All rights reserved.
@ -34,6 +34,8 @@
* @(#)tan.s 8.1 (Berkeley) 6/4/93
*/
#include <machine/asm.h>
/* This is the implementation of Peter Tang's double precision
* tangent for the VAX using Bob Corbett's argument reduction.
*
@ -46,12 +48,7 @@
* method: true range reduction to [-pi/4,pi/4], P. Tang & B. Corbett
* S. McDonald, April 4, 1985
*/
.text
.align 1
.globl _tan
.type _tan,@function
_tan: .word 0xffc # save r2-r11
ENTRY(tan, 0x0fc0) # save r6-r11
movq 4(ap),r0
bicw3 $0x807f,r0,r2
beql 1f # if x is zero or reserved operand then return x
@ -64,7 +61,7 @@ _tan: .word 0xffc # save r2-r11
* Clear the IV & FU bits.
*/
bicpsw $0x0060
jsb libm$argred
jsb _C_LABEL(__libm_argred)+2
/*
* At this point,
* r0 contains the quadrant number, 0, 1, 2, or 3;
@ -79,7 +76,7 @@ _tan: .word 0xffc # save r2-r11
* Call sine. r4 = 0 implies sine.
*/
movl $0,r4
jsb libm$sincos
jsb _C_LABEL(__libm_sincos)+2
/*
* Save sin(x) in r11/r10 .
*/
@ -90,7 +87,7 @@ _tan: .word 0xffc # save r2-r11
movq (sp)+,r0
movq (sp)+,r2
movl $1,r4
jsb libm$sincos
jsb _C_LABEL(__libm_sincos)+2
divd3 r0,r10,r0
bispsw (sp)+
1: ret

View File

@ -1,4 +1,4 @@
/* $NetBSD: n_pow.c,v 1.4 1999/07/02 15:37:37 simonb Exp $ */
/* $NetBSD: n_pow.c,v 1.5 2000/07/14 04:51:01 matt Exp $ */
/*
* Copyright (c) 1985, 1993
* The Regents of the University of California. All rights reserved.
@ -126,6 +126,12 @@ const static double zero=0.0, one=1.0, two=2.0, negone= -1.0;
static double pow_P __P((double, double));
float powf(x,y)
float x,y;
{
return pow((double) x, (double) (y));
}
double pow(x,y)
double x,y;
{