doprnt.s

<B>Digital的Unix操作系统VAX 4.2源码</B>

	ashp	$-1,$fltprec+2,12(r5),$5,$fltprec,(r5)	# Round on the 17th digit to 
						# get the 16 significant digits
	bvc	1f				# slr001 br if no overflow
/*
 *		SLR001
 *  This will fix the problem when the rounding causes the most signif. digit
 *  to cause an overflow.
 */
	ashp	$-1,$fltprec+2,12(r5),$5,$fltprec,(r5)	# Round on the 17th digit to 
						# get the 16 significant digits
	incl	dec_exp(r7)			# Adjust the exponent to
						# account for the rounding.

1:	bisb2	sig(r7),7(r5)			# Set sign
finis:
	addl2	$local_frame, sp		# restore stack pointer
	movl	r7, r1			# restore frame pointer
	rsb				# return to caller

zero:
 #	movl	string_addr(r7), r1	# get string address
 #	movc5	$0, (sp), $^a/0/, sig_digits(r7), (r1) # zero fill string
	clrl	offset(r7)		# clear offset and exponent
	movl	$1,dec_exp(r7)
	clrl	sig(r7)			# clear sign
	cvtlp	$0,$fltprec,(r5)
	movl	r7, r1			# restore frame pointer
	rsb				# return to caller


/* this is the subroutine which does the multiple
* precision multiplies. it is called with bsb or jsb
* with r2 containing a pointer to an appropriate
* entry in the power-of-ten table. binexp & binnum
*  are multiplied by the reciprocal of this entry, 
* with the results going to
* binexp & binnum, and decexp is updated with the
* power_of_ten value.
* this routine clobbers r0-r1, r3-r6, and changes r2.
*/

 #ots$$cvt_mul::
rmul:					# entry point	

		/* find the reciprocal table entry pointed to by r2.
		* r2 contains the base (t0) plus an "index". the
		* reciprocal entry has an adr of "t0-index" which
		* is calculated by 2*t0-(t0+index), or 2*t0-r2.
		*/

	moval	t0, r1			# get base adr
	addl2	r1, r1			# 2*base
	subl3	r2, r1, r2		# get adr of reciprocal entry

	clrl	prodf-4(r8)		# init product
	clrq	prodf+0(r8)
	clrq	prodf+8(r8)

	clrq	cry+0(r8)		# clear carries
	clrq	cry+8(r8)

		/* this macro has the function r=a*b, with the carries
		* going into the 4 l-words at "cry". a and b are
		* unsigned long-words. r is an unsigned double long-word.
		*  removing this macro definition (which is only used once), 
		*  and expanding the code where it is used, obscures the function.
		*/

/* macro lmul a, b, r, ?l1, ?l2, ?l3
*	movl	a, r0			# get first operand
*	jeql	l3			# skip if zero
*	movl	b, r1			# get second operand
*	jeql	l3			# skip if zero
*	emul	r0, r1, $0, r0		# form product of a and b
*	tstl	a
*	jgeq	l1
*	addl2	b, r1			# if a<0, fixup for neg sign
*l1:	tstl	b
*	jgeq	l2
*	addl2	a, r1			# if b<0, fixup for neg sign
*l2:	addl2	r0, r			# add low product into result
*	adwc	r1, 4+r			# add hi product into result
*	adwc	$0, cry+8-prodf+r	# and save carries
*l3:
*	.endm lmul
*/
		/* the following loop forms all the cross-products
		* required for a 4-long-word by 4-long-word multiply.
		* only the high 4 long-words are accumulated. the byte
		* table at "bytab" shows the indicies used for the
		* long-word operands and the resulting double-long-
		* word products.
		*/

	movaw	bytab, r3		# init byte table index
bytlup:	cvtbl	(r3)+, r4		# setup 1st index
	jlss	bytdun			# and test for end
	cvtbl	(r3)+, r5		# setup 2nd index
	cvtbl	(r3)+, r6		# setup 3rd index

/*	lmul	binnum(r8)[r4], 0(r2)[r5], prodf_4(r8)[r6]  */

		movl	binnum(r8)[r4], r0	# get first operand
		jeql	l3			# skip if zero
		movl	0(r2)[r5], r1		# get second operand
		jeql	l3			# skip if zero
		emul	r0, r1, $0, r0		# form product of a and b
		tstl	binnum(r8)[r4]
		jgeq	l1
		addl2	0(r2)[r5], r1		# if a<0, fixup for neg sign
l1:		tstl    0(r2)[r5]	
		jgeq	l2
		addl2	binnum(r8)[r4], r1	# if b<0, fixup for neg sign
l2:		addl2	r0, prodf_4(r8)[r6]	# add low product into result
		adwc	r1, 4+prodf_4(r8)[r6]	# add hi product into result
 #		adwc	$0, cry+8-prodf+prodf_4(r8)[r6]	# and save carries
/*
*	This had to be hard code because of preprocessor problems
*/
			adwc	$0,0x30(r8)[r6]
l3:

	jbr	bytlup			# loop

bytdun:
/*	incl	cry+0(r8)		; small extra fudge */
	addl2	cry+0(r8), prodf+0(r8)	# put carries into sum
	adwc	cry+4(r8), prodf+4(r8)
	adwc	cry+8(r8), prodf+8(r8)
	adwc	cry+12(r8), prodf+12(r8)

	extzv	$31, $1, prodf+12(r8), r1	# get normalize bit

		/* normalized operands cannot produce a result
		* un-normalized by more than one bit position. so
		* if norm_bit=1, shift left by 0
		* if norm_bit=0, shift left by 1 and sub 1 from exp
		*/

	jneq	nosub1			# xfer if norm_bit = 1
	decl	binexp(r8)		# norm_bit = 0, sub 1 from exponent

		# move the product from prodf to binnum, normalizing
		# it one bit position if required.

nosub1:	addl2	$31, r1				# do extv's from bit 31 or 32

	extv	r1, $32, prodf-4(r8), binnum+0(r8)	# shift left 0 or 1 bit
	extv	r1, $32, prodf+0(r8), binnum+4(r8)
	extv	r1, $32, prodf+4(r8), binnum+8(r8)
	extv	r1, $32, prodf+8(r8), binnum+12(r8)

	cvtwl	t_bexp(r2), r1		# extract binary exponent
	addl2	r1, binexp(r8)		# add exponents for mul

		/* the binary exponent moves toward zero while the
		* decimal exponent moves away from zero by an amount
		* about equal to log(bin exp).
		*/

	cvtwl	t_dexp(r2), r1		# get equivalent decimal exponent
	subl2	r1, dec_exp(r7)		# and sub it from result exp

	rsb				# return

 #	.sbttl	tables


 #.macro number a1, a2, a3, a4, a5, a6, a7
 #	.long ^x'a5+<<^x'a6@-31>&1>, ^x'a4, ^x'a3, ^x'a2
 #	.word ^d<a1>, ^d<a7>
 #.endm number

#define number(a1,a2,a3,a4,a5,a6,a7)	 \
	.long 0x/**/a5+((0x/**/a6>>31)&1), 0x/**/a4,0x/**/a3,0x/**/a2  \
 ;	.word a1,a7

		# this macro creates a table entry of the following form:

		#	.long < least sig bits>	:     0(r2)
		#	.long <   ........    >	:     4(r2)
		#	.long <   ........    >	:     8(r2)
		#	.long < most sig bits >	:    12(r2)
		#	.word < binary exp    >	:t_bexp(r2)
		#	.word < decimal exp   >	:t_dexp(r2)

 #t_bexp=16	# the binary exponent is bytes 16-17 of each table entry
 #t_dexp=18	# the decimal exponent is bytes 18-19


/*            value = fraction * 2**power_of_2 = 10**power_of_10

*	the hex fraction is stored as a 4 long-word unsigned integer,
*	left normalized, with the binary point left of bit 31
*	of the most significant long-word.

*	the fraction is guaranteed correct for the four high-order
*	long-words. about 16 bits of the fifth low-order long-word may
*	be in error. the check line at the bottom of the table is
*	the product of the first and last table entries. it would
*	equal exactly 1.0 if every bit of the 5 long-words were correct.

*	      decimal,<-------5 long-word hex fraction----------->, decimal
*	       power ,<--msb--------------------------------lsb-->,  power
*	       of 2                                                  of 10
*
*
*/

	.align	2

tsmall:
 #	number (-27213,d986c20b,686da869,5d1d4fd8,5b05f4c2,eef0fb87,-8192)
	number (-13606,a6dd04c8,d2ce9fde,2de38123,a1c3cffc,203028da,-4096)
	number ( -6803,ceae534f,34362de4,492512d4,f2ead2cb,8263aa10,-2048)
	number ( -3401,a2a682a5,da57c0bd,87a60158,6bd3f698,f53e881e,-1024)
	number ( -1700,9049ee32,db23d21c,7132d332,e3f204d4,e73177c2, -512)
	number (  -850,c0314325,637a1939,fa911155,fefb5308,a23e2b15, -256)
	number (  -425,ddd0467c,64bce4a0,ac7cb3f6,d05ddbde,e26ca3df, -128)
	number (  -212,a87fea27,a539e9a5,3f2398d7,47b36224,2a1fed70,  -64)
 # tm32:
	number (  -106,cfb11ead,453994ba,67de18ed,a5814af2,b5b1a20,  -32)
/*	number (  -102,81ceb32c,4b43fcf4,80eacf94,8770ced7,4718f05a,  -31)
*	number (   -99,a2425ff7,5e14fc31,a1258379,a94d028d,18df2c73,  -30)
*	number (   -96,cad2f7f5,359a3b3e, 96ee458,13a04330,5f16f793,  -29)
*	number (   -93,fd87b5f2,8300ca0d,8bca9d6e,188853fc,76dcb57b,  -28)
*	number (   -89,9e74d1b7,91e07e48,775ea264,cf55347d,ca49f16f,  -27)
*	number (   -86,c6120625,76589dda,95364afe,32a819d,3cdc6dcd,  -26)
*	number (   -83,f79687ae,d3eec551,3a83ddbd,83f52204,8c138944,  -25)
*	number (   -79,9abe14cd,44753b52,c4926a96,72793542,d78c35ce,  -24)
*	number (   -76,c16d9a00,95928a27,75b7053c,f178293,8d6f434a,  -23)
*	number (   -73,f1c90080,baf72cb1,5324c68b,12dd6338,70cb1420,  -22)
*	number (   -69,971da050,74da7bee,d3f6fc16,ebca5e03,467eec97,  -21)
*	number (   -66,bce50864,92111aea,88f4bb1c,a6bcf584,181ea7c0,  -20)
*	number (   -63,ec1e4a7d,b69561a5,2b31e9e3,d06c32e5,1e2651b1,  -19)
*	number (   -59,9392ee8e,921d5d07,3aff322e,62439fcf,32d7f311,  -18)
*	number (   -56,b877aa32,36a4b449,9befeb9,fad487c2,ff8defdb,  -17)
*/
tm16:	number (   -53,e69594be,c44de15b,4c2ebe68,7989a9b3,bf716bd5,  -16)
	number (   -49,901d7cf7,3ab0acd9,f9d3701,4bf60a10,57a6e369,  -15)
	number (   -46,b424dc35,95cd80f,538484c1,9ef38c94,6d909c46,  -14)
	number (   -43,e12e1342,4bb40e13,2865a5f2,6b06fb9,88f4c35a,  -13)
	number (   -39,8cbccc09,6f5088cb,f93f87b7,442e45d3,f598fa1c,  -12)
	number (   -36,afebff0b,cb24aafe,f78f69a5,1539d748,f2ff38a8,  -11)
	number (   -33,dbe6fece,bdedd5be,b573440e,5a884d1b,2fbf06d5,  -10)
	number (   -29,89705f41,36b4a597,31680a88,f8953030,fdd76447,   -9)
	number (   -26,abcc7711,8461cefc,fdc20d2b,36ba7c3d,3d4d3d5c,   -8)
	number (   -23,d6bf94d5,e57a42bc,3d329076,4691b4c,8ca08cb8,   -7)
	number (   -19,8637bd05,af6c69b5,a63f9a49,c2c1b10f,d7e457f7,   -6)
	number (   -16,a7c5ac47,1b478423,fcf80dc,33721d53,cddd6df6,   -5)
	number (   -13,d1b71758,e219652b,d3c36113,404ea4a8,c154c978,   -4)
	number (    -9,83126e97,8d4fdf3b,645a1cac,83126e9,78d4fdee,   -3)
	number (    -6,a3d70a3d,70a3d70a,3d70a3d7,a3d70a3,d70a3d6c,   -2)
	number (    -3,cccccccc,cccccccc,cccccccc,cccccccc,cccccccc,   -1)
 #ots$$a_cvt_tab::
t0:	number (     1,80000000,00000000,00000000,00000000,00000000,    0)
t1:	number (     4,a0000000,00000000,00000000,00000000,00000000,    1)
	number (     7,c8000000,00000000,00000000,00000000,00000000,    2)
	number (    10,fa000000,00000000,00000000,00000000,00000000,    3)
	number (    14,9c400000,00000000,00000000,00000000,00000000,    4)
	number (    17,c3500000,00000000,00000000,00000000,00000000,    5)
	number (    20,f4240000,00000000,00000000,00000000,00000000,    6)
	number (    24,98968000,00000000,00000000,00000000,00000000,    7)
	number (    27,bebc2000,00000000,00000000,00000000,00000000,    8)
	number (    30,ee6b2800,00000000,00000000,00000000,00000000,    9)
	number (    34,9502f900,00000000,00000000,00000000,00000000,   10)
	number (    37,ba43b740,00000000,00000000,00000000,00000000,   11)
	number (    40,e8d4a510,00000000,00000000,00000000,00000000,   12)
	number (    44,9184e72a,00000000,00000000,00000000,00000000,   13)
	number (    47,b5e620f4,80000000,00000000,00000000,00000000,   14)
	number (    50,e35fa931,a0000000,00000000,00000000,00000000,   15)
t16:	number (    54,8e1bc9bf,4000000,00000000,00000000,00000000,   16)
/*	number (    57,b1a2bc2e,c5000000,00000000,00000000,00000000,   17)
*	number (    60,de0b6b3a,76400000,00000000,00000000,00000000,   18)
*	number (    64,8ac72304,89e80000,00000000,00000000,00000000,   19)
*	number (    67,ad78ebc5,ac620000,00000000,00000000,00000000,   20)
*	number (    70,d8d726b7,177a8000,00000000,00000000,00000000,   21)
*	number (    74,87867832,6eac9000,00000000,00000000,00000000,   22)
*	number (    77,a968163f, a57b400,00000000,00000000,00000000,   23)
*	number (    80,d3c21bce,cceda100,00000000,00000000,00000000,   24)
*	number (    84,84595161,401484a0,00000000,00000000,00000000,   25)
*	number (    87,a56fa5b9,9019a5c8,00000000,00000000,00000000,   26)
*	number (    90,cecb8f27,f4200f3a,00000000,00000000,00000000,   27)
*	number (    94,813f3978,f8940984,40000000,00000000,00000000,   28)
*	number (    97,a18f07d7,36b90be5,50000000,00000000,00000000,   29)
*	number (   100,c9f2c9cd, 4674ede,a4000000,00000000,00000000,   30)
*	number (   103,fc6f7c40,45812296,4d000000,00000000,00000000,   31)
* t32:
*/
	number (   107,9dc5ada8,2b70b59d,f0200000,00000000,00000000,   32)
	number (   213,c2781f49,ffcfa6d5,3cbf6b71,c76b25fb,50f80800,   64)
	number (   426,93ba47c9,80e98cdf,c66f336c,36b10137,234f3fc,  128)
	number (   851,aa7eebfb,9df9de8d,ddbb901b,98feeab7,851e4cbb,  256)
	number (  1701,e319a0ae,a60e91c6,cc655c54,bc5058f8,9c658389,  512)
	number (  3402,c9767586,81750c17,650d3d28,f18b50ce,526b9865, 1024)
	number (  6804,9e8b3b5d,c53d5de4,a74d28ce,329ace52,6a31978c, 2048)
	number ( 13607,c4605202,8a20979a,c94c153f,804a4a92,65761f39, 4096)
 #	number ( 27214,96a3a1d1,7faf211a,c7c2892,305f4e12,72b205f, 8192)
 #                   0,ffffffff,ffffffff,ffffffff,ffffffff,ffff5eb4	; 1.0 if exact)

		# this table contains the byte indicies for the
		# multiple precision multiply cross products.
		# the 1st and 2nd entries on each line are the indicies
		# for the multiplicand and the multiplier. the third
		#  entry is the product index.

bytab:	.byte	0,3,0
	.byte	3,0,0
	.byte	2,1,0
	.byte	1,2,0
	.byte	1,3,1
	.byte	3,1,1
	.byte	2,2,1
	.byte	2,3,2
	.byte	3,2,2
	.byte	3,3,3
	.byte	-1		# end flag


# else

	clrb sign
	movd (ap)+,r5
	jeql fzero
	bgtr fpos
	mnegd r5,r5
	incb sign
fpos:
	extzv $7,$8,r5,r2		# exponent of 2
	movab -0200(r2),r2		# unbias
	mull2 $59,r2			# 59/196: 3rd convergent continued frac of log10(2)
	jlss eneg
	movab 196(r2),r2
eneg:
	movab -98(r2),r2
	divl2 $196,r2
	bsbw expten
	cmpd r0,r5
	bgtr ceil
	incl r2
ceil:	movl r2,exp
	mnegl r2,r2
	cmpl r2,$29			# 10^(29+9) is all we can handle
	bleq getman
	muld2 ten16,r5
	subl2 $16,r2
getman:	addl2 $9,r2			# -ceil(log10(x)) + 9
	jsb expten
	emodd r0,r4,r5,r0,r5		# (r0+r4)*r5; r0=int, r5=frac
fz1:	cvtlp r0,$9,16(sp)		# leading 9 digits
	ashp $8,$9,16(sp),$0,$17,4(sp)	# as top 9 of 17
	emodd ten8,$0,r5,r0,r5
	cvtlp r0,$8,16(sp)		# trailing 8 digits
		# if precision >= 17, must round here
	movl ndigit,r7			# so figure out what precision is
	pushab scien
	cmpl (sp)+,(sp)
	jleq gm1			# who called us?
	addl2 exp,r7			# float; adjust for exponent
gm1:	cmpl r7,$17
	jlss gm2
	cmpd r5,$0d0.5			# must round here; check fraction
	jlss gm2
	bisb2 $0x10,8+4(sp)		# increment l.s. digit
gm2:		# end of "round here" code
	addp4 $8,16(sp),$17,4(sp)	# combine leading and trailing
	bisb2 sign,12(sp)		# and insert sign
	rsb
fzero:	clrl r0
	movl $1,exp		# 0.000e+00 and 0.000 rather than 0.000e-01 and .000
	jbr fz1

	.align 2
lsb: .long 0x00010000		# lsb in the crazy floating-point format

	# return 10^r2 as a double float in r0||r1 and 8 extra bits of precision in r4
	# preserve r2, r5||r6
expten:
	movd $0d1.0,r0			# begin computing 10^exp10
	clrl r4				# bit counter
	movad ten1,r3			# table address
	tstl r2
	bgeq e10lp
	mnegl r2,r2			# get absolute value
	jbss $6,r2,e10lp		# flag as negative
e10lp:	jbc r4,r2,el1			# want this power?
	muld2 (r3),r0			# yes
el1:	addl2 $8,r3			# advance to next power
	aobleq $5,r4,e10lp		# through 10^32
	jbcc $6,r2,el2			# correct for negative exponent
	divd3 r0,$0d1.0,r0		# by taking reciprocal
	cmpl $28,r2
	jneq enm28
	addl2 lsb,r1			# 10**-28 needs lsb incremented
enm28:	mnegl r2,r2			# original exponent of 10
el2:	addl3 $5*8,r2,r3		# negative bit positions are illegal?
	jbc r3,xlsbh-5,eoklsb
	subl2 lsb,r1			# lsb was too high
eoklsb:
	movzbl xprec[r2],r4		# 8 extra bits
	rsb

	# powers of ten
	.align	2
ten1:	.word	0x4220,0,0,0
ten2:	.word	0x43c8,0,0,0
ten4:	.word	0x471c,0x4000,0,0
ten8:	.word	0x4dbe,0xbc20,0,0
ten16:	.word	0x5b0e,0x1bc9,0xbf04,0
ten32:	.word	0x759d,0xc5ad,0xa82b,0x70b6

	# whether lsb is too high or not
	.byte 1:0,1:0,1:0,1:0,1:1,1:0,1:1,1:0	# -40 thru -33
	.byte 1:0,1:1,1:0,1:0,1:0,1:0,1:1,1:0	# -32 thru -25
	.byte 1:0,1:0,1:1,1:1,1:1,1:1,1:0,1:0	# -24 thru -17
	.byte 1:0,1:1,1:0,1:0,1:1,1:1,1:1,1:1	# -16 thru -9
	.byte 1:1,1:1,1:1,1:0,1:0,1:0,1:0,1:1	# -8  thru -1
xlsbh:
	.byte 1:0,1:0,1:0,1:0,1:0,1:0,1:0,1:0	# 0 thru 7
	.byte 1:0,1:0,1:0,1:0,1:0,1:0,1:0,1:0	# 8 thru 15
	.byte 1:0,1:0,1:0,1:0,1:0,1:0,1:0,1:0	# 16 thru 23
	.byte 1:0,1:1,1:1,1:0,1:1,1:1,1:1,1:1	# 24 thru 31
	.byte 1:1,1:1,1:1,1:1,1:1,1:1,1:1    	# 32 thru 38

	# bytes of extra precision
	.byte           0x56,0x76,0xd3,0x88,0xb5,0x62	# -38 thru -33
	.byte 0xba,0xf5,0x32,0x3e,0x0e,0x48,0xdb,0x51	# -32 thru -25
	.byte 0x53,0x27,0xb1,0xef,0xeb,0xa5,0x07,0x49	# -24 thru -17
	.byte 0x5b,0xd9,0x0f,0x13,0xcd,0xff,0xbf,0x97	# -16 thru -9
	.byte 0xfd,0xbc,0xb6,0x23,0x2c,0x3b,0x0a,0xcd	# -8  thru -1
xprec:
	.byte 0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00	# 0  thru 7
	.byte 0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00	# 8  thru 15
	.byte 0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00	# 16 thru 23
	.byte 0x00,0xa0,0xc8,0x3a,0x84,0xe4,0xdc,0x92	# 24 thru 31
	.byte 0x9b,0x00,0xc0,0x58,0xae,0x18,0xef     	# 32 thru 38
# endif