NetBSD/sys/arch/m68k/fpsp/stanh.sa

*	MOTOROLA MICROPROCESSOR & MEMORY TECHNOLOGY GROUP
*	M68000 Hi-Performance Microprocessor Division
*	M68040 Software Package 
*
*	M68040 Software Package Copyright (c) 1993, 1994 Motorola Inc.
*	All rights reserved.
*
*	THE SOFTWARE is provided on an "AS IS" basis and without warranty.
*	To the maximum extent permitted by applicable law,
*	MOTOROLA DISCLAIMS ALL WARRANTIES WHETHER EXPRESS OR IMPLIED,
*	INCLUDING IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
*	PARTICULAR PURPOSE and any warranty against infringement with
*	regard to the SOFTWARE (INCLUDING ANY MODIFIED VERSIONS THEREOF)
*	and any accompanying written materials. 
*
*	To the maximum extent permitted by applicable law,
*	IN NO EVENT SHALL MOTOROLA BE LIABLE FOR ANY DAMAGES WHATSOEVER
*	(INCLUDING WITHOUT LIMITATION, DAMAGES FOR LOSS OF BUSINESS
*	PROFITS, BUSINESS INTERRUPTION, LOSS OF BUSINESS INFORMATION, OR
*	OTHER PECUNIARY LOSS) ARISING OF THE USE OR INABILITY TO USE THE
*	SOFTWARE.  Motorola assumes no responsibility for the maintenance
*	and support of the SOFTWARE.  
*
*	You are hereby granted a copyright license to use, modify, and
*	distribute the SOFTWARE so long as this entire notice is retained
*	without alteration in any modified and/or redistributed versions,
*	and that such modified versions are clearly identified as such.
*	No licenses are granted by implication, estoppel or otherwise
*	under any patents or trademarks of Motorola, Inc.

*
*	stanh.sa 3.1 12/10/90
*
*	The entry point sTanh computes the hyperbolic tangent of
*	an input argument; sTanhd does the same except for denormalized
*	input.
*
*	Input: Double-extended number X in location pointed to
*		by address register a0.
*
*	Output: The value tanh(X) returned in floating-point register Fp0.
*
*	Accuracy and Monotonicity: The returned result is within 3 ulps in
*		64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
*		result is subsequently rounded to double precision. The
*		result is provably monotonic in double precision.
*
*	Speed: The program stanh takes approximately 270 cycles.
*
*	Algorithm:
*
*	TANH
*	1. If |X| >= (5/2) log2 or |X| <= 2**(-40), go to 3.
*
*	2. (2**(-40) < |X| < (5/2) log2) Calculate tanh(X) by
*		sgn := sign(X), y := 2|X|, z := expm1(Y), and
*		tanh(X) = sgn*( z/(2+z) ).
*		Exit.
*
*	3. (|X| <= 2**(-40) or |X| >= (5/2) log2). If |X| < 1,
*		go to 7.
*
*	4. (|X| >= (5/2) log2) If |X| >= 50 log2, go to 6.
*
*	5. ((5/2) log2 <= |X| < 50 log2) Calculate tanh(X) by
*		sgn := sign(X), y := 2|X|, z := exp(Y),
*		tanh(X) = sgn - [ sgn*2/(1+z) ].
*		Exit.
*
*	6. (|X| >= 50 log2) Tanh(X) = +-1 (round to nearest). Thus, we
*		calculate Tanh(X) by
*		sgn := sign(X), Tiny := 2**(-126),
*		tanh(X) := sgn - sgn*Tiny.
*		Exit.
*
*	7. (|X| < 2**(-40)). Tanh(X) = X.	Exit.
*

STANH	IDNT	2,1 Motorola 040 Floating Point Software Package

	section	8
	
	include fpsp.h

X	equ	FP_SCR5
XDCARE	equ	X+2
XFRAC	equ	X+4

SGN	equ	L_SCR3

V	equ	FP_SCR6

BOUNDS1	DC.L $3FD78000,$3FFFDDCE ... 2^(-40), (5/2)LOG2

	xref	t_frcinx
	xref	t_extdnrm
	xref	setox
	xref	setoxm1

	xdef	stanhd
stanhd:
*--TANH(X) = X FOR DENORMALIZED X

	bra		t_extdnrm

	xdef	stanh
stanh:
	FMOVE.X		(a0),FP0	...LOAD INPUT

	FMOVE.X		FP0,X(a6)
	move.l		(a0),d0
	move.w		4(a0),d0
	MOVE.L		D0,X(a6)
	AND.L		#$7FFFFFFF,D0
	CMP2.L		BOUNDS1(pc),D0	...2**(-40) < |X| < (5/2)LOG2 ?
	BCS.B		TANHBORS

*--THIS IS THE USUAL CASE
*--Y = 2|X|, Z = EXPM1(Y), TANH(X) = SIGN(X) * Z / (Z+2).

	MOVE.L		X(a6),D0
	MOVE.L		D0,SGN(a6)
	AND.L		#$7FFF0000,D0
	ADD.L		#$00010000,D0	...EXPONENT OF 2|X|
	MOVE.L		D0,X(a6)
	AND.L		#$80000000,SGN(a6)
	FMOVE.X		X(a6),FP0		...FP0 IS Y = 2|X|

	move.l		d1,-(a7)
	clr.l		d1
	fmovem.x	fp0,(a0)
	bsr		setoxm1	 	...FP0 IS Z = EXPM1(Y)
	move.l		(a7)+,d1

	FMOVE.X		FP0,FP1
	FADD.S		#:40000000,FP1	...Z+2
	MOVE.L		SGN(a6),D0
	FMOVE.X		FP1,V(a6)
	EOR.L		D0,V(a6)

	FMOVE.L		d1,FPCR		;restore users exceptions
	FDIV.X		V(a6),FP0
	bra		t_frcinx

TANHBORS:
	CMP.L		#$3FFF8000,D0
	BLT.W		TANHSM

	CMP.L		#$40048AA1,D0
	BGT.W		TANHHUGE

*-- (5/2) LOG2 < |X| < 50 LOG2,
*--TANH(X) = 1 - (2/[EXP(2X)+1]). LET Y = 2|X|, SGN = SIGN(X),
*--TANH(X) = SGN -	SGN*2/[EXP(Y)+1].

	MOVE.L		X(a6),D0
	MOVE.L		D0,SGN(a6)
	AND.L		#$7FFF0000,D0
	ADD.L		#$00010000,D0	...EXPO OF 2|X|
	MOVE.L		D0,X(a6)		...Y = 2|X|
	AND.L		#$80000000,SGN(a6)
	MOVE.L		SGN(a6),D0
	FMOVE.X		X(a6),FP0		...Y = 2|X|

	move.l		d1,-(a7)
	clr.l		d1
	fmovem.x	fp0,(a0)
	bsr		setox		...FP0 IS EXP(Y)
	move.l		(a7)+,d1
	move.l		SGN(a6),d0
	FADD.S		#:3F800000,FP0	...EXP(Y)+1

	EOR.L		#$C0000000,D0	...-SIGN(X)*2
	FMOVE.S		d0,FP1		...-SIGN(X)*2 IN SGL FMT
	FDIV.X		FP0,FP1	 	...-SIGN(X)2 / [EXP(Y)+1 ]

	MOVE.L		SGN(a6),D0
	OR.L		#$3F800000,D0	...SGN
	FMOVE.S		d0,FP0		...SGN IN SGL FMT

	FMOVE.L		d1,FPCR		;restore users exceptions
	FADD.X		fp1,FP0

	bra		t_frcinx

TANHSM:
	CLR.W		XDCARE(a6)

	FMOVE.L		d1,FPCR		;restore users exceptions
	FMOVE.X		X(a6),FP0		;last inst - possible exception set

	bra		t_frcinx

TANHHUGE:
*---RETURN SGN(X) - SGN(X)EPS
	MOVE.L		X(a6),D0
	AND.L		#$80000000,D0
	OR.L		#$3F800000,D0
	FMOVE.S		d0,FP0
	AND.L		#$80000000,D0
	EOR.L		#$80800000,D0	...-SIGN(X)*EPS

	FMOVE.L		d1,FPCR		;restore users exceptions
	FADD.S		d0,FP0

	bra		t_frcinx

	end