Commit | Line | Data |
---|---|---|
1da177e4 LT |
1 | | |
2 | | ssin.sa 3.3 7/29/91 | |
3 | | | |
4 | | The entry point sSIN computes the sine of an input argument | |
5 | | sCOS computes the cosine, and sSINCOS computes both. The | |
6 | | corresponding entry points with a "d" computes the same | |
7 | | corresponding function values for denormalized inputs. | |
8 | | | |
9 | | Input: Double-extended number X in location pointed to | |
10 | | by address register a0. | |
11 | | | |
12 | | Output: The function value sin(X) or cos(X) returned in Fp0 if SIN or | |
13 | | COS is requested. Otherwise, for SINCOS, sin(X) is returned | |
14 | | in Fp0, and cos(X) is returned in Fp1. | |
15 | | | |
16 | | Modifies: Fp0 for SIN or COS; both Fp0 and Fp1 for SINCOS. | |
17 | | | |
18 | | Accuracy and Monotonicity: The returned result is within 1 ulp in | |
19 | | 64 significant bit, i.e. within 0.5001 ulp to 53 bits if the | |
20 | | result is subsequently rounded to double precision. The | |
21 | | result is provably monotonic in double precision. | |
22 | | | |
23 | | Speed: The programs sSIN and sCOS take approximately 150 cycles for | |
24 | | input argument X such that |X| < 15Pi, which is the usual | |
25 | | situation. The speed for sSINCOS is approximately 190 cycles. | |
26 | | | |
27 | | Algorithm: | |
28 | | | |
29 | | SIN and COS: | |
30 | | 1. If SIN is invoked, set AdjN := 0; otherwise, set AdjN := 1. | |
31 | | | |
32 | | 2. If |X| >= 15Pi or |X| < 2**(-40), go to 7. | |
33 | | | |
34 | | 3. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let | |
35 | | k = N mod 4, so in particular, k = 0,1,2,or 3. Overwrite | |
36 | | k by k := k + AdjN. | |
37 | | | |
38 | | 4. If k is even, go to 6. | |
39 | | | |
40 | | 5. (k is odd) Set j := (k-1)/2, sgn := (-1)**j. Return sgn*cos(r) | |
41 | | where cos(r) is approximated by an even polynomial in r, | |
42 | | 1 + r*r*(B1+s*(B2+ ... + s*B8)), s = r*r. | |
43 | | Exit. | |
44 | | | |
45 | | 6. (k is even) Set j := k/2, sgn := (-1)**j. Return sgn*sin(r) | |
46 | | where sin(r) is approximated by an odd polynomial in r | |
47 | | r + r*s*(A1+s*(A2+ ... + s*A7)), s = r*r. | |
48 | | Exit. | |
49 | | | |
50 | | 7. If |X| > 1, go to 9. | |
51 | | | |
52 | | 8. (|X|<2**(-40)) If SIN is invoked, return X; otherwise return 1. | |
53 | | | |
54 | | 9. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 3. | |
55 | | | |
56 | | SINCOS: | |
57 | | 1. If |X| >= 15Pi or |X| < 2**(-40), go to 6. | |
58 | | | |
59 | | 2. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let | |
60 | | k = N mod 4, so in particular, k = 0,1,2,or 3. | |
61 | | | |
62 | | 3. If k is even, go to 5. | |
63 | | | |
64 | | 4. (k is odd) Set j1 := (k-1)/2, j2 := j1 (EOR) (k mod 2), i.e. | |
65 | | j1 exclusive or with the l.s.b. of k. | |
66 | | sgn1 := (-1)**j1, sgn2 := (-1)**j2. | |
67 | | SIN(X) = sgn1 * cos(r) and COS(X) = sgn2*sin(r) where | |
68 | | sin(r) and cos(r) are computed as odd and even polynomials | |
69 | | in r, respectively. Exit | |
70 | | | |
71 | | 5. (k is even) Set j1 := k/2, sgn1 := (-1)**j1. | |
72 | | SIN(X) = sgn1 * sin(r) and COS(X) = sgn1*cos(r) where | |
73 | | sin(r) and cos(r) are computed as odd and even polynomials | |
74 | | in r, respectively. Exit | |
75 | | | |
76 | | 6. If |X| > 1, go to 8. | |
77 | | | |
78 | | 7. (|X|<2**(-40)) SIN(X) = X and COS(X) = 1. Exit. | |
79 | | | |
80 | | 8. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 2. | |
81 | | | |
82 | ||
83 | | Copyright (C) Motorola, Inc. 1990 | |
84 | | All Rights Reserved | |
85 | | | |
e00d82d0 MW |
86 | | For details on the license for this file, please see the |
87 | | file, README, in this same directory. | |
1da177e4 LT |
88 | |
89 | |SSIN idnt 2,1 | Motorola 040 Floating Point Software Package | |
90 | ||
91 | |section 8 | |
92 | ||
93 | #include "fpsp.h" | |
94 | ||
95 | BOUNDS1: .long 0x3FD78000,0x4004BC7E | |
96 | TWOBYPI: .long 0x3FE45F30,0x6DC9C883 | |
97 | ||
98 | SINA7: .long 0xBD6AAA77,0xCCC994F5 | |
99 | SINA6: .long 0x3DE61209,0x7AAE8DA1 | |
100 | ||
101 | SINA5: .long 0xBE5AE645,0x2A118AE4 | |
102 | SINA4: .long 0x3EC71DE3,0xA5341531 | |
103 | ||
104 | SINA3: .long 0xBF2A01A0,0x1A018B59,0x00000000,0x00000000 | |
105 | ||
106 | SINA2: .long 0x3FF80000,0x88888888,0x888859AF,0x00000000 | |
107 | ||
108 | SINA1: .long 0xBFFC0000,0xAAAAAAAA,0xAAAAAA99,0x00000000 | |
109 | ||
110 | COSB8: .long 0x3D2AC4D0,0xD6011EE3 | |
111 | COSB7: .long 0xBDA9396F,0x9F45AC19 | |
112 | ||
113 | COSB6: .long 0x3E21EED9,0x0612C972 | |
114 | COSB5: .long 0xBE927E4F,0xB79D9FCF | |
115 | ||
116 | COSB4: .long 0x3EFA01A0,0x1A01D423,0x00000000,0x00000000 | |
117 | ||
118 | COSB3: .long 0xBFF50000,0xB60B60B6,0x0B61D438,0x00000000 | |
119 | ||
120 | COSB2: .long 0x3FFA0000,0xAAAAAAAA,0xAAAAAB5E | |
121 | COSB1: .long 0xBF000000 | |
122 | ||
123 | INVTWOPI: .long 0x3FFC0000,0xA2F9836E,0x4E44152A | |
124 | ||
125 | TWOPI1: .long 0x40010000,0xC90FDAA2,0x00000000,0x00000000 | |
126 | TWOPI2: .long 0x3FDF0000,0x85A308D4,0x00000000,0x00000000 | |
127 | ||
128 | |xref PITBL | |
129 | ||
130 | .set INARG,FP_SCR4 | |
131 | ||
132 | .set X,FP_SCR5 | |
133 | .set XDCARE,X+2 | |
134 | .set XFRAC,X+4 | |
135 | ||
136 | .set RPRIME,FP_SCR1 | |
137 | .set SPRIME,FP_SCR2 | |
138 | ||
139 | .set POSNEG1,L_SCR1 | |
140 | .set TWOTO63,L_SCR1 | |
141 | ||
142 | .set ENDFLAG,L_SCR2 | |
143 | .set N,L_SCR2 | |
144 | ||
145 | .set ADJN,L_SCR3 | |
146 | ||
147 | | xref t_frcinx | |
148 | |xref t_extdnrm | |
149 | |xref sto_cos | |
150 | ||
151 | .global ssind | |
152 | ssind: | |
153 | |--SIN(X) = X FOR DENORMALIZED X | |
154 | bra t_extdnrm | |
155 | ||
156 | .global scosd | |
157 | scosd: | |
158 | |--COS(X) = 1 FOR DENORMALIZED X | |
159 | ||
160 | fmoves #0x3F800000,%fp0 | |
161 | | | |
162 | | 9D25B Fix: Sometimes the previous fmove.s sets fpsr bits | |
163 | | | |
164 | fmovel #0,%fpsr | |
165 | | | |
166 | bra t_frcinx | |
167 | ||
168 | .global ssin | |
169 | ssin: | |
170 | |--SET ADJN TO 0 | |
171 | movel #0,ADJN(%a6) | |
172 | bras SINBGN | |
173 | ||
174 | .global scos | |
175 | scos: | |
176 | |--SET ADJN TO 1 | |
177 | movel #1,ADJN(%a6) | |
178 | ||
179 | SINBGN: | |
180 | |--SAVE FPCR, FP1. CHECK IF |X| IS TOO SMALL OR LARGE | |
181 | ||
182 | fmovex (%a0),%fp0 | ...LOAD INPUT | |
183 | ||
184 | movel (%a0),%d0 | |
185 | movew 4(%a0),%d0 | |
186 | fmovex %fp0,X(%a6) | |
187 | andil #0x7FFFFFFF,%d0 | ...COMPACTIFY X | |
188 | ||
189 | cmpil #0x3FD78000,%d0 | ...|X| >= 2**(-40)? | |
190 | bges SOK1 | |
191 | bra SINSM | |
192 | ||
193 | SOK1: | |
194 | cmpil #0x4004BC7E,%d0 | ...|X| < 15 PI? | |
195 | blts SINMAIN | |
196 | bra REDUCEX | |
197 | ||
198 | SINMAIN: | |
199 | |--THIS IS THE USUAL CASE, |X| <= 15 PI. | |
200 | |--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP. | |
201 | fmovex %fp0,%fp1 | |
202 | fmuld TWOBYPI,%fp1 | ...X*2/PI | |
203 | ||
204 | |--HIDE THE NEXT THREE INSTRUCTIONS | |
205 | lea PITBL+0x200,%a1 | ...TABLE OF N*PI/2, N = -32,...,32 | |
206 | ||
207 | ||
208 | |--FP1 IS NOW READY | |
209 | fmovel %fp1,N(%a6) | ...CONVERT TO INTEGER | |
210 | ||
211 | movel N(%a6),%d0 | |
212 | asll #4,%d0 | |
213 | addal %d0,%a1 | ...A1 IS THE ADDRESS OF N*PIBY2 | |
214 | | ...WHICH IS IN TWO PIECES Y1 & Y2 | |
215 | ||
216 | fsubx (%a1)+,%fp0 | ...X-Y1 | |
217 | |--HIDE THE NEXT ONE | |
218 | fsubs (%a1),%fp0 | ...FP0 IS R = (X-Y1)-Y2 | |
219 | ||
220 | SINCONT: | |
221 | |--continuation from REDUCEX | |
222 | ||
223 | |--GET N+ADJN AND SEE IF SIN(R) OR COS(R) IS NEEDED | |
224 | movel N(%a6),%d0 | |
225 | addl ADJN(%a6),%d0 | ...SEE IF D0 IS ODD OR EVEN | |
226 | rorl #1,%d0 | ...D0 WAS ODD IFF D0 IS NEGATIVE | |
227 | cmpil #0,%d0 | |
228 | blt COSPOLY | |
229 | ||
230 | SINPOLY: | |
231 | |--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J. | |
232 | |--THEN WE RETURN SGN*SIN(R). SGN*SIN(R) IS COMPUTED BY | |
233 | |--R' + R'*S*(A1 + S(A2 + S(A3 + S(A4 + ... + SA7)))), WHERE | |
234 | |--R' = SGN*R, S=R*R. THIS CAN BE REWRITTEN AS | |
235 | |--R' + R'*S*( [A1+T(A3+T(A5+TA7))] + [S(A2+T(A4+TA6))]) | |
236 | |--WHERE T=S*S. | |
237 | |--NOTE THAT A3 THROUGH A7 ARE STORED IN DOUBLE PRECISION | |
238 | |--WHILE A1 AND A2 ARE IN DOUBLE-EXTENDED FORMAT. | |
239 | fmovex %fp0,X(%a6) | ...X IS R | |
240 | fmulx %fp0,%fp0 | ...FP0 IS S | |
241 | |---HIDE THE NEXT TWO WHILE WAITING FOR FP0 | |
242 | fmoved SINA7,%fp3 | |
243 | fmoved SINA6,%fp2 | |
244 | |--FP0 IS NOW READY | |
245 | fmovex %fp0,%fp1 | |
246 | fmulx %fp1,%fp1 | ...FP1 IS T | |
247 | |--HIDE THE NEXT TWO WHILE WAITING FOR FP1 | |
248 | ||
249 | rorl #1,%d0 | |
250 | andil #0x80000000,%d0 | |
251 | | ...LEAST SIG. BIT OF D0 IN SIGN POSITION | |
252 | eorl %d0,X(%a6) | ...X IS NOW R'= SGN*R | |
253 | ||
254 | fmulx %fp1,%fp3 | ...TA7 | |
255 | fmulx %fp1,%fp2 | ...TA6 | |
256 | ||
257 | faddd SINA5,%fp3 | ...A5+TA7 | |
258 | faddd SINA4,%fp2 | ...A4+TA6 | |
259 | ||
260 | fmulx %fp1,%fp3 | ...T(A5+TA7) | |
261 | fmulx %fp1,%fp2 | ...T(A4+TA6) | |
262 | ||
263 | faddd SINA3,%fp3 | ...A3+T(A5+TA7) | |
264 | faddx SINA2,%fp2 | ...A2+T(A4+TA6) | |
265 | ||
266 | fmulx %fp3,%fp1 | ...T(A3+T(A5+TA7)) | |
267 | ||
268 | fmulx %fp0,%fp2 | ...S(A2+T(A4+TA6)) | |
269 | faddx SINA1,%fp1 | ...A1+T(A3+T(A5+TA7)) | |
270 | fmulx X(%a6),%fp0 | ...R'*S | |
271 | ||
272 | faddx %fp2,%fp1 | ...[A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))] | |
273 | |--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING | |
274 | |--FP2 RELEASED, RESTORE NOW AND TAKE FULL ADVANTAGE OF HIDING | |
275 | ||
276 | ||
277 | fmulx %fp1,%fp0 | ...SIN(R')-R' | |
278 | |--FP1 RELEASED. | |
279 | ||
280 | fmovel %d1,%FPCR |restore users exceptions | |
281 | faddx X(%a6),%fp0 |last inst - possible exception set | |
282 | bra t_frcinx | |
283 | ||
284 | ||
285 | COSPOLY: | |
286 | |--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J. | |
287 | |--THEN WE RETURN SGN*COS(R). SGN*COS(R) IS COMPUTED BY | |
288 | |--SGN + S'*(B1 + S(B2 + S(B3 + S(B4 + ... + SB8)))), WHERE | |
289 | |--S=R*R AND S'=SGN*S. THIS CAN BE REWRITTEN AS | |
290 | |--SGN + S'*([B1+T(B3+T(B5+TB7))] + [S(B2+T(B4+T(B6+TB8)))]) | |
291 | |--WHERE T=S*S. | |
292 | |--NOTE THAT B4 THROUGH B8 ARE STORED IN DOUBLE PRECISION | |
293 | |--WHILE B2 AND B3 ARE IN DOUBLE-EXTENDED FORMAT, B1 IS -1/2 | |
294 | |--AND IS THEREFORE STORED AS SINGLE PRECISION. | |
295 | ||
296 | fmulx %fp0,%fp0 | ...FP0 IS S | |
297 | |---HIDE THE NEXT TWO WHILE WAITING FOR FP0 | |
298 | fmoved COSB8,%fp2 | |
299 | fmoved COSB7,%fp3 | |
300 | |--FP0 IS NOW READY | |
301 | fmovex %fp0,%fp1 | |
302 | fmulx %fp1,%fp1 | ...FP1 IS T | |
303 | |--HIDE THE NEXT TWO WHILE WAITING FOR FP1 | |
304 | fmovex %fp0,X(%a6) | ...X IS S | |
305 | rorl #1,%d0 | |
306 | andil #0x80000000,%d0 | |
307 | | ...LEAST SIG. BIT OF D0 IN SIGN POSITION | |
308 | ||
309 | fmulx %fp1,%fp2 | ...TB8 | |
310 | |--HIDE THE NEXT TWO WHILE WAITING FOR THE XU | |
311 | eorl %d0,X(%a6) | ...X IS NOW S'= SGN*S | |
312 | andil #0x80000000,%d0 | |
313 | ||
314 | fmulx %fp1,%fp3 | ...TB7 | |
315 | |--HIDE THE NEXT TWO WHILE WAITING FOR THE XU | |
316 | oril #0x3F800000,%d0 | ...D0 IS SGN IN SINGLE | |
317 | movel %d0,POSNEG1(%a6) | |
318 | ||
319 | faddd COSB6,%fp2 | ...B6+TB8 | |
320 | faddd COSB5,%fp3 | ...B5+TB7 | |
321 | ||
322 | fmulx %fp1,%fp2 | ...T(B6+TB8) | |
323 | fmulx %fp1,%fp3 | ...T(B5+TB7) | |
324 | ||
325 | faddd COSB4,%fp2 | ...B4+T(B6+TB8) | |
326 | faddx COSB3,%fp3 | ...B3+T(B5+TB7) | |
327 | ||
328 | fmulx %fp1,%fp2 | ...T(B4+T(B6+TB8)) | |
329 | fmulx %fp3,%fp1 | ...T(B3+T(B5+TB7)) | |
330 | ||
331 | faddx COSB2,%fp2 | ...B2+T(B4+T(B6+TB8)) | |
332 | fadds COSB1,%fp1 | ...B1+T(B3+T(B5+TB7)) | |
333 | ||
334 | fmulx %fp2,%fp0 | ...S(B2+T(B4+T(B6+TB8))) | |
335 | |--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING | |
336 | |--FP2 RELEASED. | |
337 | ||
338 | ||
339 | faddx %fp1,%fp0 | |
340 | |--FP1 RELEASED | |
341 | ||
342 | fmulx X(%a6),%fp0 | |
343 | ||
344 | fmovel %d1,%FPCR |restore users exceptions | |
345 | fadds POSNEG1(%a6),%fp0 |last inst - possible exception set | |
346 | bra t_frcinx | |
347 | ||
348 | ||
349 | SINBORS: | |
350 | |--IF |X| > 15PI, WE USE THE GENERAL ARGUMENT REDUCTION. | |
351 | |--IF |X| < 2**(-40), RETURN X OR 1. | |
352 | cmpil #0x3FFF8000,%d0 | |
353 | bgts REDUCEX | |
354 | ||
355 | ||
356 | SINSM: | |
357 | movel ADJN(%a6),%d0 | |
358 | cmpil #0,%d0 | |
359 | bgts COSTINY | |
360 | ||
361 | SINTINY: | |
362 | movew #0x0000,XDCARE(%a6) | ...JUST IN CASE | |
363 | fmovel %d1,%FPCR |restore users exceptions | |
364 | fmovex X(%a6),%fp0 |last inst - possible exception set | |
365 | bra t_frcinx | |
366 | ||
367 | ||
368 | COSTINY: | |
369 | fmoves #0x3F800000,%fp0 | |
370 | ||
371 | fmovel %d1,%FPCR |restore users exceptions | |
372 | fsubs #0x00800000,%fp0 |last inst - possible exception set | |
373 | bra t_frcinx | |
374 | ||
375 | ||
376 | REDUCEX: | |
377 | |--WHEN REDUCEX IS USED, THE CODE WILL INEVITABLY BE SLOW. | |
378 | |--THIS REDUCTION METHOD, HOWEVER, IS MUCH FASTER THAN USING | |
379 | |--THE REMAINDER INSTRUCTION WHICH IS NOW IN SOFTWARE. | |
380 | ||
381 | fmovemx %fp2-%fp5,-(%a7) | ...save FP2 through FP5 | |
382 | movel %d2,-(%a7) | |
383 | fmoves #0x00000000,%fp1 | |
384 | |--If compact form of abs(arg) in d0=$7ffeffff, argument is so large that | |
385 | |--there is a danger of unwanted overflow in first LOOP iteration. In this | |
386 | |--case, reduce argument by one remainder step to make subsequent reduction | |
387 | |--safe. | |
388 | cmpil #0x7ffeffff,%d0 |is argument dangerously large? | |
389 | bnes LOOP | |
390 | movel #0x7ffe0000,FP_SCR2(%a6) |yes | |
391 | | ;create 2**16383*PI/2 | |
392 | movel #0xc90fdaa2,FP_SCR2+4(%a6) | |
393 | clrl FP_SCR2+8(%a6) | |
394 | ftstx %fp0 |test sign of argument | |
395 | movel #0x7fdc0000,FP_SCR3(%a6) |create low half of 2**16383* | |
396 | | ;PI/2 at FP_SCR3 | |
397 | movel #0x85a308d3,FP_SCR3+4(%a6) | |
398 | clrl FP_SCR3+8(%a6) | |
399 | fblt red_neg | |
400 | orw #0x8000,FP_SCR2(%a6) |positive arg | |
401 | orw #0x8000,FP_SCR3(%a6) | |
402 | red_neg: | |
403 | faddx FP_SCR2(%a6),%fp0 |high part of reduction is exact | |
404 | fmovex %fp0,%fp1 |save high result in fp1 | |
405 | faddx FP_SCR3(%a6),%fp0 |low part of reduction | |
406 | fsubx %fp0,%fp1 |determine low component of result | |
407 | faddx FP_SCR3(%a6),%fp1 |fp0/fp1 are reduced argument. | |
408 | ||
409 | |--ON ENTRY, FP0 IS X, ON RETURN, FP0 IS X REM PI/2, |X| <= PI/4. | |
410 | |--integer quotient will be stored in N | |
411 | |--Intermediate remainder is 66-bit long; (R,r) in (FP0,FP1) | |
412 | ||
413 | LOOP: | |
414 | fmovex %fp0,INARG(%a6) | ...+-2**K * F, 1 <= F < 2 | |
415 | movew INARG(%a6),%d0 | |
416 | movel %d0,%a1 | ...save a copy of D0 | |
417 | andil #0x00007FFF,%d0 | |
418 | subil #0x00003FFF,%d0 | ...D0 IS K | |
419 | cmpil #28,%d0 | |
420 | bles LASTLOOP | |
421 | CONTLOOP: | |
422 | subil #27,%d0 | ...D0 IS L := K-27 | |
423 | movel #0,ENDFLAG(%a6) | |
424 | bras WORK | |
425 | LASTLOOP: | |
426 | clrl %d0 | ...D0 IS L := 0 | |
427 | movel #1,ENDFLAG(%a6) | |
428 | ||
429 | WORK: | |
430 | |--FIND THE REMAINDER OF (R,r) W.R.T. 2**L * (PI/2). L IS SO CHOSEN | |
431 | |--THAT INT( X * (2/PI) / 2**(L) ) < 2**29. | |
432 | ||
433 | |--CREATE 2**(-L) * (2/PI), SIGN(INARG)*2**(63), | |
434 | |--2**L * (PIby2_1), 2**L * (PIby2_2) | |
435 | ||
436 | movel #0x00003FFE,%d2 | ...BIASED EXPO OF 2/PI | |
437 | subl %d0,%d2 | ...BIASED EXPO OF 2**(-L)*(2/PI) | |
438 | ||
439 | movel #0xA2F9836E,FP_SCR1+4(%a6) | |
440 | movel #0x4E44152A,FP_SCR1+8(%a6) | |
441 | movew %d2,FP_SCR1(%a6) | ...FP_SCR1 is 2**(-L)*(2/PI) | |
442 | ||
443 | fmovex %fp0,%fp2 | |
444 | fmulx FP_SCR1(%a6),%fp2 | |
445 | |--WE MUST NOW FIND INT(FP2). SINCE WE NEED THIS VALUE IN | |
446 | |--FLOATING POINT FORMAT, THE TWO FMOVE'S FMOVE.L FP <--> N | |
447 | |--WILL BE TOO INEFFICIENT. THE WAY AROUND IT IS THAT | |
448 | |--(SIGN(INARG)*2**63 + FP2) - SIGN(INARG)*2**63 WILL GIVE | |
449 | |--US THE DESIRED VALUE IN FLOATING POINT. | |
450 | ||
451 | |--HIDE SIX CYCLES OF INSTRUCTION | |
452 | movel %a1,%d2 | |
453 | swap %d2 | |
454 | andil #0x80000000,%d2 | |
455 | oril #0x5F000000,%d2 | ...D2 IS SIGN(INARG)*2**63 IN SGL | |
456 | movel %d2,TWOTO63(%a6) | |
457 | ||
458 | movel %d0,%d2 | |
459 | addil #0x00003FFF,%d2 | ...BIASED EXPO OF 2**L * (PI/2) | |
460 | ||
461 | |--FP2 IS READY | |
462 | fadds TWOTO63(%a6),%fp2 | ...THE FRACTIONAL PART OF FP1 IS ROUNDED | |
463 | ||
464 | |--HIDE 4 CYCLES OF INSTRUCTION; creating 2**(L)*Piby2_1 and 2**(L)*Piby2_2 | |
465 | movew %d2,FP_SCR2(%a6) | |
466 | clrw FP_SCR2+2(%a6) | |
467 | movel #0xC90FDAA2,FP_SCR2+4(%a6) | |
468 | clrl FP_SCR2+8(%a6) | ...FP_SCR2 is 2**(L) * Piby2_1 | |
469 | ||
470 | |--FP2 IS READY | |
471 | fsubs TWOTO63(%a6),%fp2 | ...FP2 is N | |
472 | ||
473 | addil #0x00003FDD,%d0 | |
474 | movew %d0,FP_SCR3(%a6) | |
475 | clrw FP_SCR3+2(%a6) | |
476 | movel #0x85A308D3,FP_SCR3+4(%a6) | |
477 | clrl FP_SCR3+8(%a6) | ...FP_SCR3 is 2**(L) * Piby2_2 | |
478 | ||
479 | movel ENDFLAG(%a6),%d0 | |
480 | ||
481 | |--We are now ready to perform (R+r) - N*P1 - N*P2, P1 = 2**(L) * Piby2_1 and | |
482 | |--P2 = 2**(L) * Piby2_2 | |
483 | fmovex %fp2,%fp4 | |
484 | fmulx FP_SCR2(%a6),%fp4 | ...W = N*P1 | |
485 | fmovex %fp2,%fp5 | |
486 | fmulx FP_SCR3(%a6),%fp5 | ...w = N*P2 | |
487 | fmovex %fp4,%fp3 | |
488 | |--we want P+p = W+w but |p| <= half ulp of P | |
489 | |--Then, we need to compute A := R-P and a := r-p | |
490 | faddx %fp5,%fp3 | ...FP3 is P | |
491 | fsubx %fp3,%fp4 | ...W-P | |
492 | ||
493 | fsubx %fp3,%fp0 | ...FP0 is A := R - P | |
494 | faddx %fp5,%fp4 | ...FP4 is p = (W-P)+w | |
495 | ||
496 | fmovex %fp0,%fp3 | ...FP3 A | |
497 | fsubx %fp4,%fp1 | ...FP1 is a := r - p | |
498 | ||
499 | |--Now we need to normalize (A,a) to "new (R,r)" where R+r = A+a but | |
500 | |--|r| <= half ulp of R. | |
501 | faddx %fp1,%fp0 | ...FP0 is R := A+a | |
502 | |--No need to calculate r if this is the last loop | |
503 | cmpil #0,%d0 | |
504 | bgt RESTORE | |
505 | ||
506 | |--Need to calculate r | |
507 | fsubx %fp0,%fp3 | ...A-R | |
508 | faddx %fp3,%fp1 | ...FP1 is r := (A-R)+a | |
509 | bra LOOP | |
510 | ||
511 | RESTORE: | |
512 | fmovel %fp2,N(%a6) | |
513 | movel (%a7)+,%d2 | |
514 | fmovemx (%a7)+,%fp2-%fp5 | |
515 | ||
516 | ||
517 | movel ADJN(%a6),%d0 | |
518 | cmpil #4,%d0 | |
519 | ||
520 | blt SINCONT | |
521 | bras SCCONT | |
522 | ||
523 | .global ssincosd | |
524 | ssincosd: | |
525 | |--SIN AND COS OF X FOR DENORMALIZED X | |
526 | ||
527 | fmoves #0x3F800000,%fp1 | |
528 | bsr sto_cos |store cosine result | |
529 | bra t_extdnrm | |
530 | ||
531 | .global ssincos | |
532 | ssincos: | |
533 | |--SET ADJN TO 4 | |
534 | movel #4,ADJN(%a6) | |
535 | ||
536 | fmovex (%a0),%fp0 | ...LOAD INPUT | |
537 | ||
538 | movel (%a0),%d0 | |
539 | movew 4(%a0),%d0 | |
540 | fmovex %fp0,X(%a6) | |
541 | andil #0x7FFFFFFF,%d0 | ...COMPACTIFY X | |
542 | ||
543 | cmpil #0x3FD78000,%d0 | ...|X| >= 2**(-40)? | |
544 | bges SCOK1 | |
545 | bra SCSM | |
546 | ||
547 | SCOK1: | |
548 | cmpil #0x4004BC7E,%d0 | ...|X| < 15 PI? | |
549 | blts SCMAIN | |
550 | bra REDUCEX | |
551 | ||
552 | ||
553 | SCMAIN: | |
554 | |--THIS IS THE USUAL CASE, |X| <= 15 PI. | |
555 | |--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP. | |
556 | fmovex %fp0,%fp1 | |
557 | fmuld TWOBYPI,%fp1 | ...X*2/PI | |
558 | ||
559 | |--HIDE THE NEXT THREE INSTRUCTIONS | |
560 | lea PITBL+0x200,%a1 | ...TABLE OF N*PI/2, N = -32,...,32 | |
561 | ||
562 | ||
563 | |--FP1 IS NOW READY | |
564 | fmovel %fp1,N(%a6) | ...CONVERT TO INTEGER | |
565 | ||
566 | movel N(%a6),%d0 | |
567 | asll #4,%d0 | |
568 | addal %d0,%a1 | ...ADDRESS OF N*PIBY2, IN Y1, Y2 | |
569 | ||
570 | fsubx (%a1)+,%fp0 | ...X-Y1 | |
571 | fsubs (%a1),%fp0 | ...FP0 IS R = (X-Y1)-Y2 | |
572 | ||
573 | SCCONT: | |
574 | |--continuation point from REDUCEX | |
575 | ||
576 | |--HIDE THE NEXT TWO | |
577 | movel N(%a6),%d0 | |
578 | rorl #1,%d0 | |
579 | ||
580 | cmpil #0,%d0 | ...D0 < 0 IFF N IS ODD | |
581 | bge NEVEN | |
582 | ||
583 | NODD: | |
584 | |--REGISTERS SAVED SO FAR: D0, A0, FP2. | |
585 | ||
586 | fmovex %fp0,RPRIME(%a6) | |
587 | fmulx %fp0,%fp0 | ...FP0 IS S = R*R | |
588 | fmoved SINA7,%fp1 | ...A7 | |
589 | fmoved COSB8,%fp2 | ...B8 | |
590 | fmulx %fp0,%fp1 | ...SA7 | |
591 | movel %d2,-(%a7) | |
592 | movel %d0,%d2 | |
593 | fmulx %fp0,%fp2 | ...SB8 | |
594 | rorl #1,%d2 | |
595 | andil #0x80000000,%d2 | |
596 | ||
597 | faddd SINA6,%fp1 | ...A6+SA7 | |
598 | eorl %d0,%d2 | |
599 | andil #0x80000000,%d2 | |
600 | faddd COSB7,%fp2 | ...B7+SB8 | |
601 | ||
602 | fmulx %fp0,%fp1 | ...S(A6+SA7) | |
603 | eorl %d2,RPRIME(%a6) | |
604 | movel (%a7)+,%d2 | |
605 | fmulx %fp0,%fp2 | ...S(B7+SB8) | |
606 | rorl #1,%d0 | |
607 | andil #0x80000000,%d0 | |
608 | ||
609 | faddd SINA5,%fp1 | ...A5+S(A6+SA7) | |
610 | movel #0x3F800000,POSNEG1(%a6) | |
611 | eorl %d0,POSNEG1(%a6) | |
612 | faddd COSB6,%fp2 | ...B6+S(B7+SB8) | |
613 | ||
614 | fmulx %fp0,%fp1 | ...S(A5+S(A6+SA7)) | |
615 | fmulx %fp0,%fp2 | ...S(B6+S(B7+SB8)) | |
616 | fmovex %fp0,SPRIME(%a6) | |
617 | ||
618 | faddd SINA4,%fp1 | ...A4+S(A5+S(A6+SA7)) | |
619 | eorl %d0,SPRIME(%a6) | |
620 | faddd COSB5,%fp2 | ...B5+S(B6+S(B7+SB8)) | |
621 | ||
622 | fmulx %fp0,%fp1 | ...S(A4+...) | |
623 | fmulx %fp0,%fp2 | ...S(B5+...) | |
624 | ||
625 | faddd SINA3,%fp1 | ...A3+S(A4+...) | |
626 | faddd COSB4,%fp2 | ...B4+S(B5+...) | |
627 | ||
628 | fmulx %fp0,%fp1 | ...S(A3+...) | |
629 | fmulx %fp0,%fp2 | ...S(B4+...) | |
630 | ||
631 | faddx SINA2,%fp1 | ...A2+S(A3+...) | |
632 | faddx COSB3,%fp2 | ...B3+S(B4+...) | |
633 | ||
634 | fmulx %fp0,%fp1 | ...S(A2+...) | |
635 | fmulx %fp0,%fp2 | ...S(B3+...) | |
636 | ||
637 | faddx SINA1,%fp1 | ...A1+S(A2+...) | |
638 | faddx COSB2,%fp2 | ...B2+S(B3+...) | |
639 | ||
640 | fmulx %fp0,%fp1 | ...S(A1+...) | |
641 | fmulx %fp2,%fp0 | ...S(B2+...) | |
642 | ||
643 | ||
644 | ||
645 | fmulx RPRIME(%a6),%fp1 | ...R'S(A1+...) | |
646 | fadds COSB1,%fp0 | ...B1+S(B2...) | |
647 | fmulx SPRIME(%a6),%fp0 | ...S'(B1+S(B2+...)) | |
648 | ||
649 | movel %d1,-(%sp) |restore users mode & precision | |
650 | andil #0xff,%d1 |mask off all exceptions | |
651 | fmovel %d1,%FPCR | |
652 | faddx RPRIME(%a6),%fp1 | ...COS(X) | |
653 | bsr sto_cos |store cosine result | |
654 | fmovel (%sp)+,%FPCR |restore users exceptions | |
655 | fadds POSNEG1(%a6),%fp0 | ...SIN(X) | |
656 | ||
657 | bra t_frcinx | |
658 | ||
659 | ||
660 | NEVEN: | |
661 | |--REGISTERS SAVED SO FAR: FP2. | |
662 | ||
663 | fmovex %fp0,RPRIME(%a6) | |
664 | fmulx %fp0,%fp0 | ...FP0 IS S = R*R | |
665 | fmoved COSB8,%fp1 | ...B8 | |
666 | fmoved SINA7,%fp2 | ...A7 | |
667 | fmulx %fp0,%fp1 | ...SB8 | |
668 | fmovex %fp0,SPRIME(%a6) | |
669 | fmulx %fp0,%fp2 | ...SA7 | |
670 | rorl #1,%d0 | |
671 | andil #0x80000000,%d0 | |
672 | faddd COSB7,%fp1 | ...B7+SB8 | |
673 | faddd SINA6,%fp2 | ...A6+SA7 | |
674 | eorl %d0,RPRIME(%a6) | |
675 | eorl %d0,SPRIME(%a6) | |
676 | fmulx %fp0,%fp1 | ...S(B7+SB8) | |
677 | oril #0x3F800000,%d0 | |
678 | movel %d0,POSNEG1(%a6) | |
679 | fmulx %fp0,%fp2 | ...S(A6+SA7) | |
680 | ||
681 | faddd COSB6,%fp1 | ...B6+S(B7+SB8) | |
682 | faddd SINA5,%fp2 | ...A5+S(A6+SA7) | |
683 | ||
684 | fmulx %fp0,%fp1 | ...S(B6+S(B7+SB8)) | |
685 | fmulx %fp0,%fp2 | ...S(A5+S(A6+SA7)) | |
686 | ||
687 | faddd COSB5,%fp1 | ...B5+S(B6+S(B7+SB8)) | |
688 | faddd SINA4,%fp2 | ...A4+S(A5+S(A6+SA7)) | |
689 | ||
690 | fmulx %fp0,%fp1 | ...S(B5+...) | |
691 | fmulx %fp0,%fp2 | ...S(A4+...) | |
692 | ||
693 | faddd COSB4,%fp1 | ...B4+S(B5+...) | |
694 | faddd SINA3,%fp2 | ...A3+S(A4+...) | |
695 | ||
696 | fmulx %fp0,%fp1 | ...S(B4+...) | |
697 | fmulx %fp0,%fp2 | ...S(A3+...) | |
698 | ||
699 | faddx COSB3,%fp1 | ...B3+S(B4+...) | |
700 | faddx SINA2,%fp2 | ...A2+S(A3+...) | |
701 | ||
702 | fmulx %fp0,%fp1 | ...S(B3+...) | |
703 | fmulx %fp0,%fp2 | ...S(A2+...) | |
704 | ||
705 | faddx COSB2,%fp1 | ...B2+S(B3+...) | |
706 | faddx SINA1,%fp2 | ...A1+S(A2+...) | |
707 | ||
708 | fmulx %fp0,%fp1 | ...S(B2+...) | |
709 | fmulx %fp2,%fp0 | ...s(a1+...) | |
710 | ||
711 | ||
712 | ||
713 | fadds COSB1,%fp1 | ...B1+S(B2...) | |
714 | fmulx RPRIME(%a6),%fp0 | ...R'S(A1+...) | |
715 | fmulx SPRIME(%a6),%fp1 | ...S'(B1+S(B2+...)) | |
716 | ||
717 | movel %d1,-(%sp) |save users mode & precision | |
718 | andil #0xff,%d1 |mask off all exceptions | |
719 | fmovel %d1,%FPCR | |
720 | fadds POSNEG1(%a6),%fp1 | ...COS(X) | |
721 | bsr sto_cos |store cosine result | |
722 | fmovel (%sp)+,%FPCR |restore users exceptions | |
723 | faddx RPRIME(%a6),%fp0 | ...SIN(X) | |
724 | ||
725 | bra t_frcinx | |
726 | ||
727 | SCBORS: | |
728 | cmpil #0x3FFF8000,%d0 | |
729 | bgt REDUCEX | |
730 | ||
731 | ||
732 | SCSM: | |
733 | movew #0x0000,XDCARE(%a6) | |
734 | fmoves #0x3F800000,%fp1 | |
735 | ||
736 | movel %d1,-(%sp) |save users mode & precision | |
737 | andil #0xff,%d1 |mask off all exceptions | |
738 | fmovel %d1,%FPCR | |
739 | fsubs #0x00800000,%fp1 | |
740 | bsr sto_cos |store cosine result | |
741 | fmovel (%sp)+,%FPCR |restore users exceptions | |
742 | fmovex X(%a6),%fp0 | |
743 | bra t_frcinx | |
744 | ||
745 | |end |