slogn.s

microwindows移植到S3C44B0的源码

     lsll	%d6,%d4
     movel	%d5,%d7			| ...a copy of D5
     lsll	%d6,%d5
     negl	%d6
     addil	#32,%d6
     lsrl	%d6,%d7
     orl	%d7,%d4			| ...(D3,D4,D5) normalized

     movel	%d3,X(%a6)
     movel	%d4,XFRAC(%a6)
     movel	%d5,XFRAC+4(%a6)
     negl	%d2
     movel	%d2,ADJK(%a6)
     fmovex	X(%a6),%fp0
     moveml	(%a7)+,%d2-%d7		| ...restore registers
     lea	X(%a6),%a0
     bras	LOGBGN			| ...begin regular log(X)


	.global	slogn
slogn:
|--ENTRY POINT FOR LOG(X) FOR X FINITE, NON-ZERO, NOT NAN'S

	fmovex		(%a0),%fp0	| ...LOAD INPUT
	movel		#0x00000000,ADJK(%a6)

LOGBGN:
|--FPCR SAVED AND CLEARED, INPUT IS 2^(ADJK)*FP0, FP0 CONTAINS
|--A FINITE, NON-ZERO, NORMALIZED NUMBER.

	movel	(%a0),%d0
	movew	4(%a0),%d0

	movel	(%a0),X(%a6)
	movel	4(%a0),X+4(%a6)
	movel	8(%a0),X+8(%a6)

	cmpil	#0,%d0		| ...CHECK IF X IS NEGATIVE
	blt	LOGNEG		| ...LOG OF NEGATIVE ARGUMENT IS INVALID
	cmp2l	BOUNDS1,%d0	| ...X IS POSITIVE, CHECK IF X IS NEAR 1
	bcc	LOGNEAR1	| ...BOUNDS IS ROUGHLY [15/16, 17/16]

LOGMAIN:
|--THIS SHOULD BE THE USUAL CASE, X NOT VERY CLOSE TO 1

|--X = 2^(K) * Y, 1 <= Y < 2. THUS, Y = 1.XXXXXXXX....XX IN BINARY.
|--WE DEFINE F = 1.XXXXXX1, I.E. FIRST 7 BITS OF Y AND ATTACH A 1.
|--THE IDEA IS THAT LOG(X) = K*LOG2 + LOG(Y)
|--			 = K*LOG2 + LOG(F) + LOG(1 + (Y-F)/F).
|--NOTE THAT U = (Y-F)/F IS VERY SMALL AND THUS APPROXIMATING
|--LOG(1+U) CAN BE VERY EFFICIENT.
|--ALSO NOTE THAT THE VALUE 1/F IS STORED IN A TABLE SO THAT NO
|--DIVISION IS NEEDED TO CALCULATE (Y-F)/F. 

|--GET K, Y, F, AND ADDRESS OF 1/F.
	asrl	#8,%d0
	asrl	#8,%d0		| ...SHIFTED 16 BITS, BIASED EXPO. OF X
	subil	#0x3FFF,%d0 	| ...THIS IS K
	addl	ADJK(%a6),%d0	| ...ADJUST K, ORIGINAL INPUT MAY BE  DENORM.
	lea	LOGTBL,%a0 	| ...BASE ADDRESS OF 1/F AND LOG(F)
	fmovel	%d0,%fp1		| ...CONVERT K TO FLOATING-POINT FORMAT

|--WHILE THE CONVERSION IS GOING ON, WE GET F AND ADDRESS OF 1/F
	movel	#0x3FFF0000,X(%a6)	| ...X IS NOW Y, I.E. 2^(-K)*X
	movel	XFRAC(%a6),FFRAC(%a6)
	andil	#0xFE000000,FFRAC(%a6) | ...FIRST 7 BITS OF Y
	oril	#0x01000000,FFRAC(%a6) | ...GET F: ATTACH A 1 AT THE EIGHTH BIT
	movel	FFRAC(%a6),%d0	| ...READY TO GET ADDRESS OF 1/F
	andil	#0x7E000000,%d0	
	asrl	#8,%d0
	asrl	#8,%d0
	asrl	#4,%d0		| ...SHIFTED 20, D0 IS THE DISPLACEMENT
	addal	%d0,%a0		| ...A0 IS THE ADDRESS FOR 1/F

	fmovex	X(%a6),%fp0
	movel	#0x3fff0000,F(%a6)
	clrl	F+8(%a6)
	fsubx	F(%a6),%fp0		| ...Y-F
	fmovemx %fp2-%fp2/%fp3,-(%sp)	| ...SAVE FP2 WHILE FP0 IS NOT READY
|--SUMMARY: FP0 IS Y-F, A0 IS ADDRESS OF 1/F, FP1 IS K
|--REGISTERS SAVED: FPCR, FP1, FP2

LP1CONT1:
|--AN RE-ENTRY POINT FOR LOGNP1
	fmulx	(%a0),%fp0	| ...FP0 IS U = (Y-F)/F
	fmulx	LOGOF2,%fp1	| ...GET K*LOG2 WHILE FP0 IS NOT READY
	fmovex	%fp0,%fp2
	fmulx	%fp2,%fp2		| ...FP2 IS V=U*U
	fmovex	%fp1,KLOG2(%a6)	| ...PUT K*LOG2 IN MEMORY, FREE FP1

|--LOG(1+U) IS APPROXIMATED BY
|--U + V*(A1+U*(A2+U*(A3+U*(A4+U*(A5+U*A6))))) WHICH IS
|--[U + V*(A1+V*(A3+V*A5))]  +  [U*V*(A2+V*(A4+V*A6))]

	fmovex	%fp2,%fp3
	fmovex	%fp2,%fp1	

	fmuld	LOGA6,%fp1	| ...V*A6
	fmuld	LOGA5,%fp2	| ...V*A5

	faddd	LOGA4,%fp1	| ...A4+V*A6
	faddd	LOGA3,%fp2	| ...A3+V*A5

	fmulx	%fp3,%fp1		| ...V*(A4+V*A6)
	fmulx	%fp3,%fp2		| ...V*(A3+V*A5)

	faddd	LOGA2,%fp1	| ...A2+V*(A4+V*A6)
	faddd	LOGA1,%fp2	| ...A1+V*(A3+V*A5)

	fmulx	%fp3,%fp1		| ...V*(A2+V*(A4+V*A6))
	addal	#16,%a0		| ...ADDRESS OF LOG(F)
	fmulx	%fp3,%fp2		| ...V*(A1+V*(A3+V*A5)), FP3 RELEASED

	fmulx	%fp0,%fp1		| ...U*V*(A2+V*(A4+V*A6))
	faddx	%fp2,%fp0		| ...U+V*(A1+V*(A3+V*A5)), FP2 RELEASED

	faddx	(%a0),%fp1	| ...LOG(F)+U*V*(A2+V*(A4+V*A6))
	fmovemx  (%sp)+,%fp2-%fp2/%fp3	| ...RESTORE FP2
	faddx	%fp1,%fp0		| ...FP0 IS LOG(F) + LOG(1+U)

	fmovel	%d1,%fpcr
	faddx	KLOG2(%a6),%fp0	| ...FINAL ADD
	bra	t_frcinx


LOGNEAR1:
|--REGISTERS SAVED: FPCR, FP1. FP0 CONTAINS THE INPUT.
	fmovex	%fp0,%fp1
	fsubs	one,%fp1		| ...FP1 IS X-1
	fadds	one,%fp0		| ...FP0 IS X+1
	faddx	%fp1,%fp1		| ...FP1 IS 2(X-1)
|--LOG(X) = LOG(1+U/2)-LOG(1-U/2) WHICH IS AN ODD POLYNOMIAL
|--IN U, U = 2(X-1)/(X+1) = FP1/FP0

LP1CONT2:
|--THIS IS AN RE-ENTRY POINT FOR LOGNP1
	fdivx	%fp0,%fp1		| ...FP1 IS U
	fmovemx %fp2-%fp2/%fp3,-(%sp)	 | ...SAVE FP2
|--REGISTERS SAVED ARE NOW FPCR,FP1,FP2,FP3
|--LET V=U*U, W=V*V, CALCULATE
|--U + U*V*(B1 + V*(B2 + V*(B3 + V*(B4 + V*B5)))) BY
|--U + U*V*(  [B1 + W*(B3 + W*B5)]  +  [V*(B2 + W*B4)]  )
	fmovex	%fp1,%fp0
	fmulx	%fp0,%fp0	| ...FP0 IS V
	fmovex	%fp1,SAVEU(%a6) | ...STORE U IN MEMORY, FREE FP1
	fmovex	%fp0,%fp1	
	fmulx	%fp1,%fp1	| ...FP1 IS W

	fmoved	LOGB5,%fp3
	fmoved	LOGB4,%fp2

	fmulx	%fp1,%fp3	| ...W*B5
	fmulx	%fp1,%fp2	| ...W*B4

	faddd	LOGB3,%fp3 | ...B3+W*B5
	faddd	LOGB2,%fp2 | ...B2+W*B4

	fmulx	%fp3,%fp1	| ...W*(B3+W*B5), FP3 RELEASED

	fmulx	%fp0,%fp2	| ...V*(B2+W*B4)

	faddd	LOGB1,%fp1 | ...B1+W*(B3+W*B5)
	fmulx	SAVEU(%a6),%fp0 | ...FP0 IS U*V

	faddx	%fp2,%fp1	| ...B1+W*(B3+W*B5) + V*(B2+W*B4), FP2 RELEASED
	fmovemx (%sp)+,%fp2-%fp2/%fp3 | ...FP2 RESTORED

	fmulx	%fp1,%fp0	| ...U*V*( [B1+W*(B3+W*B5)] + [V*(B2+W*B4)] )

	fmovel	%d1,%fpcr
	faddx	SAVEU(%a6),%fp0		
	bra	t_frcinx
	rts

LOGNEG:
|--REGISTERS SAVED FPCR. LOG(-VE) IS INVALID
	bra	t_operr

	.global	slognp1d
slognp1d:
|--ENTRY POINT FOR LOG(1+Z) FOR DENORMALIZED INPUT
| Simply return the denorm

	bra	t_extdnrm

	.global	slognp1
slognp1:
|--ENTRY POINT FOR LOG(1+X) FOR X FINITE, NON-ZERO, NOT NAN'S

	fmovex	(%a0),%fp0	| ...LOAD INPUT
	fabsx	%fp0		|test magnitude
	fcmpx	LTHOLD,%fp0	|compare with min threshold
	fbgt	LP1REAL		|if greater, continue
	fmovel	#0,%fpsr		|clr N flag from compare
	fmovel	%d1,%fpcr
	fmovex	(%a0),%fp0	|return signed argument
	bra	t_frcinx

LP1REAL:
	fmovex	(%a0),%fp0	| ...LOAD INPUT
	movel	#0x00000000,ADJK(%a6)
	fmovex	%fp0,%fp1	| ...FP1 IS INPUT Z
	fadds	one,%fp0	| ...X := ROUND(1+Z)
	fmovex	%fp0,X(%a6)
	movew	XFRAC(%a6),XDCARE(%a6)
	movel	X(%a6),%d0
	cmpil	#0,%d0
	ble	LP1NEG0	| ...LOG OF ZERO OR -VE
	cmp2l	BOUNDS2,%d0
	bcs	LOGMAIN	| ...BOUNDS2 IS [1/2,3/2]
|--IF 1+Z > 3/2 OR 1+Z < 1/2, THEN X, WHICH IS ROUNDING 1+Z,
|--CONTAINS AT LEAST 63 BITS OF INFORMATION OF Z. IN THAT CASE,
|--SIMPLY INVOKE LOG(X) FOR LOG(1+Z).

LP1NEAR1:
|--NEXT SEE IF EXP(-1/16) < X < EXP(1/16)
	cmp2l	BOUNDS1,%d0
	bcss	LP1CARE

LP1ONE16:
|--EXP(-1/16) < X < EXP(1/16). LOG(1+Z) = LOG(1+U/2) - LOG(1-U/2)
|--WHERE U = 2Z/(2+Z) = 2Z/(1+X).
	faddx	%fp1,%fp1	| ...FP1 IS 2Z
	fadds	one,%fp0	| ...FP0 IS 1+X
|--U = FP1/FP0
	bra	LP1CONT2

LP1CARE:
|--HERE WE USE THE USUAL TABLE DRIVEN APPROACH. CARE HAS TO BE
|--TAKEN BECAUSE 1+Z CAN HAVE 67 BITS OF INFORMATION AND WE MUST
|--PRESERVE ALL THE INFORMATION. BECAUSE 1+Z IS IN [1/2,3/2],
|--THERE ARE ONLY TWO CASES.
|--CASE 1: 1+Z < 1, THEN K = -1 AND Y-F = (2-F) + 2Z
|--CASE 2: 1+Z > 1, THEN K = 0  AND Y-F = (1-F) + Z
|--ON RETURNING TO LP1CONT1, WE MUST HAVE K IN FP1, ADDRESS OF
|--(1/F) IN A0, Y-F IN FP0, AND FP2 SAVED.

	movel	XFRAC(%a6),FFRAC(%a6)
	andil	#0xFE000000,FFRAC(%a6)
	oril	#0x01000000,FFRAC(%a6)	| ...F OBTAINED
	cmpil	#0x3FFF8000,%d0	| ...SEE IF 1+Z > 1
	bges	KISZERO

KISNEG1:
	fmoves	TWO,%fp0
	movel	#0x3fff0000,F(%a6)
	clrl	F+8(%a6)
	fsubx	F(%a6),%fp0	| ...2-F
	movel	FFRAC(%a6),%d0
	andil	#0x7E000000,%d0
	asrl	#8,%d0
	asrl	#8,%d0
	asrl	#4,%d0		| ...D0 CONTAINS DISPLACEMENT FOR 1/F
	faddx	%fp1,%fp1		| ...GET 2Z
	fmovemx %fp2-%fp2/%fp3,-(%sp)	| ...SAVE FP2 
	faddx	%fp1,%fp0		| ...FP0 IS Y-F = (2-F)+2Z
	lea	LOGTBL,%a0	| ...A0 IS ADDRESS OF 1/F
	addal	%d0,%a0
	fmoves	negone,%fp1	| ...FP1 IS K = -1
	bra	LP1CONT1

KISZERO:
	fmoves	one,%fp0
	movel	#0x3fff0000,F(%a6)
	clrl	F+8(%a6)
	fsubx	F(%a6),%fp0		| ...1-F
	movel	FFRAC(%a6),%d0
	andil	#0x7E000000,%d0
	asrl	#8,%d0
	asrl	#8,%d0
	asrl	#4,%d0
	faddx	%fp1,%fp0		| ...FP0 IS Y-F
	fmovemx %fp2-%fp2/%fp3,-(%sp)	| ...FP2 SAVED
	lea	LOGTBL,%a0
	addal	%d0,%a0	 	| ...A0 IS ADDRESS OF 1/F
	fmoves	zero,%fp1	| ...FP1 IS K = 0
	bra	LP1CONT1

LP1NEG0:
|--FPCR SAVED. D0 IS X IN COMPACT FORM.
	cmpil	#0,%d0
	blts	LP1NEG
LP1ZERO:
	fmoves	negone,%fp0

	fmovel	%d1,%fpcr
	bra t_dz

LP1NEG:
	fmoves	zero,%fp0

	fmovel	%d1,%fpcr
	bra	t_operr

	|end