l_slogn.s

Vxworks OS source code

     movel	d5,a6@(XFRAC+4)
     negl	d2
     movel	d2,a6@(ADJK)
     fmovex	a6@(X),fp0
     moveml	A7@+,d2-d7		|...restore registers
     lea	a6@(X),a0
     jra 	LOGBGN			|...begin regular log(X)


HiX_not0:
     clrl	d6
     bfffo	d4{#0:#32},d6		|...find first 1
     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,a6@(X)
     movel	d4,a6@(XFRAC)
     movel	d5,a6@(XFRAC+4)
     negl	d2
     movel	d2,a6@(ADJK)
     fmovex	a6@(X),fp0
     moveml	A7@+,d2-d7		|...restore registers
     lea	a6@(X),a0
     jra 	LOGBGN			|...begin regular log(X)


	.globl	__l_slogn
__l_slogn:
/* |--ENTRY POINT FOR LOG(X) FOR X FINITE, NON-ZERO, NOT NAN'S */

	fmovex		A0@,fp0	|...lOAD INPUT
	movel		#0x00000000,a6@(ADJK)

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

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

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

	cmpil	#0,d0		|...CHECK IF X IS NEGATIVE
	jlt 	LOGNEG		|...lOG OF NEGATIVE ARGUMENT IS INVALID
	cmp2l	BOUNDS1,d0	|...X IS POSITIVE, CHECK IF X IS NEAR 1
	jcc 	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	a6@(ADJK),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,a6@(X)	|...X IS NOW Y, I.E. 2^(-K)*X
	movel	a6@(XFRAC),a6@(FFRAC)
	andil	#0xFE000000,a6@(FFRAC) |...FIRST 7 BITS OF Y
	oril	#0x01000000,a6@(FFRAC) |...GET F: ATTACH A 1 AT THE EIGHTH BIT
	movel	a6@(FFRAC),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	a6@(X),fp0
	movel	#0x3fff0000,a6@(F)
	clrl	a6@(F+8)
	fsubx	a6@(F),fp0		|...Y-F
	fmovemx fp2/fp3,a7@-	|...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,a6@(KLOG2)	|...PUT K*LOG2 IN MEMEORY, 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 	a7@+,fp2/fp3	|...RESTORE fp2
	faddx	fp1,fp0		|...FP0 IS LOG(F) + LOG(1+U)

	fmovel	d1,fpcr
	faddx	a6@(KLOG2),fp0	|...FINAL ADD
	jra 	__l_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/fp3,a7@-	 |...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,a6@(SAVEU) |...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	a6@(SAVEU),fp0 |...FP0 IS U*V

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

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

	fmovel	d1,fpcr
	faddx	a6@(SAVEU),fp0
	jra 	__l_t_frcinx
	rts

LOGNEG:
|--REGISTERS SAVED fpcr. LOG(-VE) IS INVALID
	jra 	__l_t_operr

	.globl	__l_slognp1d
__l_slognp1d:
|--ENTRY POINT FOR LOG(1+Z) FOR DENORMALIZED INPUT
| Simply return the __l_denorm

	jra 	__l_t_extdnrm

	.globl	__l_slognp1
__l_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
	jra 	__l_t_frcinx

LP1REAL:
	fmovex	A0@,fp0	|...lOAD INPUT
	movel	#0x00000000,a6@(ADJK)
	fmovex	fp0,fp1	|...FP1 IS INPUT Z
	fadds	one,fp0	|...X := ROUND(1+Z)
	fmovex	fp0,a6@(X)
	movew	a6@(XFRAC),a6@(XDCARE)
	movel	a6@(X),d0
	cmpil	#0,d0
	jle 	LP1NEG0	|...lOG OF ZERO OR -VE
	cmp2l	BOUNDS2,d0
	jcs 	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
	jcs 	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
	jra 	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	a6@(XFRAC),a6@(FFRAC)
	andil	#0xFE000000,a6@(FFRAC)
	oril	#0x01000000,a6@(FFRAC)	|...F OBTAINED
	cmpil	#0x3FFF8000,d0	|...SEE IF 1+Z > 1
	jge 	KISZERO

KISNEG1:
	fmoves	TWO,fp0
	movel	#0x3fff0000,a6@(F)
	clrl	a6@(F+8)
	fsubx	a6@(F),fp0	|...2-F
	movel	a6@(FFRAC),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/fp3,a7@-	|...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
	jra 	LP1CONT1

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

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

	fmovel	d1,fpcr
	jra  __l_t_dz

LP1NEG:
	fmoves	zero,fp0

	fmovel	d1,fpcr
	jra 	__l_t_operr

|	end