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
|
* $OpenBSD: sacos.sa,v 1.2 1996/05/29 21:05:34 niklas Exp $
* $NetBSD: sacos.sa,v 1.3 1994/10/26 07:49:27 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.
*
* sacos.sa 3.3 12/19/90
*
* Description: The entry point sAcos computes the inverse cosine of
* an input argument; sAcosd does the same except for denormalized
* input.
*
* Input: Double-extended number X in location pointed to
* by address register a0.
*
* Output: The value arccos(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 sCOS takes approximately 310 cycles.
*
* Algorithm:
*
* ACOS
* 1. If |X| >= 1, go to 3.
*
* 2. (|X| < 1) Calculate acos(X) by
* z := (1-X) / (1+X)
* acos(X) = 2 * atan( sqrt(z) ).
* Exit.
*
* 3. If |X| > 1, go to 5.
*
* 4. (|X| = 1) If X > 0, return 0. Otherwise, return Pi. Exit.
*
* 5. (|X| > 1) Generate an invalid operation by 0 * infinity.
* Exit.
*
SACOS IDNT 2,1 Motorola 040 Floating Point Software Package
section 8
PI DC.L $40000000,$C90FDAA2,$2168C235,$00000000
PIBY2 DC.L $3FFF0000,$C90FDAA2,$2168C235,$00000000
xref t_operr
xref t_frcinx
xref satan
xdef sacosd
sacosd:
*--ACOS(X) = PI/2 FOR DENORMALIZED X
fmove.l d1,fpcr ...load user's rounding mode/precision
FMOVE.X PIBY2,FP0
bra t_frcinx
xdef sacos
sacos:
FMOVE.X (a0),FP0 ...LOAD INPUT
move.l (a0),d0 ...pack exponent with upper 16 fraction
move.w 4(a0),d0
ANDI.L #$7FFFFFFF,D0
CMPI.L #$3FFF8000,D0
BGE.B ACOSBIG
*--THIS IS THE USUAL CASE, |X| < 1
*--ACOS(X) = 2 * ATAN( SQRT( (1-X)/(1+X) ) )
FMOVE.S #:3F800000,FP1
FADD.X FP0,FP1 ...1+X
FNEG.X FP0 ... -X
FADD.S #:3F800000,FP0 ...1-X
FDIV.X FP1,FP0 ...(1-X)/(1+X)
FSQRT.X FP0 ...SQRT((1-X)/(1+X))
fmovem.x fp0,(a0) ...overwrite input
move.l d1,-(sp) ;save original users fpcr
clr.l d1
bsr satan ...ATAN(SQRT([1-X]/[1+X]))
fMOVE.L (sp)+,fpcr ;restore users exceptions
FADD.X FP0,FP0 ...2 * ATAN( STUFF )
bra t_frcinx
ACOSBIG:
FABS.X FP0
FCMP.S #:3F800000,FP0
fbgt t_operr ;cause an operr exception
*--|X| = 1, ACOS(X) = 0 OR PI
move.l (a0),d0 ...pack exponent with upper 16 fraction
move.w 4(a0),d0
TST.L D0 ;D0 has original exponent+fraction
BGT.B ACOSP1
*--X = -1
*Returns PI and inexact exception
FMOVE.X PI,FP0
FMOVE.L d1,FPCR
FADD.S #:00800000,FP0 ;cause an inexact exception to be put
* ;into the 040 - will not trap until next
* ;fp inst.
bra t_frcinx
ACOSP1:
FMOVE.L d1,FPCR
FMOVE.S #:00000000,FP0
rts ;Facos of +1 is exact
end
|
