argred.s

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

#ifndef lint
#static char	*sccsid = "@(#)argred.s	4.1	(ULTRIX)	7/17/90";
#endif lint

/************************************************************************
 *									*
 *			Copyright (c) 1986 by				*
 *		Digital Equipment Corporation, Maynard, MA		*
 *			All rights reserved.				*
 *									*
 *   This software is furnished under a license and may be used and	*
 *   copied  only  in accordance with the terms of such license and	*
 *   with the  inclusion  of  the  above  copyright  notice.   This	*
 *   software  or  any  other copies thereof may not be provided or	*
 *   otherwise made available to any other person.  No title to and	*
 *   ownership of the software is hereby transferred.			*
 *									*
 *   This software is  derived  from  software  received  from  the	*
 *   University    of   California,   Berkeley,   and   from   Bell	*
 *   Laboratories.  Use, duplication, or disclosure is  subject  to	*
 *   restrictions  under  license  agreements  with  University  of	*
 *   California and with AT&T.						*
 *									*
 *   The information in this software is subject to change  without	*
 *   notice  and should not be construed as a commitment by Digital	*
 *   Equipment Corporation.						*
 *									*
 *   Digital assumes no responsibility for the use  or  reliability	*
 *   of its software on equipment which is not supplied by Digital.	*
 *									*
 ************************************************************************/

/************************************************************************
*
*			Modification History
*
*	David Metsky, 12/16/86
* 001	Adapted from argred.s 1.1 (Berkeley) 
*
*************************************************************************/

# @(#)argred.s	1.1 (Berkeley) 8/21/85

#  libm$argred implements Bob Corbett's argument reduction and
#  libm$sincos implements Peter Tang's double precision sin/cos.
#  
#  Note: The two entry points libm$argred and libm$sincos are meant
#        to be used only by _sin, _cos and _tan.
#
# method: true range reduction to [-pi/4,pi/4], P. Tang  &  B. Corbett
# S. McDonald, April 4,  1985
#
	.globl	libm$argred
	.globl	libm$sincos
	.text
	.align	1

libm$argred:
#
#  Compare the argument with the largest possible that can
#  be reduced by table lookup.  r3 := |x|  will be used in  table_lookup .
#
	movd	r0,r3
	bgeq	abs1
	mnegd	r3,r3
abs1:
	cmpd	r3,$0d+4.55530934770520019583e+01
	blss	small_arg
	jsb	trigred
	rsb
small_arg:
	jsb	table_lookup
	rsb
#
#  At this point,
#	   r0  contains the quadrant number, 0, 1, 2, or 3;
#	r2/r1  contains the reduced argument as a D-format number;
#  	   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:
#
#  Compensate for a cosine entry by adding one to the quadrant number.
#
	addl2	r4,r0
#
#  Polyd clobbers  r5-r0 ;  save  X  in  r7/r6 .
#  This can be avoided by rewriting  trigred .
#
	movd	r1,r6
#
#  Likewise, save  alpha  in  r8 .
#  This can be avoided by rewriting  trigred .
#
	movf	r3,r8
#
#  Odd or even quadrant?  cosine if odd, sine otherwise.
#  Save  floor(quadrant/2) in  r9  ; it determines the final sign.
#
	rotl	$-1,r0,r9
	blss	cosine
sine:
	muld2	r1,r1		# Xsq = X * X
	polyd	r1,$7,sin_coef	# Q = P(Xsq) , of deg 7
	mulf3	$0f3.0,r8,r4	# beta = 3 * alpha
	mulf2	r0,r4		# beta = Q * beta
	addf2	r8,r4		# beta = alpha + beta
	muld2	r6,r0		# S(X) = X * Q
#	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
cosine:
	muld2	r6,r6		# Xsq = X * X
	beql	zero_arg
	mulf2	r1,r8		# beta = X * alpha
	polyd	r6,$7,cos_coef	# Q = P'(Xsq) , of deg 7
	subd3	r0,r8,r0	# beta = beta - Q
	subw2	$0x80,r6	# Xsq = Xsq / 2
	addd2	r0,r6		# Xsq = Xsq + beta
zero_arg:
	subd3	r6,$0d1.0,r0	# C(X) = 1 - Xsq
done:
	blbc	r9,even
	mnegd	r0,r0
even:
	rsb

.data
.align	2

sin_coef:
	.double	0d-7.53080332264191085773e-13	# s7 = 2^-29 -1.a7f2504ffc49f8..
	.double	0d+1.60573519267703489121e-10	# s6 = 2^-21  1.611adaede473c8..
	.double	0d-2.50520965150706067211e-08	# s5 = 2^-1a -1.ae644921ed8382..
	.double	0d+2.75573191800593885716e-06	# s4 = 2^-13  1.71de3a4b884278..
	.double	0d-1.98412698411850507950e-04	# s3 = 2^-0d -1.a01a01a0125e7d..
	.double	0d+8.33333333333325688985e-03	# s2 = 2^-07  1.11111111110e50
	.double	0d-1.66666666666666664354e-01	# s1 = 2^-03 -1.55555555555554
	.double	0d+0.00000000000000000000e+00	# s0 = 0

cos_coef:
	.double	0d-1.13006966202629430300e-11	# s7 = 2^-25 -1.8D9BA04D1374BE..
	.double	0d+2.08746646574796004700e-09	# s6 = 2^-1D  1.1EE632650350BA..
	.double	0d-2.75573073031284417300e-07	# s5 = 2^-16 -1.27E4F31411719E..
	.double	0d+2.48015872682668025200e-05	# s4 = 2^-10  1.A01A0196B902E8..
	.double	0d-1.38888888888464709200e-03	# s3 = 2^-0A -1.6C16C16C11FACE..
	.double	0d+4.16666666666664761400e-02	# s2 = 2^-05  1.5555555555539E
	.double	0d+0.00000000000000000000e+00	# s1 = 0
	.double	0d+0.00000000000000000000e+00	# s0 = 0

#
#  Multiples of  pi/2  expressed as the sum of three doubles,
#
#  trailing:	n * pi/2 ,  n = 0, 1, 2, ..., 29
#			trailing[n] ,
#
#  middle:	n * pi/2 ,  n = 0, 1, 2, ..., 29
#			middle[n]   ,
#
#  leading:	n * pi/2 ,  n = 0, 1, 2, ..., 29
#			leading[n]  ,
#
#	where
#		leading[n]  := (n * pi/2)  rounded,
#		middle[n]   := (n * pi/2  -  leading[n])  rounded,
#		trailing[n] := (( n * pi/2 - leading[n]) - middle[n])  rounded .

trailing:
	.double	0d+0.00000000000000000000e+00	#  0 * pi/2  trailing
	.double	0d+4.33590506506189049611e-35	#  1 * pi/2  trailing
	.double	0d+8.67181013012378099223e-35	#  2 * pi/2  trailing
	.double	0d+1.30077151951856714215e-34	#  3 * pi/2  trailing
	.double	0d+1.73436202602475619845e-34	#  4 * pi/2  trailing
	.double	0d-1.68390735624352669192e-34	#  5 * pi/2  trailing
	.double	0d+2.60154303903713428430e-34	#  6 * pi/2  trailing
	.double	0d-8.16726343231148352150e-35	#  7 * pi/2  trailing
	.double	0d+3.46872405204951239689e-34	#  8 * pi/2  trailing
	.double	0d+3.90231455855570147991e-34	#  9 * pi/2  trailing
	.double	0d-3.36781471248705338384e-34	# 10 * pi/2  trailing
	.double	0d-1.06379439835298071785e-33	# 11 * pi/2  trailing
	.double	0d+5.20308607807426856861e-34	# 12 * pi/2  trailing
	.double	0d+5.63667658458045770509e-34	# 13 * pi/2  trailing
	.double	0d-1.63345268646229670430e-34	# 14 * pi/2  trailing
	.double	0d-1.19986217995610764801e-34	# 15 * pi/2  trailing
	.double	0d+6.93744810409902479378e-34	# 16 * pi/2  trailing
	.double	0d-8.03640094449267300110e-34	# 17 * pi/2  trailing
	.double	0d+7.80462911711140295982e-34	# 18 * pi/2  trailing
	.double	0d-7.16921993148029483506e-34	# 19 * pi/2  trailing
	.double	0d-6.73562942497410676769e-34	# 20 * pi/2  trailing
	.double	0d-6.30203891846791677593e-34	# 21 * pi/2  trailing
	.double	0d-2.12758879670596143570e-33	# 22 * pi/2  trailing
	.double	0d+2.53800212047402350390e-33	# 23 * pi/2  trailing
	.double	0d+1.04061721561485371372e-33	# 24 * pi/2  trailing
	.double	0d+6.11729905311472319056e-32	# 25 * pi/2  trailing
	.double	0d+1.12733531691609154102e-33	# 26 * pi/2  trailing
	.double	0d-3.70049587943078297272e-34	# 27 * pi/2  trailing
	.double	0d-3.26690537292459340860e-34	# 28 * pi/2  trailing
	.double	0d-1.14812616507957271361e-34	# 29 * pi/2  trailing

middle:
	.double	0d+0.00000000000000000000e+00	#  0 * pi/2  middle
	.double	0d+5.72118872610983179676e-18	#  1 * pi/2  middle
	.double	0d+1.14423774522196635935e-17	#  2 * pi/2  middle
	.double	0d-3.83475850529283316309e-17	#  3 * pi/2  middle
	.double	0d+2.28847549044393271871e-17	#  4 * pi/2  middle
	.double	0d-2.69052076007086676522e-17	#  5 * pi/2  middle
	.double	0d-7.66951701058566632618e-17	#  6 * pi/2  middle
	.double	0d-1.54628301484890040587e-17	#  7 * pi/2  middle
	.double	0d+4.57695098088786543741e-17	#  8 * pi/2  middle
	.double	0d+1.07001849766246313192e-16	#  9 * pi/2  middle
	.double	0d-5.38104152014173353044e-17	# 10 * pi/2  middle
	.double	0d-2.14622680169080983801e-16	# 11 * pi/2  middle
	.double	0d-1.53390340211713326524e-16	# 12 * pi/2  middle
	.double	0d-9.21580002543456677056e-17	# 13 * pi/2  middle
	.double	0d-3.09256602969780081173e-17	# 14 * pi/2  middle
	.double	0d+3.03066796603896507006e-17	# 15 * pi/2  middle
	.double	0d+9.15390196177573087482e-17	# 16 * pi/2  middle
	.double	0d+1.52771359575124969107e-16	# 17 * pi/2  middle
	.double	0d+2.14003699532492626384e-16	# 18 * pi/2  middle
	.double	0d-1.68853170360202329427e-16	# 19 * pi/2  middle
	.double	0d-1.07620830402834670609e-16	# 20 * pi/2  middle
	.double	0d+3.97700719404595604379e-16	# 21 * pi/2  middle
	.double	0d-4.29245360338161967602e-16	# 22 * pi/2  middle
	.double	0d-3.68013020380794313406e-16	# 23 * pi/2  middle
	.double	0d-3.06780680423426653047e-16	# 24 * pi/2  middle
	.double	0d-2.45548340466059054318e-16	# 25 * pi/2  middle
	.double	0d-1.84316000508691335411e-16	# 26 * pi/2  middle
	.double	0d-1.23083660551323675053e-16	# 27 * pi/2  middle
	.double	0d-6.18513205939560162346e-17	# 28 * pi/2  middle
	.double	0d-6.18980636588357585202e-19	# 29 * pi/2  middle

leading:
	.double	0d+0.00000000000000000000e+00	#  0 * pi/2  leading
	.double	0d+1.57079632679489661351e+00	#  1 * pi/2  leading
	.double	0d+3.14159265358979322702e+00	#  2 * pi/2  leading
	.double	0d+4.71238898038468989604e+00	#  3 * pi/2  leading
	.double	0d+6.28318530717958645404e+00	#  4 * pi/2  leading
	.double	0d+7.85398163397448312306e+00	#  5 * pi/2  leading
	.double	0d+9.42477796076937979208e+00	#  6 * pi/2  leading
	.double	0d+1.09955742875642763501e+01	#  7 * pi/2  leading
	.double	0d+1.25663706143591729081e+01	#  8 * pi/2  leading
	.double	0d+1.41371669411540694661e+01	#  9 * pi/2  leading
	.double	0d+1.57079632679489662461e+01	# 10 * pi/2  leading
	.double	0d+1.72787595947438630262e+01	# 11 * pi/2  leading
	.double	0d+1.88495559215387595842e+01	# 12 * pi/2  leading
	.double	0d+2.04203522483336561422e+01	# 13 * pi/2  leading
	.double	0d+2.19911485751285527002e+01	# 14 * pi/2  leading
	.double	0d+2.35619449019234492582e+01	# 15 * pi/2  leading
	.double	0d+2.51327412287183458162e+01	# 16 * pi/2  leading
	.double	0d+2.67035375555132423742e+01	# 17 * pi/2  leading
	.double	0d+2.82743338823081389322e+01	# 18 * pi/2  leading
	.double	0d+2.98451302091030359342e+01	# 19 * pi/2  leading
	.double	0d+3.14159265358979324922e+01	# 20 * pi/2  leading
	.double	0d+3.29867228626928286062e+01	# 21 * pi/2  leading
	.double	0d+3.45575191894877260523e+01	# 22 * pi/2  leading
	.double	0d+3.61283155162826226103e+01	# 23 * pi/2  leading
	.double	0d+3.76991118430775191683e+01	# 24 * pi/2  leading
	.double	0d+3.92699081698724157263e+01	# 25 * pi/2  leading
	.double	0d+4.08407044966673122843e+01	# 26 * pi/2  leading
	.double	0d+4.24115008234622088423e+01	# 27 * pi/2  leading
	.double	0d+4.39822971502571054003e+01	# 28 * pi/2  leading
	.double	0d+4.55530934770520019583e+01	# 29 * pi/2  leading

twoOverPi:
	.double	0d+6.36619772367581343076e-01
	.text
	.align	1

table_lookup:
	muld3	r3,twoOverPi,r0
	cvtrdl	r0,r0			# n = nearest int to ((2/pi)*|x|) rnded
	mull3	$8,r0,r5
	subd2	leading(r5),r3		# p = (|x| - leading n*pi/2) exactly
	subd3	middle(r5),r3,r1	# q = (p - middle  n*pi/2) rounded
	subd2	r1,r3			# r = (p - q)
	subd2	middle(r5),r3		# r =  r - middle  n*pi/2
	subd2	trailing(r5),r3		# r =  r - trailing n*pi/2  rounded
#
#  If the original argument was negative,
#  negate the reduce argument and
#  adjust the octant/quadrant number.
#
	tstw	4(ap)
	bgeq	abs2
	mnegf	r1,r1
	mnegf	r3,r3
#	subb3	r0,$8,r0	...used for  pi/4  reduction -S.McD
	subb3	r0,$4,r0
abs2:
#
#  Clear all unneeded octant/quadrant bits.
#
#	bicb2	$0xf8,r0	...used for  pi/4  reduction -S.McD
	bicb2	$0xfc,r0
	rsb
#
#						p.0
	.text
	.align	2
#
# Only 256 (actually 225) bits of 2/pi are needed for VAX double
# precision; this was determined by enumerating all the nearest
# machine integer multiples of pi/2 using continued fractions.
# (8a8d3673775b7ff7 required the most bits.)		-S.McD
#
	.long	0
	.long	0
	.long	0xaef1586d
	.long	0x9458eaf7
	.long	0x10e4107f
	.long	0xd8a5664f
	.long	0x4d377036
	.long	0x09d5f47d
	.long	0x91054a7f
	.long	0xbe60db93
bits2opi:
	.long	0x00000028
	.long	0
#
#  Note: wherever you see the word `octant', read `quadrant'.
#  Currently this code is set up for  pi/2  argument reduction.
#  By uncommenting/commenting the appropriate lines, it will
#  also serve as a  pi/4  argument reduction code.
#  

#						p.1
#  Trigred  preforms argument reduction
#  for the trigonometric functions.  It
#  takes one input argument, a D-format
#  number in  r1/r0 .  The magnitude of
#  the input argument must be greater
#  than or equal to  1/2 .  Trigred produces
#  three results:  the number of the octant
#  occupied by the argument, the reduced 
#  argument, and an extension of the
#  reduced argument.  The octant number is 
#  returned in  r0 .  The reduced argument 
#  is returned as a D-format number in 
#  r2/r1 .  An 8 bit extension of the 
#  reduced argument is returned as an 
#  F-format number in r3.
#						p.2
trigred:
#
#  Save the sign of the input argument.
#
	movw	r0,-(sp)
#
#  Extract the exponent field.
#
	extzv	$7,$7,r0,r2
#
#  Convert the fraction part of the input
#  argument into a quadword integer.
#
	bicw2	$0xff80,r0
	bisb2	$0x80,r0	# -S.McD
	rotl	$16,r0,r0
	rotl	$16,r1,r1
#
#  If  r1  is negative, add  1  to  r0 .  This
#  adjustment is made so that the two's
#  complement multiplications done later
#  will produce unsigned results.
#
	bgeq	posmid
	incl	r0
posmid:
#						p.3
#
#  Set  r3  to the address of the first quadword
#  used to obtain the needed portion of  2/pi .
#  The address is longword aligned to ensure
#  efficient access.
#
	ashl	$-3,r2,r3
	bicb2	$3,r3
	subl3	r3,$bits2opi,r3
#
#  Set  r2  to the size of the shift needed to 
#  obtain the correct portion of  2/pi .
#
	bicb2	$0xe0,r2
#						p.4
#
#  Move the needed  128  bits of  2/pi  into
#  r11 - r8 .  Adjust the numbers to allow
#  for unsigned multiplication.
#
	ashq	r2,(r3),r10

	subl2	$4,r3
	ashq	r2,(r3),r9
	bgeq	signoff1
	incl	r11