satan.s

RTEMS (Real-Time Executive for Multiprocessor Systems) is a free open source real-time operating sys

	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