summaryrefslogtreecommitdiff
path: root/sys/arch/m68k/fpsp/scosh.sa
blob: 64aafbfaa3d3e5c19361b55b72387e41ac43160f (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
*	$OpenBSD: scosh.sa,v 1.2 1996/05/29 21:05:36 niklas Exp $
*	$NetBSD: scosh.sa,v 1.2 1994/10/26 07:49:39 cgd Exp $

*	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.

*
*	scosh.sa 3.1 12/10/90
*
*	The entry point sCosh computes the hyperbolic cosine of
*	an input argument; sCoshd does the same except for denormalized
*	input.
*
*	Input: Double-extended number X in location pointed to
*		by address register a0.
*
*	Output: The value cosh(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 sCOSH takes approximately 250 cycles.
*
*	Algorithm:
*
*	COSH
*	1. If |X| > 16380 log2, go to 3.
*
*	2. (|X| <= 16380 log2) Cosh(X) is obtained by the formulae
*		y = |X|, z = exp(Y), and
*		cosh(X) = (1/2)*( z + 1/z ).
*		Exit.
*
*	3. (|X| > 16380 log2). If |X| > 16480 log2, go to 5.
*
*	4. (16380 log2 < |X| <= 16480 log2)
*		cosh(X) = sign(X) * exp(|X|)/2.
*		However, invoking exp(|X|) may cause premature overflow.
*		Thus, we calculate sinh(X) as follows:
*		Y	:= |X|
*		Fact	:=	2**(16380)
*		Y'	:= Y - 16381 log2
*		cosh(X) := Fact * exp(Y').
*		Exit.
*
*	5. (|X| > 16480 log2) sinh(X) must overflow. Return
*		Huge*Huge to generate overflow and an infinity with
*		the appropriate sign. Huge is the largest finite number in
*		extended format. Exit.
*

SCOSH	IDNT	2,1 Motorola 040 Floating Point Software Package

	section	8

	xref	t_ovfl
	xref	t_frcinx
	xref	setox

T1	DC.L $40C62D38,$D3D64634 ... 16381 LOG2 LEAD
T2	DC.L $3D6F90AE,$B1E75CC7 ... 16381 LOG2 TRAIL

TWO16380 DC.L $7FFB0000,$80000000,$00000000,$00000000

	xdef	scoshd
scoshd:
*--COSH(X) = 1 FOR DENORMALIZED X

	FMOVE.S		#:3F800000,FP0

	FMOVE.L		d1,FPCR
	FADD.S		#:00800000,FP0
	bra		t_frcinx

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

	move.l		(a0),d0
	move.w		4(a0),d0
	ANDI.L		#$7FFFFFFF,d0
	CMPI.L		#$400CB167,d0
	BGT.B		COSHBIG

*--THIS IS THE USUAL CASE, |X| < 16380 LOG2
*--COSH(X) = (1/2) * ( EXP(X) + 1/EXP(X) )

	FABS.X		FP0		...|X|

	move.l		d1,-(sp)
	clr.l		d1
	fmovem.x	fp0,(a0)	;pass parameter to setox
	bsr		setox		...FP0 IS EXP(|X|)
	FMUL.S		#:3F000000,FP0	...(1/2)EXP(|X|)
	move.l		(sp)+,d1

	FMOVE.S		#:3E800000,FP1	...(1/4)
	FDIV.X		FP0,FP1	 	...1/(2 EXP(|X|))

	FMOVE.L		d1,FPCR
	FADD.X		fp1,FP0

	bra		t_frcinx

COSHBIG:
	CMPI.L		#$400CB2B3,d0
	BGT.B		COSHHUGE

	FABS.X		FP0
	FSUB.D		T1(pc),FP0		...(|X|-16381LOG2_LEAD)
	FSUB.D		T2(pc),FP0		...|X| - 16381 LOG2, ACCURATE

	move.l		d1,-(sp)
	clr.l		d1
	fmovem.x	fp0,(a0)
	bsr		setox
	fmove.l		(sp)+,fpcr

	FMUL.X		TWO16380(pc),FP0
	bra		t_frcinx

COSHHUGE:
	fmove.l		#0,fpsr		;clr N bit if set by source
	bclr.b		#7,(a0)		;always return positive value
	fmovem.x	(a0),fp0
	bra		t_ovfl

	end