cos(0) debugged

utzoo!decvax!duke!chico!harpo!npois!houxi!ihnss!vax135!jfr utzoo!decvax!duke!chico!harpo!npois!houxi!ihnss!vax135!jfr
Tue Feb 16 15:22:50 AEST 1982


ecvt has nothing to do with what happened, since it is not
invoked at all in this case.  As dmr noted, the problem lies
in the misuse of floating-point arithmetic.  One look at the
three trailing zeroes on each coefficient in the sin.s file
should have been a dead give away.  I wrote doprnt.s that way
for this very reason.

Our version of atof.s (October 1979) is accurate to machine
precision; in fact, it carries 63 bits of precision until
stuffing the exponent.

There is an argument that doprnt's refusal to provide more
than 17 decimal digits is incorrect because some numbers
with common logarithms between 16 and 17 can be represented
in fewer than 56 bits.  The following diff shows how to get
18 digits from doprnt, and in fact this would have prevented
the cos(0) anomaly.  However, you aren't always entitled to
18 digits, and more digits only mask the true problem
of adopting the wrong strategy.

		John F. Reiser
--------------------------------------------------------------------
463,464c462,463
< fnarro:	subl3 $17,r7,r0			# where to round
< 	ashp r0,$17,(sp),$5,r7,16(sp)	# do it
---
> fnarro:	subl3 $18,r7,r0			# where to round
> 	ashp r0,$18,(sp),$5,r7,16(sp)	# do it
576,577c575,576
< snarro:	subl3 $17,r7,r0		# rounding position
< 	ashp r0,$17,(sp),$5,r7,16(sp) # shift and round
---
> snarro:	subl3 $18,r7,r0		# rounding position
> 	ashp r0,$18,(sp),$5,r7,16(sp) # shift and round
646c645
< 	# convert double-floating at (ap) to 17-digit packed at (sp),
---
> 	# convert double-floating at (ap) to 18-digit packed at (sp),
678,681c677,680
< 	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
---
> 	ashp $9,$9,16(sp),$0,$18,4(sp)	# as top 9 of 18
> 	emodd ten9,$0,r5,r0,r5
> 	cvtlp r0,$9,16(sp)		# trailing 9 digits
> 		# if precision >= 18, must round here
687c686
< gm1:	cmpl r7,$17
---
> gm1:	cmpl r7,$18
691c690
< 	bisb2 $0x10,8+4(sp)		# increment l.s. digit
---
> 	bisb2 $0x10,9+4(sp)		# increment l.s. digit
693,694c692,693
< 	addp4 $8,16(sp),$17,4(sp)	# combine leading and trailing
< 	bisb2 sign,12(sp)		# and insert sign
---
> 	addp4 $9,16(sp),$18,4(sp)	# combine leading and trailing
> 	bisb2 sign,9+4(sp)		# and insert sign
700c699,700
< 	.align 2
---
> 	.align 3
> ten9:	.word 0x4f6e,0x6b28,0,0



More information about the Comp.bugs.4bsd.ucb-fixes mailing list