satan.s

VXWORKS源代码

	.text

	.globl	__x_satand
__x_satand:
|--ENTRY POINT FOR ATAN(X) FOR DENORMALIZED ARGUMENT

	jra 		__x_t_extdnrm

	.globl	__x_satan
__x_satan:
/* |--ENTRY POINT FOR ATAN(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S */

	fmovex		A0@,fp0	|...lOAD INPUT

	movel		A0@,d0
	movew		A0@(4),d0
	fmovex		fp0,a6@(X)
	andil		#0x7FFFFFFF,d0

	cmpil		#0x3FFB8000,d0		|...|X| >= 1/16?
	jge 		ATANOK1
	jra 		ATANSM

ATANOK1:
	cmpil		#0x4002FFFF,d0		|...|X| < 16 ?
	jle 		ATANMAIN
	jra 		ATANBIG


|--THE MOST LIKELY CASE, |X| IN [1/16, 16). WE USE TABLE TECHNIQUE
|--THE IDEA IS ATAN(X) = ATAN(F) + ATAN( [X-F] / [1+XF] ).
|--SO IF F IS CHOSEN TO BE CLOSE TO X AND ATAN(F) IS STORED IN
|--A TABLE, ALL WE NEED IS TO APPROXIMATE ATAN(U) WHERE
|--U = (X-F)/(1+XF) IS SMALL (REMEMBER F IS CLOSE TO X). IT IS
|--TRUE THAT A DIVIDE IS NOW NEEDED, BUT THE APPROXIMATION FOR
|--ATAN(U) IS A VERY SHORT POLYNOMIAL AND THE INDEXING TO
|--FETCH F AND SAVING OF REGISTERS CAN BE ALL HIDED UNDER THE
|--DIVIDE. IN THE END THIS METHOD IS MUCH FASTER THAN A TRADITIONAL
|--ONE. NOTE ALSO THAT THE TRADITIONAL SCHEME THAT APPROXIMATE
|--ATAN(X) DIRECTLY WILL NEED TO USE A RATIONAL APPROXIMATION
|--(DIVISION NEEDED) ANYWAY BECAUSE A POLYNOMIAL APPROXIMATION
|--WILL INVOLVE A VERY LONG POLYNOMIAL.

|--NOW WE SEE X AS +-2^K * 1.BBBBBBB....B <- 1. + 63 BITS
|--WE CHOSE F TO BE +-2^K * 1.BBBB1
|--THAT IS IT MATCHES THE EXPONENT AND FIRST 5 BITS OF X, THE
|--SIXTH BITS IS SET TO BE 1. SINCE K = -4, -3, |..., 3, THERE
/* |--ARE ONLY 8 TIMES 16 = 2^7 = 128 |F|'S. SINCE ATAN(-|F|) IS */
|-- -ATAN(|F|), WE NEED TO STORE ONLY ATAN(|F|).

ATANMAIN:

	movew		#0x0000,a6@(XDCARE)	|...CLEAN UP X JUST IN CASE
	andil		#0xF8000000,a6@(XFRAC)	|...FIRST 5 BITS
	oril		#0x04000000,a6@(XFRAC)	|...SET 6-TH BIT TO 1
	movel		#0x00000000,a6@(XFRACLO)	|...lOCATION OF X IS NOW F

	fmovex		fp0,fp1			|...FP1 IS X
	fmulx		a6@(X),fp1		|...FP1 IS X*F, NOTE THAT X*F > 0
	fsubx		a6@(X),fp0		|...FP0 IS X-F
	.long 0xf23c44a2,0x3f800000		/*  fadds  &0x3F800000,fp1 */
	fdivx		fp1,fp0			|...FP0 IS U = (X-F)/(1+X*F)

|--WHILE THE DIVISION IS TAKING ITS TIME, WE FETCH ATAN(|F|)
|--CREATE ATAN(F) AND STORE IT IN ATANF, AND
|--SAVE REGISTERS FP2.

	movel		d2,a7@-	|...SAVE d2 TEMPORARILY
	movel		d0,d2		|...THE EXPO AND 16 BITS OF X
	andil		#0x00007800,d0	/* |...4 VARYING BITS OF F'S FRACTION */
	andil		#0x7FFF0000,d2	|...EXPONENT OF F
	subil		#0x3FFB0000,d2	|...K+4
	asrl		#1,d2
	addl		d2,d0		|...THE 7 BITS IDENTIFYING F
	asrl		#7,d0		|...INDEX INTO TBL OF ATAN(|F|)
	lea		ATANTBL,a1
	addal		d0,a1		|...ADDRESS OF ATAN(|F|)
	movel		a1@+,a6@(ATANF)
	movel		a1@+,a6@(ATANFHI)
	movel		a1@+,a6@(ATANFLO)	|...ATANF IS NOW ATAN(|F|)
	movel		a6@(X),d0		|...lOAD SIGN AND EXPO. AGAIN
	andil		#0x80000000,d0	|...SIGN(F)
	orl		d0,a6@(ATANF)	|...ATANF IS NOW SIGN(F)*ATAN(|F|)
	movel		a7@+,d2		|...RESTORE d2

/* |--THAT'S ALL I HAVE TO DO FOR NOW, */
|--BUT ALAS, THE DIVIDE IS STILL CRANKING!

|--U IN FP0, WE ARE NOW READY TO COMPUTE ATAN(U) AS
|--U + A1*U*V*(A2 + V*(A3 + V)), V = U*U
|--THE POLYNOMIAL MAY LOOK STRANGE, BUT IS NEVERTHELESS CORRECT.
|--THE NATURAL FORM IS U + U*V*(A1 + V*(A2 + V*A3))
|--WHAT WE HAVE HERE IS MERELY	A1 = A3, A2 = A1/A3, A3 = A2/A3.
|--THE REASON FOR THIS REARRANGEMENT IS TO MAKE THE INDEPENDENT
|--PARTS A1*U*V AND (A2 + |... STUFF) MORE LOAD-BALANCED


	fmovex		fp0,fp1
	fmulx		fp1,fp1
	fmoved		ATANA3,fp2
	faddx		fp1,fp2		|...A3+V
	fmulx		fp1,fp2		|...V*(A3+V)
	fmulx		fp0,fp1		|...U*V
	faddd		ATANA2,fp2	|...A2+V*(A3+V)
	fmuld		ATANA1,fp1	|...A1*U*V
	fmulx		fp2,fp1		|...A1*U*V*(A2+V*(A3+V))

	faddx		fp1,fp0		|...ATAN(U), fp1 RELEASED
	fmovel		d1,fpcr		| restore users exceptions
	faddx		a6@(ATANF),fp0	|...ATAN(X)
	jra 		__x_t_frcinx

ATANBORS:
|--|X| IS IN d0 IN COMPACT FORM. FP1, d0 SAVED.
|--FP0 IS X AND |X| <= 1/16 OR |X| >= 16.
	cmpil		#0x3FFF8000,d0
	jgt 		ATANBIG	|...I.E. |X| >= 16

ATANSM:
|--|X| <= 1/16
|--IF |X| < 2^(-40), RETURN X AS ANSWER. OTHERWISE, APPROXIMATE
|--ATAN(X) BY X + X*Y*(B1+Y*(B2+Y*(B3+Y*(B4+Y*(B5+Y*B6)))))
|--WHICH IS X + X*Y*( [B1+Z*(B3+Z*B5)] + [Y*(B2+Z*(B4+Z*B6)] )
|--WHERE Y = X*X, AND Z = Y*Y.

	cmpil		#0x3FD78000,d0
	jlt 		ATANTINY
|--COMPUTE POLYNOMIAL
	fmulx		fp0,fp0	|...FP0 IS Y = X*X


	movew		#0x0000,a6@(XDCARE)

	fmovex		fp0,fp1
	fmulx		fp1,fp1		|...FP1 IS Z = Y*Y

	fmoved		ATANB6,fp2
	fmoved		ATANB5,fp3

	fmulx		fp1,fp2		|...Z*B6
	fmulx		fp1,fp3		|...Z*B5

	faddd		ATANB4,fp2	|...B4+Z*B6
	faddd		ATANB3,fp3	|...B3+Z*B5

	fmulx		fp1,fp2		|...Z*(B4+Z*B6)
	fmulx		fp3,fp1		|...Z*(B3+Z*B5)

	faddd		ATANB2,fp2	|...B2+Z*(B4+Z*B6)
	faddd		ATANB1,fp1	|...B1+Z*(B3+Z*B5)

	fmulx		fp0,fp2		|...Y*(B2+Z*(B4+Z*B6))
	fmulx		a6@(X),fp0		|...X*Y

	faddx		fp2,fp1		|...[B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]


	fmulx		fp1,fp0	|...X*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))])

	fmovel		d1,fpcr		| restore users exceptions
	faddx		a6@(X),fp0

	jra 		__x_t_frcinx

ATANTINY:
|--|X| < 2^(-40), ATAN(X) = X
	movew		#0x0000,a6@(XDCARE)

	fmovel		d1,fpcr		| restore users exceptions
	fmovex		a6@(X),fp0	| last inst - possible exception set

	jra 		__x_t_frcinx

ATANBIG:
|--IF |X| > 2^(100), RETURN	SIGN(X)*(PI/2 - TINY). OTHERWISE,
|--RETURN SIGN(X)*PI/2 + ATAN(-1/X).
	cmpil		#0x40638000,d0
	jgt 		ATANHUGE

|--APPROXIMATE ATAN(-1/X) BY
/* |--X'+X'*Y*(C1+Y*(C2+Y*(C3+Y*(C4+Y*C5)))), X' = -1/X, Y = X'*X' */
|--THIS CAN BE RE-WRITTEN AS
/* |--X'+X'*Y*( [C1+Z*(C3+Z*C5)] + [Y*(C2+Z*C4)] ), Z = Y*Y. */

	.long 0xf23c4480,0xbf800000	/*  fmoves  &0xBF800000,fp1 */
	fdivx		fp0,fp1		|...FP1 IS -1/X


|--DIVIDE IS STILL CRANKING

	fmovex		fp1,fp0		/* |...FP0 IS X' */
	fmulx		fp0,fp0		/* |...FP0 IS Y = X'*X' */
	fmovex		fp1,a6@(X)		/* |...X IS REALLY X' */

	fmovex		fp0,fp1
	fmulx		fp1,fp1		|...FP1 IS Z = Y*Y

	fmoved		ATANC5,fp3
	fmoved		ATANC4,fp2

	fmulx		fp1,fp3		|...Z*C5
	fmulx		fp1,fp2		|...Z*B4

	faddd		ATANC3,fp3	|...C3+Z*C5
	faddd		ATANC2,fp2	|...C2+Z*C4

	fmulx		fp3,fp1		|...Z*(C3+Z*C5), fp3 RELEASED
	fmulx		fp0,fp2		|...Y*(C2+Z*C4)

	faddd		ATANC1,fp1	|...C1+Z*(C3+Z*C5)
	fmulx		a6@(X),fp0		/* |...X'*Y */

	faddx		fp2,fp1		|...[Y*(C2+Z*C4)]+[C1+Z*(C3+Z*C5)]


	fmulx		fp1,fp0		/* |...X'*Y*([B1+Z*(B3+Z*B5)] */
|					|...	+[Y*(B2+Z*(B4+Z*B6))])
	faddx		a6@(X),fp0

	fmovel		d1,fpcr		| restore users exceptions

	btst		#7,a0@
	jeq 		pos_big

neg_big:
	faddx		NPIBY2,fp0
	jra 		__x_t_frcinx

pos_big:
	faddx		PPIBY2,fp0
	jra 		__x_t_frcinx

ATANHUGE:
|--RETURN SIGN(X)*(PIBY2 - TINY) = SIGN(X)*PIBY2 - SIGN(X)*TINY
	btst		#7,a0@
	jeq 		pos_huge

neg_huge:
	fmovex		NPIBY2,fp0
	fmovel		d1,fpcr
	fsubx		NTINY,fp0
	jra 		__x_t_frcinx

pos_huge:
	fmovex		PPIBY2,fp0
	fmovel		d1,fpcr
	fsubx		PTINY,fp0
	jra 		__x_t_frcinx

|	end