slogn.s

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

     movel	%d6,%d2			// ...get k
     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