satan.s

讲述linux的初始化过程

	fmovex		(%a0),%fp0	| ...LOAD INPUT

	movel		(%a0),%d0
	movew		4(%a0),%d0
	fmovex		%fp0,X(%a6)
	andil		#0x7FFFFFFF,%d0

	cmpil		#0x3FFB8000,%d0		| ...|X| >= 1/16?
	bges		ATANOK1
	bra		ATANSM

ATANOK1:
	cmpil		#0x4002FFFF,%d0		| ...|X| < 16 ?
	bles		ATANMAIN
	bra		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,XDCARE(%a6)	| ...CLEAN UP X JUST IN CASE
	andil		#0xF8000000,XFRAC(%a6)	| ...FIRST 5 BITS
	oril		#0x04000000,XFRAC(%a6)	| ...SET 6-TH BIT TO 1
	movel		#0x00000000,XFRACLO(%a6)	| ...LOCATION OF X IS NOW F

	fmovex		%fp0,%fp1			| ...FP1 IS X
	fmulx		X(%a6),%fp1		| ...FP1 IS X*F, NOTE THAT X*F > 0
	fsubx		X(%a6),%fp0		| ...FP0 IS X-F
	fadds		#0x3F800000,%fp1		| ...FP1 IS 1 + X*F
	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)+,ATANF(%a6)
	movel		(%a1)+,ATANFHI(%a6)
	movel		(%a1)+,ATANFLO(%a6)	| ...ATANF IS NOW ATAN(|F|)
	movel		X(%a6),%d0		| ...LOAD SIGN AND EXPO. AGAIN
	andil		#0x80000000,%d0	| ...SIGN(F)
	orl		%d0,ATANF(%a6)	| ...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		ATANF(%a6),%fp0	| ...ATAN(X)
	bra		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
	bgt		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
	blt		ATANTINY
|--COMPUTE POLYNOMIAL
	fmulx		%fp0,%fp0	| ...FP0 IS Y = X*X

	
	movew		#0x0000,XDCARE(%a6)

	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		X(%a6),%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		X(%a6),%fp0

	bra		t_frcinx

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

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

	bra		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
	bgt		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.

	fmoves		#0xBF800000,%fp1	| ...LOAD -1
	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,X(%a6)		| ...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		X(%a6),%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		X(%a6),%fp0

	fmovel		%d1,%FPCR		|restore users exceptions
	
	btstb		#7,(%a0)
	beqs		pos_big

neg_big:
	faddx		NPIBY2,%fp0
	bra		t_frcinx

pos_big:
	faddx		PPIBY2,%fp0
	bra		t_frcinx

ATANHUGE:
|--RETURN SIGN(X)*(PIBY2 - TINY) = SIGN(X)*PIBY2 - SIGN(X)*TINY
	btstb		#7,(%a0)
	beqs		pos_huge

neg_huge:
	fmovex		NPIBY2,%fp0
	fmovel		%d1,%fpcr
	fsubx		NTINY,%fp0
	bra		t_frcinx

pos_huge:
	fmovex		PPIBY2,%fp0
	fmovel		%d1,%fpcr
	fsubx		PTINY,%fp0
	bra		t_frcinx
	
	|end