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